% This script produces the plots which require sampling. % In particular, the incidence boxplot and the sensitivity coefficient plots % The boxplot uses fixed values of the number of contacts and the size of the intervention -- 4 cases. % The sensitivity plots vary the number of contacts and the intervention size % As a result, the sample is used twice. % Location and Size of Output OutFolder=['OutSim' date]; %Relative path is from current directory FigUnit='centimeters'; FigSize=[8.2 8.2]; %How much of the work should be done? produceBox=1; recomputeBox=0; produceSens=1; recomputeSens=0; resample=0; NumSamp=10000; % Fixed parameters -- If any of these are in the sample, they will be overwritten. % Parameters in the sample need not be included here %Incoming fractions. These come from numbers for 15 year olds in Neil's work eta_f = .25; eta_m = 0.44; eta_x = 0.234; eta_y = 0.06; eta_if = 0; eta_ix = 0.014; eta_im = 0; eta_iy = 0; %Population growth and removal rates n0=40000; mu=1/15; %Activity Pattern g_m1 = .8; g_f1 = .8; g_y1 = 0.2; g_x1 = 0.2; g_m2=1-g_m1; g_f2=1-g_f1; g_y2=1-g_y1; g_x2=1-g_x1; %Initial condition from Neil's work yIC=(n0/mu)*[1378 947 2560-1378 3539-947 69 79 187-69 826-79]'./(2*[2560+187 3539+826 2560+187 3539+826 2560+187 3539+826 2560+187 3539+826]'); %These are computed to give the choice disability distribution for the entire 15-30 year old %population as seen in Neil2012. %These alpha values should make the IC above a steady state for choice (not for HIV) (In the absence of interventions which raise status) alpha_m=n0*(eta_m+eta_im)/(yIC(1)+yIC(5)) - mu; alpha_f=n0*(eta_f+eta_if)/(yIC(2)+yIC(6)) - mu; alpha_im=alpha_m; alpha_if=alpha_f; %Baseline c values for boxplots c_m1=80; c_f1=80; c_y1=80 ;c_x1=80; c_m2=80; c_f2=80; c_y2=240 ;c_x2=240; %Baseline alpha values for boxplots alpha_y1=0; alpha_x1=0; alpha_iy1=0; alpha_ix1=0; alpha_y2=2; alpha_x2=2; alpha_iy2=2; alpha_ix2=2; % Sampled Distribution specification -- Found in a text file % We are using uniform independent distributions. % name should exactly match the parameter names used elsewhere in this document. % start and finish give the max and min values for uniform distributions % dispName appears in the plots DistSpec='parapdfs_New.txt'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Parameter values are all set above this line. % For more changes to appearance of figures play with plotting commands below. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tic %Initialize output directory. This statement sometimes causes trouble on other systems. %If necessary remove it and create the output folder manually if exist(OutFolder, 'dir') ~= 7; mkdir(OutFolder); end; %Load Distribution Specification [name, start, finish, dispName] = textread(DistSpec,'%s %f %f %q','headerlines',1); %Initialize Figures %F1 is a single pane as specified above F1=figure(); set(F1, 'PaperUnits', FigUnit); set(F1,'PaperSize',FigSize); set(F1,'PaperPosition',[0,0,FigSize(2),FigSize(1)]); %Here are a few more possible figure settings set(F1, 'DefaultTextFontSize', 9); set(F1, 'DefaultAxesFontSize', 9); set(F1, 'DefaultAxesFontName', 'Times New Roman'); set(F1, 'DefaultTextFontName', 'Times New Roman'); %F2 is two panes arranged vertically F2=figure(); set(F2, 'PaperUnits', FigUnit); set(F2,'PaperSize',FigSize*diag([2 1])); set(F2,'PaperPosition',[0,0,FigSize(2),2*FigSize(1)]); set(F2, 'DefaultTextFontSize', 9); set(F2, 'DefaultAxesFontSize', 9); set(F2, 'DefaultAxesFontName', 'Times New Roman'); set(F2, 'DefaultTextFontName', 'Times New Roman'); if resample numpar=length(start); %Find Sample Sample=lhsdesign(NumSamp,numpar); %Scale and shift to give correct start and finish Sample=repmat(start',NumSamp,1)+Sample*diag(finish-start); %Save Sample save([OutFolder '/InputSample'], 'Sample') else %Load Sample load([OutFolder '/InputSample'], 'Sample') %Check against specification if min(Sample)finish' %Note to self -- There are ways the sample could fail to match that would not trigger this warning -- consider an update later. disp('Warning: Loaded Sample does not match specification') end if length(Sample(:,1))~=NumSamp disp('Warning: Loaded Sample size does not match NumSamp. Updating to prevent crashing') NumSampOld=NumSamp; NumSamp=length(Sample(:,1)); end end if produceBox if recomputeBox %Initialize output matrix OutBox=zeros(NumSamp,4); %we will be computing incidence for 4 cases. for i=1:4 for j=1:NumSamp %Import the sampled parameters for k=1:length(start) eval([char(name(k)) '=' num2str(Sample(j,k)) ';']); %This statement sets each of the parameters listed in 'name' to its sampled value -- use of eval is sketchy, end %Here we're undoing what we just did for the parameters that we want fixed. if i==1 % Low CD contact, no support c_m=c_m1; c_f=c_f1; c_y=c_y1; c_x=c_x1 ; alpha_y=alpha_y1; alpha_x=alpha_x1; alpha_iy=alpha_iy1; alpha_ix=alpha_ix1; elseif i==2 %Low CD contact, support c_m=c_m1; c_f=c_f1; c_y=c_y1; c_x=c_x1 ; alpha_y=alpha_y2; alpha_x=alpha_x2; alpha_iy=alpha_iy2; alpha_ix=alpha_ix2; elseif i==3 %High CD contact, no support c_m=c_m2; c_f=c_f2; c_y=c_y2; c_x=c_x2 ; alpha_y=alpha_y1; alpha_x=alpha_x1; alpha_iy=alpha_iy1; alpha_ix=alpha_ix1; elseif i==4 %High CD contact, support c_m=c_m2; c_f=c_f2; c_y=c_y2; c_x=c_x2 ; alpha_y=alpha_y2; alpha_x=alpha_x2; alpha_iy=alpha_iy2; alpha_ix=alpha_ix2; end %All the parameters should be set appropriately. Now compute the outcomes: F=@(t,y) system_no_recovered(t, y, g_m1, g_f1, g_x1, g_y1, g_m2, g_f2, g_y2, g_x2, beta_mf, beta_fm, beta_fy, beta_yf, beta_yx, beta_xy, beta_mx, beta_xm, eta_m, eta_f, eta_y, eta_x, eta_im, eta_if, eta_iy, eta_ix, n0, alpha_f, alpha_x, alpha_m, alpha_y,alpha_if, alpha_ix,alpha_im, alpha_iy,c_m, c_f, c_y, c_x, mu); [t,y]=ode45(F, [0,15], yIC); %Yearly instantaneous incidence is given by the system with non-disease transfer rates equal to zero. Inc=system_no_recovered(0,y(end,:)',g_m1, g_f1, g_x1, g_y1, g_m2, g_f2, g_y2, g_x2, beta_mf, beta_fm, beta_fy, beta_yf, beta_yx, beta_xy, beta_mx, beta_xm, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,c_m, c_f, c_y, c_x, 0); OutBox(j,i)=sum(Inc(5:8))/sum(y(end,:)); end end %Save BoxOutput save([OutFolder '/BoxSample'], 'OutBox') else %Load BoxOutput load([OutFolder '/BoxSample'], 'OutBox'); end figure(F1) boxplot(OutBox, 'labels', {'A', 'B', 'C', 'D'}); A=findobj(gca,'Tag','Upper Whisker'); set(A,'LineStyle','-') A=findobj(gca,'Tag','Lower Whisker'); set(A,'LineStyle','-') print('-depsc',[OutFolder '/Box1']) end if produceSens if recomputeSens %Initialize output matrix OutSens=zeros(NumSamp,2); %Output is prevalence for men and women. for i=1:NumSamp%Import the sampled parameters for k=1:length(start)%This statement sets each of the parameters listed in 'name' to its sampled value -- use of eval is sketchy, eval([char(name(k)) '=' num2str(Sample(i,k)) ';']); end %All the parameters should be set appropriately. Now compute the outcomes: F=@(t,y) system_no_recovered(t, y, g_m1, g_f1, g_x1, g_y1, g_m2, g_f2, g_y2, g_x2, beta_mf, beta_fm, beta_fy, beta_yf, beta_yx, beta_xy, beta_mx, beta_xm, eta_m, eta_f, eta_y, eta_x, eta_im, eta_if, eta_iy, eta_ix, n0, alpha_f, alpha_x, alpha_m, alpha_y,alpha_if, alpha_ix,alpha_im, alpha_iy,c_m, c_f, c_y, c_x, mu); [t,y]=ode45(F, [0,15], yIC); OutSens(i,1)=sum(y(end,[5 7]))/sum(y(end,[1 3 5 7])); OutSens(i,2)=sum(y(end,[6 8]))/sum(y(end,[2 4 6 8])); end %Save OutSens SensSample=Sample; save([OutFolder '/SensSample'], 'OutSens', 'SensSample') else %Load BoxOutput load([OutFolder '/SensSample'], 'OutSens', 'SensSample'); end rho=zeros(length(OutSens(1,:)), length(start)); for j=1:length(start); rho(:,j)=partialcorr(OutSens,SensSample(:,j),SensSample(:,[1:j-1 j+1:length(SensSample(1,:))]), 'type','Spearman'); end %%Modify for divided printing format rho=[rho(:,1:8) zeros(2,2) rho(:,9:12) zeros(2,2) rho(:,13:16) zeros(2,1)]; dispName=[dispName(1:8); {'Transmission: '; ' '}; dispName(9:12); {'Support: '; ' '}; dispName(13:16); {'Contacts: '}]; rholength=length(rho(1,:)); figure(F1) barh(rho(1,:)); axis([-1 1 0 rholength]) set(gca,'XGrid','on') set(gca,'YTick',1:rholength) set(gca,'YTickLabel',dispName) print('-depsc',[OutFolder '/SensM']) barh(rho(2,:)); axis([-1 1 0 rholength]) set(gca,'XGrid','on') set(gca,'YTick',1:length(rholength)) set(gca,'YTickLabel',dispName) print('-depsc',[OutFolder '/SensF']) figure(F2) subplot(2,1,1) barh(rho(1,:)); hold on; plot([-1 1],[9.5 9.5],'k:'); plot([-1 1],[15.5 15.5],'k:'); axis([-1 1 0 rholength]) set(gca,'XGrid','on') set(gca,'YTick',1:rholength) set(gca,'YTickLabel',dispName) title('A') subplot(2,1,2) barh(rho(2,:)); hold on; plot([-1 1],[9.5 9.5],'k:'); plot([-1 1],[15.5 15.5],'k:'); axis([-1 1 0 rholength]) set(gca,'XGrid','on') set(gca,'YTick',1:rholength) set(gca,'YTickLabel',dispName) title('B') print('-depsc',[OutFolder '/Sens2']) end toc