# NOT RUN {
### REPRODUCE THE RESULTS IN MUSEKIWA ET AL (2012), TABLES 3 AND 4
# INDEPENDENT RANDOM EFFECTS, NO WITHIN-STUDY CORRELATION (MODEL 1)
mod1 <- mixmeta(logOR~0+factor(time), S=logORvar, random=~0+factor(time)|study,
bscov="diag", data=gliomas)
print(summary(mod1), digits=2, report="var")
# COMPOUND-SYMMETRY RANDOM EFFECTS, NO WITHIN-STUDY CORRELATION (MODEL 2)
# NB: THIS REQUIRES A TWO-LEVEL MODEL WITH THE INNER-LEVEL VARIANCE FIXED TO 0
unit <- factor(seq(nrow(gliomas)))
mod2 <- mixmeta(logOR~0+factor(time), S=logORvar, random=~1|study/unit,
bscov=c("unstr","fixed"), data=gliomas, control=list(Psifix=list(unit=0)))
print(summary(mod2), digits=2, report="var")
# BUILD THE HETEROGENEOUS AR1 WITHIN-STUDY ERRORS (CORRELATION AT 0.61)
cormat <- 0.61^abs(col(matrix(1,4,4)) - row(col(matrix(1,4,4))))
addS <- lapply(split(sqrt(gliomas$logORvar), gliomas$study), inputcov, cormat)
addS <- lapply(addS, function(x) x[apply(!is.na(x),1,any), apply(!is.na(x),2,any)])
# INDEPENDENT RANDOM EFFECTS, HAR1 WITHIN-STUDY CORRELATION (MODEL 4)
mod4 <- mixmeta(logOR~0+factor(time), random=~0+factor(time)|study,
bscov="diag", data=gliomas, control=list(addSlist=addS))
print(summary(mod4), digits=2, report="var")
# UNSTRUCTURED RANDOM EFFECTS, HAR1 WITHIN-STUDY CORRELATION (MODEL 6)
mod6 <- update(mod4, bscov="unstr")
print(summary(mod6), digits=2, report="var")
# COMPARE THE FIT
AIC(mod1, mod2, mod4, mod6)
# }
# NOT RUN {
### MORE FLEXIBLE MODELLING OF RANDOM EFFECTS
# RE-RUN BEST FITTING MODEL WITH ML (ALLOWS TESTING OF FIXED EFFECTS)
mod4ml <- update(mod4, method="ml")
# RANDOM-SLOPE MODEL WITH TIME AS CONTINUOUS (CENTERED IN random)
modnew <- mixmeta(logOR~time, random=~I(time-15)|study, bscov="diag",
method="ml", data=gliomas, control=list(addSlist=addS, maxiter=200))
print(summary(modnew), digits=2, report="var")
# COMPARE
AIC(mod4ml, modnew)
# }
# NOT RUN {
### SEE help(dbs) FOR A COMPLEMENTARY EXAMPLE
# }
Run the code above in your browser using DataLab