diff --git a/.nojekyll b/.nojekyll index fd2a5cf..6995fa3 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -20822a39 \ No newline at end of file +10f7f294 \ No newline at end of file diff --git a/chapter_10.html b/chapter_10.html index 2c6b188..5ffc48f 100644 --- a/chapter_10.html +++ b/chapter_10.html @@ -854,16 +854,16 @@

sucra.uc$ranking.random 
 ADA160/80     FIL100     FIL200 GOL200/100      INF10       INF5    OZA0.92 
-0.27000000 0.16363636 0.46727273 0.63727273 0.57272727 0.75636364 0.75363636 
+0.26090909 0.18909091 0.40363636 0.59545455 0.60909091 0.74090909 0.77727273 
        PBO      TOF10      UPA45       UST6     VED300 
-0.02090909 0.38909091 0.96727273 0.37000000 0.63181818
+0.02090909 0.40636364 0.98727273 0.40727273 0.60181818
-

These results indicate that 96.7% of the evaluated treatments are worse than UPA45.

+

These results indicate that 98.7% of the evaluated treatments are worse than UPA45.

diff --git a/chapter_11.html b/chapter_11.html index 7868212..d1c2406 100644 --- a/chapter_11.html +++ b/chapter_11.html @@ -246,6 +246,8 @@

Table of contents

  • 7.2.2 Individual participant data
  • 7.2.3 Hierarchical metaregression
    • +
  • Version info
    • +
  • References
    • @@ -604,6 +606,69 @@

    +

    Version info

    +

    This chapter was rendered using the following version of R and its packages:

    +
    +
    +
    R version 4.2.3 (2023-03-15)
+Platform: x86_64-pc-linux-gnu (64-bit)
+Running under: Ubuntu 22.04.3 LTS
+
+Matrix products: default
+BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
+LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
+
+locale:
+ [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
+ [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
+ [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
+[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
+
+attached base packages:
+[1] stats     graphics  grDevices utils     datasets  methods   base     
+
+other attached packages:
+[1] ggplot2_3.4.4    meta_6.5-0       jarbes_2.0.0     dplyr_1.1.4     
+[5] table1_1.4.3     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.20-45     glue_1.6.2         
+[16] digest_0.6.33       RColorBrewer_1.1-3  promises_1.2.1     
+[19] minqa_1.2.6         rvest_1.0.3         colorspace_2.1-0   
+[22] R2jags_0.7-1        Matrix_1.6-3        mcmcplots_0.4.3    
+[25] htmltools_0.5.7     httpuv_1.6.12       plyr_1.8.9         
+[28] pkgconfig_2.0.3     purrr_1.0.2         xtable_1.8-4       
+[31] scales_1.2.1        webshot_0.5.5       svglite_2.1.2      
+[34] rjags_4-14          later_1.3.1         ggstats_0.5.1      
+[37] metadat_1.2-0       lme4_1.1-35.1       tibble_3.2.1       
+[40] farver_2.1.1        generics_0.1.3      ellipsis_0.3.2     
+[43] withr_2.5.2         cli_3.6.1           magrittr_2.0.3     
+[46] mime_0.12           evaluate_0.23       GGally_2.2.0       
+[49] fansi_1.0.5         nlme_3.1-162        MASS_7.3-58.2      
+[52] xml2_1.3.5          tools_4.2.3         lifecycle_1.0.4    
+[55] stringr_1.5.1       munsell_0.5.0       compiler_4.2.3     
+[58] systemfonts_1.0.5   rlang_1.1.2         nloptr_2.0.3       
+[61] grid_4.2.3          rstudioapi_0.15.0   CompQuadForm_1.4.3 
+[64] htmlwidgets_1.6.3   miniUI_0.1.1.1      labeling_0.4.3     
+[67] rmarkdown_2.25      boot_1.3-28.1       gtable_0.3.4       
+[70] codetools_0.2-19    abind_1.4-5         R6_2.5.1           
+[73] gridExtra_2.3       knitr_1.45          denstrip_1.5.4     
+[76] fastmap_1.1.1       utf8_1.2.4          mathjaxr_1.6-0     
+[79] ggExtra_0.10.1      stringi_1.8.2       parallel_4.2.3     
+[82] Rcpp_1.0.11         vctrs_0.6.4         tidyselect_1.2.0   
+[85] xfun_0.41           coda_0.19-4
    +
    +
    +
    +
    +

    References

    @@ -612,7 +677,6 @@

    diff --git a/chapter_12.html b/chapter_12.html index 9120d86..6c092a8 100644 --- a/chapter_12.html +++ b/chapter_12.html @@ -297,101 +297,59 @@

    8 

    8.1 Introduction

    We first load the relevant R scripts:

    -
    -
    source("resources/chapter 12/sim.r")
    -
    -
    
-Attaching package: 'dplyr'
    -
    -
    -
    The following objects are masked from 'package:stats':
-
-    filter, lag
    -
    -
    -
    The following objects are masked from 'package:base':
-
-    intersect, setdiff, setequal, union
    -
    -
    -
    
-Attaching package: 'mice'
    -
    -
    -
    The following object is masked from 'package:stats':
-
-    filter
    -
    -
    -
    The following objects are masked from 'package:base':
-
-    cbind, rbind
    -
    -
    -
    
-Attaching package: 'nlme'
    -
    -
    -
    The following object is masked from 'package:dplyr':
-
-    collapse
    -
    -
    -
    
-Attaching package: 'MASS'
    -
    -
    -
    The following object is masked from 'package:dplyr':
-
-    select
    -
    -
    source("resources/chapter 12/fig_functions.r")
-source("resources/chapter 12/mlmi.r")
    +
    +
    source("resources/chapter 12/sim.r")
+source("resources/chapter 12/fig_functions.r")
+source("resources/chapter 12/mlmi.r")

    8.2 Example dataset

    Below, we generate an example dataset that contains information on the treatment allocation x and three baseline covariates age, sex and edss (EDSS at treatment start). The discrete outcome y represents the Expanded Disability Status Scale (EDSS) score after time months of treatment exposure. Briefly, the EDSS is a semi-continuous measure that varies from 0 (no disability) to 10 (death).

    -
    -
    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
-                          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
-                          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
    +
    +
    set.seed(9843626)
+
+dataset  <- sim_data_EDSS(npatients = 500,
+                          ncenters = 10,
+                          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, 
+                          beta_t = 0.014, 
+                          beta_t2 = 0,   
+                          delta_xt = 0,
+                          delta_xt2 = 0,
+                          p_female = 0.75, 
+                          beta_female = -0.2 ,
+                          delta_xf = 0,        
+                          rho = 0.8,           
+                          corFUN = corAR1,      
+                          tx_alloc_FUN = treatment_alloc_confounding_v2)
    +
    +
    +
    +
    `summarise()` has grouped output by 'time'. You can override using the
+`.groups` argument.
    -

    -
    Distribution of the EDSS score at each time point
    +
    Longitudinal distribution of the EDSS score. The shaded area depicts the interquartile range.

    We remove the outcome y according to the informative visit process that depends on the received treatment, gender, and age.

    -
    dataset_visit <- censor_visits_a5(dataset, seed = 12345) %>% 
-  dplyr::select(-y) %>%
-  mutate(time_x = time*x)
    +
    dataset_visit <- censor_visits_a5(dataset, seed = 12345) %>% 
+  dplyr::select(-y) %>%
+  mutate(time_x = time*x)

    In the censored data, a total of 17 out of 5000 patients have a visit at time=60.

    @@ -400,73 +358,81 @@

    8.3.1 Original data

    -
    -
    origdat60 <- dataset %>% filter(time == 60)
-
-# Predict probability of treatment allocation
-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'))))
-
-# Estimate 
-fit_ref_m <- tidy(lm(y ~ x, weight = ipt, data = origdat60), conf.int = TRUE)
    +
    +
    origdat60 <- dataset %>% filter(time == 60)
+
+# Predict probability of treatment allocation
+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'))))
+
+# Estimate 
+fit_ref_m <- tidy(lm(y ~ x, weight = ipt, data = origdat60), conf.int = TRUE)

    8.3.2 Doubly-weighted marginal treatment effect

    We here implement inverse probability of response weights into the estimating equations to adjust for nonrandom missingness Coulombe, Moodie, and Platt (2020).

    -
    -
    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   
-
-obsdat60 <- obsdat60 %>% mutate(rho_i = 1/exp(gamma["(Intercept)"] +
-                                                          gamma["x"]*x +
-                                                          gamma["sex"]*sex +
-                                                          gamma["age"]*age))
-
-# Predict probability of treatment allocation
-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'))))
-
-
-fit_w <- tidy(lm(y_obs ~ x, weights = ipt*rho_i, data = obsdat60), conf.int = TRUE)
    +
    +
    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   
+
+obsdat60 <- obsdat60 %>% mutate(rho_i = 1/exp(gamma["(Intercept)"] +
+                                                          gamma["x"]*x +
+                                                          gamma["sex"]*sex +
+                                                          gamma["age"]*age))
+
+# Predict probability of treatment allocation
+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'))))
+
+
+fit_w <- tidy(lm(y_obs ~ x, weights = ipt*rho_i, data = obsdat60), 
+              conf.int = TRUE)

    8.3.3 Multilevel multiple imputation

    We adopt the imputation approach proposed by Debray et al. (2023). Briefly, we impute the entire vector of y_obs for all 61 potential visits and generate 10 imputed datasets. Note: mlmi currently does not support imputation of treatment-covariate interaction terms.

    -
    imp <- impute_y_mice_3l(dataset_visit, seed = 12345)
    +
    imp <- impute_y_mice_3l(dataset_visit, seed = 12345)

    We can now estimate the treatment effect in each imputed dataset

    -
    -
    # Predict probability of treatment allocation
-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'))))
-  
-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)
-  
-  # Estimate 
-  fit <- tidy(lm(y_imp ~ x, weight = ipt, data = dati), conf.int = TRUE) 
-  
-  Q[i] <- fit %>% filter(term == "x") %>% pull(estimate)
-  U[i] <- (fit %>% filter(term == "x") %>% pull(std.error))**2
-}
-
-fit_mlmi <- pool.scalar(Q = Q, U = U)
    +
    +
    # Predict probability of treatment allocation
+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'))))
+  
+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)
+  
+  # Estimate 
+  fit <- tidy(lm(y_imp ~ x, weight = ipt, data = dati), conf.int = TRUE) 
+  
+  Q[i] <- fit %>% filter(term == "x") %>% pull(estimate)
+  U[i] <- (fit %>% filter(term == "x") %>% pull(std.error))**2
+}
+
+fit_mlmi <- pool.scalar(Q = Q, U = U)
    @@ -475,18 +441,18 @@

    8.4.1 Doubly -weighted marginal treatment effect total

    -
    obsdatall <- dataset_visit %>% mutate(visit = ifelse(is.na(y_obs),0,1))  
-gamma <- glm(visit ~ x + sex + age + edss, family = 'binomial', data = obsdatall)$coef   
-obsdatall <- obsdatall %>% mutate(rho_i = 1/exp(gamma["(Intercept)"] +
-                                                gamma["x"]*x +
-                                                gamma["sex"]*sex +
-                                                gamma["age"]*age))
-# Predict probability of treatment allocation
-fitps <- glm(x ~ age + sex + edss, family='binomial', data = obsdatall)
-# Derive the propensity score
-obsdatall <- obsdatall %>% 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 = obsdatall), conf.int = TRUE)
    +
    obsdatall <- dataset_visit %>% mutate(visit = ifelse(is.na(y_obs),0,1))  
+gamma <- glm(visit ~ x + sex + age + edss, family = 'binomial', data = obsdatall)$coef   
+obsdatall <- obsdatall %>% mutate(rho_i = 1/exp(gamma["(Intercept)"] +
+                                                gamma["x"]*x +
+                                                gamma["sex"]*sex +
+                                                gamma["age"]*age))
+# Predict probability of treatment allocation
+fitps <- glm(x ~ age + sex + edss, family='binomial', data = obsdatall)
+# Derive the propensity score
+obsdatall <- obsdatall %>% 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 = obsdatall), conf.int = TRUE)
    @@ -525,21 +491,21 @@

    Version info

    [5] mice_3.16.0 ggplot2_3.4.4 broom_1.0.5 dplyr_1.1.4 loaded via a namespace (and not attached): - [1] shape_1.4.6 tidyselect_1.2.0 xfun_0.41 purrr_1.0.2 - [5] splines_4.2.3 lattice_0.20-45 colorspace_2.1-0 vctrs_0.6.4 - [9] generics_0.1.3 htmltools_0.5.7 yaml_2.3.7 pan_1.9 -[13] utf8_1.2.4 survival_3.5-3 rlang_1.1.2 jomo_2.7-6 -[17] pillar_1.9.0 nloptr_2.0.3 withr_2.5.2 glue_1.6.2 -[21] RColorBrewer_1.1-3 foreach_1.5.2 lifecycle_1.0.4 munsell_0.5.0 -[25] gtable_0.3.4 htmlwidgets_1.6.3 codetools_0.2-19 evaluate_0.23 -[29] labeling_0.4.3 knitr_1.45 fastmap_1.1.1 fansi_1.0.5 -[33] Rcpp_1.0.11 scales_1.2.1 backports_1.4.1 jsonlite_1.8.7 -[37] farver_2.1.1 lme4_1.1-35.1 digest_0.6.33 grid_4.2.3 -[41] cli_3.6.1 tools_4.2.3 magrittr_2.0.3 glmnet_4.1-8 -[45] tibble_3.2.1 tidyr_1.3.0 pkgconfig_2.0.3 ellipsis_0.3.2 -[49] Matrix_1.6-3 minqa_1.2.6 rmarkdown_2.25 iterators_1.0.14 -[53] mitml_0.4-5 R6_2.5.1 boot_1.3-28.1 rpart_4.1.19 -[57] nnet_7.3-18 compiler_4.2.3 + [1] shape_1.4.6 tidyselect_1.2.0 xfun_0.41 purrr_1.0.2 + [5] splines_4.2.3 lattice_0.20-45 colorspace_2.1-0 vctrs_0.6.4 + [9] generics_0.1.3 htmltools_0.5.7 yaml_2.3.7 pan_1.9 +[13] utf8_1.2.4 survival_3.5-3 rlang_1.1.2 jomo_2.7-6 +[17] pillar_1.9.0 nloptr_2.0.3 withr_2.5.2 glue_1.6.2 +[21] foreach_1.5.2 lifecycle_1.0.4 munsell_0.5.0 gtable_0.3.4 +[25] htmlwidgets_1.6.3 codetools_0.2-19 evaluate_0.23 labeling_0.4.3 +[29] knitr_1.45 fastmap_1.1.1 fansi_1.0.5 Rcpp_1.0.11 +[33] scales_1.2.1 backports_1.4.1 jsonlite_1.8.7 farver_2.1.1 +[37] lme4_1.1-35.1 digest_0.6.33 grid_4.2.3 cli_3.6.1 +[41] tools_4.2.3 magrittr_2.0.3 glmnet_4.1-8 tibble_3.2.1 +[45] tidyr_1.3.0 pkgconfig_2.0.3 ellipsis_0.3.2 Matrix_1.6-3 +[49] minqa_1.2.6 rmarkdown_2.25 iterators_1.0.14 mitml_0.4-5 +[53] R6_2.5.1 boot_1.3-28.1 rpart_4.1.19 nnet_7.3-18 +[57] compiler_4.2.3
    diff --git a/chapter_12_files/figure-html/unnamed-chunk-3-1.png b/chapter_12_files/figure-html/unnamed-chunk-3-1.png index fca3066..736d5ad 100644 Binary files a/chapter_12_files/figure-html/unnamed-chunk-3-1.png and b/chapter_12_files/figure-html/unnamed-chunk-3-1.png differ diff --git a/chapter_16.html b/chapter_16.html index a8408af..288cdd6 100644 --- a/chapter_16.html +++ b/chapter_16.html @@ -288,12 +288,12 @@ 

    head(ds)
      studyid treat          z1         z2  y
-1       1     0  0.24140876  0.8999285 11
-2       1     1 -0.08337922 -1.2523967  4
-3       1     0  0.54133588  1.1391411 11
-4       1     0 -0.72552223  0.8181371 11
-5       1     0 -1.34843084  0.5481249 11
-6       1     1  0.57662307 -0.5316352  4
    +1 1 0 -0.06226787 -0.3046572 11 +2 1 1 0.22873835 -1.3755850 8 +3 1 1 -0.39734360 0.4309557 10 +4 1 0 -0.48128849 -0.8760153 11 +5 1 1 -0.07360034 -1.8460789 7 +6 1 0 -2.00590597 -0.6318041 10

    The simulated dataset contains information on the following variables:

    @@ -320,9 +320,9 @@

    0
    -(N=301)     +(N=296)     1
    -(N=299)     +(N=304) Overall
    (N=600)     @@ -336,15 +336,15 @@

    Mean (SD) --0.0562 (0.986) -0.0454 (0.890) --0.00558 (0.940) +0.0522 (1.00) +-0.0126 (0.990) +0.0194 (0.995) Median [Min, Max] --0.0694 [-2.72, 2.39] -0.0749 [-2.57, 2.10] -0.0164 [-2.72, 2.39] +0.0611 [-2.68, 3.33] +-0.0719 [-2.29, 3.32] +-0.0371 [-2.68, 3.33] z2 @@ -354,15 +354,15 @@

    Mean (SD) --0.0214 (0.959) --0.0183 (0.953) --0.0198 (0.955) +-0.0736 (0.929) +0.00776 (0.978) +-0.0324 (0.954) Median [Min, Max] --0.164 [-2.53, 2.59] --0.0303 [-2.36, 3.00] --0.0852 [-2.53, 3.00] +-0.0719 [-2.79, 2.79] +-0.0234 [-2.52, 2.77] +-0.0617 [-2.79, 2.79] studyid @@ -372,38 +372,38 @@

    1 -41 (13.6%) -59 (19.7%) +46 (15.5%) +54 (17.8%) 100 (16.7%) 2 -52 (17.3%) -48 (16.1%) +52 (17.6%) +48 (15.8%) 100 (16.7%) 3 -53 (17.6%) -47 (15.7%) +52 (17.6%) +48 (15.8%) 100 (16.7%) 4 -54 (17.9%) -46 (15.4%) +46 (15.5%) +54 (17.8%) 100 (16.7%) 5 -57 (18.9%) -43 (14.4%) +50 (16.9%) +50 (16.4%) 100 (16.7%) 6 -44 (14.6%) -56 (18.7%) +50 (16.9%) +50 (16.4%) 100 (16.7%) @@ -456,7 +456,7 @@

    +<environment: 0x5556cce09328>

    We can fit the treatment effect model as follows:

    @@ -479,36 +479,36 @@

    +delta[2] -2.6947 -2.4129 -2.1227 -1.8953 -1.5725 +gamma[1] -0.5627 -0.5429 -0.5233 -0.4980 -0.4536 +gamma[2] 0.4830 0.5322 0.5561 0.5734 0.6010 +sd 0.5870 0.7426 0.8418 0.9520 1.2676 @@ -519,12 +519,12 @@ 

    round(treatment.effect(ipd, samples, newpatient = c(z1 = 1, z2 = 0.5)), 2)
    0.025   0.5 0.975 
--4.38 -3.01 -1.91
    +-2.92 -2.32 -1.77
    -

    This means that the predicted outcome for patient with covariate values z1 = 1 and z2 = 0.5 will differ by -3.01 units when receiving the active treatment (treat = 1) as compared to the control treatment (treat = 0).

    +

    This means that the predicted outcome for patient with covariate values z1 = 1 and z2 = 0.5 will differ by -2.32 units when receiving the active treatment (treat = 1) as compared to the control treatment (treat = 0).

    We can also predict treatment benefit for all patients in the sample, and look at the distribution of predicted benefit.

    -
    +
    library(dplyr)
 library(ggplot2)
 
@@ -542,10 +542,11 @@ 
    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) + 
-  facet_wrap(~study) + 
-  ylab("Density") +
-  xlab("Predicted treatment benefit")  + theme_bw()
    + geom_vline(data = summbenefit, aes(xintercept = meanbenefit), + linewidth = 0.5, lty = 2) + + facet_wrap(~study) + + ylab("Density") + + xlab("Predicted treatment benefit") + theme_bw()
    @@ -583,7 +584,7 @@ 

    round(treatment.effect(ipd, samples, newpatient = c(1,0.5)), 2)

    0.025   0.5 0.975 
--3.69 -2.89 -2.04
    +-3.37 -2.28 -1.30
    @@ -598,12 +599,12 @@

    head(ds2)

      studyid treat         w1          w2 y
-1       1     1  1.1959217 -0.60042578 1
-2       1     0  0.4634961  0.08469719 0
-3       1     0  0.3403400  0.01961329 0
-4       1     0 -0.6867991  0.60877432 1
-5       1     0 -0.8644221 -0.07853599 0
-6       1     1 -1.3439234 -0.38140175 1
    +1 1 0 -0.5612891 0.65769118 1 +2 1 1 1.1342349 0.85682668 0 +3 1 1 0.1005061 -2.20312003 0 +4 1 1 -0.5153136 -0.57015407 0 +5 1 1 -1.0159330 -0.78384436 0 +6 1 0 -0.2962515 -0.01832073 0

    The simulated dataset contains information on the following variables:

    @@ -646,15 +647,15 @@

    Mean (SD) -0.0169 (0.987) --0.0413 (0.985) --0.0117 (0.986) +-0.00720 (1.02) +0.0526 (0.985) +0.0222 (1.00) Median [Min, Max] --0.0168 [-3.41, 2.50] --0.120 [-2.74, 3.02] --0.0624 [-3.41, 3.02] +-0.0492 [-2.73, 2.89] +0.0449 [-3.23, 3.33] +-0.0241 [-3.23, 3.33] w2 @@ -664,15 +665,15 @@

    Mean (SD) --0.111 (0.959) --0.0378 (1.07) --0.0750 (1.02) +0.0465 (1.00) +-0.0340 (1.01) +0.00695 (1.01) Median [Min, Max] --0.159 [-3.46, 2.03] --0.0395 [-3.29, 2.40] --0.125 [-3.46, 2.40] +0.0894 [-2.37, 2.93] +0.0399 [-2.72, 4.21] +0.0748 [-2.72, 4.21] studyid @@ -682,38 +683,38 @@

    1 -48 (15.7%) -52 (17.6%) +47 (15.4%) +53 (18.0%) 100 (16.7%) 2 -52 (17.0%) -48 (16.3%) +46 (15.1%) +54 (18.3%) 100 (16.7%) 3 -48 (15.7%) -52 (17.6%) +56 (18.4%) +44 (14.9%) 100 (16.7%) 4 -54 (17.7%) -46 (15.6%) +50 (16.4%) +50 (16.9%) 100 (16.7%) 5 -53 (17.4%) -47 (15.9%) +51 (16.7%) +49 (16.6%) 100 (16.7%) 6 -50 (16.4%) -50 (16.9%) +55 (18.0%) +45 (15.3%) 100 (16.7%) @@ -771,7 +772,7 @@

    +<environment: 0x5556d0a8e470>
    @@ -790,37 +791,37 @@

    + 2.5% 25% 50% 75% 97.5% +alpha[1] -0.4415 -0.2814894 -0.14889 -0.0004158 0.29590 +alpha[2] -1.8719 -1.5500880 -1.27938 -1.1610952 -0.93129 +alpha[3] -1.6696 -1.4417289 -1.15871 -1.0087951 -0.86845 +alpha[4] -1.5900 -1.2770702 -1.11096 -0.8363314 -0.43803 +alpha[5] -1.3737 -1.0954269 -0.92882 -0.7846082 -0.43321 +alpha[6] -1.6649 -1.3707223 -1.19790 -1.0033566 -0.84776 +beta[1] -0.3137 -0.1518571 -0.09500 -0.0398302 0.09449 +beta[2] -0.4063 -0.1586791 -0.13276 -0.0716571 0.10386 +delta[1] 0.0000 0.0000000 0.00000 0.0000000 0.00000 +delta[2] -2.0235 -0.8224947 -0.52883 -0.0107413 0.79044 +gamma[1] -0.1809 0.0003475 0.12003 0.2052612 0.34881 +gamma[2] -0.1912 -0.0437462 0.04734 0.1713008 0.46497 +sd 0.8210 1.1530501 1.41743 1.7317748 1.99153

    The predicted treatment benefit for a new patient with covariates w1 = 1.6 and w2 = 1.3 is given as:

    @@ -828,10 +829,10 @@ 

    round(treatment.effect(ipd2, samples, newpatient = c(w1 = 1.6, w2 = 1.3)), 2)
    0.025   0.5 0.975 
- 0.41  0.86  1.27
    + 0.15 0.91 2.30
    -

    In other words, the aforementioned patient 0.86 (95% Credibility Interval: 0.41 to 1.27)

    +

    In other words, the aforementioned patient 0.91 (95% Credibility Interval: 0.15 to 2.3)

    @@ -846,13 +847,13 @@ 

    ds3 <- generate_ipdnma_example(type = "continuous")
 head(ds3)
    -
      studyid treat          z1         z2  y
-1       1     2 -0.97337590 -0.6695393  8
-2       1     1 -1.02809238 -0.6817386 11
-3       1     2  2.35703179  1.6272272  8
-4       1     2 -0.04532726  0.6444098  9
-5       1     2 -0.11708756  1.6158472 10
-6       1     1  1.08949679 -0.6346255 11
    +
      studyid treat         z1         z2  y
+1       1     2 -0.7399958  0.4011151  9
+2       1     1  1.4632718 -1.2571879 11
+3       1     2 -0.6717355 -0.0582061  8
+4       1     2  0.3056371 -1.3880050  6
+5       1     1 -1.5699080  1.3314832 11
+6       1     2 -0.8608130  0.6028625  9

    Let us look into the data a bit in more detail:

    @@ -874,11 +875,11 @@

    1
    -(N=346)     +(N=328)     2
    -(N=350)     +(N=392) 3
    -(N=304)     +(N=280) Overall
    (N=1000)     @@ -893,17 +894,17 @@

    Mean (SD) --0.103 (1.01) -0.0226 (1.01) --0.0151 (1.07) --0.0325 (1.03) +0.0245 (0.992) +-0.0143 (1.02) +-0.0318 (0.980) +-0.00650 (0.998) Median [Min, Max] --0.160 [-3.18, 2.42] -0.0532 [-2.86, 3.01] --0.00116 [-3.34, 3.12] --0.0408 [-3.34, 3.12] +-0.0445 [-2.58, 2.53] +-0.0218 [-3.11, 3.44] +0.0391 [-3.62, 2.49] +0.00453 [-3.62, 3.44] z2 @@ -914,17 +915,17 @@

    Mean (SD) -0.0224 (0.891) -0.148 (1.07) -0.000382 (1.04) -0.0596 (1.00) +0.00244 (0.988) +0.00503 (0.993) +0.0784 (0.945) +0.0247 (0.978) Median [Min, Max] --0.0321 [-2.29, 2.92] -0.129 [-2.74, 3.44] -0.0611 [-2.85, 2.99] -0.0600 [-2.85, 3.44] +-0.0347 [-2.70, 3.00] +0.0170 [-2.93, 2.50] +0.0460 [-3.24, 2.42] +0.0178 [-3.24, 3.00] studyid @@ -935,72 +936,72 @@

    1 -52 (15.0%) -48 (13.7%) +39 (11.9%) +61 (15.6%) 0 (0%) 100 (10.0%) 2 -40 (11.6%) -60 (17.1%) +42 (12.8%) +58 (14.8%) 0 (0%) 100 (10.0%) 3 -50 (14.5%) -50 (14.3%) +48 (14.6%) +52 (13.3%) 0 (0%) 100 (10.0%) 4 -52 (15.0%) +53 (16.2%) 0 (0%) -48 (15.8%) +47 (16.8%) 100 (10.0%) 5 -47 (13.6%) +51 (15.5%) 0 (0%) -53 (17.4%) +49 (17.5%) 100 (10.0%) 6 0 (0%) -44 (12.6%) -56 (18.4%) +57 (14.5%) +43 (15.4%) 100 (10.0%) 7 0 (0%) -55 (15.7%) -45 (14.8%) +52 (13.3%) +48 (17.1%) 100 (10.0%) 8 -35 (10.1%) -35 (10.0%) -30 (9.9%) +28 (8.5%) +41 (10.5%) +31 (11.1%) 100 (10.0%) 9 -42 (12.1%) -29 (8.3%) -29 (9.5%) +37 (11.3%) +31 (7.9%) +32 (11.4%) 100 (10.0%) 10 -28 (8.1%) -29 (8.3%) -43 (14.1%) +30 (9.1%) +40 (10.2%) +30 (10.7%) 100 (10.0%) @@ -1071,7 +1072,7 @@

    +<environment: 0x5556d0d738d8> 
    samples <- ipd.run(ipd3, n.chains = 2, n.iter = 20, 
                    pars.save = c("alpha", "beta", "delta", "sd", "gamma"))
    @@ -1098,54 +1099,54 @@

    + 2.5% 25% 50% 75% 97.5% +alpha[1] 10.87592 10.92486 10.9432 10.9871 11.0311 +alpha[2] 7.95960 8.02352 8.0541 8.0954 8.1444 +alpha[3] 10.39161 10.43707 10.4611 10.4893 10.5395 +alpha[4] 9.54116 9.57556 9.5917 9.6213 9.6779 +alpha[5] 12.82662 12.87380 12.9012 12.9424 12.9958 +alpha[6] 13.12556 13.16652 13.2028 13.2304 13.2911 +alpha[7] 7.37369 7.41851 7.4599 7.4968 7.5248 +alpha[8] 11.03027 11.07663 11.1139 11.1397 11.2133 +alpha[9] 10.09934 10.15093 10.1636 10.2006 10.2633 +alpha[10] 9.15669 9.21080 9.2401 9.2838 9.3546 +beta[1] 0.16187 0.18112 0.1972 0.2054 0.2386 +beta[2] 0.25950 0.27918 0.2930 0.3026 0.3165 +delta[1] 0.00000 0.00000 0.0000 0.0000 0.0000 +delta[2] -3.00894 -2.94303 -2.9166 -2.8915 -2.8531 +delta[3] -1.25176 -1.20392 -1.1548 -1.1198 -1.0930 +gamma[1,1] 0.00000 0.00000 0.0000 0.0000 0.0000 +gamma[2,1] -0.61841 -0.57430 -0.5570 -0.5424 -0.5117 +gamma[3,1] -0.35406 -0.32074 -0.2948 -0.2799 -0.2645 +gamma[1,2] 0.00000 0.00000 0.0000 0.0000 0.0000 +gamma[2,2] 0.51180 0.54303 0.5596 0.5952 0.6236 +gamma[3,2] 0.39076 0.42440 0.4402 0.4571 0.4813 +sd 0.06018 0.09209 0.1048 0.1211 0.1674

    As before, we can use the treatment.effect() function of bipd to estimate relative effects for new patients.

    @@ -1154,11 +1155,11 @@ 

    $`treatment 2`
     0.025       0.5     0.975 
--2.641458 -2.463162 -2.249292 
+-2.543043 -2.331263 -2.197777 
 
 $`treatment 3`
      0.025        0.5      0.975 
--0.8048313 -0.6370218 -0.4668667
    +-0.7342422 -0.5628488 -0.4501180

    This gives us the relative effects for all treatments versus the reference. To obtain relative effects between active treatments we need some more coding:

    @@ -1175,7 +1176,7 @@

    samples.all$gamma.3.2.*newpatient[2]) )
    -
    [1] -1.806833
    +
    [1] -1.764445
    quantile(samples.all$delta.2.+samples.all$gamma.2.1.*
            newpatient[1]+samples.all$gamma.2.2.*newpatient[2]
@@ -1184,7 +1185,7 @@ 
             , probs = 0.025)
         2.5% 
--1.998861
    +-1.915623 
    quantile(samples.all$delta.2.+samples.all$gamma.2.1.*
            newpatient[1]+samples.all$gamma.2.2.*newpatient[2]
@@ -1193,7 +1194,7 @@ 
             , probs = 0.975)
        97.5% 
--1.625211
    +-1.651104
    @@ -1219,9 +1220,9 @@ 

              
            s1 s2 s3
-  treat A: 56 53 34
-  treat B: 44  0 29
-  treat C:  0 47 37
    + treat A: 42 44 36 + treat B: 58 0 33 + treat C: 0 56 31 @@ -1350,18 +1351,18 @@

    study 1 --2.968 (SE = 0.060 ) +-2.888 (SE = 0.048 ) study 2 --1.007 (SE = 0.074 ) +-1.106 (SE = 0.054 ) study 3 --2.930 (SE = 0.076 ) --0.909 (SE = 0.061 ) +-2.753 (SE = 0.062 ) +-0.933 (SE = 0.065 ) @@ -1427,23 +1428,23 @@

    +delta[2] -2.9742 -2.9214 -2.8965 -2.8662 -2.7813 +delta[3] -1.1950 -1.1387 -1.1114 -1.0835 -1.0364 +gamma[1,1] -0.8720 -0.7937 -0.7615 -0.7297 -0.6023 +gamma[2,1] 0.8906 0.9780 1.0053 1.0446 1.0874 +gamma[1,2] -0.6110 -0.5440 -0.5103 -0.4851 -0.4518 +gamma[2,2] 0.2981 0.3422 0.3660 0.3900 0.4520 
    # calculate  treatment effects
 samples.all = data.frame(rbind(samps4.3[[1]], samps4.3[[2]]))
@@ -1454,21 +1455,21 @@ 
        samples.all$gamma.2.1.*newpatient[2]
 )
    -
    [1] -2.023269
    +
    [1] -1.642012
    quantile(samples.all$delta.2.+samples.all$gamma.1.1.*newpatient[1]+
            samples.all$gamma.2.1.*newpatient[2]
          , probs = 0.025)
    -
         2.5% 
--2.304903
    +
        2.5% 
+-1.89312 
    quantile(samples.all$delta.2.+samples.all$gamma.1.1.*newpatient[1]+
            samples.all$gamma.2.1.*newpatient[2]
          , probs = 0.975)
        97.5% 
--1.770899
    +-1.432919
    diff --git a/chapter_16_files/figure-html/fig-predben_continuous_outcome-1.png b/chapter_16_files/figure-html/fig-predben_continuous_outcome-1.png index c1478f1..6215e8d 100644 Binary files a/chapter_16_files/figure-html/fig-predben_continuous_outcome-1.png and b/chapter_16_files/figure-html/fig-predben_continuous_outcome-1.png differ diff --git a/search.json b/search.json index 875596f..48a2516 100644 --- a/search.json +++ b/search.json @@ -228,7 +228,7 @@ "href": "chapter_10.html#network-meta-analysis-of-clinical-trials", "title": "6  Systematic review and meta-analysis of Real-World Evidence", "section": "6.3 Network meta-analysis of clinical trials", - "text": "6.3 Network meta-analysis of clinical trials\nWe here use the R packages netmeta for conducting a frequentist network meta-analysis. A detailed tutorial on the use of netmeta is available from the book Doing Meta-Analysis with R: A Hands-On Guide.\n\n6.3.1 Interventions for coronavirus disease 2019\nWe here consider data from a study which aimed to assess the comparative effectiveness of remdesivir and tocilizumab for reducing mortality in hospitalised COVID-19 patients. 80 trials were identified from two published network meta-analyses (Selvarajan et al. 2022), (Siemieniuk et al. 2020), a living COVID-19 trial database (COVID-NMA Initiative) [Covid-NMA.com], and a clinical trial database [clinicaltrials.gov]. Trials were included in this study if the patient population included hospitalized COVID-19 patients, active treatment was remdesivir or tocilizumab, comparator treatment was placebo or standard care, short-term mortality data was available, and the trial was published. 21 trials were included. For included trials, a risk of bias score was extracted from the COVID-NMA Initiative.\n\n\n\n\n\nstudlab\ntreat1\ntreat2\nevent1\nn1\nevent2\nn2\n\n\n\n\nAder\nREM\nSTD\n34\n414\n37\n418\n\n\nBeigel (ACTT-1)\nREM\nSTD\n59\n541\n77\n521\n\n\nBroman\nTOCI\nSTD\n1\n57\n0\n29\n\n\nCriner\nREM\nSTD\n4\n384\n4\n200\n\n\nDeclerq (COV-AID)\nTOCI\nSTD\n10\n81\n9\n74\n\n\nGordon (REMAP-CAP)\nTOCI\nSTD\n83\n353\n116\n358\n\n\nHermine (CORIMUNO)\nTOCI\nSTD\n7\n63\n8\n67\n\n\nHorby (RECOVERY)\nTOCI\nSTD\n621\n2022\n729\n2094\n\n\nIslam\nREM\nSTD\n0\n30\n0\n30\n\n\nMahajan\nREM\nSTD\n5\n34\n3\n36\n\n\nPan (WHO Solidarity)\nREM\nSTD\n602\n4146\n643\n4129\n\n\nRosas (COVACTA)\nTOCI\nSTD\n58\n294\n28\n144\n\n\nRutgers\nTOCI\nSTD\n21\n174\n34\n180\n\n\nSalama (EMPACTA)\nTOCI\nSTD\n26\n249\n11\n128\n\n\nSalvarani\nTOCI\nSTD\n2\n60\n1\n63\n\n\nSoin (COVINTOC)\nTOCI\nSTD\n11\n92\n15\n88\n\n\nSpinner\nREM\nSTD\n5\n384\n4\n200\n\n\nStone (BACC-BAY)\nTOCI\nSTD\n9\n161\n4\n82\n\n\nTalaschian\nTOCI\nSTD\n5\n17\n4\n19\n\n\nVeiga (TOCIBRAS)\nTOCI\nSTD\n14\n65\n6\n64\n\n\nWang\nREM\nSTD\n22\n158\n10\n78\n\n\n\n\n\n\n\nThe corresponding network is displayed below:\n\n\n\n\n\nEvidence network of the 21 coronavirus-19 trials\n\n\n\n\nWe use the following command to calculate the log odds ratios and corresponding standard errors for each study:\n\ncovid <- pairwise(treat = treat, \n event = event, \n n = n, \n studlab = studlab, \n sm = \"OR\")\nhead(covid)\n\n\n\n\n\n\nTE\nseTE\nstudlab\ntreat1\ntreat2\nevent1\nn1\nevent2\nn2\nincr\nallstudies\n\n\n\n\n-0.0819293\n0.2483849\nAder\nREM\nSTD\n34\n414\n37\n418\n0.0\nFALSE\n\n\n-0.3483875\n0.1851030\nBeigel (ACTT-1)\nREM\nSTD\n59\n541\n77\n521\n0.0\nFALSE\n\n\n0.4487619\n1.6487159\nBroman\nTOCI\nSTD\n1\n57\n0\n29\n0.5\nFALSE\n\n\n-0.6620566\n0.7125543\nCriner\nREM\nSTD\n4\n384\n4\n200\n0.0\nFALSE\n\n\n0.0170679\n0.4904898\nDeclerq (COV-AID)\nTOCI\nSTD\n10\n81\n9\n74\n0.0\nFALSE\n\n\n-0.4442338\n0.1688337\nGordon (REMAP-CAP)\nTOCI\nSTD\n83\n353\n116\n358\n0.0\nFALSE\n\n\n\n\n\n\n\nBelow, we conduct a random effects network meta-analysis where we consider standard care (STD) as the control treatment. Note that we have one study where zero cell counts occur, this study will not contribute to the NMA as the log odds ratio and its standard error cannot be determined.\n\nNMA.covid <- netmeta(TE = TE, seTE = seTE, treat1 = treat1, treat2 = treat2,\n studlab = studlab, data = covid, sm = \"OR\", ref = \"STD\",\n comb.random = TRUE, common = FALSE, warn = FALSE)\nNMA.covid \n\nNumber of studies: k = 20\nNumber of pairwise comparisons: m = 20\nNumber of treatments: n = 3\nNumber of designs: d = 2\n\nRandom effects model\n\nTreatment estimate (sm = 'OR', comparison: other treatments vs 'STD'):\n OR 95%-CI z p-value\nREM 0.8999 [0.8067; 1.0039] -1.89 0.0588\nSTD . . . .\nTOCI 0.8301 [0.7434; 0.9268] -3.31 0.0009\n\nQuantifying heterogeneity / inconsistency:\ntau^2 = 0; tau = 0; I^2 = 0% [0.0%; 48.9%]\n\nTests of heterogeneity (within designs) and inconsistency (between designs):\n Q d.f. p-value\nTotal 16.38 18 0.5663\nWithin designs 16.38 18 0.5663\nBetween designs 0.00 0 --\n\n\nA league table of the treatment effect estimates is given below:\n\nnetleague(NMA.covid)\n\nLeague table (random effects model):\n \n REM 0.8999 [0.8067; 1.0039] .\n 0.8999 [0.8067; 1.0039] STD 1.2047 [1.0789; 1.3451]\n 1.0842 [0.9282; 1.2663] 1.2047 [1.0789; 1.3451] TOCI\n\n\nWe can also present the results in a forest plot:\n\n\n\n\n\nWe now consider a Bayesian random effects network meta-analysis that analyzes the observed event counts using a binomial link function.\n\nbdata <- data.frame(study = studlab,\n treatment = treat,\n responders = event,\n sampleSize = n)\n\nnetwork <- mtc.network(data.ab = bdata)\n\nmodel <- mtc.model(network,\n likelihood = \"binom\",\n link = \"log\",\n linearModel = \"random\",\n n.chain = 3)\n\n\n# Adaptation\nmcmc1 <- mtc.run(model, n.adapt = 1000, n.iter = 1000, thin = 10)\n\nCompiling model graph\n Resolving undeclared variables\n Allocating nodes\nGraph information:\n Observed stochastic nodes: 42\n Unobserved stochastic nodes: 45\n Total graph size: 930\n\nInitializing model\n\n# Sampling\nmcmc2 <- mtc.run(model, n.adapt = 10000, n.iter = 100000, thin = 10)\n\nCompiling model graph\n Resolving undeclared variables\n Allocating nodes\nGraph information:\n Observed stochastic nodes: 42\n Unobserved stochastic nodes: 45\n Total graph size: 930\n\nInitializing model\n\n\nWe can extract the pooled treatment effect estimates from the posterior distribution. When using STD as control group, we have:\n\nsummary(relative.effect(mcmc2, t1 = \"STD\"))\n\n\nResults on the Log Risk Ratio scale\n\nIterations = 10010:110000\nThinning interval = 10 \nNumber of chains = 3 \nSample size per chain = 10000 \n\n1. Empirical mean and standard deviation for each variable,\n plus standard error of the mean:\n\n Mean SD Naive SE Time-series SE\nd.STD.REM -0.1065 0.09808 0.0005663 0.0008046\nd.STD.TOCI -0.1121 0.08177 0.0004721 0.0008786\nsd.d 0.1117 0.08670 0.0005006 0.0017016\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nd.STD.REM -0.318606 -0.15967 -0.10277 -0.04958 0.08610\nd.STD.TOCI -0.256438 -0.16289 -0.11923 -0.06945 0.07634\nsd.d 0.004353 0.04436 0.09359 0.15819 0.32473\n\n\nThe corresponding odds ratios are as follows:\n\n\n\n\n\nComparison\n95% CrI\n\n\n\n\nREM vs. STD\n0.9 (0.73; 1.09)\n\n\nTOCI vs. STD\n0.89 (0.77; 1.08)\n\n\nREM vs. TOCI\n1.02 (0.75; 1.27)\n\n\n\n\n\n\n\nFinally, we expand the COVID-19 network with trials investigating the effectiveness of hydroxychloroquine (HCQ), lopinavir/ritonavir (LOPI), dexamethasone (DEXA) or interferon-\\(\\beta\\) (INTB) (Selvarajan et al. 2022). The corresponding network is displayed below:\n\n\n\n\n\nEvidence network of the 33 coronavirus-19 trials\n\n\n\n\nWe conducted a random effects network meta-analysis, results are depicted below:\n\n\nNumber of studies: k = 33\nNumber of pairwise comparisons: m = 33\nNumber of treatments: n = 7\nNumber of designs: d = 6\n\nRandom effects model\n\nTreatment estimate (sm = 'OR', comparison: other treatments vs 'STD'):\n OR 95%-CI z p-value 95%-PI\nDEXA 0.8557 [0.7558; 0.9688] -2.46 0.0139 [0.7463; 0.9812]\nHCQ 1.1809 [0.8934; 1.5610] 1.17 0.2428 [0.8786; 1.5872]\nINTB 1.1606 [0.9732; 1.3841] 1.66 0.0973 [0.9604; 1.4026]\nLOPI 1.0072 [0.8906; 1.1392] 0.11 0.9085 [0.8794; 1.1537]\nREM 0.8983 [0.8014; 1.0070] -1.84 0.0658 [0.7913; 1.0199]\nSTD . . . . .\nTOCI 0.8304 [0.7410; 0.9306] -3.20 0.0014 [0.7316; 0.9426]\n\nQuantifying heterogeneity / inconsistency:\ntau^2 = 0.0004; tau = 0.0205; I^2 = 0.6% [0.0%; 42.3%]\n\nTests of heterogeneity (within designs) and inconsistency (between designs):\n Q d.f. p-value\nTotal 27.18 27 0.4543\nWithin designs 27.18 27 0.4543\nBetween designs 0.00 0 --\n\n\nWe can calculate the P score for each treatment as follows:\n\nnetrank(NMA.covidf)\n\n P-score\nTOCI 0.9070\nDEXA 0.8357\nREM 0.7143\nSTD 0.4027\nLOPI 0.3899\nHCQ 0.1336\nINTB 0.1166\n\n\n\n\n6.3.2 Pharmacologic treatments for chronic obstructive pulmonary disease\nIn this example, we consider the resuls from a systematic review of randomized controlled trials on pharmacologic treatments for chronic obstructive pulmonary disease (Baker, Baker, and Coleman 2009). The primary outcome, occurrence of one or more episodes of COPD exacerbation, is binary (yes / no). For this outcome, five drug treatments (fluticasone, budesonide, salmeterol, formoterol, tiotropium) and two combinations (fluticasone + salmeterol, budesonide + formoterol) were compared to placebo. The authors considered the two combinations as separate treatments instead of evaluating the individual components.\n\ndata(Baker2009)\n\n\n\n\n\n\nstudy\nyear\nid\ntreatment\nexac\ntotal\n\n\n\n\nLlewellyn-Jones 1996\n1996\n1\nFluticasone\n0\n8\n\n\nLlewellyn-Jones 1996\n1996\n1\nPlacebo\n3\n8\n\n\nBoyd 1997\n1997\n2\nSalmeterol\n47\n229\n\n\nBoyd 1997\n1997\n2\nPlacebo\n59\n227\n\n\nPaggiaro 1998\n1998\n3\nFluticasone\n45\n142\n\n\nPaggiaro 1998\n1998\n3\nPlacebo\n51\n139\n\n\n\n\n\n\n\n\nBaker <- pairwise(treat = treatment,\n event = exac,\n n = total,\n studlab = id,\n sm = \"OR\",\n data = Baker2009)\n\nNMA.COPD <- netmeta(TE = TE, seTE = seTE, treat1 = treat1, treat2 = treat2,\n studlab = studlab, data = Baker, sm = \"OR\", ref = \"Placebo\",\n comb.random = TRUE)\n\nWarning: Comparisons with missing TE / seTE or zero seTE not considered in\nnetwork meta-analysis.\n\n\nComparisons not considered in network meta-analysis:\n studlab treat1 treat2 TE seTE\n 39 Fluticasone+Salmeterol Placebo NA NA\n 39 Fluticasone+Salmeterol Salmeterol NA NA\n 39 Salmeterol Placebo NA NA\n\nnetgraph(NMA.COPD)\n\n\n\n\n\n\n6.3.3 Advanced Therapies for Ulcerative Colitis\nIn this example, we consider a systematic literature review of Phase 3 randomized controlled trials investigating the following advanced therapies: infliximab, adalimumab, vedolizumab, golimumab, tofacitinib, ustekinumab, filgotinib, ozanimod, and upadacitinib (Panaccione et al. 2023). This review included 48 RCTs, from which 23 were found eligible for inclusion in a network meta-analysis. The included RCT populations were largely comparable in their baseline characteristics, though some heterogeneity was noted in weight, disease duration, extent of disease, and concomitant medications. A risk of bias assessment showed a low risk of bias for all included RCTs, which were all industry sponsored.\nWe here focus on the synthesis of 18 trials that contributed efficacy data for induction in bio-naive populations. The following FDA- and/or EMA-approved biologic or SMD doses were investigated:\n\nAdalimumab subcutaneous 160 mg at week 0, 80 mg at week 2, and 40 mg at week 4 (ADA160/80)\nInfliximab intravenous 5 mg/kg (INF5) at weeks 0, 2, and 6 then every 8 weeks\nInfliximab intravenous 10 mg/kg (INF10) at weeks 0, 2, and 6 then every 8 weeks\nFilgotinib oral 100 mg once daily (FIL100)\nFilgotinib oral 200 mg once daily (FIL200)\nGolimumab subcutaneous 200 mg at week 0 and 100 mg at week 2 (GOL200/100)\nOzanimod oral 0.23 mg once daily for 4 days, 0.46 mg once daily for 3 days, then 0.92 mg once daily (OZA0.92)\nTofacitinib oral 10 mg twice daily for 8 weeks (TOF10)\nUpadacitinib oral 45 mg once daily for 8 weeks (UPA45)\nUstekinumab intravenous 6 mg/kg at week 0 (UST6)\nVedolizumab intravenous 300 mg at weeks 0, 2, and 6 (VED300)\n\nThe reference treatment is placebo (PBO).\n\n\n\nEfficacy outcomes (i.e., clinical remission) data of induction bio-naïve populations \n\n\nstudlab\ntreat1\ntreat2\nevent1\nn1\nevent2\nn2\n\n\n\n\nACT-1\nINF10\nINF5\n39\n122\n47\n121\n\n\nACT-1\nINF10\nPBO\n39\n122\n18\n121\n\n\nACT-1\nINF5\nPBO\n47\n121\n18\n121\n\n\nACT-2\nINF10\nINF5\n33\n120\n41\n121\n\n\nACT-2\nINF10\nPBO\n33\n120\n7\n123\n\n\nACT-2\nINF5\nPBO\n41\n121\n7\n123\n\n\nGEMINI 1\nVED300\nPBO\n30\n130\n5\n76\n\n\nJapic CTI-060298\nINF5\nPBO\n21\n104\n11\n104\n\n\nJiang 2015\nINF5\nPBO\n22\n41\n9\n41\n\n\nM10-447\nADA160/80\nPBO\n9\n90\n11\n96\n\n\nNCT01551290\nINF5\nPBO\n11\n50\n5\n49\n\n\nNCT02039505\nVED300\nPBO\n22\n79\n6\n41\n\n\nOCTAVE 1\nTOF10\nPBO\n56\n222\n9\n57\n\n\nOCTAVE 2\nTOF10\nPBO\n43\n195\n4\n47\n\n\nPURSUIT-SC\nGOL200/100\nPBO\n45\n253\n16\n251\n\n\nSELECTION\nFIL100\nFIL200\n47\n277\n60\n245\n\n\nSELECTION\nFIL100\nPBO\n47\n277\n17\n137\n\n\nSELECTION\nFIL200\nPBO\n60\n245\n17\n137\n\n\nTRUE NORTH\nOZA0.92\nPBO\n66\n299\n10\n151\n\n\nU-ACCOMPLISH\nUPA45\nPBO\n54\n166\n3\n81\n\n\nU-ACHIEVE Study 2\nUPA45\nPBO\n41\n145\n4\n72\n\n\nULTRA-1\nADA160/80\nPBO\n24\n130\n12\n130\n\n\nULTRA-2\nADA160/80\nPBO\n32\n150\n16\n145\n\n\nUNIFI\nUST6\nPBO\n27\n147\n15\n151\n\n\n\n\n\n\n\nThe corresponding network is displayed below:\n\n\n\n\n\nEvidence network of 18 trials that contributed efficacy data for induction in bio-naive populations\n\n\n\n\nBelow, we conduct a random effects network meta-analysis of the reported study effects (expressed as odds ratio) and consider placebo (treat = \"PBO\") as the control treatment.\n\nNMA.uc <- netmeta(TE = TE, seTE = seTE, treat1 = treat1, treat2 = treat2,\n studlab = studlab, data = UlcerativeColitis, sm = \"OR\", \n ref = \"PBO\", common = FALSE, comb.random = TRUE)\nNMA.uc\n\nAll treatments except FIL100 and UST6 are significantly more efficacious than PBO at inducing clinical remission. We can now estimate the probabilities of each treatment being at each possible rank and the SUCRAs (Surface Under the Cumulative RAnking curve):\n\nsucra.uc <- rankogram(NMA.uc, nsim = 100, random = TRUE, common = FALSE, \n small.values = \"undesirable\")\n\n# Exctract the SUCRA values\nsucra.uc$ranking.random\n\n ADA160/80 FIL100 FIL200 GOL200/100 INF10 INF5 OZA0.92 \n0.27000000 0.16363636 0.46727273 0.63727273 0.57272727 0.75636364 0.75363636 \n PBO TOF10 UPA45 UST6 VED300 \n0.02090909 0.38909091 0.96727273 0.37000000 0.63181818 \n\n\nThese results indicate that 96.7% of the evaluated treatments are worse than UPA45." + "text": "6.3 Network meta-analysis of clinical trials\nWe here use the R packages netmeta for conducting a frequentist network meta-analysis. A detailed tutorial on the use of netmeta is available from the book Doing Meta-Analysis with R: A Hands-On Guide.\n\n6.3.1 Interventions for coronavirus disease 2019\nWe here consider data from a study which aimed to assess the comparative effectiveness of remdesivir and tocilizumab for reducing mortality in hospitalised COVID-19 patients. 80 trials were identified from two published network meta-analyses (Selvarajan et al. 2022), (Siemieniuk et al. 2020), a living COVID-19 trial database (COVID-NMA Initiative) [Covid-NMA.com], and a clinical trial database [clinicaltrials.gov]. Trials were included in this study if the patient population included hospitalized COVID-19 patients, active treatment was remdesivir or tocilizumab, comparator treatment was placebo or standard care, short-term mortality data was available, and the trial was published. 21 trials were included. For included trials, a risk of bias score was extracted from the COVID-NMA Initiative.\n\n\n\n\n\nstudlab\ntreat1\ntreat2\nevent1\nn1\nevent2\nn2\n\n\n\n\nAder\nREM\nSTD\n34\n414\n37\n418\n\n\nBeigel (ACTT-1)\nREM\nSTD\n59\n541\n77\n521\n\n\nBroman\nTOCI\nSTD\n1\n57\n0\n29\n\n\nCriner\nREM\nSTD\n4\n384\n4\n200\n\n\nDeclerq (COV-AID)\nTOCI\nSTD\n10\n81\n9\n74\n\n\nGordon (REMAP-CAP)\nTOCI\nSTD\n83\n353\n116\n358\n\n\nHermine (CORIMUNO)\nTOCI\nSTD\n7\n63\n8\n67\n\n\nHorby (RECOVERY)\nTOCI\nSTD\n621\n2022\n729\n2094\n\n\nIslam\nREM\nSTD\n0\n30\n0\n30\n\n\nMahajan\nREM\nSTD\n5\n34\n3\n36\n\n\nPan (WHO Solidarity)\nREM\nSTD\n602\n4146\n643\n4129\n\n\nRosas (COVACTA)\nTOCI\nSTD\n58\n294\n28\n144\n\n\nRutgers\nTOCI\nSTD\n21\n174\n34\n180\n\n\nSalama (EMPACTA)\nTOCI\nSTD\n26\n249\n11\n128\n\n\nSalvarani\nTOCI\nSTD\n2\n60\n1\n63\n\n\nSoin (COVINTOC)\nTOCI\nSTD\n11\n92\n15\n88\n\n\nSpinner\nREM\nSTD\n5\n384\n4\n200\n\n\nStone (BACC-BAY)\nTOCI\nSTD\n9\n161\n4\n82\n\n\nTalaschian\nTOCI\nSTD\n5\n17\n4\n19\n\n\nVeiga (TOCIBRAS)\nTOCI\nSTD\n14\n65\n6\n64\n\n\nWang\nREM\nSTD\n22\n158\n10\n78\n\n\n\n\n\n\n\nThe corresponding network is displayed below:\n\n\n\n\n\nEvidence network of the 21 coronavirus-19 trials\n\n\n\n\nWe use the following command to calculate the log odds ratios and corresponding standard errors for each study:\n\ncovid <- pairwise(treat = treat, \n event = event, \n n = n, \n studlab = studlab, \n sm = \"OR\")\nhead(covid)\n\n\n\n\n\n\nTE\nseTE\nstudlab\ntreat1\ntreat2\nevent1\nn1\nevent2\nn2\nincr\nallstudies\n\n\n\n\n-0.0819293\n0.2483849\nAder\nREM\nSTD\n34\n414\n37\n418\n0.0\nFALSE\n\n\n-0.3483875\n0.1851030\nBeigel (ACTT-1)\nREM\nSTD\n59\n541\n77\n521\n0.0\nFALSE\n\n\n0.4487619\n1.6487159\nBroman\nTOCI\nSTD\n1\n57\n0\n29\n0.5\nFALSE\n\n\n-0.6620566\n0.7125543\nCriner\nREM\nSTD\n4\n384\n4\n200\n0.0\nFALSE\n\n\n0.0170679\n0.4904898\nDeclerq (COV-AID)\nTOCI\nSTD\n10\n81\n9\n74\n0.0\nFALSE\n\n\n-0.4442338\n0.1688337\nGordon (REMAP-CAP)\nTOCI\nSTD\n83\n353\n116\n358\n0.0\nFALSE\n\n\n\n\n\n\n\nBelow, we conduct a random effects network meta-analysis where we consider standard care (STD) as the control treatment. Note that we have one study where zero cell counts occur, this study will not contribute to the NMA as the log odds ratio and its standard error cannot be determined.\n\nNMA.covid <- netmeta(TE = TE, seTE = seTE, treat1 = treat1, treat2 = treat2,\n studlab = studlab, data = covid, sm = \"OR\", ref = \"STD\",\n comb.random = TRUE, common = FALSE, warn = FALSE)\nNMA.covid \n\nNumber of studies: k = 20\nNumber of pairwise comparisons: m = 20\nNumber of treatments: n = 3\nNumber of designs: d = 2\n\nRandom effects model\n\nTreatment estimate (sm = 'OR', comparison: other treatments vs 'STD'):\n OR 95%-CI z p-value\nREM 0.8999 [0.8067; 1.0039] -1.89 0.0588\nSTD . . . .\nTOCI 0.8301 [0.7434; 0.9268] -3.31 0.0009\n\nQuantifying heterogeneity / inconsistency:\ntau^2 = 0; tau = 0; I^2 = 0% [0.0%; 48.9%]\n\nTests of heterogeneity (within designs) and inconsistency (between designs):\n Q d.f. p-value\nTotal 16.38 18 0.5663\nWithin designs 16.38 18 0.5663\nBetween designs 0.00 0 --\n\n\nA league table of the treatment effect estimates is given below:\n\nnetleague(NMA.covid)\n\nLeague table (random effects model):\n \n REM 0.8999 [0.8067; 1.0039] .\n 0.8999 [0.8067; 1.0039] STD 1.2047 [1.0789; 1.3451]\n 1.0842 [0.9282; 1.2663] 1.2047 [1.0789; 1.3451] TOCI\n\n\nWe can also present the results in a forest plot:\n\n\n\n\n\nWe now consider a Bayesian random effects network meta-analysis that analyzes the observed event counts using a binomial link function.\n\nbdata <- data.frame(study = studlab,\n treatment = treat,\n responders = event,\n sampleSize = n)\n\nnetwork <- mtc.network(data.ab = bdata)\n\nmodel <- mtc.model(network,\n likelihood = \"binom\",\n link = \"log\",\n linearModel = \"random\",\n n.chain = 3)\n\n\n# Adaptation\nmcmc1 <- mtc.run(model, n.adapt = 1000, n.iter = 1000, thin = 10)\n\nCompiling model graph\n Resolving undeclared variables\n Allocating nodes\nGraph information:\n Observed stochastic nodes: 42\n Unobserved stochastic nodes: 45\n Total graph size: 930\n\nInitializing model\n\n# Sampling\nmcmc2 <- mtc.run(model, n.adapt = 10000, n.iter = 100000, thin = 10)\n\nCompiling model graph\n Resolving undeclared variables\n Allocating nodes\nGraph information:\n Observed stochastic nodes: 42\n Unobserved stochastic nodes: 45\n Total graph size: 930\n\nInitializing model\n\n\nWe can extract the pooled treatment effect estimates from the posterior distribution. When using STD as control group, we have:\n\nsummary(relative.effect(mcmc2, t1 = \"STD\"))\n\n\nResults on the Log Risk Ratio scale\n\nIterations = 10010:110000\nThinning interval = 10 \nNumber of chains = 3 \nSample size per chain = 10000 \n\n1. Empirical mean and standard deviation for each variable,\n plus standard error of the mean:\n\n Mean SD Naive SE Time-series SE\nd.STD.REM -0.1070 0.09687 0.0005593 0.0007884\nd.STD.TOCI -0.1119 0.08258 0.0004768 0.0008304\nsd.d 0.1122 0.08608 0.0004970 0.0017143\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nd.STD.REM -0.311299 -0.16110 -0.10389 -0.05106 0.08548\nd.STD.TOCI -0.256925 -0.16347 -0.11929 -0.06926 0.08092\nsd.d 0.004401 0.04555 0.09467 0.15862 0.32142\n\n\nThe corresponding odds ratios are as follows:\n\n\n\n\n\nComparison\n95% CrI\n\n\n\n\nREM vs. STD\n0.9 (0.73; 1.09)\n\n\nTOCI vs. STD\n0.89 (0.77; 1.08)\n\n\nREM vs. TOCI\n1.02 (0.75; 1.27)\n\n\n\n\n\n\n\nFinally, we expand the COVID-19 network with trials investigating the effectiveness of hydroxychloroquine (HCQ), lopinavir/ritonavir (LOPI), dexamethasone (DEXA) or interferon-\\(\\beta\\) (INTB) (Selvarajan et al. 2022). The corresponding network is displayed below:\n\n\n\n\n\nEvidence network of the 33 coronavirus-19 trials\n\n\n\n\nWe conducted a random effects network meta-analysis, results are depicted below:\n\n\nNumber of studies: k = 33\nNumber of pairwise comparisons: m = 33\nNumber of treatments: n = 7\nNumber of designs: d = 6\n\nRandom effects model\n\nTreatment estimate (sm = 'OR', comparison: other treatments vs 'STD'):\n OR 95%-CI z p-value 95%-PI\nDEXA 0.8557 [0.7558; 0.9688] -2.46 0.0139 [0.7463; 0.9812]\nHCQ 1.1809 [0.8934; 1.5610] 1.17 0.2428 [0.8786; 1.5872]\nINTB 1.1606 [0.9732; 1.3841] 1.66 0.0973 [0.9604; 1.4026]\nLOPI 1.0072 [0.8906; 1.1392] 0.11 0.9085 [0.8794; 1.1537]\nREM 0.8983 [0.8014; 1.0070] -1.84 0.0658 [0.7913; 1.0199]\nSTD . . . . .\nTOCI 0.8304 [0.7410; 0.9306] -3.20 0.0014 [0.7316; 0.9426]\n\nQuantifying heterogeneity / inconsistency:\ntau^2 = 0.0004; tau = 0.0205; I^2 = 0.6% [0.0%; 42.3%]\n\nTests of heterogeneity (within designs) and inconsistency (between designs):\n Q d.f. p-value\nTotal 27.18 27 0.4543\nWithin designs 27.18 27 0.4543\nBetween designs 0.00 0 --\n\n\nWe can calculate the P score for each treatment as follows:\n\nnetrank(NMA.covidf)\n\n P-score\nTOCI 0.9070\nDEXA 0.8357\nREM 0.7143\nSTD 0.4027\nLOPI 0.3899\nHCQ 0.1336\nINTB 0.1166\n\n\n\n\n6.3.2 Pharmacologic treatments for chronic obstructive pulmonary disease\nIn this example, we consider the resuls from a systematic review of randomized controlled trials on pharmacologic treatments for chronic obstructive pulmonary disease (Baker, Baker, and Coleman 2009). The primary outcome, occurrence of one or more episodes of COPD exacerbation, is binary (yes / no). For this outcome, five drug treatments (fluticasone, budesonide, salmeterol, formoterol, tiotropium) and two combinations (fluticasone + salmeterol, budesonide + formoterol) were compared to placebo. The authors considered the two combinations as separate treatments instead of evaluating the individual components.\n\ndata(Baker2009)\n\n\n\n\n\n\nstudy\nyear\nid\ntreatment\nexac\ntotal\n\n\n\n\nLlewellyn-Jones 1996\n1996\n1\nFluticasone\n0\n8\n\n\nLlewellyn-Jones 1996\n1996\n1\nPlacebo\n3\n8\n\n\nBoyd 1997\n1997\n2\nSalmeterol\n47\n229\n\n\nBoyd 1997\n1997\n2\nPlacebo\n59\n227\n\n\nPaggiaro 1998\n1998\n3\nFluticasone\n45\n142\n\n\nPaggiaro 1998\n1998\n3\nPlacebo\n51\n139\n\n\n\n\n\n\n\n\nBaker <- pairwise(treat = treatment,\n event = exac,\n n = total,\n studlab = id,\n sm = \"OR\",\n data = Baker2009)\n\nNMA.COPD <- netmeta(TE = TE, seTE = seTE, treat1 = treat1, treat2 = treat2,\n studlab = studlab, data = Baker, sm = \"OR\", ref = \"Placebo\",\n comb.random = TRUE)\n\nWarning: Comparisons with missing TE / seTE or zero seTE not considered in\nnetwork meta-analysis.\n\n\nComparisons not considered in network meta-analysis:\n studlab treat1 treat2 TE seTE\n 39 Fluticasone+Salmeterol Placebo NA NA\n 39 Fluticasone+Salmeterol Salmeterol NA NA\n 39 Salmeterol Placebo NA NA\n\nnetgraph(NMA.COPD)\n\n\n\n\n\n\n6.3.3 Advanced Therapies for Ulcerative Colitis\nIn this example, we consider a systematic literature review of Phase 3 randomized controlled trials investigating the following advanced therapies: infliximab, adalimumab, vedolizumab, golimumab, tofacitinib, ustekinumab, filgotinib, ozanimod, and upadacitinib (Panaccione et al. 2023). This review included 48 RCTs, from which 23 were found eligible for inclusion in a network meta-analysis. The included RCT populations were largely comparable in their baseline characteristics, though some heterogeneity was noted in weight, disease duration, extent of disease, and concomitant medications. A risk of bias assessment showed a low risk of bias for all included RCTs, which were all industry sponsored.\nWe here focus on the synthesis of 18 trials that contributed efficacy data for induction in bio-naive populations. The following FDA- and/or EMA-approved biologic or SMD doses were investigated:\n\nAdalimumab subcutaneous 160 mg at week 0, 80 mg at week 2, and 40 mg at week 4 (ADA160/80)\nInfliximab intravenous 5 mg/kg (INF5) at weeks 0, 2, and 6 then every 8 weeks\nInfliximab intravenous 10 mg/kg (INF10) at weeks 0, 2, and 6 then every 8 weeks\nFilgotinib oral 100 mg once daily (FIL100)\nFilgotinib oral 200 mg once daily (FIL200)\nGolimumab subcutaneous 200 mg at week 0 and 100 mg at week 2 (GOL200/100)\nOzanimod oral 0.23 mg once daily for 4 days, 0.46 mg once daily for 3 days, then 0.92 mg once daily (OZA0.92)\nTofacitinib oral 10 mg twice daily for 8 weeks (TOF10)\nUpadacitinib oral 45 mg once daily for 8 weeks (UPA45)\nUstekinumab intravenous 6 mg/kg at week 0 (UST6)\nVedolizumab intravenous 300 mg at weeks 0, 2, and 6 (VED300)\n\nThe reference treatment is placebo (PBO).\n\n\n\nEfficacy outcomes (i.e., clinical remission) data of induction bio-naïve populations \n\n\nstudlab\ntreat1\ntreat2\nevent1\nn1\nevent2\nn2\n\n\n\n\nACT-1\nINF10\nINF5\n39\n122\n47\n121\n\n\nACT-1\nINF10\nPBO\n39\n122\n18\n121\n\n\nACT-1\nINF5\nPBO\n47\n121\n18\n121\n\n\nACT-2\nINF10\nINF5\n33\n120\n41\n121\n\n\nACT-2\nINF10\nPBO\n33\n120\n7\n123\n\n\nACT-2\nINF5\nPBO\n41\n121\n7\n123\n\n\nGEMINI 1\nVED300\nPBO\n30\n130\n5\n76\n\n\nJapic CTI-060298\nINF5\nPBO\n21\n104\n11\n104\n\n\nJiang 2015\nINF5\nPBO\n22\n41\n9\n41\n\n\nM10-447\nADA160/80\nPBO\n9\n90\n11\n96\n\n\nNCT01551290\nINF5\nPBO\n11\n50\n5\n49\n\n\nNCT02039505\nVED300\nPBO\n22\n79\n6\n41\n\n\nOCTAVE 1\nTOF10\nPBO\n56\n222\n9\n57\n\n\nOCTAVE 2\nTOF10\nPBO\n43\n195\n4\n47\n\n\nPURSUIT-SC\nGOL200/100\nPBO\n45\n253\n16\n251\n\n\nSELECTION\nFIL100\nFIL200\n47\n277\n60\n245\n\n\nSELECTION\nFIL100\nPBO\n47\n277\n17\n137\n\n\nSELECTION\nFIL200\nPBO\n60\n245\n17\n137\n\n\nTRUE NORTH\nOZA0.92\nPBO\n66\n299\n10\n151\n\n\nU-ACCOMPLISH\nUPA45\nPBO\n54\n166\n3\n81\n\n\nU-ACHIEVE Study 2\nUPA45\nPBO\n41\n145\n4\n72\n\n\nULTRA-1\nADA160/80\nPBO\n24\n130\n12\n130\n\n\nULTRA-2\nADA160/80\nPBO\n32\n150\n16\n145\n\n\nUNIFI\nUST6\nPBO\n27\n147\n15\n151\n\n\n\n\n\n\n\nThe corresponding network is displayed below:\n\n\n\n\n\nEvidence network of 18 trials that contributed efficacy data for induction in bio-naive populations\n\n\n\n\nBelow, we conduct a random effects network meta-analysis of the reported study effects (expressed as odds ratio) and consider placebo (treat = \"PBO\") as the control treatment.\n\nNMA.uc <- netmeta(TE = TE, seTE = seTE, treat1 = treat1, treat2 = treat2,\n studlab = studlab, data = UlcerativeColitis, sm = \"OR\", \n ref = \"PBO\", common = FALSE, comb.random = TRUE)\nNMA.uc\n\nAll treatments except FIL100 and UST6 are significantly more efficacious than PBO at inducing clinical remission. We can now estimate the probabilities of each treatment being at each possible rank and the SUCRAs (Surface Under the Cumulative RAnking curve):\n\nsucra.uc <- rankogram(NMA.uc, nsim = 100, random = TRUE, common = FALSE, \n small.values = \"undesirable\")\n\n# Exctract the SUCRA values\nsucra.uc$ranking.random\n\n ADA160/80 FIL100 FIL200 GOL200/100 INF10 INF5 OZA0.92 \n0.26090909 0.18909091 0.40363636 0.59545455 0.60909091 0.74090909 0.77727273 \n PBO TOF10 UPA45 UST6 VED300 \n0.02090909 0.40636364 0.98727273 0.40727273 0.60181818 \n\n\nThese results indicate that 98.7% of the evaluated treatments are worse than UPA45." }, { "objectID": "chapter_10.html#version-info", @@ -256,28 +256,42 @@ "href": "chapter_11.html#hierarchical-meta-regression", "title": "7  Individual Participant Data Meta-analysis of clinical trials and real-world data", "section": "7.2 Hierarchical Meta-Regression", - "text": "7.2 Hierarchical Meta-Regression\nWe illustrate the implementation of hierarchical meta-regression using an example that involves the following data sources:\n\nAggregate data from 35 randomized trials investigating the efficacy of adjunctive treatments in managing diabetic foot problems compared with routine care\nIndividual participant data from a prospective cohort study investigating patient and limb survival in patients with diabetic foot ulcers\n\n\n7.2.1 Aggregate data\nWe first retrieve the randomized evidence and summarize the treatment effect estimates using a random effects meta-analysis:\n\nlibrary(dplyr)\nlibrary(jarbes)\nlibrary(meta)\n\ndata(\"healing\")\n\nresults.ADJ <- metabin(event.c = y_c, n.c = n_c,\n event.e = y_t, n.e = n_t,\n studlab = Study, data = healing,\n sm = \"OR\", \n prediction = TRUE)\n\nThe corresponding forest plot is depicted below. The endpoint is healing without amputations within a period less than or equal to 1 year.\n\n\n\n\n\nThe random effects meta-analysis yielded a pooled odds ratio of 1.90. However, substantial between-study heterogeneity was found, with \\(\\tau\\) = 0.46.\n\n\n7.2.2 Individual participant data\nSubsequently, we retrieve the individual participant data:\n\ndata(\"healingipd\")\nIPD <- healingipd %>% dplyr::select(healing.without.amp, PAD, neuropathy,\n first.ever.lesion, no.continuous.care, \n male, diab.typ2, insulin, HOCHD, \n HOS, CRF, dialysis, DNOAP, smoking.ever, \n diabdur, wagner.class)\n\nBriefly, these IPD were obtained from a prospective cohort study enrolling consecutive patients with diabetic foot ulcers (DFUs) and without previous major amputation in a single diabetes center between June 1998 and December 1999 (Morbach et al. 2012). The baseline characteristics of the study population is summarized below:\n\n\n\n\n\n\n\n\n\n\n\n\nOverall\n(N=260)\n\n\n\n\nAge (years)\n\n\n\nMean (SD)\n68.9 (10.9)\n\n\nMedian [Min, Max]\n70.0 [25.0, 91.0]\n\n\nDiabetes duration (years)\n\n\n\nMean (SD)\n15.9 (10.6)\n\n\nMedian [Min, Max]\n14.0 [0, 53.0]\n\n\nSex\n\n\n\nFemale\n106 (40.8%)\n\n\nMale\n154 (59.2%)\n\n\nEver smoker\n\n\n\nYes\n154 (59.2%)\n\n\nNo\n106 (40.8%)\n\n\nDiabetes type 2\n\n\n\nYes\n229 (88.1%)\n\n\nNo\n31 (11.9%)\n\n\nPeripheral arterial disease\n\n\n\nYes\n148 (56.9%)\n\n\nNo\n112 (43.1%)\n\n\nNeuropathy\n\n\n\nYes\n224 (86.2%)\n\n\nNo\n36 (13.8%)\n\n\nFirst ever lesion\n\n\n\nYes\n114 (43.8%)\n\n\nNo\n146 (56.2%)\n\n\nNo continuous care\n\n\n\nYes\n177 (68.1%)\n\n\nNo\n83 (31.9%)\n\n\nInsulin dependent\n\n\n\nYes\n174 (66.9%)\n\n\nNo\n86 (33.1%)\n\n\nHistory of coronary events (CHD)\n\n\n\nYes\n52 (20.0%)\n\n\nNo\n208 (80.0%)\n\n\nHistory of stroke\n\n\n\nYes\n55 (21.2%)\n\n\nNo\n205 (78.8%)\n\n\nCharcot foot syndrome\n\n\n\nYes\n52 (20.0%)\n\n\nNo\n208 (80.0%)\n\n\nDialysis\n\n\n\nYes\n9 (3.5%)\n\n\nNo\n251 (96.5%)\n\n\nDNOAP\n\n\n\nYes\n29 (11.2%)\n\n\nNo\n231 (88.8%)\n\n\nWagner score\n\n\n\n1-2\n142 (54.6%)\n\n\n3-4-5\n118 (45.4%)\n\n\n\n\n\n\n\nAs depicted above, IPD are available from 260 patients. Some of these patients have similar characteristics to those enrolled in the randomized trials. However, other patients have comorbidities, where one or more risk factors prevent them to participate in the RCTs due to ethical reasons. For example, 118 patients have severe ulcer lesions (Wagner score 3 to 4), and 77 patients suffer from severe ulcer lesions and peripheral arterial disease (PAD). The question is: Can we generalize the benefit of adjuvant therapies observed in the RCTs to the subgroups of patients encountered in clinical practice?\n\n\n7.2.3 Hierarchical metaregression\nWe fitted an HMR model to the available RWD and published AD:\n\nset.seed(2022)\n\nAD <- healing %>% dplyr::select(yc = y_c, nc = n_c, \n yt = y_t, nt = n_t, Study = Study)\n\nmx2 <- hmr(data = AD, # Published aggregate data\n two.by.two = FALSE, # \n dataIPD = IPD, # Data frame of the IPD \n re = \"sm\", # Random effects model: \"sm\" scale mixtures \n link = \"logit\", # Link function of the random effects\n sd.mu.1 = 1, # Scale parameter for the prior of mu.1\n sd.mu.2 = 1, # Scale parameter for the prior of mu.2 \n sd.mu.phi = 1, # Scale parameter for the prior of mu.phi \n sigma.1.upper = 5, # Upper bound of the prior of sigma.1 \n sigma.2.upper = 5, # Upper bound of the prior of sigma.2\n sigma.beta.upper = 5, # Upper bound of the prior of sigma.beta\n sd.Fisher.rho = 1.25, # Scale parameter for the prior of rho\n df.estimate = TRUE, # If TRUE the degrees of freedom are estimated\n df.lower = 3, # Lower bound of the df's prior\n df.upper = 10, # Upper bound of the df's prior\n nr.chains = 2, # Number of MCMC chains\n nr.iterations = 10000, # Total number of iterations\n nr.adapt = 1000, # Number of iteration for burnin \n nr.thin = 1) # Thinning rate\n\nWe start our analysis by visualizing the conflict of evidence between the different types of data and study types. The figure below depicts the posterior distribution of \\(\\mu_{\\phi}\\), which is the mean bias of the IPD-NRS compared to the AD-RCTs control groups. The posterior distribution has a substantial probability mass on the right of zero, which indicates that in average the IPD-NRS patients present a better prognoses than the AD-RCTs control groups. That means that taking the IPD-NRS results at face value would be misleading if we aim to combine them with a meta-analysis of AD-RCTs.\n\n\n\n\n\nConflict of evidence analysis. The left panel shows the prior to posterior sensitivity analysis of bias mean between the RCTs and the IPD-NRS. The right panel depicts the posterior distribution of the outliers detection weights.\n\n\n\n\nThe figure below presents the posterior distribution of the weights \\(w_{i}\\) for each study included in the HMR. These posteriors are summarized using a forest plot, where posterior intervals substantially greater than one indicate outliers. One important aspect of the HMR is that those outliers are automatically down-weighted in the analysis.\n\n\n\n\n\nPosterior distribution of the weights for each study included in the HMR\n\n\n\n\n\n\n\n\nMorbach, Stephan, Heike Furchert, Ute Gröblinghoff, Heribert Hoffmeier, Kerstin Kersten, Gerd-Thomas Klauke, Ulrike Klemp, et al. 2012. “Long-Term Prognosis of Diabetic Foot Patients and Their Limbs.” Diabetes Care 35 (10): 2021–27. https://doi.org/10.2337/dc12-0200." + "text": "7.2 Hierarchical Meta-Regression\nWe illustrate the implementation of hierarchical meta-regression using an example that involves the following data sources:\n\nAggregate data from 35 randomized trials investigating the efficacy of adjunctive treatments in managing diabetic foot problems compared with routine care\nIndividual participant data from a prospective cohort study investigating patient and limb survival in patients with diabetic foot ulcers\n\n\n7.2.1 Aggregate data\nWe first retrieve the randomized evidence and summarize the treatment effect estimates using a random effects meta-analysis:\n\nlibrary(dplyr)\nlibrary(jarbes)\nlibrary(meta)\n\ndata(\"healing\")\n\nresults.ADJ <- metabin(event.c = y_c, n.c = n_c,\n event.e = y_t, n.e = n_t,\n studlab = Study, data = healing,\n sm = \"OR\", \n prediction = TRUE)\n\nThe corresponding forest plot is depicted below. The endpoint is healing without amputations within a period less than or equal to 1 year.\n\n\n\n\n\nThe random effects meta-analysis yielded a pooled odds ratio of 1.90. However, substantial between-study heterogeneity was found, with \\(\\tau\\) = 0.46.\n\n\n7.2.2 Individual participant data\nSubsequently, we retrieve the individual participant data:\n\ndata(\"healingipd\")\nIPD <- healingipd %>% dplyr::select(healing.without.amp, PAD, neuropathy,\n first.ever.lesion, no.continuous.care, \n male, diab.typ2, insulin, HOCHD, \n HOS, CRF, dialysis, DNOAP, smoking.ever, \n diabdur, wagner.class)\n\nBriefly, these IPD were obtained from a prospective cohort study enrolling consecutive patients with diabetic foot ulcers (DFUs) and without previous major amputation in a single diabetes center between June 1998 and December 1999 (Morbach et al. 2012). The baseline characteristics of the study population is summarized below:\n\n\n\n\n\n\n\n\n\n\n\n\nOverall\n(N=260)\n\n\n\n\nAge (years)\n\n\n\nMean (SD)\n68.9 (10.9)\n\n\nMedian [Min, Max]\n70.0 [25.0, 91.0]\n\n\nDiabetes duration (years)\n\n\n\nMean (SD)\n15.9 (10.6)\n\n\nMedian [Min, Max]\n14.0 [0, 53.0]\n\n\nSex\n\n\n\nFemale\n106 (40.8%)\n\n\nMale\n154 (59.2%)\n\n\nEver smoker\n\n\n\nYes\n154 (59.2%)\n\n\nNo\n106 (40.8%)\n\n\nDiabetes type 2\n\n\n\nYes\n229 (88.1%)\n\n\nNo\n31 (11.9%)\n\n\nPeripheral arterial disease\n\n\n\nYes\n148 (56.9%)\n\n\nNo\n112 (43.1%)\n\n\nNeuropathy\n\n\n\nYes\n224 (86.2%)\n\n\nNo\n36 (13.8%)\n\n\nFirst ever lesion\n\n\n\nYes\n114 (43.8%)\n\n\nNo\n146 (56.2%)\n\n\nNo continuous care\n\n\n\nYes\n177 (68.1%)\n\n\nNo\n83 (31.9%)\n\n\nInsulin dependent\n\n\n\nYes\n174 (66.9%)\n\n\nNo\n86 (33.1%)\n\n\nHistory of coronary events (CHD)\n\n\n\nYes\n52 (20.0%)\n\n\nNo\n208 (80.0%)\n\n\nHistory of stroke\n\n\n\nYes\n55 (21.2%)\n\n\nNo\n205 (78.8%)\n\n\nCharcot foot syndrome\n\n\n\nYes\n52 (20.0%)\n\n\nNo\n208 (80.0%)\n\n\nDialysis\n\n\n\nYes\n9 (3.5%)\n\n\nNo\n251 (96.5%)\n\n\nDNOAP\n\n\n\nYes\n29 (11.2%)\n\n\nNo\n231 (88.8%)\n\n\nWagner score\n\n\n\n1-2\n142 (54.6%)\n\n\n3-4-5\n118 (45.4%)\n\n\n\n\n\n\n\nAs depicted above, IPD are available from 260 patients. Some of these patients have similar characteristics to those enrolled in the randomized trials. However, other patients have comorbidities, where one or more risk factors prevent them to participate in the RCTs due to ethical reasons. For example, 118 patients have severe ulcer lesions (Wagner score 3 to 4), and 77 patients suffer from severe ulcer lesions and peripheral arterial disease (PAD). The question is: Can we generalize the benefit of adjuvant therapies observed in the RCTs to the subgroups of patients encountered in clinical practice?\n\n\n7.2.3 Hierarchical metaregression\nWe fitted an HMR model to the available RWD and published AD:\n\nset.seed(2022)\n\nAD <- healing %>% dplyr::select(yc = y_c, nc = n_c, \n yt = y_t, nt = n_t, Study = Study)\n\nmx2 <- hmr(data = AD, # Published aggregate data\n two.by.two = FALSE, # \n dataIPD = IPD, # Data frame of the IPD \n re = \"sm\", # Random effects model: \"sm\" scale mixtures \n link = \"logit\", # Link function of the random effects\n sd.mu.1 = 1, # Scale parameter for the prior of mu.1\n sd.mu.2 = 1, # Scale parameter for the prior of mu.2 \n sd.mu.phi = 1, # Scale parameter for the prior of mu.phi \n sigma.1.upper = 5, # Upper bound of the prior of sigma.1 \n sigma.2.upper = 5, # Upper bound of the prior of sigma.2\n sigma.beta.upper = 5, # Upper bound of the prior of sigma.beta\n sd.Fisher.rho = 1.25, # Scale parameter for the prior of rho\n df.estimate = TRUE, # If TRUE the degrees of freedom are estimated\n df.lower = 3, # Lower bound of the df's prior\n df.upper = 10, # Upper bound of the df's prior\n nr.chains = 2, # Number of MCMC chains\n nr.iterations = 10000, # Total number of iterations\n nr.adapt = 1000, # Number of iteration for burnin \n nr.thin = 1) # Thinning rate\n\nWe start our analysis by visualizing the conflict of evidence between the different types of data and study types. The figure below depicts the posterior distribution of \\(\\mu_{\\phi}\\), which is the mean bias of the IPD-NRS compared to the AD-RCTs control groups. The posterior distribution has a substantial probability mass on the right of zero, which indicates that in average the IPD-NRS patients present a better prognoses than the AD-RCTs control groups. That means that taking the IPD-NRS results at face value would be misleading if we aim to combine them with a meta-analysis of AD-RCTs.\n\n\n\n\n\nConflict of evidence analysis. The left panel shows the prior to posterior sensitivity analysis of bias mean between the RCTs and the IPD-NRS. The right panel depicts the posterior distribution of the outliers detection weights.\n\n\n\n\nThe figure below presents the posterior distribution of the weights \\(w_{i}\\) for each study included in the HMR. These posteriors are summarized using a forest plot, where posterior intervals substantially greater than one indicate outliers. One important aspect of the HMR is that those outliers are automatically down-weighted in the analysis.\n\n\n\n\n\nPosterior distribution of the weights for each study included in the HMR" + }, + { + "objectID": "chapter_11.html#version-info", + "href": "chapter_11.html#version-info", + "title": "7  Individual Participant Data Meta-analysis of clinical trials and real-world data", + "section": "Version info", + "text": "Version info\nThis chapter was rendered using the following version of R and its packages:\n\n\nR version 4.2.3 (2023-03-15)\nPlatform: x86_64-pc-linux-gnu (64-bit)\nRunning under: Ubuntu 22.04.3 LTS\n\nMatrix products: default\nBLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3\nLAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so\n\nlocale:\n [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 \n [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8 \n [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C \n[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C \n\nattached base packages:\n[1] stats graphics grDevices utils datasets methods base \n\nother attached packages:\n[1] ggplot2_3.4.4 meta_6.5-0 jarbes_2.0.0 dplyr_1.1.4 \n[5] table1_1.4.3 kableExtra_1.3.4\n\nloaded via a namespace (and not attached):\n [1] httr_1.4.7 tidyr_1.3.0 sfsmisc_1.1-16 \n [4] jsonlite_1.8.7 viridisLite_0.4.2 splines_4.2.3 \n [7] Formula_1.2-5 shiny_1.8.0 metafor_4.4-0 \n[10] yaml_2.3.7 numDeriv_2016.8-1.1 R2WinBUGS_2.1-21 \n[13] pillar_1.9.0 lattice_0.20-45 glue_1.6.2 \n[16] digest_0.6.33 RColorBrewer_1.1-3 promises_1.2.1 \n[19] minqa_1.2.6 rvest_1.0.3 colorspace_2.1-0 \n[22] R2jags_0.7-1 Matrix_1.6-3 mcmcplots_0.4.3 \n[25] htmltools_0.5.7 httpuv_1.6.12 plyr_1.8.9 \n[28] pkgconfig_2.0.3 purrr_1.0.2 xtable_1.8-4 \n[31] scales_1.2.1 webshot_0.5.5 svglite_2.1.2 \n[34] rjags_4-14 later_1.3.1 ggstats_0.5.1 \n[37] metadat_1.2-0 lme4_1.1-35.1 tibble_3.2.1 \n[40] farver_2.1.1 generics_0.1.3 ellipsis_0.3.2 \n[43] withr_2.5.2 cli_3.6.1 magrittr_2.0.3 \n[46] mime_0.12 evaluate_0.23 GGally_2.2.0 \n[49] fansi_1.0.5 nlme_3.1-162 MASS_7.3-58.2 \n[52] xml2_1.3.5 tools_4.2.3 lifecycle_1.0.4 \n[55] stringr_1.5.1 munsell_0.5.0 compiler_4.2.3 \n[58] systemfonts_1.0.5 rlang_1.1.2 nloptr_2.0.3 \n[61] grid_4.2.3 rstudioapi_0.15.0 CompQuadForm_1.4.3 \n[64] htmlwidgets_1.6.3 miniUI_0.1.1.1 labeling_0.4.3 \n[67] rmarkdown_2.25 boot_1.3-28.1 gtable_0.3.4 \n[70] codetools_0.2-19 abind_1.4-5 R6_2.5.1 \n[73] gridExtra_2.3 knitr_1.45 denstrip_1.5.4 \n[76] fastmap_1.1.1 utf8_1.2.4 mathjaxr_1.6-0 \n[79] ggExtra_0.10.1 stringi_1.8.2 parallel_4.2.3 \n[82] Rcpp_1.0.11 vctrs_0.6.4 tidyselect_1.2.0 \n[85] xfun_0.41 coda_0.19-4" + }, + { + "objectID": "chapter_11.html#references", + "href": "chapter_11.html#references", + "title": "7  Individual Participant Data Meta-analysis of clinical trials and real-world data", + "section": "References", + "text": "References\n\n\n\n\nMorbach, Stephan, Heike Furchert, Ute Gröblinghoff, Heribert Hoffmeier, Kerstin Kersten, Gerd-Thomas Klauke, Ulrike Klemp, et al. 2012. “Long-Term Prognosis of Diabetic Foot Patients and Their Limbs.” Diabetes Care 35 (10): 2021–27. https://doi.org/10.2337/dc12-0200." }, { "objectID": "chapter_12.html#introduction", "href": "chapter_12.html#introduction", "title": "8  Dealing with irregular and informative visits", "section": "8.1 Introduction", - "text": "8.1 Introduction\nWe first load the relevant R scripts:\n\nsource(\"resources/chapter 12/sim.r\")\n\n\nAttaching package: 'dplyr'\n\n\nThe following objects are masked from 'package:stats':\n\n filter, lag\n\n\nThe following objects are masked from 'package:base':\n\n intersect, setdiff, setequal, union\n\n\n\nAttaching package: 'mice'\n\n\nThe following object is masked from 'package:stats':\n\n filter\n\n\nThe following objects are masked from 'package:base':\n\n cbind, rbind\n\n\n\nAttaching package: 'nlme'\n\n\nThe following object is masked from 'package:dplyr':\n\n collapse\n\n\n\nAttaching package: 'MASS'\n\n\nThe following object is masked from 'package:dplyr':\n\n select\n\nsource(\"resources/chapter 12/fig_functions.r\")\nsource(\"resources/chapter 12/mlmi.r\")" + "text": "8.1 Introduction\nWe first load the relevant R scripts:\n\nsource(\"resources/chapter 12/sim.r\")\nsource(\"resources/chapter 12/fig_functions.r\")\nsource(\"resources/chapter 12/mlmi.r\")" }, { "objectID": "chapter_12.html#example-dataset", "href": "chapter_12.html#example-dataset", "title": "8  Dealing with irregular and informative visits", "section": "8.2 Example dataset", - "text": "8.2 Example dataset\nBelow, we generate an example dataset that contains information on the treatment allocation x and three baseline covariates age, sex and edss (EDSS at treatment start). The discrete outcome y represents the Expanded Disability Status Scale (EDSS) score after time months of treatment exposure. Briefly, the EDSS is a semi-continuous measure that varies from 0 (no disability) to 10 (death).\n\nset.seed(9843626)\n\ndataset <- sim_data_EDSS(npatients = 500,\n ncenters = 10,\n follow_up = 12*5, # Total follow-up (number of months)\n sd_a_t = 0.5, # DGM - Within-visit variation in EDSS scores\n baseline_EDSS = 1.3295, # DGM - Mean baseline EDDS score\n sd_alpha_ij = 1.46, # DGM - Between-subject variation in baseline EDSS\n sd_beta1_j = 0.20, # DGM - Between-site variation in baseline EDSS\n mean_age = 42.41,\n sd_age = 10.53,\n min_age = 18,\n beta_age = 0.05, # DGM - prognostic effect of age\n beta_t = 0.014, # DGM - prognostic effect of time\n beta_t2 = 0, # DGM - prognostic effect of time squared\n delta_xt = 0, # DGM - interaction treatment time\n delta_xt2 = 0, # 0.0005 # DGM - interaction treatment time2\n p_female = 0.75, \n beta_female = -0.2 , ## DGM - prognostic effect of male sex\n delta_xf = 0, ## DGM - interaction sex treatment \n rho = 0.8, # DGM - autocorrelation of between alpha_tij\n corFUN = corAR1, # DGM - correlation structure of the latent EDSS scores\n tx_alloc_FUN = treatment_alloc_confounding_v2 ) ## or treatment_alloc_randomized\n\n\n\n\n\n\nDistribution of the EDSS score at each time point\n\n\n\n\nWe remove the outcome y according to the informative visit process that depends on the received treatment, gender, and age.\n\ndataset_visit <- censor_visits_a5(dataset, seed = 12345) %>% \n dplyr::select(-y) %>%\n mutate(time_x = time*x)\n\nIn the censored data, a total of 17 out of 5000 patients have a visit at time=60." + "text": "8.2 Example dataset\nBelow, we generate an example dataset that contains information on the treatment allocation x and three baseline covariates age, sex and edss (EDSS at treatment start). The discrete outcome y represents the Expanded Disability Status Scale (EDSS) score after time months of treatment exposure. Briefly, the EDSS is a semi-continuous measure that varies from 0 (no disability) to 10 (death).\n\nset.seed(9843626)\n\ndataset <- sim_data_EDSS(npatients = 500,\n ncenters = 10,\n follow_up = 12*5,\n sd_a_t = 0.5, \n baseline_EDSS = 1.3295, \n sd_alpha_ij = 1.46, \n sd_beta1_j = 0.20, \n mean_age = 42.41,\n sd_age = 10.53,\n min_age = 18,\n beta_age = 0.05, \n beta_t = 0.014, \n beta_t2 = 0, \n delta_xt = 0,\n delta_xt2 = 0,\n p_female = 0.75, \n beta_female = -0.2 ,\n delta_xf = 0, \n rho = 0.8, \n corFUN = corAR1, \n tx_alloc_FUN = treatment_alloc_confounding_v2)\n\n\n\n`summarise()` has grouped output by 'time'. You can override using the\n`.groups` argument.\n\n\n\n\n\nLongitudinal distribution of the EDSS score. The shaded area depicts the interquartile range.\n\n\n\n\nWe remove the outcome y according to the informative visit process that depends on the received treatment, gender, and age.\n\ndataset_visit <- censor_visits_a5(dataset, seed = 12345) %>% \n dplyr::select(-y) %>%\n mutate(time_x = time*x)\n\nIn the censored data, a total of 17 out of 5000 patients have a visit at time=60." }, { "objectID": "chapter_12.html#estimation-of-treatment-effect", "href": "chapter_12.html#estimation-of-treatment-effect", "title": "8  Dealing with irregular and informative visits", "section": "8.3 Estimation of treatment effect", - "text": "8.3 Estimation of treatment effect\nWe will estimate the marginal treatment effect at time time=60.\n\n8.3.1 Original data\n\norigdat60 <- dataset %>% filter(time == 60)\n\n# Predict probability of treatment allocation\nfitps <- glm(x ~ age + sex + edss, family = 'binomial', \n data = origdat60)\n\n# Derive the propensity score\norigdat60 <- origdat60 %>% mutate(ipt = ifelse(x == 1, 1/predict(fitps, type = 'response'),\n 1/(1-predict(fitps, type = 'response'))))\n\n# Estimate \nfit_ref_m <- tidy(lm(y ~ x, weight = ipt, data = origdat60), conf.int = TRUE) \n\n\n\n8.3.2 Doubly-weighted marginal treatment effect\nWe here implement inverse probability of response weights into the estimating equations to adjust for nonrandom missingness Coulombe, Moodie, and Platt (2020).\n\nobsdat60 <- dataset_visit %>% mutate(visit = ifelse(is.na(y_obs),0,1)) %>% filter(time == 60)\n\ngamma <- glm(visit ~ x + sex + age + edss, family = 'binomial', data = obsdat60)$coef \n\nobsdat60 <- obsdat60 %>% mutate(rho_i = 1/exp(gamma[\"(Intercept)\"] +\n gamma[\"x\"]*x +\n gamma[\"sex\"]*sex +\n gamma[\"age\"]*age))\n\n# Predict probability of treatment allocation\nfitps <- glm(x ~ age + sex + edss, family='binomial', data = obsdat60)\n\n# Derive the propensity score\nobsdat60 <- obsdat60 %>% mutate(ipt = ifelse(x==1, 1/predict(fitps, type='response'),\n 1/(1-predict(fitps, type='response'))))\n\n\nfit_w <- tidy(lm(y_obs ~ x, weights = ipt*rho_i, data = obsdat60), conf.int = TRUE)\n\n\n\n8.3.3 Multilevel multiple imputation\nWe adopt the imputation approach proposed by Debray et al. (2023). Briefly, we impute the entire vector of y_obs for all 61 potential visits and generate 10 imputed datasets. Note: mlmi currently does not support imputation of treatment-covariate interaction terms.\n\nimp <- impute_y_mice_3l(dataset_visit, seed = 12345)\n\nWe can now estimate the treatment effect in each imputed dataset\n\n# Predict probability of treatment allocation\nfitps <- glm(x ~ age + sex + edss, family='binomial', data = dataset_visit)\n \n# Derive the propensity score\ndataset_visit <- dataset_visit %>% mutate(ipt = ifelse(x==1, 1/predict(fitps, type='response'),\n 1/(1-predict(fitps, type='response'))))\n \nQ <- U <- rep(NA, 10) # Error variances\n\nfor (i in seq(10)) {\n dati <- cbind(dataset_visit[,c(\"x\",\"ipt\",\"time\")], y_imp = imp[,i]) %>% filter(time == 60)\n \n # Estimate \n fit <- tidy(lm(y_imp ~ x, weight = ipt, data = dati), conf.int = TRUE) \n \n Q[i] <- fit %>% filter(term == \"x\") %>% pull(estimate)\n U[i] <- (fit %>% filter(term == \"x\") %>% pull(std.error))**2\n}\n\nfit_mlmi <- pool.scalar(Q = Q, U = U)" + "text": "8.3 Estimation of treatment effect\nWe will estimate the marginal treatment effect at time time=60.\n\n8.3.1 Original data\n\norigdat60 <- dataset %>% filter(time == 60)\n\n# Predict probability of treatment allocation\nfitps <- glm(x ~ age + sex + edss, family = 'binomial', \n data = origdat60)\n\n# Derive the propensity score\norigdat60 <- origdat60 %>% \n mutate(ipt = ifelse(x == 1, 1/predict(fitps, type = 'response'),\n 1/(1 - predict(fitps, type = 'response'))))\n\n# Estimate \nfit_ref_m <- tidy(lm(y ~ x, weight = ipt, data = origdat60), conf.int = TRUE) \n\n\n\n8.3.2 Doubly-weighted marginal treatment effect\nWe here implement inverse probability of response weights into the estimating equations to adjust for nonrandom missingness Coulombe, Moodie, and Platt (2020).\n\nobsdat60 <- dataset_visit %>% \n mutate(visit = ifelse(is.na(y_obs),0,1)) %>% \n filter(time == 60)\n\ngamma <- glm(visit ~ x + sex + age + edss, \n family = 'binomial', data = obsdat60)$coef \n\nobsdat60 <- obsdat60 %>% mutate(rho_i = 1/exp(gamma[\"(Intercept)\"] +\n gamma[\"x\"]*x +\n gamma[\"sex\"]*sex +\n gamma[\"age\"]*age))\n\n# Predict probability of treatment allocation\nfitps <- glm(x ~ age + sex + edss, family='binomial', data = obsdat60)\n\n# Derive the propensity score\nobsdat60 <- obsdat60 %>% \n mutate(ipt = ifelse(x==1, 1/predict(fitps, type='response'),\n 1/(1 - predict(fitps, type='response'))))\n\n\nfit_w <- tidy(lm(y_obs ~ x, weights = ipt*rho_i, data = obsdat60), \n conf.int = TRUE)\n\n\n\n8.3.3 Multilevel multiple imputation\nWe adopt the imputation approach proposed by Debray et al. (2023). Briefly, we impute the entire vector of y_obs for all 61 potential visits and generate 10 imputed datasets. Note: mlmi currently does not support imputation of treatment-covariate interaction terms.\n\nimp <- impute_y_mice_3l(dataset_visit, seed = 12345)\n\nWe can now estimate the treatment effect in each imputed dataset\n\n# Predict probability of treatment allocation\nfitps <- glm(x ~ age + sex + edss, family = 'binomial', data = dataset_visit)\n \n# Derive the propensity score\ndataset_visit <- dataset_visit %>% \n mutate(ipt = ifelse(x == 1, 1/predict(fitps, type = 'response'),\n 1/(1 - predict(fitps, type = 'response'))))\n \nQ <- U <- rep(NA, 10) # Error variances\n\nfor (i in seq(10)) {\n dati <- cbind(dataset_visit[,c(\"x\",\"ipt\",\"time\")], y_imp = imp[,i]) %>% \n filter(time == 60)\n \n # Estimate \n fit <- tidy(lm(y_imp ~ x, weight = ipt, data = dati), conf.int = TRUE) \n \n Q[i] <- fit %>% filter(term == \"x\") %>% pull(estimate)\n U[i] <- (fit %>% filter(term == \"x\") %>% pull(std.error))**2\n}\n\nfit_mlmi <- pool.scalar(Q = Q, U = U)" }, { "objectID": "chapter_12.html#reproduce-the-results-using-all-data-to-compute-the-marginal-effect-with-iiv-weighted", @@ -298,7 +312,7 @@ "href": "chapter_12.html#version-info", "title": "8  Dealing with irregular and informative visits", "section": "Version info", - "text": "Version info\nThis chapter was rendered using the following version of R and its packages:\n\n\nR version 4.2.3 (2023-03-15)\nPlatform: x86_64-pc-linux-gnu (64-bit)\nRunning under: Ubuntu 22.04.3 LTS\n\nMatrix products: default\nBLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3\nLAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so\n\nlocale:\n [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 \n [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8 \n [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C \n[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C \n\nattached base packages:\n[1] stats graphics grDevices utils datasets methods base \n\nother attached packages:\n[1] sparseMVN_0.2.2 truncnorm_1.0-9 MASS_7.3-58.2 nlme_3.1-162 \n[5] mice_3.16.0 ggplot2_3.4.4 broom_1.0.5 dplyr_1.1.4 \n\nloaded via a namespace (and not attached):\n [1] shape_1.4.6 tidyselect_1.2.0 xfun_0.41 purrr_1.0.2 \n [5] splines_4.2.3 lattice_0.20-45 colorspace_2.1-0 vctrs_0.6.4 \n [9] generics_0.1.3 htmltools_0.5.7 yaml_2.3.7 pan_1.9 \n[13] utf8_1.2.4 survival_3.5-3 rlang_1.1.2 jomo_2.7-6 \n[17] pillar_1.9.0 nloptr_2.0.3 withr_2.5.2 glue_1.6.2 \n[21] RColorBrewer_1.1-3 foreach_1.5.2 lifecycle_1.0.4 munsell_0.5.0 \n[25] gtable_0.3.4 htmlwidgets_1.6.3 codetools_0.2-19 evaluate_0.23 \n[29] labeling_0.4.3 knitr_1.45 fastmap_1.1.1 fansi_1.0.5 \n[33] Rcpp_1.0.11 scales_1.2.1 backports_1.4.1 jsonlite_1.8.7 \n[37] farver_2.1.1 lme4_1.1-35.1 digest_0.6.33 grid_4.2.3 \n[41] cli_3.6.1 tools_4.2.3 magrittr_2.0.3 glmnet_4.1-8 \n[45] tibble_3.2.1 tidyr_1.3.0 pkgconfig_2.0.3 ellipsis_0.3.2 \n[49] Matrix_1.6-3 minqa_1.2.6 rmarkdown_2.25 iterators_1.0.14 \n[53] mitml_0.4-5 R6_2.5.1 boot_1.3-28.1 rpart_4.1.19 \n[57] nnet_7.3-18 compiler_4.2.3" + "text": "Version info\nThis chapter was rendered using the following version of R and its packages:\n\n\nR version 4.2.3 (2023-03-15)\nPlatform: x86_64-pc-linux-gnu (64-bit)\nRunning under: Ubuntu 22.04.3 LTS\n\nMatrix products: default\nBLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3\nLAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so\n\nlocale:\n [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 \n [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8 \n [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C \n[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C \n\nattached base packages:\n[1] stats graphics grDevices utils datasets methods base \n\nother attached packages:\n[1] sparseMVN_0.2.2 truncnorm_1.0-9 MASS_7.3-58.2 nlme_3.1-162 \n[5] mice_3.16.0 ggplot2_3.4.4 broom_1.0.5 dplyr_1.1.4 \n\nloaded via a namespace (and not attached):\n [1] shape_1.4.6 tidyselect_1.2.0 xfun_0.41 purrr_1.0.2 \n [5] splines_4.2.3 lattice_0.20-45 colorspace_2.1-0 vctrs_0.6.4 \n [9] generics_0.1.3 htmltools_0.5.7 yaml_2.3.7 pan_1.9 \n[13] utf8_1.2.4 survival_3.5-3 rlang_1.1.2 jomo_2.7-6 \n[17] pillar_1.9.0 nloptr_2.0.3 withr_2.5.2 glue_1.6.2 \n[21] foreach_1.5.2 lifecycle_1.0.4 munsell_0.5.0 gtable_0.3.4 \n[25] htmlwidgets_1.6.3 codetools_0.2-19 evaluate_0.23 labeling_0.4.3 \n[29] knitr_1.45 fastmap_1.1.1 fansi_1.0.5 Rcpp_1.0.11 \n[33] scales_1.2.1 backports_1.4.1 jsonlite_1.8.7 farver_2.1.1 \n[37] lme4_1.1-35.1 digest_0.6.33 grid_4.2.3 cli_3.6.1 \n[41] tools_4.2.3 magrittr_2.0.3 glmnet_4.1-8 tibble_3.2.1 \n[45] tidyr_1.3.0 pkgconfig_2.0.3 ellipsis_0.3.2 Matrix_1.6-3 \n[49] minqa_1.2.6 rmarkdown_2.25 iterators_1.0.14 mitml_0.4-5 \n[53] R6_2.5.1 boot_1.3-28.1 rpart_4.1.19 nnet_7.3-18 \n[57] compiler_4.2.3" }, { "objectID": "chapter_12.html#references", @@ -312,14 +326,14 @@ "href": "chapter_16.html#estimating-heterogeneous-treatment-effects-in-pairwise-meta-analysis", "title": "9  Prediction of individual treatment effect using data from multiple studies", "section": "9.1 Estimating heterogeneous treatment effects in pairwise meta-analysis", - "text": "9.1 Estimating heterogeneous treatment effects in pairwise meta-analysis\nWe hereby provide code for estimating patient-level treatment effects for the case when we have patient-level data from multiple randomized trials.\n\n9.1.1 Example of a continuous outcome\n\n9.1.1.1 Setup\nWe start by simulating an artificial dataset using the R package bipd:\n\nlibrary(bipd)\nds <- generate_ipdma_example(type = \"continuous\")\n\nLet us have a look at the dataset:\n\nhead(ds)\n\n studyid treat z1 z2 y\n1 1 0 0.24140876 0.8999285 11\n2 1 1 -0.08337922 -1.2523967 4\n3 1 0 0.54133588 1.1391411 11\n4 1 0 -0.72552223 0.8181371 11\n5 1 0 -1.34843084 0.5481249 11\n6 1 1 0.57662307 -0.5316352 4\n\n\nThe simulated dataset contains information on the following variables:\n\nthe trial indicator studyid\nthe treatment indicator treat, which takes the values 0 for control and 1 for active treatment\ntwo prognostic variables z1 and z2\nthe continuous outcome y\n\n\n\n\n\n\n\nTable 9.1: The simulated dataset with a continuous outcome\n\n\n\n\n\n\n\n\n\n0\n(N=301)\n1\n(N=299)\nOverall\n(N=600)\n\n\n\n\nz1\n\n\n\n\n\nMean (SD)\n-0.0562 (0.986)\n0.0454 (0.890)\n-0.00558 (0.940)\n\n\nMedian [Min, Max]\n-0.0694 [-2.72, 2.39]\n0.0749 [-2.57, 2.10]\n0.0164 [-2.72, 2.39]\n\n\nz2\n\n\n\n\n\nMean (SD)\n-0.0214 (0.959)\n-0.0183 (0.953)\n-0.0198 (0.955)\n\n\nMedian [Min, Max]\n-0.164 [-2.53, 2.59]\n-0.0303 [-2.36, 3.00]\n-0.0852 [-2.53, 3.00]\n\n\nstudyid\n\n\n\n\n\n1\n41 (13.6%)\n59 (19.7%)\n100 (16.7%)\n\n\n2\n52 (17.3%)\n48 (16.1%)\n100 (16.7%)\n\n\n3\n53 (17.6%)\n47 (15.7%)\n100 (16.7%)\n\n\n4\n54 (17.9%)\n46 (15.4%)\n100 (16.7%)\n\n\n5\n57 (18.9%)\n43 (14.4%)\n100 (16.7%)\n\n\n6\n44 (14.6%)\n56 (18.7%)\n100 (16.7%)\n\n\n\n\n\n\n\n\n\n\n9.1.1.2 Model fitting\nWe synthesize the evidence using a Bayesian random effects meta-analysis model. The model is given in Equation 16.7 of the book. First we need set up the data and create the model:\n\nipd <- with(ds, ipdma.model.onestage(y = y, study = studyid, treat = treat,\n X = cbind(z1, z2), \n response = \"normal\", \n shrinkage = \"none\"), \n type = \"random\")\n\nThe JAGS model can be accessed as follows:\n\nipd$model.JAGS\n\nfunction () \n{\n for (i in 1:Np) {\n y[i] ~ dnorm(mu[i], sigma)\n mu[i] <- alpha[studyid[i]] + inprod(beta[], X[i, ]) + \n (1 - equals(treat[i], 1)) * inprod(gamma[], X[i, \n ]) + d[studyid[i], treat[i]]\n }\n sigma ~ dgamma(0.001, 0.001)\n for (j in 1:Nstudies) {\n d[j, 1] <- 0\n d[j, 2] ~ dnorm(delta[2], tau)\n }\n sd ~ dnorm(0, 1)\n T(0, )\n tau <- pow(sd, -2)\n delta[1] <- 0\n delta[2] ~ dnorm(0, 0.001)\n for (j in 1:Nstudies) {\n alpha[j] ~ dnorm(0, 0.001)\n }\n for (k in 1:Ncovariate) {\n beta[k] ~ dnorm(0, 0.001)\n }\n for (k in 1:Ncovariate) {\n gamma[k] ~ dnorm(0, 0.001)\n }\n}\n<environment: 0x5584451454f8>\n\n\nWe can fit the treatment effect model as follows:\n\nsamples <- ipd.run(ipd, n.chains = 2, n.iter = 20,\n pars.save = c(\"alpha\", \"beta\", \"delta\", \"sd\", \"gamma\"))\ntrtbenefit <- round(treatment.effect(ipd, samples, newpatient = c(z1 = 1, z2 = 0.5)), 2)\n\nHere are the estimated model parameters:\n\nsummary(samples)\n\n\nIterations = 2001:2020\nThinning interval = 1 \nNumber of chains = 2 \nSample size per chain = 20 \n\n1. Empirical mean and standard deviation for each variable,\n plus standard error of the mean:\n\n Mean SD Naive SE Time-series SE\nalpha[1] 11.0097 0.04718 0.007459 0.005721\nalpha[2] 8.0012 0.05172 0.008178 0.015132\nalpha[3] 10.5388 0.04920 0.007780 0.009529\nalpha[4] 9.5949 0.04693 0.007420 0.011516\nalpha[5] 13.0222 0.05621 0.008887 0.015491\nalpha[6] 15.7343 0.04603 0.007278 0.010399\nbeta[1] 0.1821 0.01954 0.003090 0.002915\nbeta[2] 0.2777 0.02037 0.003221 0.004587\ndelta[1] 0.0000 0.00000 0.000000 0.000000\ndelta[2] -2.7620 0.67509 0.106741 0.107957\ngamma[1] -0.4652 0.02908 0.004597 0.006635\ngamma[2] 0.5169 0.02678 0.004235 0.004570\nsd 1.8287 0.42482 0.067169 0.109751\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nalpha[1] 10.9227 10.9779 11.0163 11.0411 11.0878\nalpha[2] 7.9112 7.9685 8.0121 8.0458 8.0812\nalpha[3] 10.4667 10.5044 10.5322 10.5606 10.6264\nalpha[4] 9.5032 9.5740 9.5972 9.6238 9.7002\nalpha[5] 12.9421 12.9776 13.0109 13.0608 13.1283\nalpha[6] 15.6519 15.7099 15.7394 15.7631 15.8075\nbeta[1] 0.1430 0.1674 0.1859 0.1951 0.2127\nbeta[2] 0.2414 0.2664 0.2756 0.2927 0.3149\ndelta[1] 0.0000 0.0000 0.0000 0.0000 0.0000\ndelta[2] -4.1502 -3.2521 -2.8107 -2.2364 -1.6786\ngamma[1] -0.5276 -0.4799 -0.4636 -0.4498 -0.4121\ngamma[2] 0.4603 0.5034 0.5222 0.5361 0.5501\nsd 1.1511 1.4946 1.8143 1.9751 2.8172\n\n\n\n\n9.1.1.3 Prection\nWe can now predict the individualized treatment effect for a new patient with covariate values z1 = 1 and z2 = 0.5.\n\nround(treatment.effect(ipd, samples, newpatient = c(z1 = 1, z2 = 0.5)), 2)\n\n0.025 0.5 0.975 \n-4.38 -3.01 -1.91 \n\n\nThis means that the predicted outcome for patient with covariate values z1 = 1 and z2 = 0.5 will differ by -3.01 units when receiving the active treatment (treat = 1) as compared to the control treatment (treat = 0).\nWe can also predict treatment benefit for all patients in the sample, and look at the distribution of predicted benefit.\n\nlibrary(dplyr)\nlibrary(ggplot2)\n\nds <- ds %>% mutate(benefit = NA,\n study = paste(\"Trial\", studyid)) \n\nfor (i in seq(nrow(ds))) {\n newpat <- as.matrix(ds[i, c(\"z1\", \"z2\")])\n ds$benefit[i] <- treatment.effect(ipd, samples, newpatient = newpat)[\"0.5\"]\n}\n\nsummbenefit <- ds %>% group_by(study) %>% \n summarize(mediabenefit = median(benefit), meanbenefit = mean(benefit))\n\nggplot(ds, aes(x = benefit)) + \n geom_histogram(aes(y = after_stat(density)), alpha = 0.3) + \n geom_density() +\n geom_vline(data = summbenefit, aes(xintercept = meanbenefit), linewidth = 0.5, lty = 2) + \n facet_wrap(~study) + \n ylab(\"Density\") +\n xlab(\"Predicted treatment benefit\") + theme_bw()\n\n\n\n\nFigure 9.1: Distribution of predicted treatment benefit in each trial. Dashed lines represent the trial mean.\n\n\n\n\n\n\n9.1.1.4 Penalization\nLet us repeat the analysis, but this time while penalizing the treatment-covariate coefficients using a Bayesian LASSO prior.\n\nipd <- with(ds, ipdma.model.onestage(y = y, study = studyid, \n treat = treat,\n X = cbind(z1, z2), \n response = \"normal\", \n shrinkage = \"laplace\"), \n type = \"random\")\n\nsamples <- ipd.run(ipd, n.chains = 2, n.iter = 20, \n pars.save = c(\"alpha\", \"beta\", \"delta\", \"sd\", \"gamma\"))\n\nCompiling model graph\n Resolving undeclared variables\n Allocating nodes\nGraph information:\n Observed stochastic nodes: 600\n Unobserved stochastic nodes: 20\n Total graph size: 6039\n\nInitializing model\n\nround(treatment.effect(ipd, samples, newpatient = c(1,0.5)), 2)\n\n0.025 0.5 0.975 \n-3.69 -2.89 -2.04 \n\n\n\n\n\n9.1.2 Example of a binary outcome\n\n9.1.2.1 Setup\nWe now present the case of a binary outcome. We first generate a dataset as before, using the bipd package.\n\nds2 <- generate_ipdma_example(type = \"binary\")\nhead(ds2)\n\n studyid treat w1 w2 y\n1 1 1 1.1959217 -0.60042578 1\n2 1 0 0.4634961 0.08469719 0\n3 1 0 0.3403400 0.01961329 0\n4 1 0 -0.6867991 0.60877432 1\n5 1 0 -0.8644221 -0.07853599 0\n6 1 1 -1.3439234 -0.38140175 1\n\n\nThe simulated dataset contains information on the following variables:\n\nthe trial indicator studyid\nthe treatment indicator treat, which takes the values 0 for control and 1 for active treatment\ntwo prognostic variables w1 and w2\nthe binary outcome y\n\n\n\n\n\n\n\nTable 9.2: The simulated dataset with a binary outcome\n\n\n\n\n\n\n\n\n\n0\n(N=305)\n1\n(N=295)\nOverall\n(N=600)\n\n\n\n\nw1\n\n\n\n\n\nMean (SD)\n0.0169 (0.987)\n-0.0413 (0.985)\n-0.0117 (0.986)\n\n\nMedian [Min, Max]\n-0.0168 [-3.41, 2.50]\n-0.120 [-2.74, 3.02]\n-0.0624 [-3.41, 3.02]\n\n\nw2\n\n\n\n\n\nMean (SD)\n-0.111 (0.959)\n-0.0378 (1.07)\n-0.0750 (1.02)\n\n\nMedian [Min, Max]\n-0.159 [-3.46, 2.03]\n-0.0395 [-3.29, 2.40]\n-0.125 [-3.46, 2.40]\n\n\nstudyid\n\n\n\n\n\n1\n48 (15.7%)\n52 (17.6%)\n100 (16.7%)\n\n\n2\n52 (17.0%)\n48 (16.3%)\n100 (16.7%)\n\n\n3\n48 (15.7%)\n52 (17.6%)\n100 (16.7%)\n\n\n4\n54 (17.7%)\n46 (15.6%)\n100 (16.7%)\n\n\n5\n53 (17.4%)\n47 (15.9%)\n100 (16.7%)\n\n\n6\n50 (16.4%)\n50 (16.9%)\n100 (16.7%)\n\n\n\n\n\n\n\n\n\n\n9.1.2.2 Model fitting\nWe use a Bayesian random effects model with binomial likelihood. This is similar to the model 16.7 of the book, but with a Binomial likelihood, i.e. \n\\[\ny_{ij}\\sim \\text{Binomial}(\\pi_{ij}) \\\\\n\\] \\[\n\\text{logit}(\\pi_{ij})=a_j+\\delta_j t_{ij}+ \\sum_{l=1}^{L}\\beta_l x_{ij}+ \\sum_{l=1}^{L}\\gamma_l x_{ij} t_{ij}\n\\] The remaining of the model is as in the book. We can penalize the estimated parameters for effect modification (\\(\\gamma\\)’s), using a Bayesian LASSO. We can do this using again the bipd package:\n\nipd2 <- with(ds2, ipdma.model.onestage(y = y, study = studyid, treat = treat,\n X = cbind(w1, w2), \n response = \"binomial\", \n shrinkage = \"laplace\"), \n type = \"random\", hy.prior = list(\"dunif\", 0, 1))\n\nipd2$model.JAGS\n\nfunction () \n{\n for (i in 1:Np) {\n y[i] ~ dbern(p[i])\n logit(p[i]) <- alpha[studyid[i]] + inprod(beta[], X[i, \n ]) + (1 - equals(treat[i], 1)) * inprod(gamma[], \n X[i, ]) + d[studyid[i], treat[i]]\n }\n for (j in 1:Nstudies) {\n d[j, 1] <- 0\n d[j, 2] ~ dnorm(delta[2], tau)\n }\n sd ~ dnorm(0, 1)\n T(0, )\n tau <- pow(sd, -2)\n delta[1] <- 0\n delta[2] ~ dnorm(0, 0.001)\n for (j in 1:Nstudies) {\n alpha[j] ~ dnorm(0, 0.001)\n }\n for (k in 1:Ncovariate) {\n beta[k] ~ dnorm(0, 0.001)\n }\n tt <- lambda\n lambda <- pow(lambda.inv, -1)\n lambda.inv ~ dunif(0, 5)\n for (k in 1:Ncovariate) {\n gamma[k] ~ ddexp(0, tt)\n }\n}\n<environment: 0x558448d79c10>\n\n\n\nsamples <- ipd.run(ipd2, n.chains = 2, n.iter = 20, \n pars.save = c(\"alpha\", \"beta\", \"delta\", \"sd\", \"gamma\"))\nsummary(samples)\n\n\n\n\nIterations = 2001:2020\nThinning interval = 1 \nNumber of chains = 2 \nSample size per chain = 20 \n\n1. Empirical mean and standard deviation for each variable,\n plus standard error of the mean:\n\n Mean SD Naive SE Time-series SE\nalpha[1] 0.784188 0.32122 0.05079 0.06507\nalpha[2] 0.791775 0.24145 0.03818 0.04003\nalpha[3] 0.915208 0.21282 0.03365 0.04596\nalpha[4] 0.928365 0.23917 0.03782 0.05305\nalpha[5] 0.795365 0.29118 0.04604 0.06140\nalpha[6] 1.125847 0.27497 0.04348 0.05068\nbeta[1] 0.006016 0.11411 0.01804 0.02394\nbeta[2] 0.165149 0.08558 0.01353 0.01393\ndelta[1] 0.000000 0.00000 0.00000 0.00000\ndelta[2] -0.253590 0.21988 0.03477 0.03949\ngamma[1] 0.003805 0.12806 0.02025 0.03517\ngamma[2] 0.004944 0.11258 0.01780 0.01571\nsd 0.450381 0.18644 0.02948 0.06072\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nalpha[1] 0.24633 0.57703 0.8375355 1.01716 1.1619\nalpha[2] 0.28460 0.61171 0.8505306 0.96326 1.0852\nalpha[3] 0.54327 0.73884 0.9203480 1.05797 1.3027\nalpha[4] 0.35900 0.79618 0.9281926 1.08213 1.3064\nalpha[5] 0.41042 0.57199 0.6861846 1.02946 1.3361\nalpha[6] 0.70948 0.91596 1.1206712 1.35633 1.5907\nbeta[1] -0.23670 -0.06909 0.0389493 0.08078 0.1939\nbeta[2] -0.01212 0.11471 0.1716787 0.22121 0.3231\ndelta[1] 0.00000 0.00000 0.0000000 0.00000 0.0000\ndelta[2] -0.71335 -0.42277 -0.2665279 -0.08168 0.1310\ngamma[1] -0.23104 -0.07517 0.0007591 0.05893 0.2166\ngamma[2] -0.15422 -0.07898 -0.0080998 0.08461 0.2238\nsd 0.11577 0.31971 0.4541077 0.59847 0.7210\n\n\nThe predicted treatment benefit for a new patient with covariates w1 = 1.6 and w2 = 1.3 is given as:\n\nround(treatment.effect(ipd2, samples, newpatient = c(w1 = 1.6, w2 = 1.3)), 2)\n\n0.025 0.5 0.975 \n 0.41 0.86 1.27 \n\n\nIn other words, the aforementioned patient 0.86 (95% Credibility Interval: 0.41 to 1.27)" + "text": "9.1 Estimating heterogeneous treatment effects in pairwise meta-analysis\nWe hereby provide code for estimating patient-level treatment effects for the case when we have patient-level data from multiple randomized trials.\n\n9.1.1 Example of a continuous outcome\n\n9.1.1.1 Setup\nWe start by simulating an artificial dataset using the R package bipd:\n\nlibrary(bipd)\nds <- generate_ipdma_example(type = \"continuous\")\n\nLet us have a look at the dataset:\n\nhead(ds)\n\n studyid treat z1 z2 y\n1 1 0 -0.06226787 -0.3046572 11\n2 1 1 0.22873835 -1.3755850 8\n3 1 1 -0.39734360 0.4309557 10\n4 1 0 -0.48128849 -0.8760153 11\n5 1 1 -0.07360034 -1.8460789 7\n6 1 0 -2.00590597 -0.6318041 10\n\n\nThe simulated dataset contains information on the following variables:\n\nthe trial indicator studyid\nthe treatment indicator treat, which takes the values 0 for control and 1 for active treatment\ntwo prognostic variables z1 and z2\nthe continuous outcome y\n\n\n\n\n\n\n\nTable 9.1: The simulated dataset with a continuous outcome\n\n\n\n\n\n\n\n\n\n0\n(N=296)\n1\n(N=304)\nOverall\n(N=600)\n\n\n\n\nz1\n\n\n\n\n\nMean (SD)\n0.0522 (1.00)\n-0.0126 (0.990)\n0.0194 (0.995)\n\n\nMedian [Min, Max]\n0.0611 [-2.68, 3.33]\n-0.0719 [-2.29, 3.32]\n-0.0371 [-2.68, 3.33]\n\n\nz2\n\n\n\n\n\nMean (SD)\n-0.0736 (0.929)\n0.00776 (0.978)\n-0.0324 (0.954)\n\n\nMedian [Min, Max]\n-0.0719 [-2.79, 2.79]\n-0.0234 [-2.52, 2.77]\n-0.0617 [-2.79, 2.79]\n\n\nstudyid\n\n\n\n\n\n1\n46 (15.5%)\n54 (17.8%)\n100 (16.7%)\n\n\n2\n52 (17.6%)\n48 (15.8%)\n100 (16.7%)\n\n\n3\n52 (17.6%)\n48 (15.8%)\n100 (16.7%)\n\n\n4\n46 (15.5%)\n54 (17.8%)\n100 (16.7%)\n\n\n5\n50 (16.9%)\n50 (16.4%)\n100 (16.7%)\n\n\n6\n50 (16.9%)\n50 (16.4%)\n100 (16.7%)\n\n\n\n\n\n\n\n\n\n\n9.1.1.2 Model fitting\nWe synthesize the evidence using a Bayesian random effects meta-analysis model. The model is given in Equation 16.7 of the book. First we need set up the data and create the model:\n\nipd <- with(ds, ipdma.model.onestage(y = y, study = studyid, treat = treat,\n X = cbind(z1, z2), \n response = \"normal\", \n shrinkage = \"none\"), \n type = \"random\")\n\nThe JAGS model can be accessed as follows:\n\nipd$model.JAGS\n\nfunction () \n{\n for (i in 1:Np) {\n y[i] ~ dnorm(mu[i], sigma)\n mu[i] <- alpha[studyid[i]] + inprod(beta[], X[i, ]) + \n (1 - equals(treat[i], 1)) * inprod(gamma[], X[i, \n ]) + d[studyid[i], treat[i]]\n }\n sigma ~ dgamma(0.001, 0.001)\n for (j in 1:Nstudies) {\n d[j, 1] <- 0\n d[j, 2] ~ dnorm(delta[2], tau)\n }\n sd ~ dnorm(0, 1)\n T(0, )\n tau <- pow(sd, -2)\n delta[1] <- 0\n delta[2] ~ dnorm(0, 0.001)\n for (j in 1:Nstudies) {\n alpha[j] ~ dnorm(0, 0.001)\n }\n for (k in 1:Ncovariate) {\n beta[k] ~ dnorm(0, 0.001)\n }\n for (k in 1:Ncovariate) {\n gamma[k] ~ dnorm(0, 0.001)\n }\n}\n<environment: 0x5556cce09328>\n\n\nWe can fit the treatment effect model as follows:\n\nsamples <- ipd.run(ipd, n.chains = 2, n.iter = 20,\n pars.save = c(\"alpha\", \"beta\", \"delta\", \"sd\", \"gamma\"))\ntrtbenefit <- round(treatment.effect(ipd, samples, newpatient = c(z1 = 1, z2 = 0.5)), 2)\n\nHere are the estimated model parameters:\n\nsummary(samples)\n\n\nIterations = 2001:2020\nThinning interval = 1 \nNumber of chains = 2 \nSample size per chain = 20 \n\n1. Empirical mean and standard deviation for each variable,\n plus standard error of the mean:\n\n Mean SD Naive SE Time-series SE\nalpha[1] 10.9849 0.05156 0.008153 0.010193\nalpha[2] 8.0479 0.05597 0.008850 0.010776\nalpha[3] 10.5161 0.05724 0.009050 0.013192\nalpha[4] 9.5786 0.05669 0.008964 0.009557\nalpha[5] 12.7748 0.04828 0.007633 0.010383\nalpha[6] 15.7653 0.04403 0.006961 0.013683\nbeta[1] 0.2210 0.02398 0.003792 0.006286\nbeta[2] 0.2415 0.02058 0.003254 0.005248\ndelta[1] 0.0000 0.00000 0.000000 0.000000\ndelta[2] -2.1425 0.34904 0.055188 0.058311\ngamma[1] -0.5197 0.03174 0.005018 0.008381\ngamma[2] 0.5524 0.03240 0.005122 0.007422\nsd 0.8678 0.18934 0.029938 0.030147\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nalpha[1] 10.8945 10.9453 10.9882 11.0272 11.0552\nalpha[2] 7.9621 8.0092 8.0404 8.0770 8.1599\nalpha[3] 10.4196 10.4736 10.5094 10.5600 10.6018\nalpha[4] 9.4696 9.5374 9.5820 9.6229 9.6749\nalpha[5] 12.6896 12.7393 12.7767 12.8050 12.8493\nalpha[6] 15.6834 15.7370 15.7701 15.7952 15.8465\nbeta[1] 0.1774 0.2003 0.2241 0.2337 0.2647\nbeta[2] 0.2006 0.2296 0.2436 0.2542 0.2878\ndelta[1] 0.0000 0.0000 0.0000 0.0000 0.0000\ndelta[2] -2.6947 -2.4129 -2.1227 -1.8953 -1.5725\ngamma[1] -0.5627 -0.5429 -0.5233 -0.4980 -0.4536\ngamma[2] 0.4830 0.5322 0.5561 0.5734 0.6010\nsd 0.5870 0.7426 0.8418 0.9520 1.2676\n\n\n\n\n9.1.1.3 Prection\nWe can now predict the individualized treatment effect for a new patient with covariate values z1 = 1 and z2 = 0.5.\n\nround(treatment.effect(ipd, samples, newpatient = c(z1 = 1, z2 = 0.5)), 2)\n\n0.025 0.5 0.975 \n-2.92 -2.32 -1.77 \n\n\nThis means that the predicted outcome for patient with covariate values z1 = 1 and z2 = 0.5 will differ by -2.32 units when receiving the active treatment (treat = 1) as compared to the control treatment (treat = 0).\nWe can also predict treatment benefit for all patients in the sample, and look at the distribution of predicted benefit.\n\nlibrary(dplyr)\nlibrary(ggplot2)\n\nds <- ds %>% mutate(benefit = NA,\n study = paste(\"Trial\", studyid)) \n\nfor (i in seq(nrow(ds))) {\n newpat <- as.matrix(ds[i, c(\"z1\", \"z2\")])\n ds$benefit[i] <- treatment.effect(ipd, samples, newpatient = newpat)[\"0.5\"]\n}\n\nsummbenefit <- ds %>% group_by(study) %>% \n summarize(mediabenefit = median(benefit), meanbenefit = mean(benefit))\n\nggplot(ds, aes(x = benefit)) + \n geom_histogram(aes(y = after_stat(density)), alpha = 0.3) + \n geom_density() +\n geom_vline(data = summbenefit, aes(xintercept = meanbenefit), \n linewidth = 0.5, lty = 2) + \n facet_wrap(~study) + \n ylab(\"Density\") +\n xlab(\"Predicted treatment benefit\") + theme_bw()\n\n\n\n\nFigure 9.1: Distribution of predicted treatment benefit in each trial. Dashed lines represent the trial mean.\n\n\n\n\n\n\n9.1.1.4 Penalization\nLet us repeat the analysis, but this time while penalizing the treatment-covariate coefficients using a Bayesian LASSO prior.\n\nipd <- with(ds, ipdma.model.onestage(y = y, study = studyid, \n treat = treat,\n X = cbind(z1, z2), \n response = \"normal\", \n shrinkage = \"laplace\"), \n type = \"random\")\n\nsamples <- ipd.run(ipd, n.chains = 2, n.iter = 20, \n pars.save = c(\"alpha\", \"beta\", \"delta\", \"sd\", \"gamma\"))\n\nCompiling model graph\n Resolving undeclared variables\n Allocating nodes\nGraph information:\n Observed stochastic nodes: 600\n Unobserved stochastic nodes: 20\n Total graph size: 6039\n\nInitializing model\n\nround(treatment.effect(ipd, samples, newpatient = c(1,0.5)), 2)\n\n0.025 0.5 0.975 \n-3.37 -2.28 -1.30 \n\n\n\n\n\n9.1.2 Example of a binary outcome\n\n9.1.2.1 Setup\nWe now present the case of a binary outcome. We first generate a dataset as before, using the bipd package.\n\nds2 <- generate_ipdma_example(type = \"binary\")\nhead(ds2)\n\n studyid treat w1 w2 y\n1 1 0 -0.5612891 0.65769118 1\n2 1 1 1.1342349 0.85682668 0\n3 1 1 0.1005061 -2.20312003 0\n4 1 1 -0.5153136 -0.57015407 0\n5 1 1 -1.0159330 -0.78384436 0\n6 1 0 -0.2962515 -0.01832073 0\n\n\nThe simulated dataset contains information on the following variables:\n\nthe trial indicator studyid\nthe treatment indicator treat, which takes the values 0 for control and 1 for active treatment\ntwo prognostic variables w1 and w2\nthe binary outcome y\n\n\n\n\n\n\n\nTable 9.2: The simulated dataset with a binary outcome\n\n\n\n\n\n\n\n\n\n0\n(N=305)\n1\n(N=295)\nOverall\n(N=600)\n\n\n\n\nw1\n\n\n\n\n\nMean (SD)\n-0.00720 (1.02)\n0.0526 (0.985)\n0.0222 (1.00)\n\n\nMedian [Min, Max]\n-0.0492 [-2.73, 2.89]\n0.0449 [-3.23, 3.33]\n-0.0241 [-3.23, 3.33]\n\n\nw2\n\n\n\n\n\nMean (SD)\n0.0465 (1.00)\n-0.0340 (1.01)\n0.00695 (1.01)\n\n\nMedian [Min, Max]\n0.0894 [-2.37, 2.93]\n0.0399 [-2.72, 4.21]\n0.0748 [-2.72, 4.21]\n\n\nstudyid\n\n\n\n\n\n1\n47 (15.4%)\n53 (18.0%)\n100 (16.7%)\n\n\n2\n46 (15.1%)\n54 (18.3%)\n100 (16.7%)\n\n\n3\n56 (18.4%)\n44 (14.9%)\n100 (16.7%)\n\n\n4\n50 (16.4%)\n50 (16.9%)\n100 (16.7%)\n\n\n5\n51 (16.7%)\n49 (16.6%)\n100 (16.7%)\n\n\n6\n55 (18.0%)\n45 (15.3%)\n100 (16.7%)\n\n\n\n\n\n\n\n\n\n\n9.1.2.2 Model fitting\nWe use a Bayesian random effects model with binomial likelihood. This is similar to the model 16.7 of the book, but with a Binomial likelihood, i.e. \n\\[\ny_{ij}\\sim \\text{Binomial}(\\pi_{ij}) \\\\\n\\] \\[\n\\text{logit}(\\pi_{ij})=a_j+\\delta_j t_{ij}+ \\sum_{l=1}^{L}\\beta_l x_{ij}+ \\sum_{l=1}^{L}\\gamma_l x_{ij} t_{ij}\n\\] The remaining of the model is as in the book. We can penalize the estimated parameters for effect modification (\\(\\gamma\\)’s), using a Bayesian LASSO. We can do this using again the bipd package:\n\nipd2 <- with(ds2, ipdma.model.onestage(y = y, study = studyid, treat = treat,\n X = cbind(w1, w2), \n response = \"binomial\", \n shrinkage = \"laplace\"), \n type = \"random\", hy.prior = list(\"dunif\", 0, 1))\n\nipd2$model.JAGS\n\nfunction () \n{\n for (i in 1:Np) {\n y[i] ~ dbern(p[i])\n logit(p[i]) <- alpha[studyid[i]] + inprod(beta[], X[i, \n ]) + (1 - equals(treat[i], 1)) * inprod(gamma[], \n X[i, ]) + d[studyid[i], treat[i]]\n }\n for (j in 1:Nstudies) {\n d[j, 1] <- 0\n d[j, 2] ~ dnorm(delta[2], tau)\n }\n sd ~ dnorm(0, 1)\n T(0, )\n tau <- pow(sd, -2)\n delta[1] <- 0\n delta[2] ~ dnorm(0, 0.001)\n for (j in 1:Nstudies) {\n alpha[j] ~ dnorm(0, 0.001)\n }\n for (k in 1:Ncovariate) {\n beta[k] ~ dnorm(0, 0.001)\n }\n tt <- lambda\n lambda <- pow(lambda.inv, -1)\n lambda.inv ~ dunif(0, 5)\n for (k in 1:Ncovariate) {\n gamma[k] ~ ddexp(0, tt)\n }\n}\n<environment: 0x5556d0a8e470>\n\n\n\nsamples <- ipd.run(ipd2, n.chains = 2, n.iter = 20, \n pars.save = c(\"alpha\", \"beta\", \"delta\", \"sd\", \"gamma\"))\nsummary(samples)\n\n\n\n\nIterations = 2001:2020\nThinning interval = 1 \nNumber of chains = 2 \nSample size per chain = 20 \n\n1. Empirical mean and standard deviation for each variable,\n plus standard error of the mean:\n\n Mean SD Naive SE Time-series SE\nalpha[1] -0.12870 0.2232 0.03530 0.03561\nalpha[2] -1.35792 0.2735 0.04325 0.08085\nalpha[3] -1.21444 0.2672 0.04225 0.05449\nalpha[4] -1.05525 0.3042 0.04810 0.07045\nalpha[5] -0.92933 0.2450 0.03874 0.05174\nalpha[6] -1.20645 0.2515 0.03977 0.03976\nbeta[1] -0.10056 0.1082 0.01710 0.02352\nbeta[2] -0.11918 0.1155 0.01826 0.03240\ndelta[1] 0.00000 0.0000 0.00000 0.00000\ndelta[2] -0.46724 0.7439 0.11761 0.16827\ngamma[1] 0.10270 0.1523 0.02409 0.03066\ngamma[2] 0.06879 0.1718 0.02717 0.04397\nsd 1.42643 0.3439 0.05438 0.09879\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nalpha[1] -0.4415 -0.2814894 -0.14889 -0.0004158 0.29590\nalpha[2] -1.8719 -1.5500880 -1.27938 -1.1610952 -0.93129\nalpha[3] -1.6696 -1.4417289 -1.15871 -1.0087951 -0.86845\nalpha[4] -1.5900 -1.2770702 -1.11096 -0.8363314 -0.43803\nalpha[5] -1.3737 -1.0954269 -0.92882 -0.7846082 -0.43321\nalpha[6] -1.6649 -1.3707223 -1.19790 -1.0033566 -0.84776\nbeta[1] -0.3137 -0.1518571 -0.09500 -0.0398302 0.09449\nbeta[2] -0.4063 -0.1586791 -0.13276 -0.0716571 0.10386\ndelta[1] 0.0000 0.0000000 0.00000 0.0000000 0.00000\ndelta[2] -2.0235 -0.8224947 -0.52883 -0.0107413 0.79044\ngamma[1] -0.1809 0.0003475 0.12003 0.2052612 0.34881\ngamma[2] -0.1912 -0.0437462 0.04734 0.1713008 0.46497\nsd 0.8210 1.1530501 1.41743 1.7317748 1.99153\n\n\nThe predicted treatment benefit for a new patient with covariates w1 = 1.6 and w2 = 1.3 is given as:\n\nround(treatment.effect(ipd2, samples, newpatient = c(w1 = 1.6, w2 = 1.3)), 2)\n\n0.025 0.5 0.975 \n 0.15 0.91 2.30 \n\n\nIn other words, the aforementioned patient 0.91 (95% Credibility Interval: 0.15 to 2.3)" }, { "objectID": "chapter_16.html#estimating-heterogeous-treatment-effects-in-network-meta-analysis", "href": "chapter_16.html#estimating-heterogeous-treatment-effects-in-network-meta-analysis", "title": "9  Prediction of individual treatment effect using data from multiple studies", "section": "9.2 Estimating heterogeous treatment effects in network meta-analysis", - "text": "9.2 Estimating heterogeous treatment effects in network meta-analysis\n\n9.2.1 Example of a continuous outcome\n\n9.2.1.1 Setup\nWe use again the bipd package to simulate a dataset:\n\nds3 <- generate_ipdnma_example(type = \"continuous\")\nhead(ds3)\n\n studyid treat z1 z2 y\n1 1 2 -0.97337590 -0.6695393 8\n2 1 1 -1.02809238 -0.6817386 11\n3 1 2 2.35703179 1.6272272 8\n4 1 2 -0.04532726 0.6444098 9\n5 1 2 -0.11708756 1.6158472 10\n6 1 1 1.08949679 -0.6346255 11\n\n\nLet us look into the data a bit in more detail:\n\n\n\n\n\n\nTable 9.3: The simulated dataset with a continuous outcome\n\n\n\n\n\n\n\n\n\n\n1\n(N=346)\n2\n(N=350)\n3\n(N=304)\nOverall\n(N=1000)\n\n\n\n\nz1\n\n\n\n\n\n\nMean (SD)\n-0.103 (1.01)\n0.0226 (1.01)\n-0.0151 (1.07)\n-0.0325 (1.03)\n\n\nMedian [Min, Max]\n-0.160 [-3.18, 2.42]\n0.0532 [-2.86, 3.01]\n-0.00116 [-3.34, 3.12]\n-0.0408 [-3.34, 3.12]\n\n\nz2\n\n\n\n\n\n\nMean (SD)\n0.0224 (0.891)\n0.148 (1.07)\n0.000382 (1.04)\n0.0596 (1.00)\n\n\nMedian [Min, Max]\n-0.0321 [-2.29, 2.92]\n0.129 [-2.74, 3.44]\n0.0611 [-2.85, 2.99]\n0.0600 [-2.85, 3.44]\n\n\nstudyid\n\n\n\n\n\n\n1\n52 (15.0%)\n48 (13.7%)\n0 (0%)\n100 (10.0%)\n\n\n2\n40 (11.6%)\n60 (17.1%)\n0 (0%)\n100 (10.0%)\n\n\n3\n50 (14.5%)\n50 (14.3%)\n0 (0%)\n100 (10.0%)\n\n\n4\n52 (15.0%)\n0 (0%)\n48 (15.8%)\n100 (10.0%)\n\n\n5\n47 (13.6%)\n0 (0%)\n53 (17.4%)\n100 (10.0%)\n\n\n6\n0 (0%)\n44 (12.6%)\n56 (18.4%)\n100 (10.0%)\n\n\n7\n0 (0%)\n55 (15.7%)\n45 (14.8%)\n100 (10.0%)\n\n\n8\n35 (10.1%)\n35 (10.0%)\n30 (9.9%)\n100 (10.0%)\n\n\n9\n42 (12.1%)\n29 (8.3%)\n29 (9.5%)\n100 (10.0%)\n\n\n10\n28 (8.1%)\n29 (8.3%)\n43 (14.1%)\n100 (10.0%)\n\n\n\n\n\n\n\n\n\n\n9.2.1.2 Model fitting\nWe will use the model shown in Equation 16.8 in the book. In addition, we will use Bayesian LASSO to penalize the treatment-covariate interactions.\n\nipd3 <- with(ds3, ipdnma.model.onestage(y = y, study = studyid, treat = treat, \n X = cbind(z1, z2), \n response = \"normal\", \n shrinkage = \"laplace\", \n type = \"random\"))\nipd3$model.JAGS\n\nfunction () \n{\n for (i in 1:Np) {\n y[i] ~ dnorm(mu[i], sigma)\n mu[i] <- alpha[studyid[i]] + inprod(beta[], X[i, ]) + \n inprod(gamma[treat[i], ], X[i, ]) + d[studyid[i], \n treatment.arm[i]]\n }\n sigma ~ dgamma(0.001, 0.001)\n for (i in 1:Nstudies) {\n w[i, 1] <- 0\n d[i, 1] <- 0\n for (k in 2:na[i]) {\n d[i, k] ~ dnorm(mdelta[i, k], taudelta[i, k])\n mdelta[i, k] <- delta[t[i, k]] - delta[t[i, 1]] + \n sw[i, k]\n taudelta[i, k] <- tau * 2 * (k - 1)/k\n w[i, k] <- d[i, k] - delta[t[i, k]] + delta[t[i, \n 1]]\n sw[i, k] <- sum(w[i, 1:(k - 1)])/(k - 1)\n }\n }\n sd ~ dnorm(0, 1)\n T(0, )\n tau <- pow(sd, -2)\n delta[1] <- 0\n for (k in 2:Ntreat) {\n delta[k] ~ dnorm(0, 0.001)\n }\n for (j in 1:Nstudies) {\n alpha[j] ~ dnorm(0, 0.001)\n }\n for (k in 1:Ncovariate) {\n beta[k] ~ dnorm(0, 0.001)\n }\n lambda[1] <- 0\n lambda.inv[1] <- 0\n for (m in 2:Ntreat) {\n tt[m] <- lambda[m] * sigma\n lambda[m] <- pow(lambda.inv[m], -1)\n lambda.inv[m] ~ dunif(0, 5)\n }\n for (k in 1:Ncovariate) {\n gamma[1, k] <- 0\n for (m in 2:Ntreat) {\n gamma[m, k] ~ ddexp(0, tt[m])\n }\n }\n}\n<environment: 0x55844907a308>\n\nsamples <- ipd.run(ipd3, n.chains = 2, n.iter = 20, \n pars.save = c(\"alpha\", \"beta\", \"delta\", \"sd\", \"gamma\"))\n\nCompiling model graph\n Resolving undeclared variables\n Allocating nodes\nGraph information:\n Observed stochastic nodes: 1000\n Unobserved stochastic nodes: 35\n Total graph size: 10141\n\nInitializing model\n\nsummary(samples)\n\n\nIterations = 2001:2020\nThinning interval = 1 \nNumber of chains = 2 \nSample size per chain = 20 \n\n1. Empirical mean and standard deviation for each variable,\n plus standard error of the mean:\n\n Mean SD Naive SE Time-series SE\nalpha[1] 11.0984 0.05546 0.008768 0.016654\nalpha[2] 8.0310 0.04769 0.007541 0.009679\nalpha[3] 10.5087 0.04524 0.007153 0.009658\nalpha[4] 9.5946 0.04952 0.007830 0.010258\nalpha[5] 12.9679 0.04653 0.007356 0.010654\nalpha[6] 13.2682 0.04732 0.007482 0.010367\nalpha[7] 7.5078 0.04455 0.007045 0.007418\nalpha[8] 11.1698 0.04510 0.007131 0.009550\nalpha[9] 10.1576 0.03590 0.005676 0.005631\nalpha[10] 9.2693 0.05017 0.007932 0.011032\nbeta[1] 0.2103 0.01864 0.002947 0.005830\nbeta[2] 0.3254 0.02057 0.003252 0.006833\ndelta[1] 0.0000 0.00000 0.000000 0.000000\ndelta[2] -2.9253 0.06482 0.010248 0.013653\ndelta[3] -1.1187 0.05289 0.008363 0.008103\ngamma[1,1] 0.0000 0.00000 0.000000 0.000000\ngamma[2,1] -0.5765 0.02944 0.004654 0.008253\ngamma[3,1] -0.3579 0.02442 0.003861 0.005350\ngamma[1,2] 0.0000 0.00000 0.000000 0.000000\ngamma[2,2] 0.5498 0.03554 0.005619 0.012004\ngamma[3,2] 0.4372 0.02777 0.004391 0.008825\nsd 0.1432 0.03491 0.005520 0.006474\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nalpha[1] 11.00852 11.0632 11.0934 11.1230 11.2241\nalpha[2] 7.95427 8.0048 8.0303 8.0516 8.1236\nalpha[3] 10.43685 10.4756 10.5060 10.5516 10.5810\nalpha[4] 9.49399 9.5614 9.6023 9.6278 9.6697\nalpha[5] 12.89895 12.9397 12.9632 12.9851 13.0433\nalpha[6] 13.18158 13.2385 13.2684 13.3092 13.3368\nalpha[7] 7.41889 7.4823 7.5125 7.5310 7.5836\nalpha[8] 11.09612 11.1381 11.1746 11.1962 11.2503\nalpha[9] 10.09623 10.1340 10.1537 10.1831 10.2269\nalpha[10] 9.19161 9.2345 9.2589 9.2951 9.3850\nbeta[1] 0.17774 0.1943 0.2170 0.2243 0.2361\nbeta[2] 0.29120 0.3104 0.3249 0.3389 0.3700\ndelta[1] 0.00000 0.0000 0.0000 0.0000 0.0000\ndelta[2] -3.03557 -2.9760 -2.9293 -2.8730 -2.7999\ndelta[3] -1.20401 -1.1651 -1.1231 -1.0734 -1.0185\ngamma[1,1] 0.00000 0.0000 0.0000 0.0000 0.0000\ngamma[2,1] -0.62032 -0.5970 -0.5864 -0.5590 -0.5197\ngamma[3,1] -0.40027 -0.3758 -0.3609 -0.3366 -0.3165\ngamma[1,2] 0.00000 0.0000 0.0000 0.0000 0.0000\ngamma[2,2] 0.49567 0.5272 0.5455 0.5731 0.6175\ngamma[3,2] 0.38928 0.4150 0.4390 0.4591 0.4823\nsd 0.09005 0.1208 0.1478 0.1615 0.2236\n\n\nAs before, we can use the treatment.effect() function of bipd to estimate relative effects for new patients.\n\ntreatment.effect(ipd3, samples, newpatient= c(1,2))\n\n$`treatment 2`\n 0.025 0.5 0.975 \n-2.641458 -2.463162 -2.249292 \n\n$`treatment 3`\n 0.025 0.5 0.975 \n-0.8048313 -0.6370218 -0.4668667 \n\n\nThis gives us the relative effects for all treatments versus the reference. To obtain relative effects between active treatments we need some more coding:\n\nsamples.all=data.frame(rbind(samples[[1]], samples[[2]]))\nnewpatient= c(1,2)\nnewpatient <- (newpatient - ipd3$scale_mean)/ipd3$scale_sd\n\nmedian(\n samples.all$delta.2.+samples.all$gamma.2.1.*\n newpatient[1]+samples.all$gamma.2.2.*newpatient[2]\n-\n (samples.all$delta.3.+samples.all$gamma.3.1.*newpatient[1]+\n samples.all$gamma.3.2.*newpatient[2])\n)\n\n[1] -1.806833\n\nquantile(samples.all$delta.2.+samples.all$gamma.2.1.*\n newpatient[1]+samples.all$gamma.2.2.*newpatient[2]\n -(samples.all$delta.3.+samples.all$gamma.3.1.*newpatient[1]+\n samples.all$gamma.3.2.*newpatient[2])\n , probs = 0.025)\n\n 2.5% \n-1.998861 \n\nquantile(samples.all$delta.2.+samples.all$gamma.2.1.*\n newpatient[1]+samples.all$gamma.2.2.*newpatient[2]\n -(samples.all$delta.3.+samples.all$gamma.3.1.*newpatient[1]+\n samples.all$gamma.3.2.*newpatient[2])\n , probs = 0.975)\n\n 97.5% \n-1.625211 \n\n\n\n\n\n9.2.2 Modeling patient-level relative effects using randomized and observational evidence for a network of treatments\nWe will now follow Chapter 16.3.5 from the book. In this analysis we will not use penalization, and we will assume fixed effects. For an example with penalization and random effects, see part 2 of this vignettte.\n\n9.2.2.1 Setup\nWe generate a very simple dataset of three studies comparing three treatments. We will assume 2 RCTs and 1 non-randomized trial:\n\nds4 <- generate_ipdnma_example(type = \"continuous\")\nds4 <- ds4 %>% filter(studyid %in% c(1,4,10)) %>%\n mutate(studyid = factor(studyid) %>%\n recode_factor(\n \"1\" = \"1\",\n \"4\" = \"2\",\n \"10\" = \"3\"),\n design = ifelse(studyid == \"3\", \"nrs\", \"rct\"))\n\nThe sample size is as follows:\n\n\n \n s1 s2 s3\n treat A: 56 53 34\n treat B: 44 0 29\n treat C: 0 47 37\n\n\n\n\n9.2.2.2 Model fitting\nWe will use the design-adjusted model, equation 16.9 in the book. We will fit a two-stage fixed effects meta-analysis and we will use a variance inflation factor. The code below is used to specify the analysis of each individual study. Briefly, in each study we adjust the treatment effect for the prognostic factors z1 and z2, as well as their interaction with treat.\n\nlibrary(rjags)\n\nLoading required package: coda\n\n\nLinked to JAGS 4.3.0\n\n\nLoaded modules: basemod,bugs\n\nfirst.stage <- \"\nmodel{\n\nfor (i in 1:N){\n y[i] ~ dnorm(mu[i], tau) \n mu[i] <- a + inprod(b[], X[i,]) + inprod(c[,treat[i]], X[i,]) + d[treat[i]] \n}\nsigma ~ dunif(0, 5)\ntau <- pow(sigma, -2)\n\na ~ dnorm(0, 0.001)\n\nfor(k in 1:Ncovariate){\n b[k] ~ dnorm(0,0.001)\n}\n\nfor(k in 1:Ncovariate){\n c[k,1] <- 0\n}\n\ntauGamma <- pow(sdGamma,-1)\nsdGamma ~ dunif(0, 5)\n\nfor(k in 1:Ncovariate){\n for(t in 2:Ntreat){\n c[k,t] ~ ddexp(0, tauGamma)\n }\n}\n\nd[1] <- 0\nfor(t in 2:Ntreat){\n d[t] ~ dnorm(0, 0.001)\n}\n}\"\n\nSubsequently, we estimate the relative treatment effects in the first (randomized) study comparing treatments A and B:\n\nmodel1.spec <- textConnection(first.stage) \ndata1 <- with(ds4 %>% filter(studyid == 1), \n list(y = y,\n N = length(y), \n X = cbind(z1,z2), \n treat = treat,\n Ncovariate = 2, \n Ntreat = 2))\njags.m <- jags.model(model1.spec, data = data1, n.chains = 2, n.adapt = 500,\n quiet = TRUE)\nparams <- c(\"d\", \"c\") \nsamps4.1 <- coda.samples(jags.m, params, n.iter = 50)\nsamps.all.s1 <- data.frame(as.matrix(samps4.1))\n\nsamps.all.s1 <- samps.all.s1[, c(\"c.1.2.\", \"c.2.2.\", \"d.2.\")]\ndelta.1 <- colMeans(samps.all.s1)\ncov.1 <- var(samps.all.s1)\n\nWe repeat the analysis for the second (randomized) study comparing treatments A and C:\n\nmodel1.spec <- textConnection(first.stage) \ndata2 <- with(ds4 %>% filter(studyid == 2), \n list(y = y,\n N = length(y), \n X = cbind(z1,z2), \n treat = ifelse(treat == 3, 2, treat),\n Ncovariate = 2, \n Ntreat = 2))\njags.m <- jags.model(model1.spec, data = data2, n.chains = 2, n.adapt = 100,\n quiet = TRUE)\nparams <- c(\"d\", \"c\") \nsamps4.2 <- coda.samples(jags.m, params, n.iter = 50)\nsamps.all.s2 <- data.frame(as.matrix(samps4.2))\nsamps.all.s2 <- samps.all.s2[, c(\"c.1.2.\", \"c.2.2.\", \"d.2.\")]\ndelta.2 <- colMeans(samps.all.s2)\ncov.2 <- var(samps.all.s2)\n\nFinally, we analyze the third (non-randomized) study comparing treatments A, B, and C:\n\nmodel1.spec <- textConnection(first.stage) \ndata3 <- with(ds4 %>% filter(studyid == 3), \n list(y = y,\n N = length(y), \n X = cbind(z1,z2), \n treat = treat,\n Ncovariate = 2, \n Ntreat = 3))\njags.m <- jags.model(model1.spec, data = data3, n.chains = 2, n.adapt = 100,\n quiet = TRUE)\nparams <- c(\"d\", \"c\") \nsamps4.3 <- coda.samples(jags.m, params, n.iter = 50)\nsamps.all.s3 <- data.frame(as.matrix(samps4.3))\n\nsamps.all.s3 <- samps.all.s3[, c(\"c.1.2.\", \"c.2.2.\", \"d.2.\", \"c.1.3.\", \n \"c.2.3.\", \"d.3.\")]\ndelta.3 <- colMeans(samps.all.s3)\ncov.3 <- var(samps.all.s3)\n\nThe corresponding treatment effect estimates are depicted below:\n\n\n\n\nTable 9.4: Treatment effect estimates.\n\n\nstudy\nB versus A\nC versus A\n\n\n\n\nstudy 1\n-2.968 (SE = 0.060 )\n\n\n\nstudy 2\n\n-1.007 (SE = 0.074 )\n\n\nstudy 3\n-2.930 (SE = 0.076 )\n-0.909 (SE = 0.061 )\n\n\n\n\n\n\n\n\nWe can now fit the second stage of the network meta-analysis. The corresponding JAGS model is specified below:\n\nsecond.stage <-\n\"model{\n \n #likelihood\n y1 ~ dmnorm(Mu1, Omega1)\n y2 ~ dmnorm(Mu2, Omega2)\n y3 ~ dmnorm(Mu3, Omega3*W)\n\n \n Omega1 <- inverse(cov.1)\n Omega2 <- inverse(cov.2)\n Omega3 <- inverse(cov.3)\n\n Mu1 <- c(gamma[,1], delta[2])\n Mu2 <- c(gamma[,2], delta[3]) \n Mu3 <- c(gamma[,1], delta[2],gamma[,2], delta[3])\n \n #parameters\n for(i in 1:2){\n gamma[i,1] ~ dnorm(0, 0.001)\n gamma[i,2] ~ dnorm(0, 0.001)\n }\n \n delta[1] <- 0\n delta[2] ~ dnorm(0, 0.001)\n delta[3] ~ dnorm(0, 0.001)\n \n}\n\"\n\nWe can fit as follows:\n\nmodel1.spec <- textConnection(second.stage) \ndata3 <- list(y1 = delta.1, y2 = delta.2, y3 = delta.3, \n cov.1 = cov.1, cov.2 = cov.2, cov.3 = cov.3, W = 0.5)\n\njags.m <- jags.model(model1.spec, data = data3, n.chains = 2, n.adapt = 50,\n quiet = TRUE)\nparams <- c(\"delta\", \"gamma\") \nsamps4.3 <- coda.samples(jags.m, params, n.iter = 50)\n\n\nsummary(samps4.3)\n\n\nIterations = 1:50\nThinning interval = 1 \nNumber of chains = 2 \nSample size per chain = 50 \n\n1. Empirical mean and standard deviation for each variable,\n plus standard error of the mean:\n\n Mean SD Naive SE Time-series SE\ndelta[1] 0.0000 0.00000 0.000000 0.000000\ndelta[2] -2.9334 0.06130 0.006130 0.006159\ndelta[3] -0.9916 0.05605 0.005605 0.005625\ngamma[1,1] -0.7941 0.04857 0.004857 0.004882\ngamma[2,1] 0.8471 0.08325 0.008325 0.009960\ngamma[1,2] -0.5136 0.06171 0.006171 0.005653\ngamma[2,2] 0.2771 0.07953 0.007953 0.009764\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\ndelta[1] 0.0000 0.0000 0.0000 0.0000 0.0000\ndelta[2] -3.0117 -2.9806 -2.9412 -2.9017 -2.8267\ndelta[3] -1.0883 -1.0346 -0.9940 -0.9529 -0.8877\ngamma[1,1] -0.8847 -0.8270 -0.7939 -0.7650 -0.6989\ngamma[2,1] 0.6984 0.8259 0.8578 0.8865 0.9501\ngamma[1,2] -0.6175 -0.5543 -0.5108 -0.4833 -0.4025\ngamma[2,2] 0.1453 0.2304 0.2803 0.3232 0.4166\n\n# calculate treatment effects\nsamples.all = data.frame(rbind(samps4.3[[1]], samps4.3[[2]]))\nnewpatient = c(1,2)\n\nmedian(\n samples.all$delta.2. + samples.all$gamma.1.1.*newpatient[1] +\n samples.all$gamma.2.1.*newpatient[2]\n)\n\n[1] -2.023269\n\nquantile(samples.all$delta.2.+samples.all$gamma.1.1.*newpatient[1]+\n samples.all$gamma.2.1.*newpatient[2]\n , probs = 0.025)\n\n 2.5% \n-2.304903 \n\nquantile(samples.all$delta.2.+samples.all$gamma.1.1.*newpatient[1]+\n samples.all$gamma.2.1.*newpatient[2]\n , probs = 0.975)\n\n 97.5% \n-1.770899" + "text": "9.2 Estimating heterogeous treatment effects in network meta-analysis\n\n9.2.1 Example of a continuous outcome\n\n9.2.1.1 Setup\nWe use again the bipd package to simulate a dataset:\n\nds3 <- generate_ipdnma_example(type = \"continuous\")\nhead(ds3)\n\n studyid treat z1 z2 y\n1 1 2 -0.7399958 0.4011151 9\n2 1 1 1.4632718 -1.2571879 11\n3 1 2 -0.6717355 -0.0582061 8\n4 1 2 0.3056371 -1.3880050 6\n5 1 1 -1.5699080 1.3314832 11\n6 1 2 -0.8608130 0.6028625 9\n\n\nLet us look into the data a bit in more detail:\n\n\n\n\n\n\nTable 9.3: The simulated dataset with a continuous outcome\n\n\n\n\n\n\n\n\n\n\n1\n(N=328)\n2\n(N=392)\n3\n(N=280)\nOverall\n(N=1000)\n\n\n\n\nz1\n\n\n\n\n\n\nMean (SD)\n0.0245 (0.992)\n-0.0143 (1.02)\n-0.0318 (0.980)\n-0.00650 (0.998)\n\n\nMedian [Min, Max]\n-0.0445 [-2.58, 2.53]\n-0.0218 [-3.11, 3.44]\n0.0391 [-3.62, 2.49]\n0.00453 [-3.62, 3.44]\n\n\nz2\n\n\n\n\n\n\nMean (SD)\n0.00244 (0.988)\n0.00503 (0.993)\n0.0784 (0.945)\n0.0247 (0.978)\n\n\nMedian [Min, Max]\n-0.0347 [-2.70, 3.00]\n0.0170 [-2.93, 2.50]\n0.0460 [-3.24, 2.42]\n0.0178 [-3.24, 3.00]\n\n\nstudyid\n\n\n\n\n\n\n1\n39 (11.9%)\n61 (15.6%)\n0 (0%)\n100 (10.0%)\n\n\n2\n42 (12.8%)\n58 (14.8%)\n0 (0%)\n100 (10.0%)\n\n\n3\n48 (14.6%)\n52 (13.3%)\n0 (0%)\n100 (10.0%)\n\n\n4\n53 (16.2%)\n0 (0%)\n47 (16.8%)\n100 (10.0%)\n\n\n5\n51 (15.5%)\n0 (0%)\n49 (17.5%)\n100 (10.0%)\n\n\n6\n0 (0%)\n57 (14.5%)\n43 (15.4%)\n100 (10.0%)\n\n\n7\n0 (0%)\n52 (13.3%)\n48 (17.1%)\n100 (10.0%)\n\n\n8\n28 (8.5%)\n41 (10.5%)\n31 (11.1%)\n100 (10.0%)\n\n\n9\n37 (11.3%)\n31 (7.9%)\n32 (11.4%)\n100 (10.0%)\n\n\n10\n30 (9.1%)\n40 (10.2%)\n30 (10.7%)\n100 (10.0%)\n\n\n\n\n\n\n\n\n\n\n9.2.1.2 Model fitting\nWe will use the model shown in Equation 16.8 in the book. In addition, we will use Bayesian LASSO to penalize the treatment-covariate interactions.\n\nipd3 <- with(ds3, ipdnma.model.onestage(y = y, study = studyid, treat = treat, \n X = cbind(z1, z2), \n response = \"normal\", \n shrinkage = \"laplace\", \n type = \"random\"))\nipd3$model.JAGS\n\nfunction () \n{\n for (i in 1:Np) {\n y[i] ~ dnorm(mu[i], sigma)\n mu[i] <- alpha[studyid[i]] + inprod(beta[], X[i, ]) + \n inprod(gamma[treat[i], ], X[i, ]) + d[studyid[i], \n treatment.arm[i]]\n }\n sigma ~ dgamma(0.001, 0.001)\n for (i in 1:Nstudies) {\n w[i, 1] <- 0\n d[i, 1] <- 0\n for (k in 2:na[i]) {\n d[i, k] ~ dnorm(mdelta[i, k], taudelta[i, k])\n mdelta[i, k] <- delta[t[i, k]] - delta[t[i, 1]] + \n sw[i, k]\n taudelta[i, k] <- tau * 2 * (k - 1)/k\n w[i, k] <- d[i, k] - delta[t[i, k]] + delta[t[i, \n 1]]\n sw[i, k] <- sum(w[i, 1:(k - 1)])/(k - 1)\n }\n }\n sd ~ dnorm(0, 1)\n T(0, )\n tau <- pow(sd, -2)\n delta[1] <- 0\n for (k in 2:Ntreat) {\n delta[k] ~ dnorm(0, 0.001)\n }\n for (j in 1:Nstudies) {\n alpha[j] ~ dnorm(0, 0.001)\n }\n for (k in 1:Ncovariate) {\n beta[k] ~ dnorm(0, 0.001)\n }\n lambda[1] <- 0\n lambda.inv[1] <- 0\n for (m in 2:Ntreat) {\n tt[m] <- lambda[m] * sigma\n lambda[m] <- pow(lambda.inv[m], -1)\n lambda.inv[m] ~ dunif(0, 5)\n }\n for (k in 1:Ncovariate) {\n gamma[1, k] <- 0\n for (m in 2:Ntreat) {\n gamma[m, k] ~ ddexp(0, tt[m])\n }\n }\n}\n<environment: 0x5556d0d738d8>\n\nsamples <- ipd.run(ipd3, n.chains = 2, n.iter = 20, \n pars.save = c(\"alpha\", \"beta\", \"delta\", \"sd\", \"gamma\"))\n\nCompiling model graph\n Resolving undeclared variables\n Allocating nodes\nGraph information:\n Observed stochastic nodes: 1000\n Unobserved stochastic nodes: 35\n Total graph size: 10141\n\nInitializing model\n\nsummary(samples)\n\n\nIterations = 2001:2020\nThinning interval = 1 \nNumber of chains = 2 \nSample size per chain = 20 \n\n1. Empirical mean and standard deviation for each variable,\n plus standard error of the mean:\n\n Mean SD Naive SE Time-series SE\nalpha[1] 10.9513 0.04430 0.007005 0.011211\nalpha[2] 8.0546 0.05247 0.008296 0.014855\nalpha[3] 10.4641 0.04407 0.006968 0.006754\nalpha[4] 9.6001 0.03953 0.006251 0.009281\nalpha[5] 12.9067 0.05125 0.008103 0.011440\nalpha[6] 13.2034 0.04994 0.007896 0.013482\nalpha[7] 7.4560 0.05069 0.008016 0.009881\nalpha[8] 11.1103 0.05053 0.007989 0.006599\nalpha[9] 10.1725 0.04351 0.006879 0.009881\nalpha[10] 9.2435 0.05389 0.008521 0.010676\nbeta[1] 0.1960 0.02048 0.003238 0.006558\nbeta[2] 0.2901 0.01747 0.002763 0.004416\ndelta[1] 0.0000 0.00000 0.000000 0.000000\ndelta[2] -2.9222 0.04655 0.007360 0.013167\ndelta[3] -1.1632 0.05075 0.008025 0.009775\ngamma[1,1] 0.0000 0.00000 0.000000 0.000000\ngamma[2,1] -0.5592 0.02856 0.004516 0.007406\ngamma[3,1] -0.3009 0.02584 0.004085 0.005820\ngamma[1,2] 0.0000 0.00000 0.000000 0.000000\ngamma[2,2] 0.5665 0.03096 0.004895 0.011632\ngamma[3,2] 0.4402 0.02537 0.004011 0.006529\nsd 0.1091 0.02682 0.004240 0.004268\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nalpha[1] 10.87592 10.92486 10.9432 10.9871 11.0311\nalpha[2] 7.95960 8.02352 8.0541 8.0954 8.1444\nalpha[3] 10.39161 10.43707 10.4611 10.4893 10.5395\nalpha[4] 9.54116 9.57556 9.5917 9.6213 9.6779\nalpha[5] 12.82662 12.87380 12.9012 12.9424 12.9958\nalpha[6] 13.12556 13.16652 13.2028 13.2304 13.2911\nalpha[7] 7.37369 7.41851 7.4599 7.4968 7.5248\nalpha[8] 11.03027 11.07663 11.1139 11.1397 11.2133\nalpha[9] 10.09934 10.15093 10.1636 10.2006 10.2633\nalpha[10] 9.15669 9.21080 9.2401 9.2838 9.3546\nbeta[1] 0.16187 0.18112 0.1972 0.2054 0.2386\nbeta[2] 0.25950 0.27918 0.2930 0.3026 0.3165\ndelta[1] 0.00000 0.00000 0.0000 0.0000 0.0000\ndelta[2] -3.00894 -2.94303 -2.9166 -2.8915 -2.8531\ndelta[3] -1.25176 -1.20392 -1.1548 -1.1198 -1.0930\ngamma[1,1] 0.00000 0.00000 0.0000 0.0000 0.0000\ngamma[2,1] -0.61841 -0.57430 -0.5570 -0.5424 -0.5117\ngamma[3,1] -0.35406 -0.32074 -0.2948 -0.2799 -0.2645\ngamma[1,2] 0.00000 0.00000 0.0000 0.0000 0.0000\ngamma[2,2] 0.51180 0.54303 0.5596 0.5952 0.6236\ngamma[3,2] 0.39076 0.42440 0.4402 0.4571 0.4813\nsd 0.06018 0.09209 0.1048 0.1211 0.1674\n\n\nAs before, we can use the treatment.effect() function of bipd to estimate relative effects for new patients.\n\ntreatment.effect(ipd3, samples, newpatient= c(1,2))\n\n$`treatment 2`\n 0.025 0.5 0.975 \n-2.543043 -2.331263 -2.197777 \n\n$`treatment 3`\n 0.025 0.5 0.975 \n-0.7342422 -0.5628488 -0.4501180 \n\n\nThis gives us the relative effects for all treatments versus the reference. To obtain relative effects between active treatments we need some more coding:\n\nsamples.all=data.frame(rbind(samples[[1]], samples[[2]]))\nnewpatient= c(1,2)\nnewpatient <- (newpatient - ipd3$scale_mean)/ipd3$scale_sd\n\nmedian(\n samples.all$delta.2.+samples.all$gamma.2.1.*\n newpatient[1]+samples.all$gamma.2.2.*newpatient[2]\n-\n (samples.all$delta.3.+samples.all$gamma.3.1.*newpatient[1]+\n samples.all$gamma.3.2.*newpatient[2])\n)\n\n[1] -1.764445\n\nquantile(samples.all$delta.2.+samples.all$gamma.2.1.*\n newpatient[1]+samples.all$gamma.2.2.*newpatient[2]\n -(samples.all$delta.3.+samples.all$gamma.3.1.*newpatient[1]+\n samples.all$gamma.3.2.*newpatient[2])\n , probs = 0.025)\n\n 2.5% \n-1.915623 \n\nquantile(samples.all$delta.2.+samples.all$gamma.2.1.*\n newpatient[1]+samples.all$gamma.2.2.*newpatient[2]\n -(samples.all$delta.3.+samples.all$gamma.3.1.*newpatient[1]+\n samples.all$gamma.3.2.*newpatient[2])\n , probs = 0.975)\n\n 97.5% \n-1.651104 \n\n\n\n\n\n9.2.2 Modeling patient-level relative effects using randomized and observational evidence for a network of treatments\nWe will now follow Chapter 16.3.5 from the book. In this analysis we will not use penalization, and we will assume fixed effects. For an example with penalization and random effects, see part 2 of this vignettte.\n\n9.2.2.1 Setup\nWe generate a very simple dataset of three studies comparing three treatments. We will assume 2 RCTs and 1 non-randomized trial:\n\nds4 <- generate_ipdnma_example(type = \"continuous\")\nds4 <- ds4 %>% filter(studyid %in% c(1,4,10)) %>%\n mutate(studyid = factor(studyid) %>%\n recode_factor(\n \"1\" = \"1\",\n \"4\" = \"2\",\n \"10\" = \"3\"),\n design = ifelse(studyid == \"3\", \"nrs\", \"rct\"))\n\nThe sample size is as follows:\n\n\n \n s1 s2 s3\n treat A: 42 44 36\n treat B: 58 0 33\n treat C: 0 56 31\n\n\n\n\n9.2.2.2 Model fitting\nWe will use the design-adjusted model, equation 16.9 in the book. We will fit a two-stage fixed effects meta-analysis and we will use a variance inflation factor. The code below is used to specify the analysis of each individual study. Briefly, in each study we adjust the treatment effect for the prognostic factors z1 and z2, as well as their interaction with treat.\n\nlibrary(rjags)\n\nLoading required package: coda\n\n\nLinked to JAGS 4.3.0\n\n\nLoaded modules: basemod,bugs\n\nfirst.stage <- \"\nmodel{\n\nfor (i in 1:N){\n y[i] ~ dnorm(mu[i], tau) \n mu[i] <- a + inprod(b[], X[i,]) + inprod(c[,treat[i]], X[i,]) + d[treat[i]] \n}\nsigma ~ dunif(0, 5)\ntau <- pow(sigma, -2)\n\na ~ dnorm(0, 0.001)\n\nfor(k in 1:Ncovariate){\n b[k] ~ dnorm(0,0.001)\n}\n\nfor(k in 1:Ncovariate){\n c[k,1] <- 0\n}\n\ntauGamma <- pow(sdGamma,-1)\nsdGamma ~ dunif(0, 5)\n\nfor(k in 1:Ncovariate){\n for(t in 2:Ntreat){\n c[k,t] ~ ddexp(0, tauGamma)\n }\n}\n\nd[1] <- 0\nfor(t in 2:Ntreat){\n d[t] ~ dnorm(0, 0.001)\n}\n}\"\n\nSubsequently, we estimate the relative treatment effects in the first (randomized) study comparing treatments A and B:\n\nmodel1.spec <- textConnection(first.stage) \ndata1 <- with(ds4 %>% filter(studyid == 1), \n list(y = y,\n N = length(y), \n X = cbind(z1,z2), \n treat = treat,\n Ncovariate = 2, \n Ntreat = 2))\njags.m <- jags.model(model1.spec, data = data1, n.chains = 2, n.adapt = 500,\n quiet = TRUE)\nparams <- c(\"d\", \"c\") \nsamps4.1 <- coda.samples(jags.m, params, n.iter = 50)\nsamps.all.s1 <- data.frame(as.matrix(samps4.1))\n\nsamps.all.s1 <- samps.all.s1[, c(\"c.1.2.\", \"c.2.2.\", \"d.2.\")]\ndelta.1 <- colMeans(samps.all.s1)\ncov.1 <- var(samps.all.s1)\n\nWe repeat the analysis for the second (randomized) study comparing treatments A and C:\n\nmodel1.spec <- textConnection(first.stage) \ndata2 <- with(ds4 %>% filter(studyid == 2), \n list(y = y,\n N = length(y), \n X = cbind(z1,z2), \n treat = ifelse(treat == 3, 2, treat),\n Ncovariate = 2, \n Ntreat = 2))\njags.m <- jags.model(model1.spec, data = data2, n.chains = 2, n.adapt = 100,\n quiet = TRUE)\nparams <- c(\"d\", \"c\") \nsamps4.2 <- coda.samples(jags.m, params, n.iter = 50)\nsamps.all.s2 <- data.frame(as.matrix(samps4.2))\nsamps.all.s2 <- samps.all.s2[, c(\"c.1.2.\", \"c.2.2.\", \"d.2.\")]\ndelta.2 <- colMeans(samps.all.s2)\ncov.2 <- var(samps.all.s2)\n\nFinally, we analyze the third (non-randomized) study comparing treatments A, B, and C:\n\nmodel1.spec <- textConnection(first.stage) \ndata3 <- with(ds4 %>% filter(studyid == 3), \n list(y = y,\n N = length(y), \n X = cbind(z1,z2), \n treat = treat,\n Ncovariate = 2, \n Ntreat = 3))\njags.m <- jags.model(model1.spec, data = data3, n.chains = 2, n.adapt = 100,\n quiet = TRUE)\nparams <- c(\"d\", \"c\") \nsamps4.3 <- coda.samples(jags.m, params, n.iter = 50)\nsamps.all.s3 <- data.frame(as.matrix(samps4.3))\n\nsamps.all.s3 <- samps.all.s3[, c(\"c.1.2.\", \"c.2.2.\", \"d.2.\", \"c.1.3.\", \n \"c.2.3.\", \"d.3.\")]\ndelta.3 <- colMeans(samps.all.s3)\ncov.3 <- var(samps.all.s3)\n\nThe corresponding treatment effect estimates are depicted below:\n\n\n\n\nTable 9.4: Treatment effect estimates.\n\n\nstudy\nB versus A\nC versus A\n\n\n\n\nstudy 1\n-2.888 (SE = 0.048 )\n\n\n\nstudy 2\n\n-1.106 (SE = 0.054 )\n\n\nstudy 3\n-2.753 (SE = 0.062 )\n-0.933 (SE = 0.065 )\n\n\n\n\n\n\n\n\nWe can now fit the second stage of the network meta-analysis. The corresponding JAGS model is specified below:\n\nsecond.stage <-\n\"model{\n \n #likelihood\n y1 ~ dmnorm(Mu1, Omega1)\n y2 ~ dmnorm(Mu2, Omega2)\n y3 ~ dmnorm(Mu3, Omega3*W)\n\n \n Omega1 <- inverse(cov.1)\n Omega2 <- inverse(cov.2)\n Omega3 <- inverse(cov.3)\n\n Mu1 <- c(gamma[,1], delta[2])\n Mu2 <- c(gamma[,2], delta[3]) \n Mu3 <- c(gamma[,1], delta[2],gamma[,2], delta[3])\n \n #parameters\n for(i in 1:2){\n gamma[i,1] ~ dnorm(0, 0.001)\n gamma[i,2] ~ dnorm(0, 0.001)\n }\n \n delta[1] <- 0\n delta[2] ~ dnorm(0, 0.001)\n delta[3] ~ dnorm(0, 0.001)\n \n}\n\"\n\nWe can fit as follows:\n\nmodel1.spec <- textConnection(second.stage) \ndata3 <- list(y1 = delta.1, y2 = delta.2, y3 = delta.3, \n cov.1 = cov.1, cov.2 = cov.2, cov.3 = cov.3, W = 0.5)\n\njags.m <- jags.model(model1.spec, data = data3, n.chains = 2, n.adapt = 50,\n quiet = TRUE)\nparams <- c(\"delta\", \"gamma\") \nsamps4.3 <- coda.samples(jags.m, params, n.iter = 50)\n\n\nsummary(samps4.3)\n\n\nIterations = 1:50\nThinning interval = 1 \nNumber of chains = 2 \nSample size per chain = 50 \n\n1. Empirical mean and standard deviation for each variable,\n plus standard error of the mean:\n\n Mean SD Naive SE Time-series SE\ndelta[1] 0.0000 0.00000 0.000000 0.000000\ndelta[2] -2.8927 0.04639 0.004639 0.005083\ndelta[3] -1.1114 0.04134 0.004134 0.004096\ngamma[1,1] -0.7580 0.06400 0.006400 0.006432\ngamma[2,1] 1.0019 0.06612 0.006612 0.006604\ngamma[1,2] -0.5170 0.04438 0.004438 0.005018\ngamma[2,2] 0.3689 0.03892 0.003892 0.003897\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\ndelta[1] 0.0000 0.0000 0.0000 0.0000 0.0000\ndelta[2] -2.9742 -2.9214 -2.8965 -2.8662 -2.7813\ndelta[3] -1.1950 -1.1387 -1.1114 -1.0835 -1.0364\ngamma[1,1] -0.8720 -0.7937 -0.7615 -0.7297 -0.6023\ngamma[2,1] 0.8906 0.9780 1.0053 1.0446 1.0874\ngamma[1,2] -0.6110 -0.5440 -0.5103 -0.4851 -0.4518\ngamma[2,2] 0.2981 0.3422 0.3660 0.3900 0.4520\n\n# calculate treatment effects\nsamples.all = data.frame(rbind(samps4.3[[1]], samps4.3[[2]]))\nnewpatient = c(1,2)\n\nmedian(\n samples.all$delta.2. + samples.all$gamma.1.1.*newpatient[1] +\n samples.all$gamma.2.1.*newpatient[2]\n)\n\n[1] -1.642012\n\nquantile(samples.all$delta.2.+samples.all$gamma.1.1.*newpatient[1]+\n samples.all$gamma.2.1.*newpatient[2]\n , probs = 0.025)\n\n 2.5% \n-1.89312 \n\nquantile(samples.all$delta.2.+samples.all$gamma.1.1.*newpatient[1]+\n samples.all$gamma.2.1.*newpatient[2]\n , probs = 0.975)\n\n 97.5% \n-1.432919" }, { "objectID": "chapter_16.html#version-info",