Degrees of Freedom Std. Error estimates in mixed model of a parallel group design with repeated measurements are not as I expected

Question

I have multilevel data and I would like to use mixed effects linear model. For each subject (sampleName variable), I have repeated measurements (peptide_Id variable). I first did model this data in two steps: step 1 - get an estimate for each subject, step 2 - fit a fixed effects model to estimate treatment group (dilution. variable) coefficients.

library(lmerTest)
mb <- readRDS(url("https://github.com/wolski/pairseqsim2/files/4072167/mb.zip"))

lm_m <- lm(transformedIntensity ~ sampleName, mb)

protAnn <- mb %>% select(-transformedIntensity, -peptide_Id) %>% distinct()

est <- coef(lm_m) 
est <- c(est[1] , est[1] - est[-1])
names(est)[1] <- "sampleNamea~10"
names(est) <- gsub("sampleName","", names(est))
protInt <- enframe(est, name = "sampleName", "IntEstimate")
protLevel <- inner_join(protAnn, protInt)

protmodel <- lm(IntEstimate ~ dilution.,data = protLevel)
df.residual(protmodel)
coef(summary(protmodel ))

> df.residual(protmodel)
[1] 15
> coef(summary(protmodel))
               Estimate Std. Error    t value    Pr(>|t|)
(Intercept) -0.21742659 0.05958290 -3.6491445 0.002374004
dilution.b   0.06030848 0.08426294  0.7157178 0.485158848
dilution.c   0.06765230 0.08426294  0.8028714 0.434586826
dilution.d   0.14464911 0.08426294  1.7166397 0.106616968
dilution.e   0.19043001 0.08426294  2.2599497 0.039132403

Now, I thought to use mixed effect models, hoping that they might have some advantages over the two-stage modelling. I am modelling the treatment group (dilution.) as a fixed effect and the repeated measure as a random effect.

lm.mixed <- lme4::lmer(transformedIntensity ~ dilution. + (1 | peptide_Id), data = mb)
coef(summary(lm.mixed))

What confuses me about the result, and I am struggling to understand, is why the Std. Error estimates for the fixed effects coefficient (except of the intercept) are so small compared with the two-stage modelling above.

> coef(summary(lm.mixed))
               Estimate Std. Error    t value
(Intercept) -0.26464840 0.29013559 -0.9121542
dilution.b  -0.07686173 0.07268624 -1.0574454
dilution.c  -0.06208256 0.07195788 -0.8627625
dilution.d  -0.13610607 0.07134742 -1.9076523
dilution.e  -0.18188697 0.07134742 -2.5493138

They are actually the same (except the intercept) as if fitting a fixed effect model i.e.

lm.fixed.pep <- lm(transformedIntensity ~ dilution. + peptide_Id, data = mb)
coef(summary(lm.fixed.pep))

Furthermore, if I would like to obtain p-values using the package lmerTest. However, the degrees of freedom are also in my opinion way too high for all coefficient except for the intercept. On the other hand for the intercept coefficient the estimates of the standard error are (I think) to high and degrees of freedom are too low.

library(lmerTest)
lm.mixed <- lmerTest::lmer(transformedIntensity ~ dilution. + (1|peptide_Id), data = mb)
coef(summary(lm.mixed))

> coef(summary(lm.mixed))
               Estimate Std. Error         df    t value   Pr(>|t|)
(Intercept) -0.26464840 0.29013559   6.309638 -0.9121542 0.39520126
dilution.b  -0.07686173 0.07268624 125.004465 -1.0574454 0.29234760
dilution.c  -0.06208256 0.07195788 125.002425 -0.8627625 0.38992031
dilution.d  -0.13610607 0.07134742 125.004251 -1.9076523 0.05873013
dilution.e  -0.18188697 0.07134742 125.004251 -2.5493138 0.01200242

I would like to use the model to compute group differences using a contrast vector and the lmerTest::contest method but fear that because of the too small Std.error estimate and too large degrees of freedom the obtained p-values will be completely off.

Here are my questions?

• Is the mixed model I came up with correct?

If yes:

• why does the intercept Std. Error estimate is so much different from the other treatment groups std.error estimates?
• why are the denominator degrees of freedom differ even more?

If no:

• What would be a better model?

> packageVersion("lmerTest")
[1] ‘3.1.0’

Answer 1

In the two-stage approach, in the first step, you summarize the repeated measurements per peptide_Id into a single number. This inevitably leads to some information loss. To give another example, say that you measure the blood pressure of a patient ten times. Even though these measurements are correlated, they contain more information than their average. The only case in which they would contain the same amount of information as the average is when they are perfectly correlated. That is if the ten measurements are exactly the same.

In the mixed model, you are using all the data into one model, and therefore you have more information leading to smaller standard errors.

More specifically, regarding your other questions,

• For the reason mentioned above but also for other reasons (e.g., if you have missing data), the mixed model approach is better than the two-stage approach.
• The difference in the degrees of freedom is because again different pieces of information/data are used in the two approaches. Moreover, the issue of calculating the denominator degrees of freedom in mixed models is a complex one.
• What would be a better model will depend on your aims, but in any case, you should look for a model that describes the correlations in the repeated measurements of each peptide_Id sufficiently well.

Answer 2

A better model would include (1|sampleName).

The problem of the large number of degrees of freedom in the model discussed is alleviated by adding a further random effect (1|sampleName). The df are now in the same range as when doing the two stage modelling. Also the standard errors are almost the same as for the model based on summarized values.

library(lmerTest)
library(lmerTest)
mb <- readRDS(url("https://github.com/wolski/pairseqsim2/files/4072167/mb.zip"))


lm.mixed <- lmerTest::lmer(
         transformedIntensity ~ dilution. + (1|peptide_Id) + (1|sampleName),
         data = mb)
coef(summary(lm.mixed))

> coef(summary(lm.mixed))
               Estimate Std. Error        df    t value   Pr(>|t|)
(Intercept) -0.26479628 0.29171751  6.445259 -0.9077147 0.39669757
dilution.b  -0.07474148 0.08416964 15.436659 -0.8879862 0.38817524
dilution.c  -0.06362291 0.08356263 15.011593 -0.7613799 0.45822509
dilution.d  -0.13595819 0.08304649 14.671464 -1.6371336 0.12286642
dilution.e  -0.18173909 0.08304649 14.671464 -2.1884019 0.04525707

However, I do not fully understand, why it is necessary to include this random effect and is there a different way in doing it?