[R] meta analysis with repeated measure-designs?
Mike Cheung
mikewlcheung at gmail.com
Fri Jun 18 11:21:37 CEST 2010
Dear Gerrit,
Sorry. There was an error in my previous code. As a record, the
followings are the revised code.
#### Robust SE based on Hedges et al., (2010) Eq. 6 on Research
Synthesis Methods
#### rma.obj: object fitted by metafor()
#### cluster: indicator for clusters of studies
robustSE <- function(rma.obj, cluster=NULL, CI=.95) {
# m: no. of clusters; assumed independent if not specified
# rma.obj$not.na: complete cases
if (is.null(cluster)) {
m=length(rma.obj$X[rma.obj$not.na,1])
} else {
m=nlevels(unique(as.factor(cluster[rma.obj$not.na])))
}
res2 <- diag(residuals(rma.obj)^2)
X <- rma.obj$X
b <- rma.obj$b
W <- diag(1/(rma.obj$vi+rma.obj$tau2)) # Use vi+tau2
meat <- t(X) %*% W %*% res2 %*% W %*% X # W is symmetric
bread <- solve( t(X) %*% W %*% X)
V.R <- bread %*% meat %*% bread # Robust sampling covariance matrix
p <- length(b) # no. of predictors
including intercept
se <- sqrt( diag(V.R)*m/(m-p) ) # small sample adjustment (Eq.7)
tval <- b/se
pval <- 2*(1-pt(abs(tval),df=(m-p)))
crit <- qt( (1-CI)/2, df=(m-p), lower.tail=FALSE )
ci.lb <- b-crit*se
ci.ub <- b+crit*se
data.frame(estimate=b, se=se, tval=tval, pval=pval, ci.lb=ci.lb, ci.ub=ci.ub)
}
library(metafor)
data(dat.bcg)
### calculate log relative risks and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
dat <- cbind(dat.bcg, dat)
### random-effects model
fit1 <- rma(yi, vi, data=dat, method="DL")
summary(fit1)
estimate se zval pval ci.lb ci.ub
-0.7141 0.1787 -3.9952 <.0001 -1.0644 -0.3638 ***
robustSE(fit1)
estimate se tval pval ci.lb ci.ub
intrcpt -0.7141172 0.1791445 -3.986265 0.001805797 -1.104439 -0.323795
### mixed-effects model with two moderators (absolute latitude and
publication year)
fit2 <- rma(yi, vi, mods=cbind(ablat, year), data=dat, method="DL")
summary(fit2)
estimate se zval pval ci.lb ci.ub
intrcpt -1.2798 25.7550 -0.0497 0.9604 -51.7586 49.1990
ablat -0.0288 0.0090 -3.2035 0.0014 -0.0464 -0.0112 **
year 0.0008 0.0130 0.0594 0.9526 -0.0247 0.0262
robustSE(fit2)
estimate se tval pval
ci.lb ci.ub
intrcpt -1.2797914381 22.860353022 -0.05598301 0.956458098
-52.21583218 49.65624930
ablat -0.0287644840 0.007212163 -3.98832970 0.002566210
-0.04483418 -0.01269478
year 0.0007720798 0.011550188 0.06684565 0.948022174
-0.02496334 0.02650750
--
---------------------------------------------------------------------
Mike W.L. Cheung Phone: (65) 6516-3702
Department of Psychology Fax: (65) 6773-1843
National University of Singapore
http://courses.nus.edu.sg/course/psycwlm/internet/
---------------------------------------------------------------------
On Wed, Jun 16, 2010 at 4:47 PM, Mike Cheung <mikewlcheung at gmail.com> wrote:
> Dear Gerrit,
>
> If the correlations of the dependent effect sizes are unknown, one
> approach is to conduct the meta-analysis by assuming that the effect
> sizes are independent. A robust standard error is then calculated to
> adjust for the dependence. You may refer to Hedges et. al., (2010) for
> more information. I have coded it here for reference.
>
> Hedges, L. V., Tipton, E., & Johnson, M. C. (2010). Robust variance
> estimation in meta-regression with dependent effect size estimates.
> Research Synthesis Methods, 1(1), 39-65. doi:10.1002/jrsm.5
>
> Regards,
> Mike
> --
> ---------------------------------------------------------------------
> Mike W.L. Cheung Phone: (65) 6516-3702
> Department of Psychology Fax: (65) 6773-1843
> National University of Singapore
> http://courses.nus.edu.sg/course/psycwlm/internet/
> ---------------------------------------------------------------------
>
> library(metafor)
>
> #### Robust SE based on Hedges et al., (2010) Eq. 6
> #### rma.obj: object fitted by metafor()
> #### cluster: indicator for clusters of studies
> robustSE <- function(rma.obj, cluster=NULL, CI=.95) {
> # m: no. of clusters; assumed independent if not specified
> if (is.null(cluster)) {
> m=nrow(rma.obj$X)
> } else {
> m=nlevels(unique(as.factor(cluster)))
> }
> res2 <- diag(residuals(rma.obj)^2)
> X <- rma.obj$X
> b <- rma.obj$b
> W <- diag(1/(rma.obj$vi+rma.obj$tau2)) # Use vi+tau2
> meat <- t(X) %*% W %*% res2 %*% W %*% X # W is symmetric
> bread <- solve( t(X) %*% W %*% X)
> V.R <- bread %*% meat %*% bread # Robust sampling covariance matrix
> p <- length(b) # no. of predictors
> including intercept
> se.R <- sqrt( diag(V.R)*m/(m-p) ) # small sample adjustment (Eq.7)
> t.R <- b/se.R
> p.R <- 2*(1-pt(abs(t.R),df=(m-p)))
> crit <- qt( 1-CI/2, df=(m-p) )
> ci.lb <- b-crit*se.R
> ci.ub <- b+crit*se.R
> data.frame(estimate=b, se.R=se.R, t.R=t.R, p.R=p.R, ci.lb=ci.lb, ci.ub=ci.ub)
> }
>
>
> ## Example: calculate log relative risks and corresponding sampling variances
> data(dat.bcg)
> dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
> dat <- cbind(dat.bcg, dat)
>
> ## random-effects model
> fit1 <- rma(yi, vi, data=dat, method="DL")
> summary(fit1)
>
> estimate se zval pval ci.lb ci.ub
> -0.7141 0.1787 -3.9952 <.0001 -1.0644 -0.3638 ***
>
> ## Robust SE
> robustSE(fit1)
> estimate se.R t.R p.R ci.lb ci.ub
> intrcpt -0.7141172 0.1791445 -3.986265 0.001805797 -0.725588 -0.7026465
>
> ## mixed-effects model with two moderators (absolute latitude and
> publication year)
> fit2 <- rma(yi, vi, mods=cbind(ablat, year), data=dat, method="DL")
> summary(fit2)
>
> estimate se zval pval ci.lb ci.ub
> intrcpt -1.2798 25.7550 -0.0497 0.9604 -51.7586 49.1990
> ablat -0.0288 0.0090 -3.2035 0.0014 -0.0464 -0.0112 **
> year 0.0008 0.0130 0.0594 0.9526 -0.0247 0.0262
>
> ### Robust SE
> robustSE(fit2)
> estimate se.R t.R p.R ci.lb
> intrcpt -1.2797914381 22.860353022 -0.05598301 0.956458098 -2.749670e+00
> ablat -0.0287644840 0.007212163 -3.98832970 0.002566210 -2.922821e-02
> year 0.0007720798 0.011550188 0.06684565 0.948022174 2.942412e-05
> ci.ub
> intrcpt 0.190086881
> ablat -0.028300755
> year 0.001514735
>
More information about the R-help
mailing list