How to get an effect size in nlme?

Question

I have recently been asked to provide an effect size for a mixed-effects linear model which I had fit using nlme. After doing a lot of hunting, it's become clear that effect size with mixed models is not a straightforward or well-defined measure. This thorough explanation is specific to the lme4 package. I'm wondering if there is a comparable solution for finding R^2 for models fit using nlme.

My model looks something like this:

model <-lme(depvar ~ (var1 + var2 + var3)^2, data, random = list( sub = ~ 1 + var3))

Further, is there a way to get effect sizes for main effects and the interactions?

Answer 1

You can report the likelihood ratio test as an effect size measure. I'm not sure what the exact design of your overall model is, but say you're interested in a two-way repeated measures design where you want to assess the main effects of var1, var2, and the var1*var2 interaction.

To get the likelihood ratio, you can take a multilevel approach using nlme whereby you compare models against each other in a nested design.

Here's how you set the null model (where the dependent variable is predicted by its overall mean).

nullModel <- lme(depvar ~ 1, random = ~1 | id/var1/var2, data = data, method = "ML") 

The var1 model

var1Model <- lme(depvar ~ var1, random = ~1 | id/var1/var2, data = data,
method = "ML")

The var2 model

var2Model <- lme(depvar ~ var1 + var2, random = ~1 | id/var1/var2, data = data, method = "ML")

And finally, the interaction or 'full' model, which includes the main effects and the interaction

IModel <- lme(depvar ~ var1 * var2, random = ~1 | id/var1/var2, data = data, method = "ML")

Then, you compare the models using the anova() function

anova(nullModel, var1Model, var2Model, IModel)

The L.ratios in your output are ratios of how much more likely the data is under a given model compared to another. You get p-values.

Here's example anova output from a dataset I'm currently working on (note the different model names). Here we can see that that the data are only 1.12 times more likely under model 2 (cond_model_dp) than the model 1 (baseline_dp), the null model. The interaction model (fullModel_dp) isn't much better than the model with both main effects (em_model_dp). According the p-values below, none of the models are significantly better fits of the data.

              Model df      AIC      BIC    logLik   Test  L.Ratio p-value
baseline_dp       1  5 1073.524 1086.549 -531.7618                        
cond_model_dp     2  7 1076.408 1094.644 -531.2038 1 vs 2 1.115851  0.5724
em_model_dp       3  6 1072.636 1088.267 -530.3180 2 vs 3 1.771603  0.1832
fullModel_dp      4 10 1076.956 1103.008 -528.4780 3 vs 4 3.680001  0.4510

While many people (particularly in the biobehavioral sciences) would expect an eta-squared statistic for such an analysis, the likelihood ratio is arguably a "better" statistic as it's easier to interpret the magnitude of the effect. That is, Cohen's guidelines of small, medium, and, and large effects for eta squared are just rules of thumb, whereas the likelihood ratio is more immediately intuitive as it operates like an odds ratio.

Finally, if you have any missing data, just add "na.action = na.exclude" to your model, like this

var1Model <- lme(depvar ~ var1, random = ~1 | id/var1/var2, data = data,
method = "ML", na.action = na.exclude)