knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
We present here some examples to highlight the differences between the gANOVA()
function and the built-in spherical distribution of lmer()
. Both functions estimate a spherical correlation structure. However the parametrization using gANOVA()
increases the parameter space and widen the range of the models that can be estimated without having estimates at the boundary of the parameter space (e.g. variances estimated at 0).
We need the following packages:
library(lme4) library(gANOVA) library(tidyverse)
Then, we construct the design of a simple experiment. It contains 12 participants, a within-participant factor with 4 levels and 3 replications for each participant within each level. The design is simply:
set.seed(42) df <- expand.grid(pt = paste0("P",str_pad(1:12,width = 2,pad = "0")), A = paste0("A",1:4), r = 1:3, stringsAsFactors = F)%>% select(-r)
Then we simulate 2 different response variables. There are y_diff
and y_same
and they differ only by the variance of the random intercepts. The response y_diff
has a lower variance of the random intercepts compared to the one of the random interactions.
# add +(1|pt) random intercept df<- df%>% nest(data=-pt)%>% mutate(ranef_pt = rnorm(n(), sd = .8))%>% unnest(data) # add +(1|pt:A) random interaction df<- df%>% nest(data=-c(pt,A))%>% mutate(ranef_ptA = rnorm(n(), sd = 3))%>% unnest(data) # add an error term df<- df%>% mutate(err = rnorm(n(),sd = 4)) # create 2 dependant variables df <- df%>% transmute(pt = pt, A= A, y_diff = ranef_pt + ranef_ptA + err, y_same = ranef_pt*2 +ranef_ptA + err)
Finally, we estimate the model using y_same
as dependant variable:
lmer_same <- lmer(y_same ~ A+ (1|pt) + (1|pt:A), data=df, contrasts = list(A = contr.sum)) summary(lmer_same)
and then gANOVA:
gANOVA_same <- gANOVA(y_same ~ A + (1|pt|A), data=df, contrasts = list(A = contr.sum)) summary(gANOVA_same)
Using y_same
the 2 estimations are the same (likelihood, fixed parameters, etc). Here the reparametrizations of gANOVA()
does not improve the results as the variance of the random intercept is big enough with respect to the variance of the random interaction.
However, using y_diff
, the variability of the random intercepts is smaller. The estimation using the parametrization of lmer()
produces a singular fit:
lmer_diff <- lmer(y_diff~ A + (1|pt)+(1|pt:A),data=df,contrasts = list(A = contr.sum)) summary(lmer_diff)
Estimating the model using gANOVA()
reduces the constraint on the parameter space (of the variance of the intercept). The gANOVA
parametrization produces a non-zero estimate of the variance of the random intercept. It follows that the likelihood is smaller:
mod_gANOVA <- gANOVA(y_diff~ A+ (1|pt|A),data=df,contrasts = list(A = contr.sum)) summary(mod_gANOVA)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.