I have multilevel data and I would like to use mixed effects linear model.
For each subject (
sampleName variable), I have repeated measurements (
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 , est - est[-1]) names(est) <- "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)  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?
- 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?
- What would be a better model?
> packageVersion("lmerTest")  ‘3.1.0’