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’