Re-ran all code to check outputs and updated the numerical algebra ch…

…eck script
kgostic · Mar 29, 2019 · 2ea2f8e · 2ea2f8e
1 parent cb10b6b
commit 2ea2f8e
Show file tree

Hide file tree

Showing 7 changed files with 30 additions and 9 deletions.
diff --git a/Numerical_algebra_check.R b/Numerical_algebra_check.R
@@ -7,16 +7,36 @@
 ## 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(I|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
+##    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 ##
+  ## 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 ##
+  ## 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/code_outputs/sessionInfo.txt b/code_outputs/sessionInfo.txt
@@ -19,10 +19,11 @@ other attached packages:
 [22] 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      yaml_2.2.0         
- [8] utf8_1.1.4          rlang_0.3.1         pillar_1.3.1        glue_1.3.0          withr_2.1.2         BiocGenerics_0.26.0 modelr_0.1.2       
-[15] readxl_1.2.0        audio_0.1-5.1       bindr_0.1.1         plyr_1.8.4          munsell_0.5.0       gtable_0.2.0        cellranger_1.1.0   
-[22] rvest_0.3.2         Biobase_2.40.0      labeling_0.3        fansi_0.4.0         broom_0.5.1         Rcpp_1.0.0          scales_1.0.0       
-[29] backports_1.1.3     jsonlite_1.6        hms_0.4.2           digest_0.6.18       stringi_1.2.4       grid_3.5.0          cli_1.0.1          
-[36] tools_3.5.0         beepr_1.3           lazyeval_0.2.1      cluster_2.0.7-1     crayon_1.3.4        pkgconfig_2.0.2     xml2_1.2.0         
+ [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     
diff --git a/manuscript_figures/Fig2_par_ests.pdf b/manuscript_figures/Fig2_par_ests.pdf
diff --git a/manuscript_figures/Fig3_model_fits.pdf b/manuscript_figures/Fig3_model_fits.pdf
diff --git a/manuscript_figures/Fig4_prob_inf_given_dose_and_route.pdf b/manuscript_figures/Fig4_prob_inf_given_dose_and_route.pdf
diff --git a/scratch_figures/Profiles_basic_model.pdf b/scratch_figures/Profiles_basic_model.pdf
diff --git a/scratch_figures/pI_profiles.pdf b/scratch_figures/pI_profiles.pdf