library(lme4)
## Loading required package: Matrix
mydata <- Penicillin
head(mydata)
## diameter plate sample
## 1 27 a A
## 2 23 a B
## 3 26 a C
## 4 23 a D
## 5 23 a E
## 6 21 a F
attach(mydata)
library(lattice)
xyplot(diameter ~ sample|plate, data = mydata)
mod1 <- lmer(diameter ~ 1 + (1|sample) + (1|plate))
summary(mod1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: diameter ~ 1 + (1 | sample) + (1 | plate)
##
## REML criterion at convergence: 330.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.07923 -0.67140 0.06292 0.58377 2.97959
##
## Random effects:
## Groups Name Variance Std.Dev.
## plate (Intercept) 0.7169 0.8467
## sample (Intercept) 3.7311 1.9316
## Residual 0.3024 0.5499
## Number of obs: 144, groups: plate, 24; sample, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 22.9722 0.8086 28.41
myint <- confint(mod1, level = .90)
## Computing profile confidence intervals ...
myint
## 5 % 95 %
## .sig01 0.6620839 1.115722
## .sig02 1.1716761 3.114974
## .sigma 0.4952770 0.615432
## (Intercept) 21.6087224 24.335722
a. The diameter of the zone is the response.
b. The random effects would be defined by both the plate variable and sample variable.
c. The random effects appear to be crossed in this case.
d. Small dip at b, then back up, then negative slope throughout. F generally has the lowest diameter
e. Random effects: Groups Name Variance Std.Dev. plate (Intercept) 0.7169 0.8467
sample (Intercept) 3.7311 1.9316
Residual 0.3024 0.5499
f. Estimate Std. Error t value (Intercept) 22.9722 0.8086 28.41 2.5 % 97.5 % .sig01 0.6335660 1.1821040 .sig02 1.0957893 3.5562910 .sigma 0.4858454 0.6294535 (Intercept) 21.2666274 24.6778176
g. Sample has the highest variance
library(lme4)
mydata2 <- Pastes
head(mydata2)
## strength batch cask sample
## 1 62.8 A a A:a
## 2 62.6 A a A:a
## 3 60.1 A b A:b
## 4 62.3 A b A:b
## 5 62.7 A c A:c
## 6 63.1 A c A:c
attach(mydata2)
## The following object is masked from mydata:
##
## sample
xyplot(strength ~ batch|cask, data = mydata2)
mod2 <- lmer(strength ~ 1 + (1|batch/cask))
summary(mod2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: strength ~ 1 + (1 | batch/cask)
##
## REML criterion at convergence: 247
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4798 -0.5156 0.0095 0.4720 1.3897
##
## Random effects:
## Groups Name Variance Std.Dev.
## cask:batch (Intercept) 8.434 2.9041
## batch (Intercept) 1.657 1.2874
## Residual 0.678 0.8234
## Number of obs: 60, groups: cask:batch, 30; batch, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 60.0533 0.6769 88.72
library(nlme)
##
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
##
## lmList
randslope <- lmer(strength ~ cask + (1|batch/cask))
(teststat <- 2*as.numeric(logLik(randslope)-logLik(mod2)))
## [1] 5.702696
pchisq(teststat, df = 1, lower.tail = FALSE)
## [1] 0.01693888
a. The strength of the sample of paste is the response
b. The variables of batch, cask, and sample are all random effects in this scenario.
c. These random effects are nested because there is a natural progression to these random effects.
d. Plot looks random
e. we need to use random slope because the pvalue is below .05 which is significant. Cask variance is significant.