[R-sig-ME] Post model fitting checks in Metafor (rma.mv)

[Viechtbauer Wolfgang (STAT)]

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