In jaromilfrossard/gANOVA: Generalization of the ANOVA covariance structure to mixed models

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction

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)

Simulating data

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)

Comparing estimation

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)

jaromilfrossard/gANOVA documentation built on July 28, 2020, 5:32 p.m.