diff --git a/.nojekyll b/.nojekyll index 6995fa3..cc3f4c3 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -10f7f294 \ No newline at end of file +1c8ad7bd \ No newline at end of file diff --git a/authors.html b/authors.html new file mode 100644 index 0000000..1adc0a3 --- /dev/null +++ b/authors.html @@ -0,0 +1,481 @@ + + + + + + + + + +Comparative Effectiveness and Personalized Medicine Research Using Real-World Data - 11  Book Authors + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + +
+ +
+ + +
+ + + +
+ +
+
+

11  Book Authors

+
+ + +
+
Author
+
Affiliation
+ +
+

Thomas Debray

+
+
+

+ Smart Data Analysis and Statistics B.V. +

+
+
+ +
+ + + + +
+ + +
+ +
+

11.1 About this book

+

We gratefully acknowledge the contribution from the following authors:

+ + +
+ +
+ + +
+ + + + \ No newline at end of file diff --git a/chapter_03.html b/chapter_03.html index 5053883..b7b8d74 100644 --- a/chapter_03.html +++ b/chapter_03.html @@ -205,6 +205,12 @@ 10  Visualization and interpretation of individualized treatment rule results + +
  • +
    + + 11  Book Authors +
    • diff --git a/chapter_06.html b/chapter_06.html index c8c9c9d..9c0f2b9 100644 --- a/chapter_06.html +++ b/chapter_06.html @@ -205,6 +205,12 @@ 10  Visualization and interpretation of individualized treatment rule results + +
  • +
    + + 11  Book Authors +
    • diff --git a/chapter_07.html b/chapter_07.html index 3283eb6..7b31140 100644 --- a/chapter_07.html +++ b/chapter_07.html @@ -210,6 +210,12 @@ 10  Visualization and interpretation of individualized treatment rule results + +
  • +
    + + 11  Book Authors +
    • diff --git a/chapter_09.html b/chapter_09.html index 1b79f8e..519ce4d 100644 --- a/chapter_09.html +++ b/chapter_09.html @@ -209,6 +209,12 @@ 10  Visualization and interpretation of individualized treatment rule results + +
  • +
    + + 11  Book Authors +
    • diff --git a/chapter_10.html b/chapter_10.html index 5ffc48f..849a25d 100644 --- a/chapter_10.html +++ b/chapter_10.html @@ -229,6 +229,12 @@ 10  Visualization and interpretation of individualized treatment rule results + +
  • +
    + + 11  Book Authors +
    • @@ -854,16 +860,16 @@

    REM vs. STD -0.9 (0.73; 1.09) +0.9 (0.73; 1.08) TOCI vs. STD @@ -887,7 +893,7 @@

    REM vs. TOCI -1.02 (0.75; 1.27) +1.02 (0.75; 1.26) @@ -1336,12 +1342,12 @@

    sucra.uc$ranking.random 
     ADA160/80     FIL100     FIL200 GOL200/100      INF10       INF5    OZA0.92 
-0.26090909 0.18909091 0.40363636 0.59545455 0.60909091 0.74090909 0.77727273 
+0.24909091 0.16909091 0.41636364 0.64272727 0.58363636 0.75000000 0.77363636 
        PBO      TOF10      UPA45       UST6     VED300 
-0.02090909 0.40636364 0.98727273 0.40727273 0.60181818
    +0.01363636 0.41181818 0.97454545 0.38545455 0.63000000
    -

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

    +

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

    diff --git a/chapter_11.html b/chapter_11.html index d1c2406..6510eb2 100644 --- a/chapter_11.html +++ b/chapter_11.html @@ -228,6 +228,12 @@ 10  Visualization and interpretation of individualized treatment rule results + +
  • +
    + + 11  Book Authors +
    • @@ -337,19 +343,25 @@

    HOS, CRF, dialysis, DNOAP, smoking.ever, diabdur, wagner.class) -

    Briefly, 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:

    -
    +

    Briefly, 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 are summarized below:

    +
    --++++ + + @@ -358,193 +370,289 @@

    + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -553,10 +661,17 @@

    7.2.3 Hierarchical metaregression

    +

    We first investigate the event rate of patients receiving routine care:

    +
    +
    +

    +
    +
    +

    The forest plot above indicates that the baseline risk in the observational study from Morbach et al. is much higher than most trials.

    We fitted an HMR model to the available RWD and published AD:

    set.seed(2022)
@@ -584,23 +699,23 @@ 
               nr.adapt = 1000,         # Number of iteration for burnin 
            nr.thin = 1)             # Thinning rate
    -

    We 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.

    -
    +

    We 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 below 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.

    +
    -

    +

    Conflict 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.

    The 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.

    -
    +
    -

    +

    Posterior distribution of the weights for each study included in the HMR
    @@ -611,7 +726,7 @@

    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)
@@ -631,7 +746,7 @@ 
    Version info
    
 [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     
+[1] meta_6.5-0       jarbes_2.0.0     dplyr_1.1.4      ggplot2_3.4.4   
 [5] table1_1.4.3     kableExtra_1.3.4
 
 loaded via a namespace (and not attached):
diff --git a/chapter_11_files/figure-html/unnamed-chunk-6-1.png b/chapter_11_files/figure-html/unnamed-chunk-6-1.png
index 9ef195d..db5891c 100644
Binary files a/chapter_11_files/figure-html/unnamed-chunk-6-1.png and b/chapter_11_files/figure-html/unnamed-chunk-6-1.png differ
diff --git a/chapter_11_files/figure-html/unnamed-chunk-7-1.png b/chapter_11_files/figure-html/unnamed-chunk-7-1.png
index 9736f61..9ef195d 100644
Binary files a/chapter_11_files/figure-html/unnamed-chunk-7-1.png and b/chapter_11_files/figure-html/unnamed-chunk-7-1.png differ
diff --git a/chapter_11_files/figure-html/unnamed-chunk-8-1.png b/chapter_11_files/figure-html/unnamed-chunk-8-1.png
new file mode 100644
index 0000000..770b120
Binary files /dev/null and b/chapter_11_files/figure-html/unnamed-chunk-8-1.png differ
diff --git a/chapter_12.html b/chapter_12.html
index 6c092a8..a65de78 100644
--- a/chapter_12.html
+++ b/chapter_12.html
@@ -224,6 +224,12 @@
   
  10  Visualization and interpretation of individualized treatment rule results
    + +
  • +
    + + 11  Book Authors +
    • diff --git a/chapter_16.html b/chapter_16.html index 288cdd6..650e2d4 100644 --- a/chapter_16.html +++ b/chapter_16.html @@ -211,6 +211,12 @@ 10  Visualization and interpretation of individualized treatment rule results
    + +
  • +
    + + 11  Book Authors +
    • @@ -287,13 +293,13 @@ 

    head(ds)
    -
      studyid treat          z1         z2  y
-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
    +
      studyid treat          z1          z2  y
+1       1     0 -0.28302549 -1.48570136 11
+2       1     1  0.43147209  1.16605306 12
+3       1     0  0.10327470  0.66403984 11
+4       1     0 -0.06095517  0.75164484 11
+5       1     0  0.75349704 -1.38128389 11
+6       1     0 -0.62034698 -0.05138872 11

    The simulated dataset contains information on the following variables:

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

    - + @@ -336,15 +342,15 @@

    - - - + + + - - - + + + @@ -354,15 +360,15 @@

    - - - + + + - - - + + + @@ -372,38 +378,38 @@

    - - + + - + - - + + - - + + - - + + - + @@ -456,7 +462,7 @@

    +<environment: 0x55b0d68b35d8>

    We can fit the treatment effect model as follows:

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

    +delta[2] -3.0989 -2.6662 -2.1267 -1.6326 -0.3009 +gamma[1] -0.5142 -0.4825 -0.4631 -0.4352 -0.4037 +gamma[2] 0.4620 0.4877 0.4970 0.5150 0.5567 +sd 0.9874 1.3414 1.6396 1.9723 3.4098 @@ -519,10 +525,10 @@ 

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

    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).

    +

    This means that the predicted outcome for patient with covariate values z1 = 1 and z2 = 0.5 will differ by -2.37 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)
@@ -584,7 +590,7 @@ 
    round(treatment.effect(ipd, samples, newpatient = c(1,0.5)), 2)
    0.025   0.5 0.975 
--3.37 -2.28 -1.30
    +-3.58 -2.50 -1.73
    @@ -598,13 +604,13 @@ 

    ds2 <- generate_ipdma_example(type = "binary")
 head(ds2)
    -
      studyid treat         w1          w2 y
-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
    +
      studyid treat          w1          w2 y
+1       1     0 -0.02528346 -1.61296121 0
+2       1     0 -0.19520161  0.42643141 1
+3       1     0 -0.30588973  0.93049704 1
+4       1     0  1.31045922 -0.84908004 0
+5       1     0  0.11606923  0.01219508 1
+6       1     0  0.68655915 -0.53556333 0

    The simulated dataset contains information on the following variables:

    @@ -631,9 +637,9 @@

    +(N=279) +(N=321) @@ -647,15 +653,15 @@

    - - - + + + - - - + + + @@ -665,15 +671,15 @@

    - - - + + + - - - + + + @@ -683,38 +689,38 @@

    - - + + - - + + - - + + - - + + - - + + - - + + @@ -772,7 +778,7 @@

    +<environment: 0x55b0da740e10>
    @@ -791,37 +797,37 @@

    + 2.5% 25% 50% 75% 97.5% +alpha[1] -0.630580 -0.38064 -0.186586 0.05060 0.24979 +alpha[2] -0.194640 -0.02546 0.159278 0.28995 0.54055 +alpha[3] 0.085465 0.30191 0.593925 0.84833 1.24441 +alpha[4] -0.122869 0.26379 0.505358 0.77610 1.03754 +alpha[5] -0.226719 -0.01484 0.153880 0.41479 0.67717 +alpha[6] -0.590947 -0.35707 -0.095570 0.09083 0.41248 +beta[1] -0.215472 -0.04105 0.008329 0.08673 0.16396 +beta[2] -0.295459 -0.06998 -0.009676 0.06950 0.17173 +delta[1] 0.000000 0.00000 0.000000 0.00000 0.00000 +delta[2] -0.877437 -0.15798 0.044860 0.20244 0.82915 +gamma[1] -0.453274 -0.31043 -0.208522 -0.10696 0.09105 +gamma[2] 0.004015 0.10249 0.188586 0.29745 0.61430 +sd 0.398060 0.45208 0.766515 1.00511 1.83004

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

    @@ -829,10 +835,10 @@ 

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

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

    +

    In other words, the aforementioned patient 0.98 (95% Credibility Interval: 0.34 to 2.02)

    @@ -848,12 +854,12 @@

    head(ds3) 
      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
    +1 1 2 -0.5727157 -1.5014736 7 +2 1 1 0.3106872 -2.7650867 11 +3 1 2 0.7890523 0.1047426 8 +4 1 2 -0.1457122 -0.4751427 8 +5 1 1 1.1084918 -2.0759722 11 +6 1 2 -1.0116378 -0.4520275 8

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

    @@ -875,11 +881,11 @@

    +(N=316) +(N=380) +(N=304) @@ -894,17 +900,17 @@

    - - - - + + + + - - - - + + + + @@ -915,17 +921,17 @@

    - - - - + + + + - - - - + + + + @@ -936,72 +942,72 @@

    - - + + - - + + - - + + - + - + - + - + - - + + - - + + - - - + + + - - - + + + - - - + + + @@ -1072,7 +1078,7 @@

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

    + 2.5% 25% 50% 75% 97.5% +alpha[1] 11.00303 11.0337 11.0577 11.0871 11.1608 +alpha[2] 7.91623 7.9585 7.9846 8.0182 8.0567 +alpha[3] 10.48634 10.5245 10.5664 10.6016 10.6821 +alpha[4] 9.54121 9.6059 9.6497 9.6733 9.7368 +alpha[5] 12.82394 12.9052 12.9294 12.9553 13.0078 +alpha[6] 13.05422 13.1470 13.1754 13.2026 13.2400 +alpha[7] 7.22379 7.2851 7.3240 7.3497 7.3896 +alpha[8] 11.08466 11.1296 11.1653 11.1995 11.2475 +alpha[9] 10.13501 10.1735 10.2025 10.2481 10.3250 +alpha[10] 9.06006 9.1164 9.1446 9.1657 9.2147 +beta[1] 0.14794 0.1693 0.1781 0.1861 0.2142 +beta[2] 0.22820 0.2435 0.2648 0.2763 0.2992 +delta[1] 0.00000 0.0000 0.0000 0.0000 0.0000 +delta[2] -3.13502 -3.0701 -3.0388 -2.9886 -2.9411 +delta[3] -1.30947 -1.2295 -1.1741 -1.1484 -1.1043 +gamma[1,1] 0.00000 0.0000 0.0000 0.0000 0.0000 +gamma[2,1] -0.60764 -0.5915 -0.5862 -0.5629 -0.5379 +gamma[3,1] -0.33034 -0.3092 -0.2877 -0.2605 -0.2330 +gamma[1,2] 0.00000 0.0000 0.0000 0.0000 0.0000 +gamma[2,2] 0.61124 0.6375 0.6540 0.6713 0.7157 +gamma[3,2] 0.47080 0.4904 0.5105 0.5379 0.5727 +sd 0.07905 0.1279 0.1474 0.1741 0.2080

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

    @@ -1155,11 +1161,11 @@ 

    $`treatment 2`
     0.025       0.5     0.975 
--2.543043 -2.331263 -2.197777 
+-2.504054 -2.345491 -2.264910 
 
 $`treatment 3`
      0.025        0.5      0.975 
--0.7342422 -0.5628488 -0.4501180
    +-0.6621586 -0.5024203 -0.3532586

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

    @@ -1176,7 +1182,7 @@

    samples.all$gamma.3.2.*newpatient[2]) )
    -
    [1] -1.764445
    +
    [1] -1.868877
    quantile(samples.all$delta.2.+samples.all$gamma.2.1.*
            newpatient[1]+samples.all$gamma.2.2.*newpatient[2]
@@ -1185,7 +1191,7 @@ 
             , probs = 0.025)
         2.5% 
--1.915623
    +-2.008437 
    quantile(samples.all$delta.2.+samples.all$gamma.2.1.*
            newpatient[1]+samples.all$gamma.2.2.*newpatient[2]
@@ -1194,7 +1200,7 @@ 
             , probs = 0.975)
        97.5% 
--1.651104
    +-1.686934
    @@ -1220,9 +1226,9 @@ 

              
            s1 s2 s3
-  treat A: 42 44 36
-  treat B: 58  0 33
-  treat C:  0 56 31
    + treat A: 50 54 26 + treat B: 50 0 44 + treat C: 0 46 30 @@ -1351,18 +1357,18 @@

    - + - + - - + +
    Healing without amputation
    +(N=165)    		No healing without amputation
    +(N=95)    		 Overall
    (N=260)
    Age (years)
    Mean (SD)69.1 (10.9)68.5 (11.0) 68.9 (10.9)
    Median [Min, Max]70.0 [25.0, 90.0]69.0 [36.0, 91.0] 70.0 [25.0, 91.0]
    Diabetes duration (years)
    Mean (SD)15.9 (10.3)15.9 (11.2) 15.9 (10.6)
    Median [Min, Max]14.0 [1.00, 53.0]14.0 [0, 50.0] 14.0 [0, 53.0]
    Sex
    Female71 (43.0%)35 (36.8%) 106 (40.8%)
    Male94 (57.0%)60 (63.2%) 154 (59.2%)
    Ever smoker
    Yes97 (58.8%)57 (60.0%) 154 (59.2%)
    No68 (41.2%)38 (40.0%) 106 (40.8%)
    Diabetes type 2
    Yes150 (90.9%)79 (83.2%) 229 (88.1%)
    No15 (9.1%)16 (16.8%) 31 (11.9%)
    Peripheral arterial disease
    Yes82 (49.7%)66 (69.5%) 148 (56.9%)
    No83 (50.3%)29 (30.5%) 112 (43.1%)
    Neuropathy
    Yes144 (87.3%)80 (84.2%) 224 (86.2%)
    No21 (12.7%)15 (15.8%) 36 (13.8%)
    First ever lesion
    Yes70 (42.4%)44 (46.3%) 114 (43.8%)
    No95 (57.6%)51 (53.7%) 146 (56.2%)
    No continuous care
    Yes115 (69.7%)62 (65.3%) 177 (68.1%)
    No50 (30.3%)33 (34.7%) 83 (31.9%)
    Insulin dependent
    Yes109 (66.1%)65 (68.4%) 174 (66.9%)
    No56 (33.9%)30 (31.6%) 86 (33.1%)
    History of coronary events (CHD)
    Yes31 (18.8%)21 (22.1%) 52 (20.0%)
    No134 (81.2%)74 (77.9%) 208 (80.0%)
    History of stroke
    Yes36 (21.8%)19 (20.0%) 55 (21.2%)
    No129 (78.2%)76 (80.0%) 205 (78.8%)
    Charcot foot syndrome
    Yes28 (17.0%)24 (25.3%) 52 (20.0%)
    No137 (83.0%)71 (74.7%) 208 (80.0%)
    Dialysis
    Yes3 (1.8%)6 (6.3%) 9 (3.5%)
    No162 (98.2%)89 (93.7%) 251 (96.5%)
    DNOAP
    Yes19 (11.5%)10 (10.5%) 29 (11.2%)
    No146 (88.5%)85 (89.5%) 231 (88.8%)
    Wagner score
    1-2115 (69.7%)27 (28.4%) 142 (54.6%)
    3-4-550 (30.3%)68 (71.6%) 118 (45.4%)
    0
    -(N=296)    		1
    (N=304)    		1
    +(N=296)    		 Overall
    (N=600)
    Mean (SD)0.0522 (1.00)-0.0126 (0.990)0.0194 (0.995)-0.00620 (0.934)0.0216 (0.993)0.00751 (0.963)
    Median [Min, Max]0.0611 [-2.68, 3.33]-0.0719 [-2.29, 3.32]-0.0371 [-2.68, 3.33]-0.0769 [-2.23, 2.94]0.0183 [-2.66, 2.95]-0.0101 [-2.66, 2.95]
    z2Mean (SD)-0.0736 (0.929)0.00776 (0.978)-0.0324 (0.954)0.00678 (0.962)0.0300 (0.985)0.0182 (0.973)
    Median [Min, Max]-0.0719 [-2.79, 2.79]-0.0234 [-2.52, 2.77]-0.0617 [-2.79, 2.79]0.0547 [-2.77, 2.82]-0.0211 [-3.47, 2.65]0.00964 [-3.47, 2.82]
    studyid146 (15.5%)54 (17.8%)55 (18.1%)45 (15.2%) 100 (16.7%)
    252 (17.6%) 48 (15.8%)52 (17.6%) 100 (16.7%)
    352 (17.6%)48 (15.8%)54 (17.8%)46 (15.5%) 100 (16.7%)
    446 (15.5%)54 (17.8%)41 (13.5%)59 (19.9%) 100 (16.7%)
    550 (16.9%)50 (16.4%)56 (18.4%)44 (14.9%) 100 (16.7%)
    650 (16.9%) 50 (16.4%)50 (16.9%) 100 (16.7%)
    0
    -(N=305)    		 1
    -(N=295)    		 Overall
    (N=600)
    Mean (SD)-0.00720 (1.02)0.0526 (0.985)0.0222 (1.00)0.00838 (1.04)0.0736 (0.985)0.0432 (1.01)
    Median [Min, Max]-0.0492 [-2.73, 2.89]0.0449 [-3.23, 3.33]-0.0241 [-3.23, 3.33]0.00312 [-2.78, 2.65]-0.0240 [-2.58, 3.38]-0.0171 [-2.78, 3.38]
    w2
    Mean (SD)0.0465 (1.00)-0.0340 (1.01)0.00695 (1.01)0.00653 (1.03)-0.0209 (1.03)-0.00813 (1.03)
    Median [Min, Max]0.0894 [-2.37, 2.93]0.0399 [-2.72, 4.21]0.0748 [-2.72, 4.21]0.00908 [-3.03, 2.80]-0.00717 [-2.66, 2.78]-0.000170 [-3.03, 2.80]
    studyid
    147 (15.4%)53 (18.0%)46 (16.5%)54 (16.8%) 100 (16.7%)
    246 (15.1%)54 (18.3%)48 (17.2%)52 (16.2%) 100 (16.7%)
    356 (18.4%)44 (14.9%)43 (15.4%)57 (17.8%) 100 (16.7%)
    450 (16.4%)50 (16.9%)47 (16.8%)53 (16.5%) 100 (16.7%)
    551 (16.7%)49 (16.6%)46 (16.5%)54 (16.8%) 100 (16.7%)
    655 (18.0%)45 (15.3%)49 (17.6%)51 (15.9%) 100 (16.7%)
    1
    -(N=328)    		 2
    -(N=392)    		 3
    -(N=280)    		 Overall
    (N=1000)
    Mean (SD)0.0245 (0.992)-0.0143 (1.02)-0.0318 (0.980)-0.00650 (0.998)-0.0407 (1.01)0.0833 (0.961)0.0152 (0.954)0.0234 (0.976)
    Median [Min, Max]-0.0445 [-2.58, 2.53]-0.0218 [-3.11, 3.44]0.0391 [-3.62, 2.49]0.00453 [-3.62, 3.44]-0.0384 [-2.72, 3.02]0.000192 [-2.45, 3.06]-0.000770 [-2.15, 3.25]-0.00994 [-2.72, 3.25]
    z2
    Mean (SD)0.00244 (0.988)0.00503 (0.993)0.0784 (0.945)0.0247 (0.978)0.0794 (1.01)-0.0205 (1.07)-0.0459 (1.04)0.00337 (1.04)
    Median [Min, Max]-0.0347 [-2.70, 3.00]0.0170 [-2.93, 2.50]0.0460 [-3.24, 2.42]0.0178 [-3.24, 3.00]0.0562 [-2.77, 3.52]-0.131 [-4.00, 2.99]-0.00942 [-2.92, 3.50]0.00781 [-4.00, 3.52]
    studyid
    139 (11.9%)61 (15.6%)50 (15.8%)50 (13.2%) 0 (0%) 100 (10.0%)
    242 (12.8%)58 (14.8%)40 (12.7%)60 (15.8%) 0 (0%) 100 (10.0%)
    348 (14.6%)52 (13.3%)43 (13.6%)57 (15.0%) 0 (0%) 100 (10.0%)
    453 (16.2%)44 (13.9%) 0 (0%)47 (16.8%)56 (18.4%) 100 (10.0%)
    551 (15.5%)47 (14.9%) 0 (0%)49 (17.5%)53 (17.4%) 100 (10.0%)
    6 0 (0%)57 (14.5%)43 (15.4%)49 (12.9%)51 (16.8%) 100 (10.0%)
    7 0 (0%)52 (13.3%)48 (17.1%)52 (13.7%)48 (15.8%) 100 (10.0%)
    828 (8.5%)41 (10.5%)31 (11.1%)31 (9.8%)40 (10.5%)29 (9.5%) 100 (10.0%)
    937 (11.3%)31 (7.9%)32 (11.4%)34 (10.8%)33 (8.7%)33 (10.9%) 100 (10.0%)
    1030 (9.1%)40 (10.2%)30 (10.7%)27 (8.5%)39 (10.3%)34 (11.2%) 100 (10.0%)
    study 1-2.888 (SE = 0.048 )-2.940 (SE = 0.050 )
    study 2 -1.106 (SE = 0.054 )-1.116 (SE = 0.052 )
    study 3-2.753 (SE = 0.062 )-0.933 (SE = 0.065 )-2.893 (SE = 0.080 )-0.925 (SE = 0.079 )
    @@ -1428,23 +1434,23 @@

    +delta[2] -2.9906 -2.9470 -2.9228 -2.8901 -2.7640 +delta[3] -1.1445 -1.0895 -1.0531 -1.0312 -0.9631 +gamma[1,1] -0.9692 -0.9152 -0.8857 -0.8516 -0.7863 +gamma[2,1] 0.8398 0.8902 0.9270 0.9568 1.0009 +gamma[1,2] -0.6224 -0.5538 -0.5268 -0.4796 -0.4323 +gamma[2,2] 0.2464 0.3251 0.3673 0.4104 0.4743 

    # calculate  treatment effects
 samples.all = data.frame(rbind(samps4.3[[1]], samps4.3[[2]]))
@@ -1455,21 +1461,21 @@ 
        samples.all$gamma.2.1.*newpatient[2]
 )
    -
    [1] -1.642012
    +
    [1] -1.935452
    quantile(samples.all$delta.2.+samples.all$gamma.1.1.*newpatient[1]+
            samples.all$gamma.2.1.*newpatient[2]
          , probs = 0.025)
    -
        2.5% 
--1.89312
    +
         2.5% 
+-2.171893 
    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.432919
    +-1.708357

    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 6215e8d..31c7481 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/chapter_18.html b/chapter_18.html index 40c38fa..c9a8290 100644 --- a/chapter_18.html +++ b/chapter_18.html @@ -83,6 +83,7 @@ + @@ -229,6 +230,12 @@ 10  Visualization and interpretation of individualized treatment rule results + +
  • +
    + + 11  Book Authors +
    • @@ -1569,6 +1576,9 @@

    References

    + + 11  Book Authors +
    diff --git a/index.html b/index.html index c6643fe..905173d 100644 --- a/index.html +++ b/index.html @@ -171,6 +171,12 @@ 10  Visualization and interpretation of individualized treatment rule results + +
  • +
    + + 11  Book Authors +
    • diff --git a/search.json b/search.json index 48a2516..3c2e82e 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.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." + "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.1080 0.09763 0.0005637 0.0008465\nd.STD.TOCI -0.1122 0.08290 0.0004787 0.0008302\nsd.d 0.1129 0.08646 0.0004992 0.0016835\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nd.STD.REM -0.317179 -0.16198 -0.10373 -0.05056 0.08034\nd.STD.TOCI -0.259671 -0.16388 -0.12013 -0.06886 0.07900\nsd.d 0.004162 0.04682 0.09526 0.15878 0.32491\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.08)\n\n\nTOCI vs. STD\n0.89 (0.77; 1.08)\n\n\nREM vs. TOCI\n1.02 (0.75; 1.26)\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.24909091 0.16909091 0.41636364 0.64272727 0.58363636 0.75000000 0.77363636 \n PBO TOF10 UPA45 UST6 VED300 \n0.01363636 0.41181818 0.97454545 0.38545455 0.63000000 \n\n\nThese results indicate that 97.5% of the evaluated treatments are worse than UPA45." }, { "objectID": "chapter_10.html#version-info", @@ -256,14 +256,14 @@ "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" + "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 are summarized below:\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nHealing without amputation\n(N=165)\nNo healing without amputation\n(N=95)\nOverall\n(N=260)\n\n\n\n\nAge (years)\n\n\n\n\n\nMean (SD)\n69.1 (10.9)\n68.5 (11.0)\n68.9 (10.9)\n\n\nMedian [Min, Max]\n70.0 [25.0, 90.0]\n69.0 [36.0, 91.0]\n70.0 [25.0, 91.0]\n\n\nDiabetes duration (years)\n\n\n\n\n\nMean (SD)\n15.9 (10.3)\n15.9 (11.2)\n15.9 (10.6)\n\n\nMedian [Min, Max]\n14.0 [1.00, 53.0]\n14.0 [0, 50.0]\n14.0 [0, 53.0]\n\n\nSex\n\n\n\n\n\nFemale\n71 (43.0%)\n35 (36.8%)\n106 (40.8%)\n\n\nMale\n94 (57.0%)\n60 (63.2%)\n154 (59.2%)\n\n\nEver smoker\n\n\n\n\n\nYes\n97 (58.8%)\n57 (60.0%)\n154 (59.2%)\n\n\nNo\n68 (41.2%)\n38 (40.0%)\n106 (40.8%)\n\n\nDiabetes type 2\n\n\n\n\n\nYes\n150 (90.9%)\n79 (83.2%)\n229 (88.1%)\n\n\nNo\n15 (9.1%)\n16 (16.8%)\n31 (11.9%)\n\n\nPeripheral arterial disease\n\n\n\n\n\nYes\n82 (49.7%)\n66 (69.5%)\n148 (56.9%)\n\n\nNo\n83 (50.3%)\n29 (30.5%)\n112 (43.1%)\n\n\nNeuropathy\n\n\n\n\n\nYes\n144 (87.3%)\n80 (84.2%)\n224 (86.2%)\n\n\nNo\n21 (12.7%)\n15 (15.8%)\n36 (13.8%)\n\n\nFirst ever lesion\n\n\n\n\n\nYes\n70 (42.4%)\n44 (46.3%)\n114 (43.8%)\n\n\nNo\n95 (57.6%)\n51 (53.7%)\n146 (56.2%)\n\n\nNo continuous care\n\n\n\n\n\nYes\n115 (69.7%)\n62 (65.3%)\n177 (68.1%)\n\n\nNo\n50 (30.3%)\n33 (34.7%)\n83 (31.9%)\n\n\nInsulin dependent\n\n\n\n\n\nYes\n109 (66.1%)\n65 (68.4%)\n174 (66.9%)\n\n\nNo\n56 (33.9%)\n30 (31.6%)\n86 (33.1%)\n\n\nHistory of coronary events (CHD)\n\n\n\n\n\nYes\n31 (18.8%)\n21 (22.1%)\n52 (20.0%)\n\n\nNo\n134 (81.2%)\n74 (77.9%)\n208 (80.0%)\n\n\nHistory of stroke\n\n\n\n\n\nYes\n36 (21.8%)\n19 (20.0%)\n55 (21.2%)\n\n\nNo\n129 (78.2%)\n76 (80.0%)\n205 (78.8%)\n\n\nCharcot foot syndrome\n\n\n\n\n\nYes\n28 (17.0%)\n24 (25.3%)\n52 (20.0%)\n\n\nNo\n137 (83.0%)\n71 (74.7%)\n208 (80.0%)\n\n\nDialysis\n\n\n\n\n\nYes\n3 (1.8%)\n6 (6.3%)\n9 (3.5%)\n\n\nNo\n162 (98.2%)\n89 (93.7%)\n251 (96.5%)\n\n\nDNOAP\n\n\n\n\n\nYes\n19 (11.5%)\n10 (10.5%)\n29 (11.2%)\n\n\nNo\n146 (88.5%)\n85 (89.5%)\n231 (88.8%)\n\n\nWagner score\n\n\n\n\n\n1-2\n115 (69.7%)\n27 (28.4%)\n142 (54.6%)\n\n\n3-4-5\n50 (30.3%)\n68 (71.6%)\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 5), 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 first investigate the event rate of patients receiving routine care:\n\n\n\n\n\nThe forest plot above indicates that the baseline risk in the observational study from Morbach et al. is much higher than most trials.\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 below 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" + "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] meta_6.5-0 jarbes_2.0.0 dplyr_1.1.4 ggplot2_3.4.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", @@ -326,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.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)" + "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.28302549 -1.48570136 11\n2 1 1 0.43147209 1.16605306 12\n3 1 0 0.10327470 0.66403984 11\n4 1 0 -0.06095517 0.75164484 11\n5 1 0 0.75349704 -1.38128389 11\n6 1 0 -0.62034698 -0.05138872 11\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=304)\n1\n(N=296)\nOverall\n(N=600)\n\n\n\n\nz1\n\n\n\n\n\nMean (SD)\n-0.00620 (0.934)\n0.0216 (0.993)\n0.00751 (0.963)\n\n\nMedian [Min, Max]\n-0.0769 [-2.23, 2.94]\n0.0183 [-2.66, 2.95]\n-0.0101 [-2.66, 2.95]\n\n\nz2\n\n\n\n\n\nMean (SD)\n0.00678 (0.962)\n0.0300 (0.985)\n0.0182 (0.973)\n\n\nMedian [Min, Max]\n0.0547 [-2.77, 2.82]\n-0.0211 [-3.47, 2.65]\n0.00964 [-3.47, 2.82]\n\n\nstudyid\n\n\n\n\n\n1\n55 (18.1%)\n45 (15.2%)\n100 (16.7%)\n\n\n2\n48 (15.8%)\n52 (17.6%)\n100 (16.7%)\n\n\n3\n54 (17.8%)\n46 (15.5%)\n100 (16.7%)\n\n\n4\n41 (13.5%)\n59 (19.9%)\n100 (16.7%)\n\n\n5\n56 (18.4%)\n44 (14.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.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: 0x55b0d68b35d8>\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.9913 0.04797 0.007584 0.010169\nalpha[2] 8.0283 0.05584 0.008829 0.017094\nalpha[3] 10.5216 0.04548 0.007191 0.007981\nalpha[4] 9.6472 0.05175 0.008182 0.014972\nalpha[5] 12.8914 0.04693 0.007420 0.012130\nalpha[6] 15.7731 0.04680 0.007399 0.010816\nbeta[1] 0.1819 0.01775 0.002806 0.003999\nbeta[2] 0.2791 0.01416 0.002238 0.001900\ndelta[1] 0.0000 0.00000 0.000000 0.000000\ndelta[2] -1.9900 0.80661 0.127536 0.111701\ngamma[1] -0.4598 0.03163 0.005001 0.004587\ngamma[2] 0.5006 0.02435 0.003850 0.003581\nsd 1.7861 0.63814 0.100899 0.157533\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nalpha[1] 10.9147 10.9601 10.9761 11.0287 11.1054\nalpha[2] 7.9388 7.9786 8.0307 8.0716 8.1128\nalpha[3] 10.4509 10.4930 10.5139 10.5533 10.6089\nalpha[4] 9.5691 9.6156 9.6467 9.6748 9.7803\nalpha[5] 12.8063 12.8631 12.8965 12.9220 12.9808\nalpha[6] 15.7019 15.7422 15.7728 15.8139 15.8417\nbeta[1] 0.1564 0.1696 0.1791 0.1932 0.2213\nbeta[2] 0.2566 0.2696 0.2773 0.2858 0.3057\ndelta[1] 0.0000 0.0000 0.0000 0.0000 0.0000\ndelta[2] -3.0989 -2.6662 -2.1267 -1.6326 -0.3009\ngamma[1] -0.5142 -0.4825 -0.4631 -0.4352 -0.4037\ngamma[2] 0.4620 0.4877 0.4970 0.5150 0.5567\nsd 0.9874 1.3414 1.6396 1.9723 3.4098\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-3.38 -2.37 -0.59 \n\n\nThis means that the predicted outcome for patient with covariate values z1 = 1 and z2 = 0.5 will differ by -2.37 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.58 -2.50 -1.73 \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.02528346 -1.61296121 0\n2 1 0 -0.19520161 0.42643141 1\n3 1 0 -0.30588973 0.93049704 1\n4 1 0 1.31045922 -0.84908004 0\n5 1 0 0.11606923 0.01219508 1\n6 1 0 0.68655915 -0.53556333 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=279)\n1\n(N=321)\nOverall\n(N=600)\n\n\n\n\nw1\n\n\n\n\n\nMean (SD)\n0.00838 (1.04)\n0.0736 (0.985)\n0.0432 (1.01)\n\n\nMedian [Min, Max]\n0.00312 [-2.78, 2.65]\n-0.0240 [-2.58, 3.38]\n-0.0171 [-2.78, 3.38]\n\n\nw2\n\n\n\n\n\nMean (SD)\n0.00653 (1.03)\n-0.0209 (1.03)\n-0.00813 (1.03)\n\n\nMedian [Min, Max]\n0.00908 [-3.03, 2.80]\n-0.00717 [-2.66, 2.78]\n-0.000170 [-3.03, 2.80]\n\n\nstudyid\n\n\n\n\n\n1\n46 (16.5%)\n54 (16.8%)\n100 (16.7%)\n\n\n2\n48 (17.2%)\n52 (16.2%)\n100 (16.7%)\n\n\n3\n43 (15.4%)\n57 (17.8%)\n100 (16.7%)\n\n\n4\n47 (16.8%)\n53 (16.5%)\n100 (16.7%)\n\n\n5\n46 (16.5%)\n54 (16.8%)\n100 (16.7%)\n\n\n6\n49 (17.6%)\n51 (15.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: 0x55b0da740e10>\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.199134 0.2653 0.04195 0.04666\nalpha[2] 0.149541 0.2204 0.03484 0.03412\nalpha[3] 0.597879 0.3615 0.05716 0.12176\nalpha[4] 0.475872 0.3589 0.05675 0.15486\nalpha[5] 0.195731 0.2644 0.04180 0.07596\nalpha[6] -0.105562 0.3038 0.04803 0.06863\nbeta[1] 0.015167 0.1008 0.01594 0.02339\nbeta[2] -0.015164 0.1270 0.02008 0.04709\ndelta[1] 0.000000 0.0000 0.00000 0.00000\ndelta[2] 0.009803 0.3823 0.06045 0.04933\ngamma[1] -0.200628 0.1432 0.02264 0.02971\ngamma[2] 0.218744 0.1689 0.02671 0.06228\nsd 0.882181 0.4627 0.07316 0.15832\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nalpha[1] -0.630580 -0.38064 -0.186586 0.05060 0.24979\nalpha[2] -0.194640 -0.02546 0.159278 0.28995 0.54055\nalpha[3] 0.085465 0.30191 0.593925 0.84833 1.24441\nalpha[4] -0.122869 0.26379 0.505358 0.77610 1.03754\nalpha[5] -0.226719 -0.01484 0.153880 0.41479 0.67717\nalpha[6] -0.590947 -0.35707 -0.095570 0.09083 0.41248\nbeta[1] -0.215472 -0.04105 0.008329 0.08673 0.16396\nbeta[2] -0.295459 -0.06998 -0.009676 0.06950 0.17173\ndelta[1] 0.000000 0.00000 0.000000 0.00000 0.00000\ndelta[2] -0.877437 -0.15798 0.044860 0.20244 0.82915\ngamma[1] -0.453274 -0.31043 -0.208522 -0.10696 0.09105\ngamma[2] 0.004015 0.10249 0.188586 0.29745 0.61430\nsd 0.398060 0.45208 0.766515 1.00511 1.83004\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.34 0.98 2.02 \n\n\nIn other words, the aforementioned patient 0.98 (95% Credibility Interval: 0.34 to 2.02)" }, { "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.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" + "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.5727157 -1.5014736 7\n2 1 1 0.3106872 -2.7650867 11\n3 1 2 0.7890523 0.1047426 8\n4 1 2 -0.1457122 -0.4751427 8\n5 1 1 1.1084918 -2.0759722 11\n6 1 2 -1.0116378 -0.4520275 8\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=316)\n2\n(N=380)\n3\n(N=304)\nOverall\n(N=1000)\n\n\n\n\nz1\n\n\n\n\n\n\nMean (SD)\n-0.0407 (1.01)\n0.0833 (0.961)\n0.0152 (0.954)\n0.0234 (0.976)\n\n\nMedian [Min, Max]\n-0.0384 [-2.72, 3.02]\n0.000192 [-2.45, 3.06]\n-0.000770 [-2.15, 3.25]\n-0.00994 [-2.72, 3.25]\n\n\nz2\n\n\n\n\n\n\nMean (SD)\n0.0794 (1.01)\n-0.0205 (1.07)\n-0.0459 (1.04)\n0.00337 (1.04)\n\n\nMedian [Min, Max]\n0.0562 [-2.77, 3.52]\n-0.131 [-4.00, 2.99]\n-0.00942 [-2.92, 3.50]\n0.00781 [-4.00, 3.52]\n\n\nstudyid\n\n\n\n\n\n\n1\n50 (15.8%)\n50 (13.2%)\n0 (0%)\n100 (10.0%)\n\n\n2\n40 (12.7%)\n60 (15.8%)\n0 (0%)\n100 (10.0%)\n\n\n3\n43 (13.6%)\n57 (15.0%)\n0 (0%)\n100 (10.0%)\n\n\n4\n44 (13.9%)\n0 (0%)\n56 (18.4%)\n100 (10.0%)\n\n\n5\n47 (14.9%)\n0 (0%)\n53 (17.4%)\n100 (10.0%)\n\n\n6\n0 (0%)\n49 (12.9%)\n51 (16.8%)\n100 (10.0%)\n\n\n7\n0 (0%)\n52 (13.7%)\n48 (15.8%)\n100 (10.0%)\n\n\n8\n31 (9.8%)\n40 (10.5%)\n29 (9.5%)\n100 (10.0%)\n\n\n9\n34 (10.8%)\n33 (8.7%)\n33 (10.9%)\n100 (10.0%)\n\n\n10\n27 (8.5%)\n39 (10.3%)\n34 (11.2%)\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: 0x55b0daa41f08>\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.0640 0.04419 0.006988 0.009942\nalpha[2] 7.9870 0.04487 0.007094 0.007783\nalpha[3] 10.5690 0.05369 0.008489 0.017340\nalpha[4] 9.6455 0.05225 0.008262 0.018824\nalpha[5] 12.9266 0.04759 0.007525 0.009543\nalpha[6] 13.1716 0.04975 0.007866 0.008900\nalpha[7] 7.3176 0.04695 0.007424 0.011618\nalpha[8] 11.1647 0.04940 0.007810 0.015052\nalpha[9] 10.2113 0.05361 0.008477 0.012564\nalpha[10] 9.1426 0.04001 0.006326 0.008012\nbeta[1] 0.1782 0.01628 0.002574 0.004144\nbeta[2] 0.2623 0.02125 0.003359 0.006661\ndelta[1] 0.0000 0.00000 0.000000 0.000000\ndelta[2] -3.0360 0.05504 0.008703 0.012784\ndelta[3] -1.1866 0.05923 0.009365 0.012837\ngamma[1,1] 0.0000 0.00000 0.000000 0.000000\ngamma[2,1] -0.5789 0.02066 0.003267 0.003929\ngamma[3,1] -0.2839 0.02922 0.004620 0.007762\ngamma[1,2] 0.0000 0.00000 0.000000 0.000000\ngamma[2,2] 0.6560 0.02739 0.004331 0.007318\ngamma[3,2] 0.5148 0.03191 0.005046 0.008789\nsd 0.1470 0.03636 0.005748 0.009147\n\n2. Quantiles for each variable:\n\n 2.5% 25% 50% 75% 97.5%\nalpha[1] 11.00303 11.0337 11.0577 11.0871 11.1608\nalpha[2] 7.91623 7.9585 7.9846 8.0182 8.0567\nalpha[3] 10.48634 10.5245 10.5664 10.6016 10.6821\nalpha[4] 9.54121 9.6059 9.6497 9.6733 9.7368\nalpha[5] 12.82394 12.9052 12.9294 12.9553 13.0078\nalpha[6] 13.05422 13.1470 13.1754 13.2026 13.2400\nalpha[7] 7.22379 7.2851 7.3240 7.3497 7.3896\nalpha[8] 11.08466 11.1296 11.1653 11.1995 11.2475\nalpha[9] 10.13501 10.1735 10.2025 10.2481 10.3250\nalpha[10] 9.06006 9.1164 9.1446 9.1657 9.2147\nbeta[1] 0.14794 0.1693 0.1781 0.1861 0.2142\nbeta[2] 0.22820 0.2435 0.2648 0.2763 0.2992\ndelta[1] 0.00000 0.0000 0.0000 0.0000 0.0000\ndelta[2] -3.13502 -3.0701 -3.0388 -2.9886 -2.9411\ndelta[3] -1.30947 -1.2295 -1.1741 -1.1484 -1.1043\ngamma[1,1] 0.00000 0.0000 0.0000 0.0000 0.0000\ngamma[2,1] -0.60764 -0.5915 -0.5862 -0.5629 -0.5379\ngamma[3,1] -0.33034 -0.3092 -0.2877 -0.2605 -0.2330\ngamma[1,2] 0.00000 0.0000 0.0000 0.0000 0.0000\ngamma[2,2] 0.61124 0.6375 0.6540 0.6713 0.7157\ngamma[3,2] 0.47080 0.4904 0.5105 0.5379 0.5727\nsd 0.07905 0.1279 0.1474 0.1741 0.2080\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.504054 -2.345491 -2.264910 \n\n$`treatment 3`\n 0.025 0.5 0.975 \n-0.6621586 -0.5024203 -0.3532586 \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.868877\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-2.008437 \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.686934 \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: 50 54 26\n treat B: 50 0 44\n treat C: 0 46 30\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.940 (SE = 0.050 )\n\n\n\nstudy 2\n\n-1.116 (SE = 0.052 )\n\n\nstudy 3\n-2.893 (SE = 0.080 )\n-0.925 (SE = 0.079 )\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.9153 0.05449 0.005449 0.005475\ndelta[3] -1.0574 0.04682 0.004682 0.006849\ngamma[1,1] -0.8822 0.04775 0.004775 0.004785\ngamma[2,1] 0.9280 0.05845 0.005845 0.005868\ngamma[1,2] -0.5210 0.05102 0.005102 0.006205\ngamma[2,2] 0.3644 0.07441 0.007441 0.007478\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.9906 -2.9470 -2.9228 -2.8901 -2.7640\ndelta[3] -1.1445 -1.0895 -1.0531 -1.0312 -0.9631\ngamma[1,1] -0.9692 -0.9152 -0.8857 -0.8516 -0.7863\ngamma[2,1] 0.8398 0.8902 0.9270 0.9568 1.0009\ngamma[1,2] -0.6224 -0.5538 -0.5268 -0.4796 -0.4323\ngamma[2,2] 0.2464 0.3251 0.3673 0.4104 0.4743\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.935452\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.171893 \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.708357" }, { "objectID": "chapter_16.html#version-info", @@ -397,5 +397,12 @@ "title": "10  Visualization and interpretation of individualized treatment rule results", "section": "References", "text": "References\n\n\n\n\nYadlowsky, Steve, Fabio Pellegrini, Federica Lionetto, Stefan Braune, and Lu Tian. 2020. “Estimation and Validation of Ratio-Based Conditional Average Treatment Effects Using Observational Data.” Journal of the American Statistical Association 116 (533): 335–52. https://doi.org/10.1080/01621459.2020.1772080.\n\n\nZhao, Lihui, Lu Tian, Tianxi Cai, Brian Claggett, and L. J. Wei. 2013. “Effectively Selecting a Target Population for a Future Comparative Study.” Journal of the American Statistical Association 108 (502): 527–39. https://doi.org/10.1080/01621459.2013.770705." + }, + { + "objectID": "authors.html#about-this-book", + "href": "authors.html#about-this-book", + "title": "11  Book Authors", + "section": "11.1 About this book", + "text": "11.1 About this book\nWe gratefully acknowledge the contribution from the following authors:" } ] \ No newline at end of file