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 <mikewlche...@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 > ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.