R help - Aug 2012 - 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.

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

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