\title{studies to illustrate model checking methods}
\description{results from 20 hypothetical randomized clinical trials examining the effectiveness of a medication for treating some disease.}
\bold{nti}    \tab \code{numeric}   \tab number of patients in the treatment group \cr
\bold{nci}    \tab \code{numeric}   \tab number of patients in the control group \cr
\bold{xti}    \tab \code{numeric}   \tab number of patients in the treatment group with remission \cr
\bold{xci}    \tab \code{numeric}   \tab number of patients in the control group with remission \cr
\bold{dose}   \tab \code{numeric}   \tab dosage of the medication provided to patients in the treatment group (in milligrams per day)
the dataset was constructed for the purposes of illustrating the model checking and diagnostic methods described in viechtbauer (2021). the code below provides the results for many of the analyses and plots discussed in the book chapter.
viechtbauer, w. (2021). model checking in meta-analysis. in c. h. schmid, t. stijnen, & i. r. white (eds.), \emph{handbook of meta-analysis} (pp. 219-254). boca raton, fl: crc press. \verb{https://doi.org/10.1201/9781315119403}
dat <- escalc(measure="or", ai=xti, n1i=nti, ci=xci, n2i=nci, add=1/2, to="all", data=dat)
res.ce <- rma(yi, vi, data=dat, method="ce") # same as method="ee"
res.re <- rma(yi, vi, data=dat, method="dl")
res.mr <- rma(yi, vi, mods = ~ dose, data=dat, method="fe")
res.me <- rma(yi, vi, mods = ~ dose, data=dat, method="dl")
forest(dat$yi, dat$vi, psize=0.8, efac=0, xlim=c(-4,6), ylim=c(-3,23),
cex=1, width=c(5,5,5), xlab="log odds ratio (lnor)",
tmp <- regplot(res.me, xlim=c(0,250), ylim=c(-1,1.5), predlim=c(0,250), shade=false, digits=1,
xlab="dosage (mg per day)", psize="seinv", plim=c(na,5), bty="l", las=1,
lty=c("solid", "dashed"), label=true, labsize=0.8, offset=c(1,0.7))
res.sub <- rma(yi, vi, mods = ~ dose, data=dat, method="dl", subset=-6)
points(tmp$xi, tmp$yi, pch=21, cex=tmp$psize, col="black", bg="darkgray")
### number of standardized deleted residuals larger than +-1.96 in each model
### plot of the standardized deleted residuals for the re and me models
plot(na, na, xlim=c(1,20), ylim=c(-4,4), xlab="study", ylab="standardized (deleted) residual",
points(1:20, rstandard(res.re)$z, type="o", pch=19, col="gray70")
legend=c("standardized residuals","standardized deleted residuals"), bty="n")
points(1:20, rstandard(res.me)$z, type="o", pch=19, col="gray70")
baujat(res.ce, main="common-effects model", xlab="squared pearson residual", ylim=c(0,5), las=1)
baujat(res.me, main="mixed-effects model", ylim=c(0,2), las=1)
### gosh plots (skipped because this takes quite some time to run)
plot(res.gosh.me, het="tau2", out=6, breaks=50, adjust=0.6, las=1)
### plot of treatment dosage against the standardized residuals
plot(dat$dose, rstandard(res.me)$z, pch=19, xlab="dosage (mg per day)",
ylab="standardized residual", xlim=c(0,250), ylim=c(-2.5,2.5), las=1)
text(dat$dose[6], rstandard(res.me)$z[6], "6", pos=4, offset=0.4)
rma(yi, vi, mods = ~ dose + i(dose^2), data=dat, method="dl")
reslof <- rma(yi, vi, mods = ~ dose + factor(dose), data=dat, method="dl", btt=3:9)
### scatter plot to illustrate the lack-of-fit model
regplot(res.me, xlim=c(0,250), ylim=c(-1.0,1.5), xlab="dosage (mg per day)", ci=false,
predlim=c(0,250), psize=1, pch=19, col="gray60", digits=1, lwd=1, bty="l", las=1)
lines(dosages, fitted(reslof)[match(dosages, dat$dose)], type="o", pch=19, cex=2, lwd=2)
legend("bottomright", legend=c("linear model", "lack-of-fit model"), pch=c(na,19), col="black",
lty="solid", lwd=c(1,2), pt.cex=c(1,2), seg.len=4, bty="n")
### checking normality of the standardized deleted residuals
qqnorm(res.me, type="rstudent", main="standardized deleted residuals", pch=19, label="out",
lwd=2, pos=24, ylim=c(-4,3), lty=c("solid", "dotted"), las=1)
### checking normality of the random effects
sav <- qqnorm(ranef(res.me)$pred, main="blups of the random effects", cex=1, pch=19,
plot(na, na, xlim=c(1,20), ylim=c(0,0.21), xaxt="n", las=1, xlab="study", ylab="hat value")
title("hat values for the ce/re models")
legend=c("common-effects model", "random-effects model"), bty="n")
### heatmap of the hat matrix for the me model
cols <- colorramppalette(c("blue", "white", "red"))(101)
xlab="influence of the observed effect of study ...", ylab="on the fitted value of study ...",
axis(2, 1:20, labels=20:1, las=1, tick=false)
title("heatmap for the mixed-effects model")
### plot of leverages versus standardized residuals for the me model
las=1, xlab="leverage (hat value)", ylab="standardized deleted residual",
text(hatvalues(res.me)[ids] + c(0,0.013,0.010), rstudent(res.me)$z[ids] - c(0.18,0,0), ids)
title("leverage vs. standardized deleted residuals")
### plot of the cook's distances for the me model
plot(1:20, cooks.distance(res.me), ylim=c(0,1.6), type="o", pch=19, las=1, xaxt="n", yaxt="n",
### plot of the leave-one-out estimates of tau^2 for the me model
plot(1:20, x$inf$tau2.del,  ylim=c(0,0.15), type="o", pch=19, las=1, xaxt="n", xlab="study",
ylab=expression(paste("estimate of ", tau^2, " without the ", italic(i), "th study")))
### plot of the covariance ratios for the me model
### fit mixed-effects model without studies 3 and/or 6
rma(yi, vi, mods = ~ dose, data=dat, method="dl", subset=-3)
rma(yi, vi, mods = ~ dose, data=dat, method="dl", subset=-6)
rma(yi, vi, mods = ~ dose, data=dat, method="dl", subset=-c(3,6))
medicine, odds ratios, outliers, model checks