Post model fitting checks in Metafor (rma.mv)
Dear Wolfgang, Sorry to contact you again after so long, I had completed one part of my meta-analysis and moved on to writing another chapter of my thesis. Now that I begin the second part of my meta-analysis I have come across a problem. For the response variable I can either use the absolute change in growth (whether studies increased or decreased growth) or I can use the growth change relative to oxidative stress (this was positive if both growth and oxidative stress were increased/decreased; but negative if they changed in opposite directions). Ideally I'd like to use the latter but when I do so, the profile plot for rho is flat (it looks fine for the abolute values). The only difference between the two are the signs and when plotting the data as histograms, both show a similar spread. Oddly, if I split the data into 3 groups based on the marker of oxidative stress, the plots look better too. Any ideas as to why this might be happening are extremely welcome. Kind regards, Shona Shona Smith PhD Researcher Room 321 Institute of Biodiversity, Animal Health and Comparative Medicine Graham Kerr Building University of Glasgow Glasgow G12 8QQ
From: Viechtbauer Wolfgang (STAT) [wolfgang.viechtbauer at maastrichtuniversity.nl]
Sent: 03 September 2014 18:40
To: Shona Smith; r-sig-mixed-models at r-project.org
Subject: RE: Post model fitting checks in Metafor (rma.mv)
Sent: 03 September 2014 18:40
To: Shona Smith; r-sig-mixed-models at r-project.org
Subject: RE: Post model fitting checks in Metafor (rma.mv)
By default (check options("contrasts")), R should use treatment contrasts (see help(contr.treatment)). So, one level of each factor is chosen as the reference level and you get coefficients that indicate the contrasts with this reference level. So, if you use 'mods = ~ Age + Treatment + Biomarker' and Age has 3 levels, Treatment has 5 levels, and Biomarker has 3 levels, then you should get 2+4+2=8 coefficients. So:
predict(res, newmods=c(1,0, 0,0,0,0, 0,0))
would give you the predicted effect for the second level of Age, reference level for Treatment, and reference level for Biomarker. And
predict(res, newmods=c(0,1, 0,0,1,0, 1,0))
would be for the third level of Age, fourth level of Treatment, and second level of Biomarker. And so on ...
You may also want to take a look at this tutorial:
http://www.metafor-project.org/doku.php/tips:testing_factors_lincoms
It doesn't cover multiple factors (I'll add one to the website soon), but should help to clarify things a bit. Also, some things won't work for 'rma.mv' models at this point (e.g., the anova() function or permutation tests). But you can use predict(), linearHypothesis() from the 'car' package, and glht() from the 'multcomp' package.
Best,
Wolfgang
> -----Original Message-----
> From: Shona Smith [mailto:s.smith.7 at research.gla.ac.uk]
> Sent: Wednesday, September 03, 2014 19:24
> To: Viechtbauer Wolfgang (STAT); r-sig-mixed-models at r-project.org
> Subject: RE: Post model fitting checks in Metafor (rma.mv)
>
> Hi Wolfgang,
>
> This definitely makes sense thanks - to look at an 'overall effect'
> without the moderators sort of defeats the purpose of putting them there
> in the first place. I was able to download your paper so I can have a
> wee look at that too, thanks. I think I shall stick with extracting
> predicted values for specific levels of moderators. I tried the predict
> function and realise I need to code the factor levels myself, but am a
> little confused about how to do so. To see how the factors were coded in
> the model, I used:
>
> predict(res, addx=TRUE)
>
> Then I know I can use 'newmods' to put in the codes for the factor levels
> I would like a predicted value for. There were a lot of '1's and '0's
> (one for each observation and level of factor) so how can I use this to
> code my factor levels?
>
> As for the profile plot for sigma2, I suppose it peaks at ~0.05 but it is
> just a straight line, declining from an REML value of ~ -71 (at Sigma2 =
> ~0.05) to an REML value of ~ -72.5 (at Sigma2 = 1).
>
> Of course missing values were the problem with the residual against
> fitted plot, it was silly of me to miss this. Thanks!
>
> Kind regards,
> Shona
>
>
> Shona Smith
> PhD Researcher
>
> Room 321
> Institute of Biodiversity, Animal Health and Comparative Medicine
> Graham Kerr Building
> University of Glasgow
> Glasgow G12 8QQ
> ________________________________________
> From: Viechtbauer Wolfgang (STAT)
> [wolfgang.viechtbauer at maastrichtuniversity.nl]
> Sent: 03 September 2014 14:32
> To: Shona Smith; r-sig-mixed-models at r-project.org
> Subject: RE: Post model fitting checks in Metafor (rma.mv)
>
> Regarding the profile plot for sigma2: I am not quite sure I understand.
> Does it 'peak' at zero? So, is the estimate (essentially) zero then? That
> would be okay (essentially means that the variability due to 'Study' is
> no larger than what would be expected due to the other variance
> components and/or sampling variability). The issue of zero variance
> components was also recently (and also in past) discussed on this list
> (not with respect to meta-analysis, but it's the same issue).
>
> Regarding the resid-fitted plot: So, looks like you have some missings,
> so the two vectors end up being of different length. This should do it:
>
> options(na.action = "na.pass")
> plot(fitted(res), rstandard(res)$z, pch=19)
>
> Also explained here: http://www.metafor-
> project.org/doku.php/tips:handling_missing_data
>
> Regarding overall effects: I personally don't think an 'overall' effect
> makes much sense when the effect size appears to be related to a number
> of moderators/covariates. Take the simplest situation, where the true
> effect is of size theta_1 for level 1 of a dichotomous covariate and
> theta_2 for level 2. Now suppose we ignore that covariate and fit a
> random-effects model. Essentially, that is a misfitted model, because the
> model assumes normally distributed true effects (and not two point
> massess). Also, where the 'average' effect then falls depends on how many
> studies in the sample are at level 1 and at level 2 of that covariate.
> That doesn't seem that sensible to me. So, instead, we can fit the model
> with the covariate and then compute predicted values, for example, for
> level 1 or level 2. Or, if you really want an 'overall' effect, we could
> use a sort of 'lsmeans' approach and say: Let's assume that in the
> population of studies, 50% are at level 1 and 50% at level 2, so let's
> compute the predicted effect for such a population (essentially fill in
> 0.5 for the 'dummy' variable when computing the predicted value). But I
> would then describe explicitly that this is what was done (as it makes
> clear what the assumption is about the population of studies).
>
> I discuss this issue a bit in this article:
>
> Viechtbauer, W. (2007). Accounting for heterogeneity via random-effects
> models and moderator analyses in meta-analysis. Zeitschrift f?r
> Psychologie / Journal of Psychology, 215(2), 104-121.
>
> (not the 'lsmeans' idea, but the problem of fitting random-effects models
> when there is a relevant moderators/covariate). If you can't access the
> article, let me know and I'll send you a copy if you are interested.
>
> Best,
> Wolfgang
>
> > -----Original Message-----
> > From: Shona Smith [mailto:s.smith.7 at research.gla.ac.uk]
> > Sent: Wednesday, September 03, 2014 13:39
> > To: Viechtbauer Wolfgang (STAT); r-sig-mixed-models at r-project.org
> > Subject: RE: Post model fitting checks in Metafor (rma.mv)
> >
> > Hi Wolfgang,
> >
> > Thanks for such a useful and quick reply. The sigma2 profile plot does
> > not have a peak - but it's not flat, it seems to just decline (even
> when
> > I expand the x axis values) - I wondered what this meant? The other
> two
> > look fine, although tau2 drops off much more gradually than rho (which
> is
> > quite a steep drop) after the parameter estimate. What would be the
> > ideal plot?
> >
> > I can now plot the standardised residuals but fitted against residuals
> > gives me an error:
> >
> > plot(fitted(res), rstandard(res)$z, pch=19)
> > Error in xy.coords(x, y, xlabel, ylabel, log) :
> > 'x' and 'y' lengths differ
> >
> > The fitted values seem to have the row number alongside each value,
> > whiles the standardised residuals don't - so perhaps this is the
> problem?
> >
> > Finally, I also wondered if it is possible to get an overall effect
> size
> > estimate for my response variable? I notice other studies have carried
> > out a separate meta-analysis to obtain this, before looking at
> > moderators, but I wasn't sure if this was correct.
> >
> > Thanks again,
> > Shona
> >
> >
> > Shona Smith
> > PhD Researcher
> >
> > Room 321
> > Institute of Biodiversity, Animal Health and Comparative Medicine
> > Graham Kerr Building
> > University of Glasgow
> > Glasgow G12 8QQ
> > ________________________________________
> > From: Viechtbauer Wolfgang (STAT)
> > [wolfgang.viechtbauer at maastrichtuniversity.nl]
> > Sent: 03 September 2014 10:28
> > To: Shona Smith; r-sig-mixed-models at r-project.org
> > Subject: RE: Post model fitting checks in Metafor (rma.mv)
> >
> > Dear Shona,
> >
> > The profile() function in metafor allows you to examine the profiled
> > (restricted) log-likelihood for a particular parameter. So, ideally,
> you
> > should do this for each variance component and correlation in the model
> > (in your case, sigma2, tau2 and rho). I am not sure what you mean by: "
> > how I know which value to specify for each?" After you have fitted your
> > model and stored the results in, let's say, 'res', then just do:
> >
> > par(mfrow=c(3,1))
> > profile(res, sigma2=1)
> > profile(res, tau2=1)
> > profile(res, rho=1)
> >
> > to get all three profile plots. And yes, the functions should peak at
> the
> > parameter estimates. If a function is flat, then this suggests that the
> > model is overparameterized.
> >
> > You could also look at how quickly the log-likelihood drops off as you
> > move away from the parameter estimate. The bounds of a 95% profile
> > likelihood CI for a particular parameter would be those two values from
> > the x-axis where the log-likelihood has dropped by 3.84/2. You could
> add
> >
> > abline(h = logLik(res) - qchisq(.95, df=1)/2, lty="dotted")
> >
> > to the figures to see that cutoff. Depending on how much data you have,
> > you may find that those CIs are quite wide. You may have to increase
> the
> > x-axis range, in case the cutoff isn't reached within the range chosen
> by
> > default by the profile function (use the 'xlim' argument).
> >
> > Indeed, you could also look at the (standardized) residuals. Use
> >
> > rstandard(res)$z
> >
> > to get those values. Or:
> >
> > plot(fitted(res), rstandard(res)$z, pch=19)
> >
> > to create a fitted values versus standardized residuals plot.
> >
> > As for the interpretation of the results when you exclude the intercept
> > (i.e., mods = ~ Age + Treatment + Biomarker - 1), you will get the
> > estimated (average) effect for each level of 'Age', but the
> coefficients
> > for 'Treatment' and 'Biomarker' are still going to be contrasts that
> > indicate how much higher/lower the (average) effect is for the levels
> > indicated, relative to the reference level.
> >
> > I hope this helps!
> >
> > Best,
> > Wolfgang
> >
> > > -----Original Message-----
> > > From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-
> > > models-bounces at r-project.org] On Behalf Of Shona Smith
> > > Sent: Tuesday, September 02, 2014 18:37
> > > To: r-sig-mixed-models at r-project.org
> > > Subject: [R-sig-ME] Post model fitting checks in Metafor (rma.mv)
> > >
> > > Hi all,
> > >
> > > I am currently conducting a meta-analysis using rma.mv in metafor.
> My
> > > model uses Hedges' d (converted to g) and includes 3 moderators (age-
> 3
> > > levels; treatment-5 levels; biomarker-3 levels). I have included 3
> > > random effects: species nested within taxonomic class (since I have
> > more
> > > than one study for some species, and species are spread over 7
> > taxonomic
> > > classes) and study separately. So my code is as follows:
> > >
> > > rma.mv(yi, vi, mods = ~ Age + Treatment + Biomarker, random = list(~
> 1
> > |
> > > Study, ~ Species | Taxonomic.class), data=mydata)
> > >
> > > I was wondering what the best method for post model fitting checks
> was
> > in
> > > rma.mv? I know in the reference manual it mentions profile.rma to
> > create
> > > a plot of the restricted log likelihood and I have done so. However,
> I
> > > wondered if I need to plot all 3 variables (sigma2, tau2 and rho) and
> > > also how I know which value to specify for each? Am I correct in
> that
> > I
> > > should see a clear peak in each graph? Is there anything else I
> should
> > > be looking for?
> > >
> > > For post model fitting checks should I also look at residual
> normality
> > > and residual against fitted values, as would be done for a typical
> > mixed
> > > model? I think the standardised residuals are best for this - I can
> > get
> > > them with rstandard.rma.mv, but it does not allow me to plot them.
> > >
> > > Finally, when I include the intercept in the model, I can see if
> there
> > > are significant differences among moderator levels. However, I was
> > > wondering what the output includes when the intercept is not
> included:
> > is
> > > this the overall effect size estimates for each moderator level?
> > >
> > > Kind regards,
> > > Shona
> > >
> > >
> > > Shona Smith
> > > PhD Student
> > >
> > > Room 321
> > > Institute of Biodiversity, Animal Health and Comparative Medicine
> > > Graham Kerr Building
> > > University of Glasgow
> > > Glasgow G12 8QQ
> > > _______________________________________________
> > > R-sig-mixed-models at r-project.org mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models