Scatterplots with lattice
lattice::xyplot(Reaction ~ Days | Subject, data=sleepstudy,
xlab="Days of sleep deprivation",
ylab="Average reaction time (ms)",
type=c("g","p","r"))
ggplot time series by subject
ggplot(sleepstudy, aes(Days, Reaction)) +
stat_smooth(method = "lm",
formula = y ~ x,
se = FALSE,
lwd = rel(.5)) +
geom_point (size = rel(.5)) +
facet_wrap ( ~ Subject, nrow = 1) +
labs(x = "Deprivation (in Days)",
y = "Reaction time (in ms)")
ggplot average mean and se
ggplot2::ggplot(sleepstudy, aes(Days, Reaction)) +
geom_point() +
stat_smooth(method = "lm", aes(Days, Reaction), se = TRUE) +
labs(x = "Deprivation (in Days)",
y = "Reaction time (in ms)") +
theme_bw()
## `geom_smooth()` using formula 'y ~ x'
#linear model for each subject
sleeplmList <- lmList(Reaction ~ Days |Subject, data = sleepstudy)
# show the each subject intercept and slope
sleeplmList
## Call: lmList(formula = Reaction ~ Days | Subject, data = sleepstudy)
## Coefficients:
## (Intercept) Days
## 308 244.1927 21.764702
## 309 205.0549 2.261785
## 310 203.4842 6.114899
## 330 289.6851 3.008073
## 331 285.7390 5.266019
## 332 264.2516 9.566768
## 333 275.0191 9.142045
## 334 240.1629 12.253141
## 335 263.0347 -2.881034
## 337 290.1041 19.025974
## 349 215.1118 13.493933
## 350 225.8346 19.504017
## 351 261.1470 6.433498
## 352 276.3721 13.566549
## 369 254.9681 11.348109
## 370 210.4491 18.056151
## 371 253.6360 9.188445
## 372 267.0448 11.298073
##
## Degrees of freedom: 180 total; 144 residual
## Residual standard error: 25.59182
#mean int and slope match lmer Fixed effects results p.67
mean(coef(sleeplmList)[,1])
## [1] 251.4051
mean(coef(sleeplmList)[,2])
## [1] 10.46729
#these are too big as they should be,
#compare with variance random effects p.67
# mle subtracts off wobble in estimated indiv regressions Y on t
# - SSR/SST
var(coef(sleeplmList)[,1])
## [1] 838.3423
var(coef(sleeplmList)[,2])
## [1] 43.01034
# quantile of intercept and slope
quantile(coef(sleeplmList)[,1])
## 0% 25% 50% 75% 100%
## 203.4842 229.4167 258.0576 273.0255 290.1041
quantile(coef(sleeplmList)[,2])
## 0% 25% 50% 75% 100%
## -2.881034 6.194548 10.432421 13.548395 21.764702
Stem-and-leaf plot
stem(coef(sleeplmList)[,2])
##
## The decimal point is 1 digit(s) to the right of the |
##
## -0 | 3
## 0 | 23
## 0 | 56699
## 1 | 011234
## 1 | 89
## 2 | 02
Correlation between intercept and slope
cor(coef(sleeplmList)[,1], coef(sleeplmList)[,2])
## [1] -0.1375534
LMER in doing REML and MLE
# first 2 do REML, second pair match
# Bates p.67 in doing MLE (REML FALSE)
# LMER, method=REML
sleeplmer <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
summary(sleeplmer)
## 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.4634 5.1793
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 612.10 24.741
## Days 35.07 5.922 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.825 36.838
## Days 10.467 1.546 6.771
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.138
# LMER, method=MLE
sleeplmer2 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, REML=FALSE)
summary(sleeplmer2)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Reaction ~ Days + (Days | Subject)
## Data: sleepstudy
##
## AIC BIC logLik deviance df.resid
## 1763.9 1783.1 -876.0 1751.9 174
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9416 -0.4656 0.0289 0.4636 5.1793
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 565.48 23.780
## Days 32.68 5.717 0.08
## Residual 654.95 25.592
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 251.405 6.632 37.907
## Days 10.467 1.502 6.968
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.138
The end