We demonstrate how to run one-stage IPD meta-analysis using this package. First, let’s generate sample IPD for illustration.
#devtools::install_github("MikeJSeo/bipd")
library(bipd)
##load in data
ds <- generate_ipdma_example(type = "continuous")
ds2 <- generate_ipdma_example(type = "binary")
head(ds)
#> studyid treat z1 z2 y
#> 1 1 0 0.3184625 -0.1616561 11
#> 2 1 1 -1.1532349 0.9836978 10
#> 3 1 1 1.3808782 0.5557304 8
#> 4 1 1 0.6179567 0.5831149 9
#> 5 1 0 -1.7102353 -0.7682330 10
#> 6 1 1 -1.0360358 -0.9662960 8The main function to set up the function for one-stage IPD meta-analysis is ipdma.model.onestage function. Refer to help(ipdma.model.onestage) for more details. Briefly to describe, “y” is the outcome of the study; “study” is a vector indicating which study the patient belongs to in a numerical sequence (i.e. 1, 2, 3, etc); “treat” is a vector indicating which treatment the patient was assigned to (i.e. 1 for treatment, 0 for placebo); “x” is a matrix of covariates for each patients; “response” is the outcome type, either “normal” or “binomial”.
Another important parameter is the “shrinkage” parameter. To specify IPD meta-analysis without shrinkage, we set shrinkage = “none”.
# continuous outcome
ipd <- with(ds, ipdma.model.onestage(y = y, study = studyid, treat = treat, X = cbind(z1, z2), response = "normal", shrinkage = "none"))To view the JAGS code that was used to run the model, we can run the following command. Note that “alpha” is the study intercept, “beta” is the coefficient for main effects of the covariates, “gamma” is the coefficient for effect modifier, and “delta” is the average treatment effect.
cat(ipd$code)
#> model {
#>
#> ########## IPD-MA model
#> for (i in 1:Np) {
#> y[i] ~ dnorm(mu[i], sigma)
#> mu[i] <- alpha[studyid[i]] + inprod(beta[], X[i,]) +
#> (1 - equals(treat[i],1)) * inprod(gamma[], X[i,]) + d[studyid[i],treat[i]]
#> }
#> sigma ~ dgamma(0.001, 0.001)
#>
#> #####treatment effect
#> for(j in 1:Nstudies){
#> d[j,1] <- 0
#> d[j,2] ~ dnorm(delta[2], tau)
#> }
#> sd ~ dnorm(0, 1)T(0,)
#> tau <- pow(sd, -2)
#>
#> ## prior distribution for the average treatment effect
#> delta[1] <- 0
#> delta[2] ~ dnorm(0, 0.001)
#>
#> ## prior distribution for the study intercept
#> for (j in 1:Nstudies){
#> alpha[j] ~ dnorm(0, 0.001)
#> }
#>
#> ## prior distribution for the main effect of the covariates
#> for(k in 1:Ncovariate){
#> beta[k] ~ dnorm(0, 0.001)
#> }
#> ## prior distribution for the effect modifiers under no shrinkage
#> for(k in 1:Ncovariate){
#> gamma[k] ~ dnorm(0, 0.001)
#> }
#> }Once the model is set up using ipdma.model.onestage function, we use ipd.run function to run the model. help(ipd.run) describes possible parameters to specify.
samples <- ipd.run(ipd, n.chains = 3, n.burnin = 500, n.iter = 5000)
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 600
#> Unobserved stochastic nodes: 19
#> Total graph size: 6034
#>
#> Initializing model
samples <- samples[,-3] #remove delta[1] which is 0
summary(samples)
#>
#> Iterations = 1501:6500
#> Thinning interval = 1
#> Number of chains = 3
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> alpha[1] 10.9504 0.05801 0.0004737 0.0008171
#> alpha[2] 8.0049 0.05795 0.0004732 0.0008772
#> alpha[4] 9.6025 0.05647 0.0004611 0.0008300
#> alpha[5] 12.9338 0.05905 0.0004821 0.0009143
#> alpha[6] 15.7748 0.05400 0.0004409 0.0007045
#> beta[1] 0.2149 0.02416 0.0001972 0.0003501
#> beta[2] 0.3201 0.02310 0.0001886 0.0003257
#> delta[1] 0.0000 0.00000 0.0000000 0.0000000
#> delta[2] -2.4939 0.19149 0.0015635 0.0015740
#> gamma[1] -0.5038 0.03292 0.0002688 0.0004952
#> gamma[2] 0.5971 0.03269 0.0002669 0.0004667
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> alpha[1] 10.8389 10.9109 10.9500 10.9890 11.0650
#> alpha[2] 7.8919 7.9659 8.0048 8.0432 8.1206
#> alpha[4] 9.4912 9.5640 9.6025 9.6409 9.7131
#> alpha[5] 12.8179 12.8945 12.9341 12.9737 13.0478
#> alpha[6] 15.6689 15.7384 15.7747 15.8111 15.8798
#> beta[1] 0.1680 0.1986 0.2150 0.2311 0.2629
#> beta[2] 0.2748 0.3048 0.3200 0.3354 0.3661
#> delta[1] 0.0000 0.0000 0.0000 0.0000 0.0000
#> delta[2] -2.8749 -2.6020 -2.4941 -2.3877 -2.1054
#> gamma[1] -0.5690 -0.5256 -0.5031 -0.4817 -0.4398
#> gamma[2] 0.5315 0.5759 0.5974 0.6188 0.6609
#plot(samples) #traceplot and posterior of parameters
#coda::gelman.plot(samples) #gelman diagnostic plotWe can find patient-specific treatment effect using the treatment.effect function. To do this we need to specify the covariate values for the patient that we want to predict patient-specific treatment effect.
For the second example, let’s use the same data, but include shrinkage (i.e. Bayesian LASSO) in the effect modifiers (i.e. treatment-covariate interactions). We can specify Bayesian LASSO by setting shrinkage = “laplace”. Lambda is the shrinkage parameter and we can set the prior for lambda using lambda.prior parameter. The default lambda prior for Bayesian LASSO is λ−1 ∼ dunif(0, 5).
ipd <- with(ds, ipdma.model.onestage(y = y, study = studyid, treat = treat, X = cbind(z1, z2), response = "normal", shrinkage = "laplace"))
samples <- ipd.run(ipd, pars.save = c("beta", "gamma", "delta", "lambda", "tt"), n.chains = 3, n.burnin = 500, n.iter = 5000)
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 600
#> Unobserved stochastic nodes: 20
#> Total graph size: 6039
#>
#> Initializing model
summary(samples)
#>
#> Iterations = 1501:6500
#> Thinning interval = 1
#> Number of chains = 3
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> beta[1] 0.2133 0.02429 0.0001984 0.0003945
#> beta[2] 0.3217 0.02309 0.0001886 0.0003453
#> delta[1] 0.0000 0.00000 0.0000000 0.0000000
#> delta[2] -2.4943 0.18845 0.0015387 0.0016286
#> gamma[1] -0.5006 0.03309 0.0002702 0.0005776
#> gamma[2] 0.5939 0.03304 0.0002697 0.0005613
#> lambda 0.3474 0.14910 0.0012174 0.0015013
#> tt 2.1767 0.93154 0.0076060 0.0095421
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> beta[1] 0.1655 0.1969 0.2134 0.2296 0.2613
#> beta[2] 0.2758 0.3067 0.3219 0.3372 0.3670
#> delta[1] 0.0000 0.0000 0.0000 0.0000 0.0000
#> delta[2] -2.8666 -2.6010 -2.4951 -2.3877 -2.1143
#> gamma[1] -0.5639 -0.5229 -0.5007 -0.4784 -0.4351
#> gamma[2] 0.5295 0.5719 0.5937 0.6157 0.6589
#> lambda 0.2039 0.2433 0.3019 0.4029 0.7548
#> tt 1.2478 1.5278 1.8958 2.5207 4.6985We can also use SSVS (stochastic search variable selection) by setting shrinkage = “SSVS”. This time let’s use the binomial dataset. “Ind” is the indicator for assigning a slab prior (instead of a spike prior) i.e. indicator for including a covariate. “eta” is the standard deviation of the slab prior.
ipd <- with(ds2, ipdma.model.onestage(y = y, study = studyid, treat = treat, X = cbind(w1, w2), response = "binomial", shrinkage = "SSVS"))
samples <- ipd.run(ipd, pars.save = c("beta", "gamma", "delta", "Ind", "eta"), n.chains = 3, n.burnin = 500, n.iter = 5000)
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 600
#> Unobserved stochastic nodes: 21
#> Total graph size: 6649
#>
#> Initializing model
summary(samples)
#>
#> Iterations = 1501:6500
#> Thinning interval = 1
#> Number of chains = 3
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> Ind[1] 0.17547 0.3804 0.0031058 0.010118
#> Ind[2] 0.18000 0.3842 0.0031370 0.009635
#> beta[1] 0.02624 0.1134 0.0009262 0.001736
#> beta[2] 0.06001 0.1155 0.0009434 0.001765
#> delta[1] 0.00000 0.0000 0.0000000 0.000000
#> delta[2] 0.15122 0.2770 0.0022621 0.008186
#> eta 1.98549 1.4831 0.0121091 0.048473
#> gamma[1] -0.01027 0.0963 0.0007863 0.001954
#> gamma[2] -0.01208 0.1004 0.0008195 0.001996
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> Ind[1] 0.0000 0.00000 0.000000 0.00000 1.0000
#> Ind[2] 0.0000 0.00000 0.000000 0.00000 1.0000
#> beta[1] -0.2013 -0.04812 0.025566 0.10074 0.2511
#> beta[2] -0.1651 -0.01693 0.058149 0.13558 0.2938
#> delta[1] 0.0000 0.00000 0.000000 0.00000 0.0000
#> delta[2] -0.3984 -0.02049 0.149945 0.32372 0.6946
#> eta 0.0420 0.60156 1.782812 3.20969 4.8020
#> gamma[1] -0.2364 -0.04667 -0.002276 0.02936 0.1859
#> gamma[2] -0.2465 -0.04778 -0.002343 0.03024 0.1814
treatment.effect(ipd, samples, newpatient = c(1,0.5)) # binary outcome reports odds ratio
#> 0.025 0.5 0.975
#> 0.6380697 1.1436818 2.0303106We now demonstrate how to run IPD network meta-analysis using this package.
##load in data
ds <- generate_ipdnma_example(type = "continuous")
ds2 <- generate_ipdnma_example(type = "binary")
head(ds)
#> studyid treat z1 z2 y
#> 1 1 1 2.3629286 -0.016977060 12
#> 2 1 1 -0.6077092 -0.109212517 11
#> 3 1 1 1.1667807 1.032182162 11
#> 4 1 1 0.6557452 0.249371112 11
#> 5 1 1 1.3332085 -0.335346796 11
#> 6 1 2 -0.3454660 0.004405287 8The main function to set up the function for one-stage IPD network meta-analysis is ipdnma.model.onestage function. The function is very similar to ipdma.model.onestage except that now we have number of treatments greater than 2. Consequently, “treat” parameter is defined differently i.e. 1 assigns baseline treatment and other treatments should be be assigned 2, 3, 4, etc.
# continuous outcome
ipd <- with(ds, ipdnma.model.onestage(y = y, study = studyid, treat = treat, X = cbind(z1, z2), response = "normal", shrinkage = "none"))
cat(ipd$code)
#> model {
#>
#> ########## IPD-NMA model
#> for (i in 1:Np) {
#> y[i] ~ dnorm(mu[i], sigma)
#> mu[i] <- alpha[studyid[i]] + inprod(beta[], X[i,]) +
#> inprod(gamma[treat[i],], X[i,]) + d[studyid[i],treatment.arm[i]]
#> }
#> sigma ~ dgamma(0.001, 0.001)
#>
#> #####treatment effect
#> for(i in 1:Nstudies){
#> w[i,1] <- 0
#> d[i,1] <- 0
#> for(k in 2:na[i]){
#> d[i,k] ~ dnorm(mdelta[i,k], taudelta[i,k])
#> mdelta[i,k] <- delta[t[i,k]] - delta[t[i,1]] + sw[i,k]
#> taudelta[i,k] <- tau * 2 * (k-1)/k
#> w[i,k] <- d[i,k] - delta[t[i,k]] + delta[t[i,1]]
#> sw[i,k] <- sum(w[i, 1:(k-1)]) / (k-1)
#> }
#> }
#> sd ~ dnorm(0, 1)T(0,)
#> tau <- pow(sd, -2)
#>
#> ## prior distribution for the average treatment effect
#> delta[1] <- 0
#> for(k in 2:Ntreat){
#> delta[k] ~ dnorm(0, 0.001)
#> }
#>
#> ## prior distribution for the study intercept
#> for (j in 1:Nstudies){
#> alpha[j] ~ dnorm(0, 0.001)
#> }
#>
#> ## prior distribution for the main effect of the covariates
#> for(k in 1:Ncovariate){
#> beta[k] ~ dnorm(0, 0.001)
#> }
#> ## prior distribution for the effect modifiers under no shrinkage
#> for(k in 1:Ncovariate){
#> gamma[1,k] <- 0
#> for(m in 2:Ntreat){
#> gamma[m,k] ~ dnorm(0, 0.001)
#> }
#> }
#> }
samples <- ipd.run(ipd, pars.save = c("beta", "gamma", "delta"), n.chains = 3, n.burnin = 500, n.iter = 5000)
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 1000
#> Unobserved stochastic nodes: 33
#> Total graph size: 10133
#>
#> Initializing model
summary(samples)
#>
#> Iterations = 1501:6500
#> Thinning interval = 1
#> Number of chains = 3
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> beta[1] 0.2180 0.01925 0.0001572 0.0003301
#> beta[2] 0.2985 0.02051 0.0001675 0.0003830
#> delta[1] 0.0000 0.00000 0.0000000 0.0000000
#> delta[2] -2.9178 0.05987 0.0004888 0.0008847
#> delta[3] -1.1356 0.06304 0.0005147 0.0009455
#> gamma[1,1] 0.0000 0.00000 0.0000000 0.0000000
#> gamma[2,1] -0.5866 0.02779 0.0002269 0.0004236
#> gamma[3,1] -0.3045 0.02968 0.0002424 0.0004253
#> gamma[1,2] 0.0000 0.00000 0.0000000 0.0000000
#> gamma[2,2] 0.5775 0.02882 0.0002354 0.0004893
#> gamma[3,2] 0.4081 0.02921 0.0002385 0.0004961
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> beta[1] 0.1801 0.2051 0.2181 0.2308 0.2554
#> beta[2] 0.2575 0.2848 0.2986 0.3124 0.3380
#> delta[1] 0.0000 0.0000 0.0000 0.0000 0.0000
#> delta[2] -3.0408 -2.9547 -2.9169 -2.8802 -2.7995
#> delta[3] -1.2614 -1.1750 -1.1352 -1.0962 -1.0098
#> gamma[1,1] 0.0000 0.0000 0.0000 0.0000 0.0000
#> gamma[2,1] -0.6403 -0.6054 -0.5867 -0.5684 -0.5314
#> gamma[3,1] -0.3622 -0.3248 -0.3046 -0.2842 -0.2465
#> gamma[1,2] 0.0000 0.0000 0.0000 0.0000 0.0000
#> gamma[2,2] 0.5217 0.5580 0.5777 0.5962 0.6348
#> gamma[3,2] 0.3514 0.3883 0.4081 0.4275 0.4658
treatment.effect(ipd, samples, newpatient = c(1,0.5))
#> $`treatment 2`
#> 0.025 0.5 0.975
#> -3.370114 -3.233811 -3.098889
#>
#> $`treatment 3`
#> 0.025 0.5 0.975
#> -1.389892 -1.250880 -1.112656We can apply shrinkage on the effect modifiers (treatment-covariate interactions) as before.
# SSVS
ipd <- with(ds, ipdnma.model.onestage(y = y, study = studyid, treat = treat, X = cbind(z1, z2), response = "normal", shrinkage = "SSVS"))
samples <- ipd.run(ipd, pars.save = c("beta", "gamma", "delta", "Ind", "eta"), n.chains = 3, n.burnin = 500, n.iter = 5000)
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 1000
#> Unobserved stochastic nodes: 38
#> Total graph size: 10155
#>
#> Initializing model
summary(samples)
#>
#> Iterations = 1501:6500
#> Thinning interval = 1
#> Number of chains = 3
#> Sample size per chain = 5000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> Ind[2,1] 0.9995 0.02160 0.0001763 0.0002697
#> Ind[3,1] 0.9733 0.16111 0.0013155 0.0063904
#> Ind[2,2] 0.9996 0.02000 0.0001633 0.0001822
#> Ind[3,2] 0.9897 0.10113 0.0008257 0.0028941
#> beta[1] 0.2177 0.01923 0.0001570 0.0003339
#> beta[2] 0.3000 0.02066 0.0001687 0.0003831
#> delta[1] 0.0000 0.00000 0.0000000 0.0000000
#> delta[2] -2.9157 0.05967 0.0004872 0.0008890
#> delta[3] -1.1344 0.06324 0.0005164 0.0009444
#> eta 0.8755 0.77748 0.0063481 0.0358468
#> gamma[1,1] 0.0000 0.00000 0.0000000 0.0000000
#> gamma[2,1] -0.5857 0.02758 0.0002252 0.0004343
#> gamma[3,1] -0.3038 0.02956 0.0002413 0.0004804
#> gamma[1,2] 0.0000 0.00000 0.0000000 0.0000000
#> gamma[2,2] 0.5752 0.02908 0.0002374 0.0004908
#> gamma[3,2] 0.4064 0.02903 0.0002370 0.0004799
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> Ind[2,1] 1.0000 1.0000 1.0000 1.0000 1.0000
#> Ind[3,1] 0.0000 1.0000 1.0000 1.0000 1.0000
#> Ind[2,2] 1.0000 1.0000 1.0000 1.0000 1.0000
#> Ind[3,2] 1.0000 1.0000 1.0000 1.0000 1.0000
#> beta[1] 0.1802 0.2047 0.2177 0.2305 0.2554
#> beta[2] 0.2596 0.2863 0.3001 0.3137 0.3406
#> delta[1] 0.0000 0.0000 0.0000 0.0000 0.0000
#> delta[2] -3.0317 -2.9528 -2.9157 -2.8779 -2.7958
#> delta[3] -1.2590 -1.1738 -1.1353 -1.0955 -1.0067
#> eta 0.3163 0.4800 0.6416 0.9152 3.8718
#> gamma[1,1] 0.0000 0.0000 0.0000 0.0000 0.0000
#> gamma[2,1] -0.6399 -0.6042 -0.5856 -0.5673 -0.5320
#> gamma[3,1] -0.3609 -0.3237 -0.3038 -0.2842 -0.2449
#> gamma[1,2] 0.0000 0.0000 0.0000 0.0000 0.0000
#> gamma[2,2] 0.5183 0.5556 0.5753 0.5945 0.6330
#> gamma[3,2] 0.3501 0.3867 0.4061 0.4258 0.4636
treatment.effect(ipd, samples, newpatient = c(1,0.5))
#> $`treatment 2`
#> 0.025 0.5 0.975
#> -3.363245 -3.232510 -3.098953
#>
#> $`treatment 3`
#> 0.025 0.5 0.975
#> -1.387533 -1.249575 -1.107381
# Bayesian LASSO
# ipd <- with(ds, ipdnma.model.onestage(y = y, study = studyid, treat = treat, X = cbind(z1, z2), response = "normal", shrinkage = "laplace", lambda.prior = list("dgamma",2,0.1)))
#samples <- ipd.run(ipd, pars.save = c("beta", "gamma", "delta", "lambda", "tt"), n.chains = 3, n.burnin = 500, n.iter = 5000)p.s. Note that in the network meta-analysis literature, “d” usually refers to average treatment effect and “delta” refers to study-specific treatment effect. In this R package, we have flipped around the two notations.