cpanderson
2012-Aug-02 22:43 UTC
[R] metafor- interpretation of moderators test for raw proportions
Hello metafor users, I'm using metafor to perform a single-effect summary estimate of the raw proportion of patients experiencing a post-operative complication, and I'm interested in seeing if this proportion differs between the three most commonly used surgical techniques. The software is working as expected, but I would like to double check on the interpretation of my mixed-effect model summary. The levels of my moderator variable surgery.type are LateralWedge, Dome, and Complex, which I've dummy coded like this: dat$LateralWedge<-ifelse(dat$Technique=="LateralWedge",1,0) datDome<-ifelse(dat$Technique=="Dome",1,0) dat$Complex<-ifelse(dat$Technique=="Complex",1,0) When i fit this random effect model: rma(yi,vi,data=dat, mods=cbind(LateralWedge,Dome,Complex)) I get this output: /Mixed-Effects Model (k = 33; tau^2 estimator: REML) logLik Deviance AIC BIC -11.2426 22.4852 30.4852 36.0900 tau^2 (estimate of residual amount of heterogeneity): 0.0925 (SE = 0.0269) tau (sqrt of the estimate of residual heterogeneity): 0.3041 Test for Residual Heterogeneity: QE(df = 30) = 654.3038, p-val < .0001 *Test of Moderators (coefficient(s) 1,2,3): QM(df = 3) = 128.7528, p-val < .0001* Model Results: estimate se zval pval ci.lb ci.ub LateralWedge *0.6462* 0.0722 8.9450 <.0001 *0.5046 0.7878* *** Dome *0.6659 * 0.1471 4.5283 <.0001 *0.3777 0.9541 * *** Complex *0.6938* 0.1306 5.3136 <.0001 * 0.4379 0.9498* ***/ My suspicion is that the test of moderators is addressing the hypothesis that the proportion of my outcome is not zero in any of the three groups, because all the estimated proportions appear to be well within the confidence bounds of the other two groups. Is this correct? If it is, I would appreciate any suggestions on how I would change that hypothesis test, so that I am testing whether the proportion in group 1 = the proportion in group 2 = the proportion in group 3? Thanks for your input, Christopher Anderson -- View this message in context: http://r.789695.n4.nabble.com/metafor-interpretation-of-moderators-test-for-raw-proportions-tp4638972.html Sent from the R help mailing list archive at Nabble.com.
Viechtbauer Wolfgang (STAT)
2012-Aug-03 08:17 UTC
[R] metafor- interpretation of moderators test for raw proportions
Are you sure that output was produced by: rma(yi, vi, data=dat, mods=cbind(LateralWedge,Dome,Complex))? Because your model does not have an intercept, which suggests that you actually used: rma(yi, vi, data=dat, mods=cbind(LateralWedge,Dome,Complex), intercept=FALSE) In that case, the coefficients are the estimated (average) proportions for each level of that moderator. If you want to actually see whether levels differ from each other, you make one of the levels your reference level, leave out the corresponding dummy, and include the model intercept: rma(yi, vi, data=dat, mods=cbind(Dome,Complex)) Now the coefficients estimate the (average) difference between the proportions for the Dome vs. LateralWedge levels and the Complex vs. LateralWedge levels. If you want to estimate the difference between Dome and Complex, you can use: rma(yi, vi, data=dat, mods=cbind(LateralWedge,Complex)) Those are all pairwise comparisons. The QM test will tell you whether that factor is significant overall. Note that you do not need to dummy code those levels manually. You could just use: rma(yi, vi, data=dat, mods = ~ Technique) With the relevel() function, you can change which level is used as the reference level. For example: rma(yi, vi, data=dat, mods = ~ relevel(Technique, ref="LateralWedge")) Best, Wolfgang -- Wolfgang Viechtbauer, Ph.D., Statistician Department of Psychiatry and Psychology School for Mental Health and Neuroscience Faculty of Health, Medicine, and Life Sciences Maastricht University, P.O. Box 616 (VIJV1) 6200 MD Maastricht, The Netherlands +31 (43) 388-4170 | http://www.wvbauer.com> -----Original Message----- > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] > On Behalf Of cpanderson > Sent: Friday, August 03, 2012 00:43 > To: r-help at r-project.org > Subject: [R] metafor- interpretation of moderators test for raw > proportions > > Hello metafor users, > I'm using metafor to perform a single-effect summary estimate of the raw > proportion of patients experiencing a post-operative complication, and I'm > interested in seeing if this proportion differs between the three most > commonly used surgical techniques. The software is working as expected, > but > I would like to double check on the interpretation of my mixed-effect > model > summary. > > The levels of my moderator variable surgery.type are LateralWedge, Dome, > and > Complex, which I've dummy coded like this: > > dat$LateralWedge<-ifelse(dat$Technique=="LateralWedge",1,0) > datDome<-ifelse(dat$Technique=="Dome",1,0) > dat$Complex<-ifelse(dat$Technique=="Complex",1,0) > > When i fit this random effect model: > rma(yi,vi,data=dat, mods=cbind(LateralWedge,Dome,Complex)) > > I get this output: > /Mixed-Effects Model (k = 33; tau^2 estimator: REML) > > logLik Deviance AIC BIC > -11.2426 22.4852 30.4852 36.0900 > > tau^2 (estimate of residual amount of heterogeneity): 0.0925 (SE = 0.0269) > tau (sqrt of the estimate of residual heterogeneity): 0.3041 > > Test for Residual Heterogeneity: > QE(df = 30) = 654.3038, p-val < .0001 > > *Test of Moderators (coefficient(s) 1,2,3): > QM(df = 3) = 128.7528, p-val < .0001* > Model Results: > > estimate se zval pval ci.lb ci.ub > LateralWedge *0.6462* 0.0722 8.9450 <.0001 *0.5046 0.7878* *** > Dome *0.6659 * 0.1471 4.5283 <.0001 *0.3777 0.9541 * *** > Complex *0.6938* 0.1306 5.3136 <.0001 * 0.4379 0.9498* ***/ > > My suspicion is that the test of moderators is addressing the hypothesis > that the proportion of my outcome is not zero in any of the three groups, > because all the estimated proportions appear to be well within the > confidence bounds of the other two groups. Is this correct? If it is, I > would appreciate any suggestions on how I would change that hypothesis > test, > so that I am testing whether the proportion in group 1 = the proportion in > group 2 = the proportion in group 3? > > Thanks for your input, > Christopher Anderson
Viechtbauer Wolfgang (STAT)
2012-Aug-04 16:31 UTC
[R] metafor- interpretation of moderators test for raw proportions
Just to follow up on what Michael wrote: I cannot reproduce that error. For example, this all works as intended: data(dat.bcg) dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg, append=TRUE) rma(yi, vi, mods = ~ alloc, data=dat) ### 'alloc' automatically converted to a factor dat$alloc <- factor(dat$alloc) ### explicitly make 'alloc' a factor rma(yi, vi, mods = ~ alloc, data=dat) ### works as before rma(yi, vi, mods = ~ factor(alloc), data=dat) ### factor() not necessary, but works If you can provide a reproducible example, I'll be glad to look into the issue. Aside from that, the error you got occurs when the design matrix is not of full rank. For example: X <- model.matrix(~ factor(alloc) - 1, data=dat) rma(yi, vi, mods = X, data=dat) will fail, because the model now has an intercept plus the 3 dummy variables for the 3 levels (setting intercept=FALSE will make this work). It seems to me that this is what happened (since in your previous post, you showed that you coded the three levels of your factor manually). But this is something different than what you describe below, so I don't know for sure. Best, Wolfgang ________________________________________ From: r-help-bounces at r-project.org [r-help-bounces at r-project.org] On Behalf Of cpanderson [christopher.p.anderson at healthpartners.com] Sent: Friday, August 03, 2012 6:03 PM To: r-help at r-project.org Subject: Re: [R] metafor- interpretation of moderators test for raw proportions Wolfgang, Thanks for your quick response. You are correct- indeed I had inadvertently left out that i had set intercept = FALSE in rma. Following your suggestions, I get the following results: Test of Moderators (coefficient(s) 2,3): QM(df = 2) = 0.2207, p-val = 0.8955 Model Results: estimate se zval pval ci.lb ci.ub intrcpt 0.6498 0.0492 13.2160 <.0001 0.5534 0.7462 *** Complex 0.0457 0.1007 0.4538 0.6500 -0.1517 0.2430 Dome 0.0244 0.1135 0.2148 0.8299 -0.1980 0.2468 This may be an even dumber question than my first one, but if you have time I'd appreciate knowing how this works. Originally I had tried the syntax mods = ~ Technique. I ended up getting the following error message: Error in qr.solve(wX, diag(k)) : singular matrix 'a' in solve if I supply the argument that mods =~ factor(Technique), I don't get the singular matrix 'a' message. But I don't understand why that should make a difference, because dat$Technique is already an object of class factor: class(dat$Technique) [1] "factor" Thanks again. Christopher