update figures

chapter_11.qmd

@@ -250,4 +250,15 @@ ggplot(dat.weights, aes(x = x, y = y,
 ```
 
 
+## Version info {.unnumbered}
+This chapter was rendered using the following version of R and its packages:
+
+```{r}
+#| echo: false
+#| message: false
+#| warning: false
+sessionInfo()
+```
+
+## References {.unnumbered}
 

chapter_12.qmd

@@ -30,6 +30,8 @@ bibliography: 'https://api.citedrive.com/bib/0d25b38b-db8f-43c4-b934-f4e2f3bd655
 We first load the relevant R scripts:
 
 ```{r}
+#| message: false
+#| warning: false
 source("resources/chapter 12/sim.r")
 source("resources/chapter 12/fig_functions.r")
 source("resources/chapter 12/mlmi.r")
@@ -44,31 +46,31 @@ set.seed(9843626)
 
 dataset  <- sim_data_EDSS(npatients = 500,
                           ncenters = 10,
-                          follow_up = 12*5, # Total follow-up (number of months)
-                          sd_a_t = 0.5,   # DGM - Within-visit variation in EDSS scores
-                          baseline_EDSS = 1.3295,    # DGM - Mean baseline EDDS score
-                          sd_alpha_ij = 1.46,    # DGM - Between-subject variation in baseline EDSS
-                          sd_beta1_j = 0.20,    # DGM - Between-site variation in baseline EDSS
+                          follow_up = 12*5,
+                          sd_a_t = 0.5,  
+                          baseline_EDSS = 1.3295, 
+                          sd_alpha_ij = 1.46, 
+                          sd_beta1_j = 0.20,  
                           mean_age = 42.41,
                           sd_age = 10.53,
                           min_age = 18,
-                          beta_age = 0.05, # DGM - prognostic effect of age
-                          beta_t = 0.014,  # DGM - prognostic effect of time
-                          beta_t2 = 0,    # DGM - prognostic effect of time squared
-                          delta_xt = 0, # DGM - interaction treatment time
-                          delta_xt2 = 0, # 0.0005    # DGM - interaction treatment time2
+                          beta_age = 0.05, 
+                          beta_t = 0.014, 
+                          beta_t2 = 0,   
+                          delta_xt = 0,
+                          delta_xt2 = 0,
                           p_female = 0.75, 
-                          beta_female = -0.2 ,  ## DGM - prognostic effect of male sex
-                          delta_xf = 0,      ## DGM - interaction sex treatment       
-                          rho = 0.8,             # DGM - autocorrelation of between alpha_tij
-                          corFUN = corAR1,       # DGM - correlation structure of the latent EDSS scores
-                          tx_alloc_FUN = treatment_alloc_confounding_v2 ) ## or treatment_alloc_randomized
+                          beta_female = -0.2 ,
+                          delta_xf = 0,        
+                          rho = 0.8,           
+                          corFUN = corAR1,      
+                          tx_alloc_FUN = treatment_alloc_confounding_v2)
 ```
 
 ```{r}
 #| echo: FALSE
-#| fig.cap: Distribution of the EDSS score at each time point
-plot_distribution_edss(dataset)
+#| fig.cap: Longitudinal distribution of the EDSS score. The shaded area depicts the interquartile range.
+ggplot_distribution_edss(dataset)
 ```
 
 We remove the outcome `y` according to the informative visit process that depends on the received treatment, gender, and age.
@@ -100,8 +102,9 @@ fitps <- glm(x ~ age + sex + edss, family = 'binomial',
              data = origdat60)
 
 # Derive the propensity score
-origdat60 <- origdat60 %>% mutate(ipt = ifelse(x == 1, 1/predict(fitps, type = 'response'),
-                                               1/(1-predict(fitps, type = 'response'))))
+origdat60 <- origdat60 %>% 
+  mutate(ipt = ifelse(x == 1, 1/predict(fitps, type = 'response'),
+                      1/(1 - predict(fitps, type = 'response'))))
 
 # Estimate 
 fit_ref_m <- tidy(lm(y ~ x, weight = ipt, data = origdat60), conf.int = TRUE) 
@@ -125,9 +128,12 @@ results <- results %>% add_row(data.frame(method = "Reference"),
 We here implement inverse probability of response weights into the estimating equations to adjust for nonrandom missingness [@Coulombe_2022, @coulombe_weighted_2021].
 
 ```{r}
-obsdat60 <- dataset_visit %>% mutate(visit = ifelse(is.na(y_obs),0,1)) %>% filter(time == 60)
+obsdat60 <- dataset_visit %>% 
+  mutate(visit = ifelse(is.na(y_obs),0,1)) %>% 
+  filter(time == 60)
 
-gamma <- glm(visit ~ x + sex + age + edss, family = 'binomial', data = obsdat60)$coef   
+gamma <- glm(visit ~ x + sex + age + edss, 
+             family = 'binomial', data = obsdat60)$coef   
 
 obsdat60 <- obsdat60 %>% mutate(rho_i = 1/exp(gamma["(Intercept)"] +
                                                           gamma["x"]*x +
@@ -138,11 +144,13 @@ obsdat60 <- obsdat60 %>% mutate(rho_i = 1/exp(gamma["(Intercept)"] +
 fitps <- glm(x ~ age + sex + edss, family='binomial', data = obsdat60)
 
 # Derive the propensity score
-obsdat60 <- obsdat60 %>% mutate(ipt = ifelse(x==1, 1/predict(fitps, type='response'),
-                                            1/(1-predict(fitps, type='response'))))
+obsdat60 <- obsdat60 %>% 
+  mutate(ipt = ifelse(x==1, 1/predict(fitps, type='response'),
+                      1/(1 - predict(fitps, type='response'))))
 
 
-fit_w <- tidy(lm(y_obs ~ x, weights = ipt*rho_i, data = obsdat60), conf.int = TRUE)
+fit_w <- tidy(lm(y_obs ~ x, weights = ipt*rho_i, data = obsdat60), 
+              conf.int = TRUE)
 ```
 
 ```{r echo = F}
@@ -169,16 +177,18 @@ We can now estimate the treatment effect in each imputed dataset
 
 ```{r}
 # Predict probability of treatment allocation
-fitps <- glm(x ~ age + sex + edss, family='binomial', data = dataset_visit)
+fitps <- glm(x ~ age + sex + edss, family = 'binomial', data = dataset_visit)
   
 # Derive the propensity score
-dataset_visit <- dataset_visit %>% mutate(ipt = ifelse(x==1, 1/predict(fitps, type='response'),
-                                                       1/(1-predict(fitps, type='response'))))
+dataset_visit <- dataset_visit %>% 
+  mutate(ipt = ifelse(x == 1, 1/predict(fitps, type = 'response'),
+                      1/(1 - predict(fitps, type = 'response'))))
   
 Q <- U <- rep(NA, 10) # Error variances
 
 for (i in seq(10)) {
-  dati <- cbind(dataset_visit[,c("x","ipt","time")], y_imp = imp[,i]) %>% filter(time == 60)
+  dati <- cbind(dataset_visit[,c("x","ipt","time")], y_imp = imp[,i]) %>% 
+    filter(time == 60)
   
   # Estimate 
   fit <- tidy(lm(y_imp ~ x, weight = ipt, data = dati), conf.int = TRUE) 

chapter_16.qmd

@@ -148,7 +148,8 @@ summbenefit <- ds %>% group_by(study) %>%
 ggplot(ds, aes(x = benefit)) + 
   geom_histogram(aes(y = after_stat(density)), alpha = 0.3) + 
   geom_density() +
-  geom_vline(data = summbenefit, aes(xintercept = meanbenefit), linewidth = 0.5, lty = 2) + 
+  geom_vline(data = summbenefit, aes(xintercept = meanbenefit), 
+             linewidth = 0.5, lty = 2) + 
   facet_wrap(~study) + 
   ylab("Density") +
   xlab("Predicted treatment benefit")  + theme_bw()

docs/chapter_11.html

@@ -246,6 +246,8 @@ <h2 id="toc-title">Table of contents</h2>
   <li><a href="#individual-participant-data" id="toc-individual-participant-data" class="nav-link" data-scroll-target="#individual-participant-data"><span class="header-section-number">7.2.2</span> Individual participant data</a></li>
   <li><a href="#hierarchical-metaregression" id="toc-hierarchical-metaregression" class="nav-link" data-scroll-target="#hierarchical-metaregression"><span class="header-section-number">7.2.3</span> Hierarchical metaregression</a></li>
   </ul></li>
+  <li><a href="#version-info" id="toc-version-info" class="nav-link" data-scroll-target="#version-info">Version info</a></li>
+  <li><a href="#references" id="toc-references" class="nav-link" data-scroll-target="#references">References</a></li>
   </ul>
 </nav>
     </div>
@@ -604,6 +606,66 @@ <h3 data-number="7.2.3" class="anchored" data-anchor-id="hierarchical-metaregres
 </div>
 </div>
 </div>
+</section>
+</section>
+<section id="version-info" class="level2 unnumbered">
+<h2 class="unnumbered anchored" data-anchor-id="version-info">Version info</h2>
+<p>This chapter was rendered using the following version of R and its packages:</p>
+<div class="cell" data-hash="chapter_11_cache/html/unnamed-chunk-8_097e9efc1fec5df7c7a25885223b5cef">
+<div class="cell-output cell-output-stdout">
+<pre><code>R version 4.2.3 (2023-03-15 ucrt)
+Platform: x86_64-w64-mingw32/x64 (64-bit)
+Running under: Windows 10 x64 (build 19045)
+
+Matrix products: default
+
+locale:
+[1] LC_COLLATE=Dutch_Netherlands.utf8  LC_CTYPE=Dutch_Netherlands.utf8   
+[3] LC_MONETARY=Dutch_Netherlands.utf8 LC_NUMERIC=C                      
+[5] LC_TIME=Dutch_Netherlands.utf8    
+
+attached base packages:
+[1] stats     graphics  grDevices utils     datasets  methods   base     
+
+other attached packages:
+ [1] meta_6.5-0       table1_1.4.3     tableone_0.13.2  dplyr_1.1.2     
+ [5] jarbes_2.0.0     GGally_2.1.2     R2jags_0.7-1     rjags_4-14      
+ [9] mcmcplots_0.4.3  coda_0.19-4      gridExtra_2.3    ggplot2_3.4.4   
+[13] kableExtra_1.3.4
+
+loaded via a namespace (and not attached):
+ [1] httr_1.4.7          tidyr_1.3.0         sfsmisc_1.1-16     
+ [4] jsonlite_1.8.7      viridisLite_0.4.2   splines_4.2.3      
+ [7] Formula_1.2-5       shiny_1.8.0         metafor_4.4-0      
+[10] yaml_2.3.7          numDeriv_2016.8-1.1 R2WinBUGS_2.1-21   
+[13] pillar_1.9.0        lattice_0.21-8      glue_1.6.2         
+[16] digest_0.6.31       RColorBrewer_1.1-3  promises_1.2.1     
+[19] minqa_1.2.6         rvest_1.0.3         colorspace_2.1-0   
+[22] htmltools_0.5.5     httpuv_1.6.12       Matrix_1.5-4.1     
+[25] survey_4.2-1        plyr_1.8.8          pkgconfig_2.0.3    
+[28] purrr_1.0.1         xtable_1.8-4        scales_1.2.1       
+[31] webshot_0.5.5       svglite_2.1.1       later_1.3.1        
+[34] metadat_1.2-0       lme4_1.1-35.1       tibble_3.2.1       
+[37] generics_0.1.3      ellipsis_0.3.2      withr_2.5.2        
+[40] cli_3.6.1           survival_3.5-5      magrittr_2.0.3     
+[43] mime_0.12           evaluate_0.23       fansi_1.0.4        
+[46] nlme_3.1-162        MASS_7.3-60         xml2_1.3.4         
+[49] tools_4.2.3         mitools_2.4         lifecycle_1.0.4    
+[52] stringr_1.5.1       munsell_0.5.0       compiler_4.2.3     
+[55] systemfonts_1.0.4   rlang_1.1.1         nloptr_2.0.3       
+[58] grid_4.2.3          rstudioapi_0.15.0   CompQuadForm_1.4.3 
+[61] htmlwidgets_1.6.2   miniUI_0.1.1.1      rmarkdown_2.25     
+[64] boot_1.3-28.1       gtable_0.3.4        abind_1.4-5        
+[67] DBI_1.1.3           reshape_0.8.9       R6_2.5.1           
+[70] knitr_1.45          denstrip_1.5.4      fastmap_1.1.1      
+[73] utf8_1.2.3          mathjaxr_1.6-0      ggExtra_0.10.1     
+[76] stringi_1.7.12      parallel_4.2.3      Rcpp_1.0.10        
+[79] vctrs_0.6.3         tidyselect_1.2.0    xfun_0.39          </code></pre>
+</div>
+</div>
+</section>
+<section id="references" class="level2 unnumbered">
+<h2 class="unnumbered anchored" data-anchor-id="references">References</h2>
 
 
 <div id="refs" class="references csl-bib-body hanging-indent" role="list">
@@ -612,7 +674,6 @@ <h3 data-number="7.2.3" class="anchored" data-anchor-id="hierarchical-metaregres
 </div>
 </div>
 </section>
-</section>
 
 </main> <!-- /main -->
 <script id="quarto-html-after-body" type="application/javascript">