P values from linear and linear mixed models
As its name implies, P values are provided by the lmerTest package, which will be illustrated with the documentation example here.
require(lmerTest)
boxplot(Reaction~Days, data=sleepstudy, main="Reaction by Days",
xlab="Days", ylab="Reaction", col="blue", border="black")
We see a trend of Reaction by Days, so it is reasonable to fit a linear regression to quantify the relationship observed,
l <- lm(Reaction ~ Days, sleepstudy)
s <- summary(l)
s##
## Call:
## lm(formula = Reaction ~ Days, data = sleepstudy)
##
## Residuals:
## Min 1Q Median 3Q Max
## -110.848 -27.483 1.546 26.142 139.953
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 251.405 6.610 38.033 < 2e-16 ***
## Days 10.467 1.238 8.454 9.89e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 47.71 on 178 degrees of freedom
## Multiple R-squared: 0.2865, Adjusted R-squared: 0.2825
## F-statistic: 71.46 on 1 and 178 DF, p-value: 9.894e-15names(s)## [1] "call" "terms" "residuals" "coefficients"
## [5] "aliased" "sigma" "df" "r.squared"
## [9] "adj.r.squared" "fstatistic" "cov.unscaled"showing significant association between Reaction and Days. We now turn to the following question: how does the association alter after accounting for individual differences? The impact of Subject effect can be revealed as follows,
boxplot(Reaction~Subject, data=sleepstudy, main="Reaction by Subject",
xlab="Subject", ylab="Reaction", col="orange", border="brown")
suggesting it is more approriate to fit a random effect model, in the sense that the Subject effect randomly plays into the association of interest:
r <- lme4::lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
s <- summary(r)
s## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (Days | Subject)
## Data: sleepstudy
##
## REML criterion at convergence: 1743.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9536 -0.4634 0.0231 0.4633 5.1793
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 611.90 24.737
## Days 35.08 5.923 0.07
## Residual 654.94 25.592
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 251.405 6.824 36.843
## Days 10.467 1.546 6.771
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.138names(s)## [1] "methTitle" "objClass" "devcomp" "isLmer"
## [5] "useScale" "logLik" "family" "link"
## [9] "ngrps" "coefficients" "sigma" "vcov"
## [13] "varcor" "AICtab" "call" "residuals"
## [17] "fitMsgs" "optinfo"We see the same estimate of effect but a larger standard error for Days in the linear mixed model compared to that in the linear regression model. We then use lmer from lmerTest.
m <- lmerTest::lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
class(m)## [1] "lmerModLmerTest"
## attr(,"package")
## [1] "lmerTest"s <-summary(m)
s## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Reaction ~ Days + (Days | Subject)
## Data: sleepstudy
##
## REML criterion at convergence: 1743.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9536 -0.4634 0.0231 0.4633 5.1793
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 611.90 24.737
## Days 35.08 5.923 0.07
## Residual 654.94 25.592
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 251.405 6.824 17.005 36.843 < 2e-16 ***
## Days 10.467 1.546 16.995 6.771 3.27e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.138names(s)## [1] "methTitle" "objClass" "devcomp" "isLmer"
## [5] "useScale" "logLik" "family" "link"
## [9] "ngrps" "coefficients" "sigma" "vcov"
## [13] "varcor" "AICtab" "call" "residuals"
## [17] "fitMsgs" "optinfo"with(s,coefficients)[2,5]## [1] 3.273014e-06The P value for Days differs by orders of magnitude from that from a linear regression.