[R-sig-ME] standard error and statistical significance in lmer versus lm

[Vincenzo Lagani]

Dear all,

I am modelling gene expression data with mixed models using lme4. My 
goal is to assess whether gene expression globally decreases or 
increases with time.

Specifically, the data consist in whole expression profiles measured at 
five different time points in two different strains of mice. For each 
time point and each strain there are three replicates. The data are not 
longitudinal, i.e., different mice are used at each time point. We had 
to remove one profile from a time point because it was not matching our 
quality criteria, so the data became not properly balanced.

On these data I am fitting the following model:

lme4Model.full <- lmer(values ~ week * genotype + (week | probesets), 
data = dataset, REML = FALSE)

where 'values' stands for the gene expression, 'week' is a scaled 
numeric reporting the age of the mice in weeks, and 'genotype' is a 
factor with two levels representing the two different genetic 
backgrounds. Each single gene is let having its own random intercept and 
slope.

What puzzles me is that when I compare this model with its simpler, 
non-mixed version:

lmModel.full <- lm(values ~ week * genotype, data = dataset)

I obtain the same coefficients but different standard errors (see 
below). Furthermore, while the interaction coefficient is not 
significant in the simple linear model, it becomes highly significant in 
the mixed model, at least according to these ANOVA tests:

lmModel <- lm(values ~ week + genotype, data = dataset)
lme4Model <- lmer(values ~ week + genotype + (week | probesets), data = 
dataset, REML = FALSE)

>anova(lmModel.full, lmModel)
Analysis of Variance Table

Model 1: values ~ week * genotype
Model 2: values ~ week + genotype
   Res.Df     RSS Df Sum of Sq      F Pr(>F)
1 414341 2162664
2 414342 2162666 -1   -2.4697 0.4732 0.4915

>anova(lme4Model.full, lme4Model)
Data: dataset
Models:
lme4Model: values ~ week + genotype + (week | probesets)
lme4Model.full: values ~ week * genotype + (week | probesets)
                Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
lme4Model       7 72376 72452 -36181    72362
lme4Model.full  8 72325 72413 -36155    72309 52.472      1  4.364e-13 ***



Assessing whether the interaction coefficient is significant is actually 
the aim of my study, and having two totally different answers confuses 
me. My understanding is that the mixed model better catches the variance 
structure of the data and thus it is able to better estimate the 
standard errors and p-values of the coefficients. Is this correct? In 
other words, can I confidently claim that the p-values obtained from the 
mixed models are "the correct ones" and that the interaction term is 
actually significant?

My apologies if this issue has been already posted on this list. Despite 
having seen multiple posts here on similar topics, I have not been able 
to find an answer to these questions.

Thanks in advance for your help. Any suggestion is very welcome.

Regards,

Vincenzo


>summary(lme4Model.full)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: values ~ week * genotype + (week | probesets)
    Data: dataset

      AIC      BIC   logLik deviance df.resid
  72325.3  72412.8 -36154.7  72309.3   414337

Scaled residuals:
      Min       1Q   Median       3Q      Max
-13.2622  -0.5246   0.0128   0.5270  23.3456

Random effects:
  Groups    Name        Variance Std.Dev. Corr
  probesets (Intercept) 5.167338 2.27318
            week        0.005029 0.07092  -0.23
  Residual              0.047063 0.21694
Number of obs: 414345, groups:  probesets, 18015

Fixed effects:
                    Estimate Std. Error t value
(Intercept)       6.5727983  0.0169414   388.0
week             -0.0151336  0.0006823   -22.2
genotypeCSB       0.2405300  0.0007121   337.8
week:genotypeCSB  0.0050430  0.0006962     7.2

Correlation of Fixed Effects:
             (Intr) week   gntCSB
week        -0.177
genotypeCSB -0.015 -0.028
wk:gntypCSB -0.001 -0.392 -0.054


>summary(lmModel.full)

Call:
lm(formula = values ~ week * genotype, data = dataset)

Residuals:
     Min      1Q  Median      3Q     Max
-5.7560 -1.9881  0.0083  1.6823  7.6407

Coefficients:
                   Estimate Std. Error  t value Pr(>|t|)
(Intercept)       6.572798   0.004407 1491.384  < 2e-16 ***
week             -0.015134   0.004547   -3.329 0.000873 ***
genotypeCSB       0.240530   0.007500   32.072  < 2e-16 ***
week:genotypeCSB  0.005043   0.007331    0.688 0.491537
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.285 on 414341 degrees of freedom
Multiple R-squared:  0.002491,	Adjusted R-squared:  0.002484
F-statistic: 344.9 on 3 and 414341 DF,  p-value: < 2.2e-16



	[[alternative HTML version deleted]]