From 03e663aa7d79c883abbb996e27e03b6ebd560341 Mon Sep 17 00:00:00 2001 From: kgostic Date: Mon, 1 Apr 2019 11:11:40 -0700 Subject: [PATCH] Updated bug in AIC calculation, and loaded missing package. --- .Rhistory | 512 ++++++++++++++++++ Full_analysis_script_revised_submission.R | 3 +- code_outputs/sessionInfo.txt | 28 +- manuscript_figures/Fig2_par_ests.pdf | Bin 13849 -> 13849 bytes manuscript_figures/Fig3_model_fits.pdf | Bin 31420 -> 31420 bytes .../Fig4_prob_inf_given_dose_and_route.pdf | Bin 10659 -> 10651 bytes scratch_figures/Profiles_basic_model.pdf | Bin 29126 -> 29126 bytes scratch_figures/pI_profiles.pdf | Bin 8034 -> 8034 bytes 8 files changed, 529 insertions(+), 14 deletions(-) create mode 100644 .Rhistory diff --git a/.Rhistory b/.Rhistory new file mode 100644 index 0000000..07a71c9 --- /dev/null +++ b/.Rhistory @@ -0,0 +1,512 @@ +## Base plot +plotdf = data.frame(alpha = aas.rat, beta = BBs.rat, nll = prof2d.rat) ## 2d profiles +plotdf$line_alphas = opt.pc.abraded.rat.skin$minimum*plotdf$beta/(1-opt.pc.abraded.rat.skin$minimum) ## Add coords for line representing alpha/(alpha+beta) = 0.02 +linepts = which(abs(plotdf$alpha/(plotdf$alpha+plotdf$beta) - opt.pc.abraded.rat.skin$minimum) < 5e-5)## extract (alpha, beta) pairs that fall on the line +linepts = linepts[c(1, 15, 34)] +plotdf[linepts, ] +pp = ggplot(subset(plotdf, alpha <=1000))+geom_point(aes(x = alpha, y = beta, color = nll)) + +scale_color_distiller(palette = "Spectral") + +geom_contour(aes(x = alpha, y = beta, z = nll), binwidth = 2, color = 'gray20') + +geom_line(aes(x = line_alphas, y = beta), color = 'white', lwd = 2) + +geom_point(data = plotdf[linepts[1],], aes(x = alpha , y = beta), shape = 7, size = 2)+ +geom_point(data = plotdf[linepts[2],], aes(x = alpha , y = beta), shape = 8, size = 2)+ +geom_point(data = plotdf[linepts[3],], aes(x = alpha , y = beta), shape = 5, size = 2)+ +geom_text(data = plotdf[linepts[1],], aes(x = alpha +90, y = beta, label = round(nll, 4))) + +geom_text(data = plotdf[linepts[2],], aes(x = alpha +90, y = beta, label = round(nll, 4))) + +geom_text(data = plotdf[linepts[3],], aes(x = alpha +90, y = beta, label = round(nll, 4))) + +ggtitle('2D likelihood surface', subtitle = 'mixture model fitted to data from rats, abraded skin')+ +xlab(expression(alpha)) + +ylab(expression(beta)) + +theme_bw() +pp +## Plot beta CDF for the three focal points +xx = seq(0, .05, by = 0.0001) +bdlist = data.frame(xx = xx, +pt1 = pbeta(q = xx, shape1 = plotdf[linepts[3],'alpha'], shape2 = plotdf[linepts[3],'beta']), +pt2 = pbeta(q = xx, shape1 = plotdf[linepts[2],'alpha'], shape2 = plotdf[linepts[2],'beta']), +pt3 = pbeta(q = xx, shape1 = plotdf[linepts[1],'alpha'], shape2 = plotdf[linepts[1],'beta'])) +bdlist = melt(bdlist, id.var = 'xx') +bdplot = ggplot(data = bdlist)+ +geom_point(aes(x = xx, y = value, shape = variable, color = variable), size = 2, alpha = .5)+ +geom_line(aes(x = xx, y = value, color = variable), alpha = .7)+ +scale_color_discrete(name = "", +breaks=c("pt1", "pt2", "pt3"), +labels=c("upper point", "middle point", "lower point")) + +scale_shape_manual(name = "", +values = c(5,8,7), +breaks=c("pt1", "pt2", "pt3"), +labels=c("upper point", "middle point", "lower point")) + +xlab('pc')+ +ylab('cumulative density')+ +ggtitle('Fitted Beta(\u03B1,\u03B2) distribution', subtitle = 'selected points in top panel specify (\u03B1, \u03B2) values')+ +theme_bw() +bdplot +outplot = arrangeGrob(pp, bdplot, nrow = 2) +ggsave(filename = figS1, plot = outplot, width = 5, height = 7, dpi = 'print', device = png()) +## Calculate delta AIC for fits to abraded skin data +## Mixture model vs basic model in hamsters +## Not possible to do this comparison in rats, because there was no real MLE in rat mixture model +AICcalc = function(n.pars, nll){ +2*n.pars+2*nll +} +AIC.basic.ham = AICcalc(n.pars = 2, nll = opt.pc.abraded.skin$value) +AIC.mixture.ham = AICcalc(n.pars = 3, nll = opt.aa.BB.abraded.skin$value) +del.AIC.hamster = c(AIC.basic.ham, AIC.mixture.ham); del.AIC.hamster = del.AIC.hamster-min(del.AIC.hamster) +del.AIC.hamster +#################################################################### +########################### SECTION 4 ########################## +#### +#### 4. Plot results! +#### +#################################################################### +#################################################################### +## Plot point estimates and CIs of each free paramter, as well as estimated Beta ditribution of pc values +{ +pdf(fig2) +par(mgp = c(2, 1, 0), mar = c(3,7,2,2)) +layout(mat = matrix(c(1,1,2,3), byrow = T, nrow = 2, ncol = 2)) +ham.col = 'dodgerblue3' +rat.col = 'gold' +## Plot1 - point estimates and CIs ---------- +xtext = c(-1.9, -8.5, -7.1, -5.4) +ytext = c(1.001, 1.001, 1.001, 1.001)+.0005 +#cols = c('firebrick1', 'deepskyblue', 'limegreen', 'chocolate1') +labs = c(expression('p'[c]*' abraded'), +expression('p'[c]*' intact'), +expression('p'[c]*' shaved'), +expression('p'[c]*' conjunctival'), +expression('p'[p]), +expression('p'[c]*' abraded'), +expression('p'[p])) +ests = c(opt.pc.abraded.skin$minimum, opt.pc.intact.skin$minimum, opt.pc.shaved.skin$minimum, opt.pc.conjunctival.skin$minimum, pP, opt.pc.abraded.rat.skin$minimum, pP.rat) +CI.lows = c(CI.abraded[1], CI.intact[1], CI.shaved[1], CI.conjunctival[1], CI.pP[1], CI.abraded.rat[1], CI.pP.rat[1]) +CI.highs = c(CI.abraded[2], CI.intact[2], CI.shaved[2], CI.conjunctival[2], CI.pP[2], CI.abraded.rat[2], CI.pP.rat[2]) +yvals = length(ests):1 +plot(log10(ests), yvals, pch = c(rep(21, 5), 22,22), bg = c(rep(ham.col, 5), rat.col, rat.col), col = c(rep('black', 5), 'gray', 'gray'), xlab = 'log10(value)', ylab = '', yaxt = 'n', bty = 'n', xlim = c(-10, 0)) +segments(x0 = log10(CI.lows), x1 = log10(CI.highs), y0 = yvals, bty = 'n', col = c(rep(ham.col, 5), rat.col, rat.col)) +axis(2, at = yvals, labels = labs, las = 2) +legend(-10, 2.5, c('hamsters', 'rats'), pch = c(21,22), yjust = .5, col = c('black', 'gray'), pt.bg = c(ham.col, rat.col), bty = 'n') +abline(h = 2.5, lty = 2) +mtext(3, text = 'A', at = -12, line = 1, font = 2, xpd = NA) +xlims = par('usr')[1:2] +par(mar = c(3, 4,3,2)) +xx = seq(0, 1, by = 0.001) +plot((xx), (dbeta(xx, shape1 = opt.aa.BB.abraded.skin$par['aa'], shape2 = opt.aa.BB.abraded.skin$par['bb'])), lty = 1, pch = 21, type = 'b', lwd = .5, ylab = 'density', xlab = expression('p'[c]), bg = ham.col) +text(x = -.3, 36, 'B', font = 2, xpd = NA) +plot((xx), (dbeta(xx, shape1 = opt.aa.BB.abraded.rat.skin$par['aa'], shape2 = opt.aa.BB.abraded.rat.skin$par['bb'])), lty = 1, bg = rat.col, pch = 22, type = 'b', lwd = .5, ylab = 'density', xlab = expression('p'[c]), col = 'gray') +text(x = -.3, 655, 'C', font = 2, xpd = NA) +dev.off() +} +########################## Plot data vs. model fits ############################ +##### Fig. 3: +## Reformat observed data for plotting +dat = reformat.route(intact) +intact.CIs = binom.confint(dat[,'died'], dat[,'n'], conf.level = .95, methods = 'wilson')[,c('lower', 'upper')] +dat = reformat.route(shaved) +shaved.CIs = binom.confint(dat[,'died'], dat[,'n'], conf.level = .95, methods = 'wilson')[,c('lower', 'upper')] +dat = reformat.route(abraded) +abraded.CIs = binom.confint(dat[,'died'], dat[,'n'], conf.level = .95, methods = 'wilson')[,c('lower', 'upper')] +dat = reformat.route(conjunctival) +conjunctival.CIs = binom.confint(dat[,'died'], dat[,'n'], conf.level = .95, methods = 'wilson')[,c('lower', 'upper')] +IP.CIs = binom.confint(IP[,'died'], IP[,'n'], conf.level = .95, methods = 'wilson')[,c('lower', 'upper')] +ratdat = reformat.route(abraded.rat) +rat.CIs = binom.confint(ratdat[,'died'], ratdat[,'n'], conf.level = .95, methods = 'wilson')[,c('lower', 'upper')] +IP.CIs.rat = binom.confint(IP.rat[,'died'], IP.rat[,'n'], conf.level = .95, methods = 'wilson')[,c('lower', 'upper')] +#Define a function to plot probability of infection at different D0s, given parameter estimates +p.inf = function(doses, best.pc, best.pp = pP){ +1-exp(-doses*best.pc*best.pp) +} +## Repeat for prob infection, given IP inoculation +p.inf.IP = function(dose, best.pp){ +1-exp(-dose*best.pp) +} +## Repeat for prob infection, given mixture model +p.inf.mixture = function(dd.in, pp.in, aa.in, bb.in, sp = 'hamsters'){ +sapply(dd.in, function(DD) {p_inf_fun(DD, pp.in, aa.in, bb.in, sp)}) +} +# Convert color to color + transparency +col2alpha <- function(col, alpha) { +col_rgb <- col2rgb(col)/255 +rgb(col_rgb[1], col_rgb[2], col_rgb[3], alpha = alpha) +} +{ +pdf(fig3, height = 5, width = 7) +par(mfrow = c(2,3), mgp = c(2,1,0), mar = c(3,3,3,1)) +dd = round(10^seq(0, 9, by = .01)) +plot(dd, p.inf(dd, opt.pc.intact.skin$minimum), ylab = 'Probability of Infection', xlab = 'Dose', type = 'l', lwd = 2, ylim = c(0,1), log = 'x', cex.lab = 1.1, xlim = c(1, 1e8)) +#Plot CI bars +segments(x0 = intact[,'dose'], y0 = intact.CIs[,1], y1 = intact.CIs[,2]) +segments(x0 = intact[,'dose']*.7, y0 = intact.CIs[,1], x1 = intact[,'dose']*1.3) +segments(x0 = intact[,'dose']*.7, y0 = intact.CIs[,2], x1 = intact[,'dose']*1.3) +# data +points(intact[,'dose'], intact[,'died']/intact[,'n'], pch = 21, cex = 1.5, bg = ham.col) +text(.1, 1.15, 'A', font = 2, xpd = NA) +#Shaved +plot(dd, p.inf(dd, opt.pc.shaved.skin$minimum, pP), ylab = 'Probability of Infection', xlab = 'Dose', type = 'l', lwd = 2, ylim = c(0,1), log = 'x', cex.lab = 1.1, xlim = c(1, 1e8)) +#Plot CI bars +segments(x0 = shaved[,'dose'], y0 = shaved.CIs[,1], y1 = shaved.CIs[,2]) +segments(x0 = shaved[,'dose']*.7, y0 = shaved.CIs[,1], x1 = shaved[,'dose']*1.3) +segments(x0 = shaved[,'dose']*.7, y0 = shaved.CIs[,2], x1 = shaved[,'dose']*1.3) +# data +points(shaved[,'dose'], shaved[,'died']/shaved[,'n'], pch = 21, cex = 1.5, bg = ham.col) +text(.1, 1.15, 'B', font = 2, xpd = NA) +#Abraded +## Basic model fit and data +plot(dd, p.inf(dd, opt.pc.abraded.skin$minimum, pP), ylab = 'Probability of Infection', xlab = 'Dose', type = 'l', lwd = 2, ylim = c(0,1), log = 'x', cex.lab = 1.1, xlim = c(1, 1e8)) +yy.mixture.plot = p.inf.mixture(dd.in = dd, aa.in = opt.aa.BB.abraded.skin$par['aa'], bb.in = opt.aa.BB.abraded.skin$par['bb'], pp.in = pP) +lines(dd, yy.mixture.plot, lty = 3, lwd = 2) +#Plot data CI bars +segments(x0 = abraded[,'dose'], y0 = abraded.CIs[,1], y1 = abraded.CIs[,2]) +segments(x0 = abraded[,'dose']*.7, y0 = abraded.CIs[,1], x1 = abraded[,'dose']*1.3) +segments(x0 = abraded[,'dose']*.7, y0 = abraded.CIs[,2], x1 = abraded[,'dose']*1.3) +# data +points(abraded[,'dose'], abraded[,'died']/abraded[,'n'], pch = 21, cex = 1.5, bg = ham.col) +legend('bottomright', legend = c('basic model', 'mixture model'), lty = c(1, 3), bty = 'n', lwd = 2) +text(.1, 1.15, 'C', font = 2, xpd = NA) +#Conjunctval +plot(dd, p.inf(dd, opt.pc.conjunctival.skin$minimum), ylab = 'Probability of Infection', xlab = 'Dose', type = 'l', lwd = 2, ylim = c(0,1), log = 'x', cex.lab = 1.1, xlim = c(1, 1e8)) +#Plot CI bars +segments(x0 = conjunctival[,'dose'], y0 = conjunctival.CIs[,1], y1 = conjunctival.CIs[,2]) +segments(x0 = conjunctival[,'dose']*.7, y0 = conjunctival.CIs[,1], x1 = conjunctival[,'dose']*1.3) +segments(x0 = conjunctival[,'dose']*.7, y0 = conjunctival.CIs[,2], x1 = conjunctival[,'dose']*1.3) +# data +points(conjunctival[,'dose'], conjunctival[,'died']/conjunctival[,'n'], pch = 21, cex = 1.5, bg = ham.col) +text(.1, 1.15, 'D', font = 2, xpd = NA) +#Abraded, rats +plot(dd, p.inf(dd, best.pc = opt.pc.abraded.rat.skin$minimum, best.pp = pP), ylab = 'Probability of Infection', xlab = 'Dose', type = 'l', lwd = 2, ylim = c(0,1), log = 'x', cex.lab = 1.1, xlim = c(1, 1e8), col = 'gray55') +# Mixture model fits +yy.mixture.plot.rat = p.inf.mixture(dd.in = dd, aa.in = opt.aa.BB.abraded.rat.skin$par['aa'], bb.in = opt.aa.BB.abraded.rat.skin$par['bb'], pp.in =pP.rat, sp = 'rats') +lines(dd, yy.mixture.plot.rat, lty = 3, lwd = 2, col = 'gray55') +#Plot data CI bars +segments(x0 = abraded.rat[,'dose'], y0 = rat.CIs[,1], y1 = rat.CIs[,2], col = 'gray55') +segments(x0 = abraded.rat[,'dose']*.7, y0 = rat.CIs[,1], x1 = abraded.rat[,'dose']*1.3, col = 'gray55') +segments(x0 = abraded.rat[,'dose']*.7, y0 = rat.CIs[,2], x1 = abraded.rat[,'dose']*1.3, col = 'gray55') +# Data +points(abraded.rat[,'dose'], abraded.rat[,'died']/abraded.rat[,'n'], pch = 22, cex = 1.5, bg = rat.col, col = 'gray') +legend('bottomright', legend = c('basic model', 'mixture model'), lty = c(1, 3), bty = 'n', lwd = c(1,2), col = 'gray55') +text(.1, 1.15, 'E', font = 2, xpd = NA) +#IP, hamsters +plot(dd, p.inf.IP(dose = dd, best.pp = pP), ylab = 'Probability of Infection', xlab = 'Dose', type = 'l', ylim = c(0,1), log = 'x', cex.lab = 1.1, xlim = c(1, 1e8), lwd = 2) +segments(x0 = IP[,'dose'], y0 = IP.CIs[,1], y1 = IP.CIs[,2], lwd = 1) +segments(x0 = IP[,'dose']*.7, y0 = IP.CIs[,1], x1 = IP[,'dose']*1.3, lwd = 1) +segments(x0 = IP[,'dose']*.7, y0 = IP.CIs[,2], x1 = IP[,'dose']*1.3, lwd = 1) +points(IP[,'dose'], IP[,'died']/IP[,'n'], pch = 21, cex = 1.7, bg = ham.col) +lines(dd, p.inf.IP(dd, pP.rat), lty = 1, lwd = 1.3, col = 'gray55') +xx.jitter = 10^(log10(IP.rat[,'dose'])+.1) # Jitter rat CIs so you can see CIs from both species +segments(x0 = xx.jitter, y0 = IP.CIs.rat[,1], y1 = IP.CIs.rat[,2], col = 'gray50', lwd = .7) +segments(x0 = IP.rat[,'dose'], y0 = IP.CIs.rat[,1], x1 = IP.rat[,'dose']*1.5, col = 'gray50', lwd = .7) +segments(x0 = IP.rat[,'dose'], y0 = IP.CIs.rat[,2], x1 = IP.rat[,'dose']*1.5, col = 'gray50', lwd = .7) +points(xx.jitter, IP.rat[,'died']/IP.rat[,'n'], pch = 22, cex = 1.2, bg = rat.col, col = 'gray') +legend('bottomright', legend = c('hamsters', 'rats'), lty = 1, bty = 'n', lwd = c(2,2), pch = c(21, 22), pt.bg = c(ham.col, rat.col), col = c('black', 'gray50')) +text(.1, 1.15, 'F', font = 2, xpd = NA) +dev.off() +} +{ +pdf(fig4, height = 7) +# Define a function to make transparent colors +tns.col = function(col.in, tns.pct){ +rgbvals = col2rgb(col.in) +tnscol = rgb(rgbvals[1], rgbvals[2], rgbvals[3], tns.pct/100*255,maxColorValue = 255) +tnscol +} +par(mfrow = c(2,1), mgp = c(2,1,0), mar = c(3,3,3,3)) +dd = c(1:8, 10^seq(1.7, 11, by = .15)) +p.inf.i = p.inf(dd, opt.pc.intact.skin$minimum); p.inf.i[which(!is.finite(p.inf.i))] = 1 +plot(dd, p.inf.i, ylab = 'P(Infection), hamsters', xlab = 'Dose', type = 'l', col = cols[1], ylim = c(0,1.3), log = 'x', cex.lab = 1.1, lwd = 3, bty= 'n', yaxt = 'n', xaxt = 'n', lty = 3) +axis(1, at = 10^(1:9), labels = 10^(1:9)) +lower = p.inf(dd, CI.intact[1]); lower[which(!is.finite(lower))] = 1 +upper = p.inf(dd, CI.intact[2]); upper[which(!is.finite(upper))] = 1 +polygon(x = c(dd, rev(dd)), y = c(lower, rev(upper)), border = NA, col = col2alpha(cols[1], .5)) +axis(side = 2, at = seq(0, 1, by = .1)) +## Initialize storage and calculate bounds -- doses at which infection becomes possible, or certain +bounds = matrix(NA, nrow = 2, ncol = 5, dimnames = list(c('min', 'always'), c('intact', 'shaved', 'conj', 'abraded', 'IP'))) +bounds[,1] = dd[c(max(which(p.inf.i < 0.1)), min(which(p.inf.i > 0.975)))] +#Shaved +p.inf.s = p.inf(dd, opt.pc.shaved.skin$minimum); p.inf.s[which(!is.finite(p.inf.s))] = 1 +lines(dd, p.inf.s, col = cols[2], lwd = 3, lty = 3) +lower = p.inf(dd, CI.shaved[1], best.pp = pP); lower[which(!is.finite(lower))] = 1 +upper = p.inf(dd, CI.shaved[2], best.pp = pP); upper[which(!is.finite(upper))] = 1 +polygon(x = c(dd, rev(dd)), y = c(lower, rev(upper)), border = NA, col = col2alpha(cols[2], .5)) +bounds[,2] = dd[c(max(which(p.inf.s < 0.1)), min(which(p.inf.s > 0.975)))] +#Abraded +yy.mixture.plot = p.inf.mixture(dd.in = dd, pp.in = pP, aa.in = opt.aa.BB.abraded.skin$par[1], bb.in = opt.aa.BB.abraded.skin$par[2]) +lines(dd, yy.mixture.plot, col = cols[3], lty = 4, lwd = 3) +############ Mixture model CIs for plotting ############### +# In the mixture model, the confidence envelope is defined by a vague elliptical curve of (alpha, beta) pairs +# Pairs on this curve have likelihoods near the threshold of the best nll + a threshold determined by the likelihood ratio method +# Above, we defined prof.acc as the set of all (alpha, beta) pairs whose neg log lilelihoods were under the threshold value +# Here, we'll downsample those values to extract pairs whose likelihood is very close to the threshold +# Randomly sample points on or near this curve from the 2D profile +# Sample 50 points per unit on the BB axis +# Only sample points whose likelihoods are close to the threshold. THese should give the worst acceptable fits. +aa.coords.reduced = bb.coords.reduced = prof.vals.reduced = NULL +threshold = LR.Threshold(opt.aa.BB.abraded.skin$value, 2) +envelope = NULL +for(ii in 1:8){ +## The first condition verifies points near the threshold +## The second condition verifies that we're sampling relatively evenly from across the full range of beta values +valid = which(prof.acc$prof>threshold-.01 & prof.acc$bb > (-1+ii) & prof.acc$bb<= (0+ii)) +smp = sample(x = valid, size = min(25, length(valid)), replace = FALSE) +envelope = rbind(envelope, prof.acc[smp, ]) +} +# plot(envelope$aa, envelope$bb) ## The reduced set of points fall on the boundary of the confidence envelope +# ## Now, calculate model-predicted probabilities of infection for each dose, using each pair of (alpha, beta) values in the envelope +wrapper = function(AA, BB){ +tryCatch({ +out = p.inf.mixture(dd.in = dd, pp.in = pP, aa.in =AA, bb.in = BB) }, +error = function(e){ +print(e) +out = rep(NA, 71) +}, +warning = function(w){ +print(w) +}, +finally = function() {return(out)} +) +} +yy.CIs = mapply(FUN = wrapper, AA = envelope$aa, BB = envelope$bb) +#### Define CI envelope for plotting as the highest and lowest y value resulting from estimates using all sampled points from the CI envelope +yy.low = apply(X = yy.CIs, MARGIN = 1, FUN = min) # Estimates correspond to doses in dd +yy.high = apply(X = yy.CIs, MARGIN = 1, FUN = max) +polygon(x = c(dd, rev(dd)), +y = c(yy.low, rev(yy.high)), +col = tns.col(cols[3], 20), border = NA) +bounds[,4] = dd[c(1, min(which(yy.mixture.plot > 0.975)))]; #c(dd[min(1, which.max(yy.high <0.1))], 1e11) +#Conjunctval +p.inf.c = p.inf(dd, opt.pc.conjunctival.skin$minimum); p.inf.c[which(!is.finite(p.inf.c))] = 1 +lines(dd, p.inf.c, col = cols[4], lwd = 3) +lower = p.inf(dd, CI.conjunctival[1]); lower[which(!is.finite(lower))] = 1 +upper = p.inf(dd, CI.conjunctival[2]); upper[which(!is.finite(upper))] = 1 +polygon(x = c(dd, rev(dd)), y = c(lower, rev(upper)), border = NA, col = col2alpha(cols[4], .5)) +bounds[,3] = dd[c(max(which(p.inf.c < 0.1)), min(which(p.inf.c > 0.975)))] +## IP inoculation +dd = c(1:9, 10^seq(.95, 11, by = .15)) +p.IP = p.inf.IP(dose = dd, best.pp = pP); p.IP[which(!is.finite(p.IP))] = 1 +lines(dd, p.IP, lwd = 3) +lower = p.inf.IP(dose = dd, best.pp = CI.pP[1]) +upper = p.inf.IP(dose = dd, best.pp = CI.pP[2]) +polygon(x = c(dd, rev(dd)), y = c(lower, rev(upper)), border = NA, col = col2alpha('black', .5)) +bounds[,5] = dd[c(1, min(which(p.IP > 0.975)))] +#cols = c('firebrick1', 'deepskyblue', 'chocolate1', 'limegreen', 'darkgreen') +ys = seq(1.03, 1.27, length = 5)+.05 +segments(x0 = bounds[1,], x1 = bounds[2,], y0 = ys, col = c(cols[c(1,2,4,3)], 'black'), lty = 2, lwd = 2) +arrows(x0 = bounds[2,], x1 = 1e11, y0 = ys, col = c(cols[c(1,2,4,3)], 'black'), lty = 1, length = .1, lwd = 2) +segments(x0 = bounds[1,], y0 = ys, y1 = ys+.02, col = c(cols[c(1,2,4,3)], 'black'), lwd =2) +segments(x0 = bounds[2,], y0 = ys, y1 = ys+.02, col = c(cols[c(1,2,4,3)], 'black'), lwd = 2) +legend(4e9, y = .6, (c('intact', 'shaved', 'conjunctival', 'abraded', 'IP, no barrier')), col = rev(c('black', cols[c(3,4,2,1)])), fill = rev(c('black', cols[c(3,4,2,1)])), cex = .8, border = NA, lty = rev(c(1, 4, 1, 2, 3)), xpd = NA, lwd = rev(c(1,1,2,2,2,2)), bty = 'n') +legend(1e3, y = 1.67, c('infection possible at this dose', 'infection almost certain at this dose'), lty = c(2, 1), cex = .8, xpd = NA, bty = 'n') +text(.1, 1.3, 'A', font = 2, xpd = NA) +#### Repeat for rats +p.inf.a = p.inf(dd, opt.pc.abraded.rat.skin$minimum, best.pp = pP.rat); p.inf.a[which(!is.finite(p.inf.a))] = 1 +plot(dd, p.inf.a, ylab = 'P(Infection), rats', xlab = 'Dose', type = 'l', col = cols[3], ylim = c(0,1.3), log = 'x', cex.lab = 1.1, lwd = 3, bty= 'n', yaxt = 'n', xaxt = 'n', lty = 4) +axis(1, at = 10^(1:9), labels = 10^(1:9)) +lower = p.inf(dd, CI.abraded.rat[1]); lower[which(!is.finite(lower))] = 1 +upper = p.inf(dd, CI.abraded.rat[2]); upper[which(!is.finite(upper))] = 1 +polygon(x = c(dd, rev(dd)), y = c(lower, rev(upper)), border = NA, col = col2alpha(cols[3], .2)) +axis(side = 2, at = seq(0, 1, by = .1)) +## Initialize storage and calculate bounds -- doses at which infection becomes possible, or certain +bounds = matrix(NA, nrow = 2, ncol = 2, dimnames = list(c('min', 'always'), c('abraded', 'IP'))) +bounds[,1] = dd[c(max(which(p.inf.a < 0.1)), min(which(p.inf.a > 0.975)))] +## IP inoculation +lines(dd, p.IP, lwd = 3) +lower = p.inf.IP(dose = dd, best.pp = CI.pP.rat[1]) +upper = p.inf.IP(dose = dd, best.pp = CI.pP.rat[2]) +polygon(x = c(dd, rev(dd)), y = c(lower, rev(upper)), border = NA, col = col2alpha('black', .5)) +bounds[,2] = dd[c(1, min(which(p.IP > 0.975)))] +#cols = c('firebrick1', 'deepskyblue', 'chocolate1', 'limegreen', 'darkgreen') +ys = seq(1.05, 1.25, length = 4)+.05 +ys = ys[1:2] +segments(x0 = bounds[1,], x1 = bounds[2,], y0 = ys, col = c(cols[3], 'black'), lty = 2, xpd = NA, lwd = 2) +arrows(x0 = bounds[2,], x1 = 1e11, y0 = ys, col = c(cols[3], 'black'), lty = 1, length = .1, xpd = NA, lwd = 2) +segments(x0 = bounds[1,], y0 = ys, y1 = ys+.02, col = c(cols[3], 'black'),xpd = NA, lwd = 2) +segments(x0 = bounds[2,], y0 = ys, y1 = ys+.02, col = c(cols[3], 'black'),xpd = NA, lwd = 2) +legend(4e9, y = .6, (c('abraded', 'IP, no barrier')), col = rev(c('black', cols[3])), fill = rev(c('black', cols[3])), cex = .8, border = NA, lty = (c(4,1)), xpd = NA, lwd = rev(c(1,3)), bty = 'n') +legend(1e3, y = 1.5, c('infection possible at this dose', 'infection almost certain at this dose'), lty = c(2, 1), cex = .8, xpd = NA, bty = 'n') +text(.1, 1.3, 'B', font = 2, xpd = NA) +dev.off() +} +### Write a table template for experimental results +ham.data = rbind(intact, shaved, abraded, conjunctival, IP) +rat.data = rbind(abraded.rat, IP.rat) +## Write output files +write.csv(ham.data, file = hamdat_summary, row.names = FALSE) +write.csv(rat.data, file = ratdat_summary, row.names = FALSE) +xx = opt.aa.BB.abraded.rat.skin$par +xx[1]/sum(xx) +opt.pc.abraded.rat.skin +## This script checks numerically that the probability of infection in the basic model can be written as: P(Infection|Dose) = 1-exp(D0*p_p*p_c) +## Long-form equation: +## Let d represent the expected inoculum dose +## Let D0 represent the actual inoculum size +## Let Dw represent the reaches the within-host environment +## Let DI represent the infecting dose +## P(DI|d) = P(D0|d)P(Dw|D0)P(DI|Dw) ## The probability of a given infectious dose, given the expected inoculum size +## P(Infection|d) = sum_DI=1:Inf sum_Dw=DI:Inf sum_D0=Dw:Inf P(D0|d)P(Dw|D0)P(DI|Dw)) ## Long form probability of infection, , given the expected inoculum size +## Calculate the summand after algebraic rearrangement into three different poisson densities +compare = function(d, D0, Dw, DI, pp, pc){ +## This is the long-form equation for P(DI|d), as given in equation 4 of the main text methods +original = dpois(x = D0, lambda = d)*dbinom(x = Dw, size = D0, prob = pc)*dbinom(x = DI, size = Dw, prob = pp) +## This is an alternate representation of the summand, as shown in equation 5 of the main text methods +rearranged = dpois(x = DI, lambda = d*pc*pp)*dpois(x = D0-Dw, lambda = d*(1-pc))*dpois(x = Dw-DI, lambda = d*pc*(1-pp)) +return(c(original, rearranged)) +} +## Test the comparison function, which will return the value of the summand using the original long-form, and the algebraic rearrangement given in equation 5 +compare(1000, 998, 100, 20, .15, .001) +## Test for 1000 randomly chosen values of d, D0, Dw, DI, pp and pc +d.test = runif(n = 1000, 1, 10^9) +D0.test = sapply(d.test, function(dd) rpois(n = 1, dd)) +Dw.test = sapply(D0.test, function(D0) runif(1, 1, D0)) +DI.test = sapply(Dw.test, function(Dw) runif(1, 1, Dw)) +pp = runif(1000, 0, 1) +pc = runif(1000, 0, 1) +## Output comparison +test.outputs = mapply(FUN = compare, d= d.test, D0 = D0.test, Dw = Dw.test, DI = DI.test, pp = pp.test, pc = pc.test) +## This script checks numerically that the probability of infection in the basic model can be written as: P(Infection|Dose) = 1-exp(D0*p_p*p_c) +## Long-form equation: +## Let d represent the expected inoculum dose +## Let D0 represent the actual inoculum size +## Let Dw represent the reaches the within-host environment +## Let DI represent the infecting dose +## P(DI|d) = P(D0|d)P(Dw|D0)P(DI|Dw) ## The probability of a given infectious dose, given the expected inoculum size +## P(Infection|d) = sum_DI=1:Inf sum_Dw=DI:Inf sum_D0=Dw:Inf P(D0|d)P(Dw|D0)P(DI|Dw)) ## Long form probability of infection, , given the expected inoculum size +## Calculate the summand after algebraic rearrangement into three different poisson densities +compare = function(d, D0, Dw, DI, pp, pc){ +## This is the long-form equation for P(DI|d), as given in equation 4 of the main text methods +original = dpois(x = D0, lambda = d)*dbinom(x = Dw, size = D0, prob = pc)*dbinom(x = DI, size = Dw, prob = pp) +## This is an alternate representation of the summand, as shown in equation 5 of the main text methods +rearranged = dpois(x = DI, lambda = d*pc*pp)*dpois(x = D0-Dw, lambda = d*(1-pc))*dpois(x = Dw-DI, lambda = d*pc*(1-pp)) +return(c(original, rearranged)) +} +## Test the comparison function, which will return the value of the summand using the original long-form, and the algebraic rearrangement given in equation 5 +compare(1000, 998, 100, 20, .15, .001) +## Test for 1000 randomly chosen values of d, D0, Dw, DI, pp and pc +d.test = runif(n = 1000, 1, 10^9) +D0.test = sapply(d.test, function(dd) rpois(n = 1, dd)) +Dw.test = sapply(D0.test, function(D0) runif(1, 1, D0)) +DI.test = sapply(Dw.test, function(Dw) runif(1, 1, Dw)) +pp.test = runif(1000, 0, 1) +pc.test = runif(1000, 0, 1) +## Output comparison +test.outputs = mapply(FUN = compare, d= d.test, D0 = D0.test, Dw = Dw.test, DI = DI.test, pp = pp.test, pc = pc.test) +warnings() +## This script checks numerically that the probability of infection in the basic model can be written as: P(Infection|Dose) = 1-exp(D0*p_p*p_c) +## Long-form equation: +## Let d represent the expected inoculum dose +## Let D0 represent the actual inoculum size +## Let Dw represent the reaches the within-host environment +## Let DI represent the infecting dose +## P(DI|d) = P(D0|d)P(Dw|D0)P(DI|Dw) ## The probability of a given infectious dose, given the expected inoculum size +## P(Infection|d) = sum_DI=1:Inf sum_Dw=DI:Inf sum_D0=Dw:Inf P(D0|d)P(Dw|D0)P(DI|Dw)) ## Long form probability of infection, , given the expected inoculum size +## Calculate the summand after algebraic rearrangement into three different poisson densities +compare = function(d, D0, Dw, DI, pp, pc){ +## This is the long-form equation for P(DI|d), as given in equation 4 of the main text methods +original = dpois(x = D0, lambda = d)*dbinom(x = Dw, size = D0, prob = pc)*dbinom(x = DI, size = Dw, prob = pp) +## This is an alternate representation of the summand, as shown in equation 5 of the main text methods +rearranged = dpois(x = DI, lambda = d*pc*pp)*dpois(x = D0-Dw, lambda = d*(1-pc))*dpois(x = Dw-DI, lambda = d*pc*(1-pp)) +return(c(original, rearranged)) +} +## Test the comparison function, which will return the value of the summand using the original long-form, and the algebraic rearrangement given in equation 5 +compare(1000, 998, 100, 20, .15, .001) +## Test for 1000 randomly chosen values of d, D0, Dw, DI, pp and pc +d.test = round(runif(n = 1000, 1, 10^9)) +D0.test = sapply(d.test, function(dd) rpois(n = 1, dd)) +Dw.test = sapply(D0.test, function(D0) runif(1, 1, D0)) +DI.test = sapply(Dw.test, function(Dw) runif(1, 1, Dw)) +pp.test = runif(1000, 0, 1) +pc.test = runif(1000, 0, 1) +## Output comparison +test.outputs = mapply(FUN = compare, d= d.test, D0 = D0.test, Dw = Dw.test, DI = DI.test, pp = pp.test, pc = pc.test) +warnings() +d.test +D0.test +Dw.test +D0.test +## This script checks numerically that the probability of infection in the basic model can be written as: P(Infection|Dose) = 1-exp(D0*p_p*p_c) +## Long-form equation: +## Let d represent the expected inoculum dose +## Let D0 represent the actual inoculum size +## Let Dw represent the reaches the within-host environment +## Let DI represent the infecting dose +## P(DI|d) = P(D0|d)P(Dw|D0)P(DI|Dw) ## The probability of a given infectious dose, given the expected inoculum size +## P(Infection|d) = sum_DI=1:Inf sum_Dw=DI:Inf sum_D0=Dw:Inf P(D0|d)P(Dw|D0)P(DI|Dw)) ## Long form probability of infection, , given the expected inoculum size +## Calculate the summand after algebraic rearrangement into three different poisson densities +compare = function(d, D0, Dw, DI, pp, pc){ +## This is the long-form equation for P(DI|d), as given in equation 4 of the main text methods +original = dpois(x = D0, lambda = d)*dbinom(x = Dw, size = D0, prob = pc)*dbinom(x = DI, size = Dw, prob = pp) +## This is an alternate representation of the summand, as shown in equation 5 of the main text methods +rearranged = dpois(x = DI, lambda = d*pc*pp)*dpois(x = D0-Dw, lambda = d*(1-pc))*dpois(x = Dw-DI, lambda = d*pc*(1-pp)) +return(c(original, rearranged)) +} +## Test the comparison function, which will return the value of the summand using the original long-form, and the algebraic rearrangement given in equation 5 +compare(1000, 998, 100, 20, .15, .001) +## Test for 1000 randomly chosen values of d, D0, Dw, DI, pp and pc +d.test = round(runif(n = 1000, 1, 10^9)) +D0.test = sapply(d.test, function(dd) rpois(n = 1, dd)) +Dw.test = round(sapply(D0.test, function(D0) runif(1, 1, D0))) +DI.test = round(sapply(Dw.test, function(Dw) runif(1, 1, Dw))) +pp.test = runif(1000, 0, 1) +pc.test = runif(1000, 0, 1) +## Output comparison +test.outputs = mapply(FUN = compare, d= d.test, D0 = D0.test, Dw = Dw.test, DI = DI.test, pp = pp.test, pc = pc.test) +test.outputs +test.outputs[1,]-test.outputs[2,] +## This script checks numerically that the probability of infection in the basic model can be written as: P(Infection|Dose) = 1-exp(D0*p_p*p_c) +## Long-form equation: +## Let d represent the expected inoculum dose +## Let D0 represent the actual inoculum size +## Let Dw represent the reaches the within-host environment +## Let DI represent the infecting dose +## P(DI|d) = P(D0|d)P(Dw|D0)P(DI|Dw) ## The probability of a given infectious dose, given the expected inoculum size +## P(Infection|d) = sum_DI=1:Inf sum_Dw=DI:Inf sum_D0=Dw:Inf P(D0|d)P(Dw|D0)P(DI|Dw)) ## Long form probability of infection, , given the expected inoculum size +## Calculate the summand after algebraic rearrangement into three different poisson densities +compare = function(d, D0, Dw, DI, pp, pc){ +## This is the long-form equation for P(DI|d), as given in equation 4 of the main text methods +original = dpois(x = D0, lambda = d)*dbinom(x = Dw, size = D0, prob = pc)*dbinom(x = DI, size = Dw, prob = pp) +## This is an alternate representation of the summand, as shown in equation 5 of the main text methods +rearranged = dpois(x = DI, lambda = d*pc*pp)*dpois(x = D0-Dw, lambda = d*(1-pc))*dpois(x = Dw-DI, lambda = d*pc*(1-pp)) +return(c(original, rearranged)) +} +## Test the comparison function, which will return the value of the summand using the original long-form, and the algebraic rearrangement given in equation 5 +compare(1000, 998, 100, 20, .15, .001) +## Test for 1000 randomly chosen values of d, D0, Dw, DI, pp and pc +d.test = round(runif(n = 1000, 1, 10^9)) +D0.test = sapply(d.test, function(dd) rpois(n = 1, dd)) +Dw.test = round(sapply(D0.test, function(D0) runif(1, 1, D0))) +DI.test = round(sapply(Dw.test, function(Dw) runif(1, 1, Dw))) +pp.test = runif(1000, 0, 1) +pc.test = runif(1000, 0, 1) +## Output comparison +test.outputs = mapply(FUN = compare, d= d.test, D0 = D0.test, Dw = Dw.test, DI = DI.test, pp = pp.test, pc = pc.test) +## Verify that both versions of the summand give exactly the same answer +test[1,]-test[2,] +## Q.E.D. +## This script checks numerically that the probability of infection in the basic model can be written as: P(Infection|Dose) = 1-exp(D0*p_p*p_c) +## Long-form equation: +## Let d represent the expected inoculum dose +## Let D0 represent the actual inoculum size +## Let Dw represent the reaches the within-host environment +## Let DI represent the infecting dose +## P(DI|d) = P(D0|d)P(Dw|D0)P(DI|Dw) ## The probability of a given infectious dose, given the expected inoculum size +## P(Infection|d) = sum_DI=1:Inf sum_Dw=DI:Inf sum_D0=Dw:Inf P(D0|d)P(Dw|D0)P(DI|Dw)) ## Long form probability of infection, , given the expected inoculum size +## Calculate the summand after algebraic rearrangement into three different poisson densities +compare = function(d, D0, Dw, DI, pp, pc){ +## This is the long-form equation for P(DI|d), as given in equation 4 of the main text methods +original = dpois(x = D0, lambda = d)*dbinom(x = Dw, size = D0, prob = pc)*dbinom(x = DI, size = Dw, prob = pp) +## This is an alternate representation of the summand, as shown in equation 5 of the main text methods +rearranged = dpois(x = DI, lambda = d*pc*pp)*dpois(x = D0-Dw, lambda = d*(1-pc))*dpois(x = Dw-DI, lambda = d*pc*(1-pp)) +return(c(original, rearranged)) +} +## Test the comparison function, which will return the value of the summand using the original long-form, and the algebraic rearrangement given in equation 5 +compare(1000, 998, 100, 20, .15, .001) +## Test for 1000 randomly chosen values of d, D0, Dw, DI, pp and pc +d.test = round(runif(n = 1000, 1, 10^9)) +D0.test = sapply(d.test, function(dd) rpois(n = 1, dd)) +Dw.test = round(sapply(D0.test, function(D0) runif(1, 1, D0))) +DI.test = round(sapply(Dw.test, function(Dw) runif(1, 1, Dw))) +pp.test = runif(1000, 0, 1) +pc.test = runif(1000, 0, 1) +## Output comparison +test.outputs = mapply(FUN = compare, d= d.test, D0 = D0.test, Dw = Dw.test, DI = DI.test, pp = pp.test, pc = pc.test) +## Verify that both versions of the summand give exactly the same answer +test.outputs[1,]-test.outputs[2,] +## Q.E.D. diff --git a/Full_analysis_script_revised_submission.R b/Full_analysis_script_revised_submission.R index 87d88b1..7f4353f 100644 --- a/Full_analysis_script_revised_submission.R +++ b/Full_analysis_script_revised_submission.R @@ -18,6 +18,7 @@ setwd('~/Dropbox/R/Hamsters/') # Set working directory library(parallel) library(ggplot2) library(binom) +library(RColorBrewer) library(gridExtra) library(reshape2) @@ -641,7 +642,7 @@ ggsave(filename = figS1, plot = outplot, width = 5, height = 7, dpi = 'print', d AICcalc = function(n.pars, nll){ 2*n.pars+2*nll } -AIC.basic.ham = AICcalc(n.pars = 2, nll = opt.pc.abraded.skin$value) +AIC.basic.ham = AICcalc(n.pars = 2, nll = opt.pc.abraded.skin$objective) AIC.mixture.ham = AICcalc(n.pars = 3, nll = opt.aa.BB.abraded.skin$value) del.AIC.hamster = c(AIC.basic.ham, AIC.mixture.ham); del.AIC.hamster = del.AIC.hamster-min(del.AIC.hamster) del.AIC.hamster diff --git a/code_outputs/sessionInfo.txt b/code_outputs/sessionInfo.txt index 034c2c6..d094f90 100644 --- a/code_outputs/sessionInfo.txt +++ b/code_outputs/sessionInfo.txt @@ -10,20 +10,22 @@ locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: -[1] parallel stats graphics grDevices utils datasets methods base +[1] parallel stats graphics grDevices utils datasets +[7] methods base other attached packages: - [1] gridExtra_2.3 RColorBrewer_1.1-2 TailRank_3.2.1 oompaBase_3.2.6 binom_1.1-1 forcats_0.3.0 stringr_1.3.1 - [8] dplyr_0.7.8 purrr_0.2.5 readr_1.3.1 tidyr_0.8.2 tibble_2.0.1 tidyverse_1.2.1 viridis_0.5.1 -[15] viridisLite_0.3.0 bindrcpp_0.2.2 ggpubr_0.2 magrittr_1.5 reshape2_1.4.3 lubridate_1.7.4 reshape_0.8.8 -[22] ggplot2_3.1.0 +[1] RColorBrewer_1.1-2 reshape2_1.4.3 gridExtra_2.3 +[4] binom_1.1-1 ggplot2_3.1.0 loaded via a namespace (and not attached): - [1] tidyselect_0.2.5 haven_2.0.0 lattice_0.20-38 oompaData_3.1.1 colorspace_1.4-0 generics_0.0.2 - [7] yaml_2.2.0 utf8_1.1.4 rlang_0.3.1 pillar_1.3.1 glue_1.3.0 withr_2.1.2 -[13] BiocGenerics_0.26.0 modelr_0.1.2 readxl_1.2.0 audio_0.1-5.1 bindr_0.1.1 plyr_1.8.4 -[19] munsell_0.5.0 gtable_0.2.0 cellranger_1.1.0 rvest_0.3.2 Biobase_2.40.0 labeling_0.3 -[25] fansi_0.4.0 broom_0.5.1 Rcpp_1.0.0 scales_1.0.0 backports_1.1.3 jsonlite_1.6 -[31] hms_0.4.2 digest_0.6.18 stringi_1.2.4 grid_3.5.0 cli_1.0.1 tools_3.5.0 -[37] beepr_1.3 lazyeval_0.2.1 cluster_2.0.7-1 crayon_1.3.4 pkgconfig_2.0.2 xml2_1.2.0 -[43] assertthat_0.2.0 httr_1.4.0 rstudioapi_0.9.0 R6_2.3.0 nlme_3.1-137 compiler_3.5.0 + [1] Rcpp_1.0.0 rstudioapi_0.9.0 bindr_0.1.1 + [4] magrittr_1.5 tidyselect_0.2.5 munsell_0.5.0 + [7] colorspace_1.4-0 R6_2.3.0 rlang_0.3.1 +[10] stringr_1.3.1 plyr_1.8.4 dplyr_0.7.8 +[13] tools_3.5.0 grid_3.5.0 gtable_0.2.0 +[16] withr_2.1.2 digest_0.6.18 yaml_2.2.0 +[19] lazyeval_0.2.1 assertthat_0.2.0 tibble_2.0.1 +[22] crayon_1.3.4 bindrcpp_0.2.2 purrr_0.2.5 +[25] glue_1.3.0 labeling_0.3 stringi_1.2.4 +[28] compiler_3.5.0 pillar_1.3.1 scales_1.0.0 +[31] pkgconfig_2.0.2 diff --git a/manuscript_figures/Fig2_par_ests.pdf b/manuscript_figures/Fig2_par_ests.pdf index 6c2c5105f8ed2d524061a4a62d9d29c06e696776..b9b597883977eeaee8c6d12760924756124e243c 100644 GIT binary patch delta 48 wcmbQ4Gc#v`nW~9_p`n4Hk)g3Bm%eX)ic4Zis)B}#m63s=r2$-SWAr{#085b#F#rGn delta 48 wcmbQ4Gc#v`nX0jorJ=E*fr+Umm%eX)ic4Zis)B}#m63s=r2$-SWAr{#08UK}PXGV_ diff --git a/manuscript_figures/Fig3_model_fits.pdf b/manuscript_figures/Fig3_model_fits.pdf index e615003f14452983fa117823316c3d9e2f11ef17..7760fa58ea1cd51eff0d89f3b25a2851dfb001cc 100644 GIT binary patch delta 50 ycmdnK z5h$y!Y|RswEO-~6k#0x{Ild%+oJ%eck4zisD$-oS!I^1~j@Mx{^wVhT$&e`ua^G_r zh0ea?J^xtd=m8U^bkn-z!r$X_K zS-+>&TYb7U4qiYVyj` zUq5L%13TTkJO9dg&zg)+(d0~Lb?Pev;g7v-h|e7AIOjsD&=?0#Q6wms6{F4+%#C)C zSxks%jhU`YhwP&u@G6stLpqc>uUb9OopV4D?mr+Nq4Xj9caPCUcxs{?x){CmJefYrlY z2FFq>qXi_krlQjG+WY5*M~q~XTd8sZmQ93bN%7mQn1`b;QGR&{<-t0 z*Bh@X+QEx5^~}c%#LwJfRQ5k*Y7bqCVFy}3e(TjW(1FPt3U#pDq{;ksE#P%4w_WIl zD=YCssmWpJ8Taidjo~^57&i29Tg;~y4%lF1IX}M1r;s6Z6wQ!QsB4=spKBX8eD~Q> z3f{Ii$T|3;w42PLv+h&fn^n3eWgW%Y2A0Hig0{O%Ma4FT?DV%KtudID#iy?C%lfczldVGI{WKFWaMH;(c@*^pGCdv$@+px#e{V&OM>T zn}i&0MixWT>cuK-Ly)JGSul#?MBeG0UXQQPc%*qcc5q_dqSS&Car9J8g2Dv^S9dKy zI)!?xCuoQt*45oqp*CU{kXYRZHYFPt9N!G@Q9JkCE9D5*2Z)DR_5+djb$8B zvw~Q>_vWJJs$T!eA8k<;f#mE2!#n-LVu)|*+M=p}KypsPz@7e;d;O4H?TGNNm!UDm z2@d)4;UI<&5{M*rs8oza|2gUV{jcUHKPONjEd^-0N*)T(pz)s*E)~yxPKc6p*X@7R zj@G6<4VKfU^)i9=pDX-xdd47v0>#?2{Day$VekZsnE<=}orePKi_Vd(&>iJf`EWRO z`*k2A433R#kzwF&GJda3i_Z!xpvycz_U_m&9Mm%MkG)h+M1qRtB%{EX<&l*@GC(EZ z%icUsCX)z3!7bZ=x^N6iVjC^nHC>~7=uYpd0J@Rap;NpuUV3s8UI{x14#*TWj1j-K zb)UKCA9cr;Iw|e@eL^gItV4m})7wWt;!d_h)`n9{aiFBXMu#TBe$B`Q3PxxwYKqCt z-*OC}9mddrpRQdWb*kJ|S$;NzH2Ikd0(U&CS!HN2u0c{SD^wuIy0Q9^%h!K{yPbP? zuFeC09Q=`nrw*iAl{R%g@B_?5?@kcjgZ7IOwe_-0^{yD6^d3134gywBGl9c^Iq(bM z*J__eT&9NDPEZ+6__R~Zf@zg;;Buo&O2 zKx5G8N4bJuRqS&4fuL|h^VL7cM}w-~SL`T9uR2B|>AK#~P#v*1EGPX%KAKFp${3ub zVYHmP*2)X<=xS!&a%$oIj68haN_c z06^{#BJg(;0=j7$fl5r;v`&&Bq*Of^wzcA2ogI9gV?oWeMk|5B1AHz30g)u(r3*>U0J&=^Q+ zG;O15QNqx2Tz2Lfgtp-h-`sjs4_D;TuYjwsj}Gl4@#j5IHkyGfzH zPxQPgMcKPbbKsy$M3MTP06eK1Ic!IH2j%9eBSI{&5fKaYcho%?e$>As#cFAPy3J$z#vmo0yFxfKZb^*{dAa$TXA>Zm}ucmh&q8?W^iz~P#~^f z$pyN$&s*2!5E-HqNFUC697`V7(8rvE?GDAT(!Dx#>oYCRSHW0vqC}oRw#WrhqZbux znnEp%G_3mqBUFqX!f;elIwele#b_?_k+XF&elJur+a{8n2jtObSKHEQC#W#hOzdRv8Q|U7 zKMHs}^AYz^`>a8S5p;^h@WHdX!d%D1Ad%)7m&GtbPK>=8=Va4gSH|DMmB}sXK-lk= z$E&14IeM~hpfmtCr83pGI+YrC6ta&y){uy?>dQO%&#Y#(;CyOGujGxEW0xDrX5c%! z6GwsE?K(8_2d+9o@=SxCHnXlLFC%DOHwq(QuGu+Yo{v>9ex71)OAmRL{q!)`Q?&>9 zpP;@|kAj*m9`h(O4fMA)tgcqtx3DOP`VTSoPTpiPyD0D`k&CHYms<)Aw_^#dMf7834?8DnzgngK z)E#G9uxAHguYz){zq2hk+R#xxMMSZH&t^dLrUmoC`A=(~Ad#QcJIa6PGTy&)o4Jmz zVd_oeCtSYU?gF8d=W{Gdi4BIm%$7rSNAYGE&&yv;xx>h);ZHbs8_x4M_=gL=O&;Rw zxL2VLf;|Je&TnfZGIZcL^&gq4_t}P?J9Q8<;3+_y4Z}~SIMIOn6<7W)HUF-1%-Q$P zVrD+h;7BN!vh)FyqpX+Z##Jn0p>k3+#h6p15UXD$S@`kPO3L=`~)>~jJ_>)2hhR(%AnAtX51St zTAA8&F4?pE#hK~JHv3LfSHvwhRAGE zSYBL~&1WK>jh2#1gc8N4m;GgNHf%UX{D5)?CrFcX&yeM1poWXg`@13U)|wJO{_yd5Cfu50mWj2l|np!Vv`W> z_r?yJeP$fWT?u4ggjnBroWe*~0Cbewih>1WYi!w0&1n6n)5Ut3@dTi$=feJ!wNDAk3s=qh1`BINkn^XH^rzcMFVaA-Pc*E-T-AL=&G-H%A!@mKeW)Hk2P^=_Z@# z<``R$bkj`yNV*v&DUG^WCbhE$JtW;6la-*urP&dYO0p5Hzq90mf==pJ6k}F&HH#o1 zpUyPA%P_PWnm!P==7oze^i`;SdL^q{AxznIGoe@YqYC!go|bRJcM}KNF}^kSjQX&( zYuLP#j=aVmv;I64D^FQ@^!t*;rd|uSM=qP%zMNTPj?ozrx0A;OeY?A>3=i8?2bFT; z)$D>ac98_q=#By_l>={oIw~^>>je z4$}F*JjMnOX`=ZOZA*QEzVjsIsyfH`bXYcJ?)>dwNqA!n#19XRLdG`IyeYu)3^E{;l?dHNc zkK-NgiuC#MvTx}|^R2Vr>u_yh=stBel`B%T!NIE`vQy?O982UE0Fg-Q!IRhG@^T35 ztDH5hzTW^s+g4o&iqu1D#a(sK9umXKz=FSIGc|58RG_}Ac_SbDVk|b;o748Zf(^5O zHhY?u7usuh$I<&@vT(1x##(SZ2<`J_6#w!Jf5_3)qOx_m7tg&;=YaUFD$?ti*&V_FES&js7MdlCr zjvld`KbGqr%&(Tb7ym5qbaOd69X!Hc>`%CAPYDGIlP^^O$Q{3E^1%wQP*JTOw127` zyS;d5bcytO^YG7UGm1i<9gMG74UqxnkX@&hd15=Kw#dUpiQ2B%w!gw-fzA@MzXRGM zRTUb>J?jIu&|9|1SE|4Twd0ljGUnysW4xlI6-}!xB8!jJ`H@Fzsj}b0EY2A4=v9

^+peJMuXT|*NPKzI*GprW$erLg4gK7or!M~7C$XDAq zi8GmI^QQ>M!_9I;V|@ur%1kfL)p4_ONNg>5T`F8}R0Y}WKO)muI5#!?tmX`N(Th`9 zKKi`z$Q|2h^YwzoLH0smBLGbm>MjWX$fxY{0@0JYttdNOS%a-=VoP26_~s?;XXQtS z+++1N_4Z!<)0ieny&kIg5(ulDMuOy$Z~>Fy;XCi1tCe4REtDgC4D^KbH}bKzQc*u? zc$|Mk>7;&dL7Zoqj_K@3GJ|A5nb>_{UG`W@IF6e!RAO0H)>#*{3Ak509xgdUoS}nm zGSOC{>uxg6NYIzE)@)hmktNZw(`mbE@;dkR(CNn&h!&7XbSoQ(YlbseWZB8fZX-2j zwnxf~QTi;MRlecNfe-BKLetbB^~}PPx^%fLT5xo0601VPjf2~;@rC*`AFELEnU9~? zT1-W!W7Py75CsQjUR?WAXcE^%c3D9^hy3|!YB=XDUH#%DeiN}W;vK_{3d9XN@*NF& z<;Ir}#T4@03chw%qQ5X_(LG%*0W&Qt zPfTK&@SpWmn0givGnYw3DJE?ybE=FgQ+l2x5_{tTO^%?(#78SQ2I&VMnV`ruFgs^fuGmQ}j)!%Xk{X`sR zN!y$Q&exX^0fgUzL&63;#`~kv)jE5l!{a@lBHy^bbyh#W4Jli@m03NID&Eu&2;bT> z{VMVT_@Of@>yLRF`_Kut7qc-_%yW~?Q~rwVd1>+YJnOj)2&)gHU*1RS()ZegWg6sL z9GTut)wZfvDZ4BgYks6dzq;xtpK^cP{fpoOHtD-R1;Z)W52MB#AzNDn6oI&NaD2G4 zwSL@vI1x+PJ0-S9{3mZN<6AyPv=ZIFEY z!%+2{%&+d@bBfM#a#^2mOebe@5B;zOh>f@yr>B$e)zhrux(!^$BzRcn0#FP;IG{GRO_%m zyQFEM$yUbC)1Mg)WSpMr(Diyequb_jdP?Jyil2G~a;jIn)|vU4k!mdo4@^P)RFX#$ zR<(_F(inUht)zPt5slk(MXL>1mJdU($0VPm_+K0@1Kpj<4NNr)J@9M_-}UHj=GgLY zv2H^6$vEIdVz`gebaFlCmrV_>Nld?1`vrSrF&z0;%{9Bx;Y(Q(pTnoJM8_*~SNSti zFMj>Ca~J3Pf{T5OFAr+O$_(I{B!naXwx-Fed`Wk-)|sp;Hxp-aa@&enO4>YNdl#zo z0b8TxT|MyFG@FlFP5cgL4yZSCXq)q0v*}Fu5(O;2nN-?{aFfo=Qa$M{c6J*k%WZA7 ztQ-owX68m3Eb1_Ula3NZ%6FfIsYSY5&tQuZvm$2u zy7i9gt7D8eFLP)}ImP5t)rx)o}Kk#gAc}{e{10Va%+tExD z=5fZeX8AFf@`ByUwn|JlPt(5rMm_gx{yw|}ujLqcqglm<1{-lhyY|0I+KPh^cv)!an0iz{l5W3QBg(VznHSpf9I({6q57ARRPFA4=T1qyVKC5z#Y&(fIk8i%w*8LG(_L z=)Dh*+|OF?hv$CxhkgC`|GIvA@3qcZd!4mUnr@CR@)?^DL>w#x78HKQMIxZ#YO7%7 zW%ra*K>{KO78m?~Ih4u#9g%18lDMpZMm$}Ed-G6vAq&k&?Y)JzVV)k^X(CD6cnTS1 z3Acu7V;48&D52V$79mv1bLN5qht%J_zlzA2k0Z*8cq#V|?3jbSXX?^sf~K$E9q-Oq z)KQZKEx^>y47^s6rvs}YKw$m+wP%3uu~yy9{=meHtGVfR_%gPZ(P>8X>c$xe*iQo* z4+zn8N_{FF#pcs{=`6F<3U7JXcc%nx}- z5)Qxi68DZ9A@_N5qZzymlj07bi{5lZto-D4-Sj=z*2$Tx3Q-@^u;vKtPIz<3W<{ON z4w zQ<82E10;in9ZTa8S`0U;dt4mWuP-myKr@>&JXiZIe`Q@Q@OH6Lgqg7O&btC__m2pz z7$F`BS`>{ds{zcqIjru%kJ_4Te6Q*U;vyD+4*^e6TtZ6b2 z99(VVZBy7>buirx79V+sJuW}-`f=aMfztEi!4WO~PpX~y`>C^waL9VvW0}Nlfw~J_ zKO-yq&@>94*JqAag%>7Pt&@@=>us7_-m3N4Dbj9>`wy&OxndYKI=M8A zZt<10f74jEToQ)mf!DrM6VTxdz1s6b2x|C<%zu<+m1bN-=ja}G>1UElq4&@qQ(Zq- z&6cg(v}20SNxuF)f+74&L>BPgv2_u;%@m&BanTh)xB12JEC|;f=zw8kmrDCa5SY>_ z8b*EJ;^g`#%=rESRVQX@&l>{Eti250Y;Jpf9tFKp5x!<;I8+)3EQ|nqY{qf;xP>v` zge^VBj9btGF0;)Pc(1uj&P>jPUzsiW%z|Idjj}Lemb`S>*4H7%?LZeR>iI z^WIJe_nl6$Rnvzdbhq*KCHs5U93)?ik=>Z<-3Eq6ttbCn(wNWF#6>5|?ub~Q4@LAc zTK{STHUSyOub}j83Ne;JpRl3~39qk1*Douh&=F#hIn>(ZJ(%4$L{78?+k48(`QW{5uq6Na&|yay6goL_ zr1PXs*mdr`2WULTP}F7)>88q9z*G)qQJA`9@*n zqJt1vy?@e@{M`ytMo6&T8IppxGrZrDh9B}>g}M`B*X@=B=IpGF6X)!#A(f-qNl=%g zDeJte=W5i$o3LcmBbcz{*!z{gb2Gu+EQ}{*!p=_{qHTMIgx^jR-6S9s+o#Q=TN+&h zKu)8E0EPiu!;NpY5KYd`A!U_VzMfp0-brqU;*!Jz7feFKs~yL@+{xFf8z=;VRGmFX zPeTQf|20xd!mhBJRbN9Tm?ev5O0S%YYZ>_prIXJk#au8A8j+~lKn+ARG3z1}?6;A{ z*{tX@|w+F+VS12d+UIFYd6q{}@ z?4h_P6ATZgBg{xJ=tX4VKUhM9_Fs$b@g1;*L@uetLF2XJmrUaP;$w4`-YlDRtq<%4 zg0G*_0;p(N70gx4vC?Zf6L{&BneeqB!(npkXZ0`zXGzk;uUSm>AL9Gu+ozMjXmIDV z20Kw3DDwmKpoZrn_A=ZKu*nKgO(}nP=RPQCSXc16GMuTRZq02 zlWk4CeU#5s=Iz1h0ka6SM7&gg5-V2!hU8Jc2KqVOd?%}{I0M*`^8h#E3uy=p-`Y6C zit$4VGB;4&s8{a+pj!*OC^Wm>sgLHR*s$9&UhzIkzxUN-__&=nq5Jyc(nlVL^)4l~ zl<2cP86YuseP2@eIUD+!uv}inRr^?m>Vg9lsVd_HO;DXW&$8ln{GH{_9lENDcYM7p zSW3UnO`p+Tzo74xq-;Y1K3J!kHcjbvCmm7``CIucg=ikAZUm_nsKVDHVej}xR=lo<_$dwfPQtK*ANbt$~BMhpsqCKE*QXT7|_2y`_3 zb=!_$C$;Yba8vzk+Eh5Aao(`Ynby=13i*`Y^;k~nb;5OFG-d~dR;l5dWkEeSj(gQZ43&5y z8dGiApGi6z&o#6R!)K_2RgB`lSl{p z|MU!ECAXS!Qs9KlPO!pai@w*|vTjl`p{LEI(w8&shpJ*itfh`3-fpLPA`XByOfZOt z^6?MvYP$?@N@kYJpNlLl2uc;dMnqtLAIk7ze_v9VMM3grG*`KM6%zQk>1DAR2AXyk zEPB?);=XJK1+wyaN~*H^u$4baSSD_r(pjCJjX$IXF}l5siC*h0Cm3pMUCEXAo2#I$ zA@@(5z7~e{71W%l$>_X0ae$ct87u2g1K+V15EgenaMXiO4_>}>2)2kOh>p$?k(ddP zNxhD?np3S%w6^f@kr>>~axC__yC zAcCbDRjTaN;#IPyJCyn9t@7A~!)Xa4YuYlK$y_cDx=J7}hS*e;8JJ4#p0e#(V@1GN zXF6Z`+J1cUs6OgC_#AW`h8vyWm637Ugd2QJ{%LIFw|uVqh*}8{_p2-7T1AbOYC{vo zaZ%vT$_a~^;;XR?LoJnUL?UI2j0HDgiP5dv_w30pX*Xj_TC1-6FBGd?!L)H(@FnF5NLpoLTI=<2Ly_dk`x;5?*xHVGJ_PX zSz@(rtFP$b1Kv2+*_E14)rH3?KEYWFSZ@W>?Q3@r1elHY1T~;OnJ|NUq?B2cSl3li zx6h`8C#dpI$(M6urNKJo-xleqyIh`TIEODti#>>6U39|k*iXDX*XyMuVzjc&Qi~Cc zn-oDOiX49B7{n)LR19KV>TC1H)twS~1V8q^c}eqb5T(`PqG&{=Ta*|B(a0}jHnw}S z0qF)f0pM43>qtu15;~z(?FZvV|LVu@!l1){fu zo?8jN4fp`pXumQ~4NH={Cvf#dkTBbr;0+=UHA_&x5h5Le*LM6eZ%2@Jzo%u?&qZA{ zoyb#iRA+6dAvOVEEE^1EL$%Gx=1T$$p$6Z8gy#Xvi7rZ?`Nc}Tr0)c_#}QGB@p*|n zB*t4GGgna{O*aT&LxJ6lw&>x=;7#jt)m*=LD3cTjY~J$)uFtGQ-~%dF(9Sr@bu-@9 zhqI?eNTe24!FkjERKW$)IECQCY0$>xh-$b*T6gtM*Eh79irLTQHf;5X2n{Sf`G6KS@Y_c~%N!m|e-=00*tNPlOt%J# zsQt9|wD-?>@C(o2Y~%~irul5Bi-C^KE|tEDTUNup0W;{%FS8DXj*;JRPx?I!I*EXzhcz?A~7v{7b;^Q z499Ctk~PqJp|)@I3uzPpo_%_e5 z41v9u{_VoBtK4)uV0F32flkh{ z0bf?ep7CJK!vdgI6C_v*B$ky)G62|pJNc+}=3@3%ST2KHkzlmFTa?bHPVq>w# z!uK1}HnmcGS$NkaqdR}Di9@}JxI8gDtwZ|zJwrtH~nonM%5IY+cC z&TMXv0dqd5D`_*^0a6tQe>2#~xtmQF|9muWZuSQP>Y7`AXh!Yr`j*c(N$xkf9awg> zkkali>Pq1LM(;8%GOcGxZC@c`XCC@ovBpm3NT@u)D0aD19^n)TF{_%+zgT*xIN|9_ z{cWS+gphN{%RSDnoFSV{W;}Yi?yGJk7lE$|6`qjxUGK>odA^0DaL31PK?B=7oD4OM0eo0jJrKhe(<5*{asXqykkeo($eosJY*E{gC=SPQi zNCPWiu$=;>i*31OZ=t$6&EwEfCjq!o=%JH<9gK3YxH$p*s=nkT&QiCY_r&9e>e9E# zJ{p10O=|%{Jd(OD^TB$qS*ieW8&NQ^=|9vC|It54NALDdDGtZr7Ja3tyHs@x>HLEn z^d6N(HK(KW>-kPR6n(9Gfiv!Os(>+6H#U}#TtVjYcj2fK3oz2&otGUa<<4K_cO~fw zIan9_5v?i3z78R~ViOxuO0t*pqD+>cH%iUpwyhmC>dwjd>{JksRFZ6;;wPS{po=P* z{!96LSktZckfwWOo2bY~GeWI#tM1A6H%Z0a{*zc*6^_A|K7n)q@Vm16V);<(rFogm zjkc}#xl#uF_K3Mgb#ERkkv*3zMf(9o7h&YCdISTY$Z>i1Vfi_ocm_G);|Cjs4et@q z9aVUYQ-m1LT%mfh?syCtumr1ps_G+J#2=@Ar(jrte4AI&kbGNDgZEwcewg2gJWiCc z2w`vTh{%XM{DuXxqo>Y4=nq zXPG^FTuFH#c~2-wj%he6mpPQ#)su=Vfje(6oHDZSFO;*iYT9(n4Bn`m6-9X=nJjl& zdEVTN4$Q2kx0uEZ za-WvTbKAAFjM&q);*24m7898YnAQ3BgP2v8K7*to-UTC>uQ12+Uk9c4HJtxA$(hhz z%-@S~__I|bOJgXv8$4%S9FJkF zwx6Q$JZz%C*7NF(hns7K)|r7~P54Ig3@vUt43|Jf3|Y`qZ)8Rp|6L5V*xwJ1R5Cfi zx`Ob*@|7d?Y;1k{w_hcHG) zWA(lLdHUgdE0?vyP|4KNslTHs!~?JYID0=6%9=6;a%56ce!g`!LieddYri#b*rd=D z`Lu5S3nz<5mHnjp0i+j(_F+Zzu?T~I7a86wP;~Rfd92;t)??xVNN34)hW>FAnyD`(HH_C`?MYQS+_~Fhty!fnD3M7n4KG~(2hlbeodFH&SLuE3 zoV(i#zV|S5VctCB+lF;sJ%7wQKB-X)?cs$>$u?n_yn|Wa5+3Ao5yx>cb{yosS_`Y^ z_w)!UoJd$9Pzz3NXcr>1jQtM+W@GC$?%&2#r;^#L|r z&z@iDdyg5>U&MEYN?NPq7Q*nK*;X)x-WPwXY%Nb!H&{n;KsjGqAA7F2&T;6a<^@oZ zV02S%iu>rCmuPZ53#41BRHrN-7OCMPR^kdw+!G#G&yT;C=p1^)7ful5$MaC{QzWO5 zyY4x15d=?*Uk%IW9JyOD;%zyF_|06BGOI zy4!yt@t-&)0kQvEEiU#<G)w0908Fng9R* delta 48 wcmaE4_sDL7nX0jorJ=Ebp{c1Rm%eX)ic4Zis)B}#m63s=r2$-SV>G)w09Q2(x&QzG