3
$\begingroup$

Let's say I have an experiment where I repeatedly measure something over Time (1:10), say 10 days.

Here I simulate some data...

set.seed(101)

N = 10 # number of repeats
n = 10 # number of individuals

# Data
Time = rep(1:N,n) # measured over time

Note that I set the fixed effect of time to 0.

# Fixed effects
X = matrix(c(Time), ncol = 1)
B = c(0) # Set fixed effect of Time's coefficient to 0 

# Random effects
id = factor(c(0L, cumsum(diff(Time) < 0))) # make id for each block of Time
Z = model.matrix(~id+id:Time) # create design matrix

u = c(rnorm(n),rnorm(n,2,1)) # random slopes and intercepts
e = rnorm(N*n) # error


y = X%*%B + Z%*%u + e # Model  

Now let's say I fit two models:

library(lme4)
fit1 = lmer(y ~ 1 + Time + (1 + Time|id))
fit2 = lmer(y ~ 1 + (1+ Time|id))

Comparing the two with anova()

anova(fit1,fit2)
Models:
fit2: y ~ 1 + (1 + Time | id)
fit1: y ~ 1 + Time + (1 + Time | id)
 Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
fit2  5 380.80 393.83 -185.40   370.80                             
fit1  6 329.41 345.05 -158.71   317.41 53.389      1  2.736e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here we see that fit1, which includes Time as a fixed effect, performs much better, but why is this true when its coefficient is set to 0?

I have a feeling it has something to do with this line.

Z = model.matrix(~id+id:Time) # create design matrix

but my intuition is failing me.

$\endgroup$
3
$\begingroup$

If you look at the summary of your models, you'll see that the estimated time effect is actually greater than zero (about 2 given your seed). I think your error is in the

u = c(rnorm(n),rnorm(n,2,1)) # random slopes and intercepts

It should actually be flipped based on your design matrix, so it reads random intercepts and slopes.

u = c(rnorm(n,2,1), rnorm(n)) # random intercepts and slopes

Before you had an intercept of zero and a slope of 2, when you actually wanted it the other way around.

This yields what you're looking for I think. Now if you use summary on mod1, you'll see that the fixed effect of time is close to zero (it would be closer with more than N = 10), and it shows that mod2 provides more parsimonious model fit according to AIC.

Data: 
Models:
fit2: y ~ 1 + (1 + Time | id)
fit1: y ~ 1 + Time + (1 + Time | id)
     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
fit2  5 352.85 365.88 -171.43   342.85                         
fit1  6 354.69 370.32 -171.35   342.69 0.1584      1     0.6907
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.