Post model fitting checks in Metafor (rma.mv)

________________________________________
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