Linear Mixed Models
Here is an example of a linear mixed model using the \(abrasion\) data set
library(faraway)
## Warning: package 'faraway' was built under R version 3.6.3
data(abrasion)
head(abrasion)
## run position material wear
## 1 1 1 C 235
## 2 1 2 D 236
## 3 1 3 B 218
## 4 1 4 A 268
## 5 2 1 A 251
## 6 2 2 B 241
This data set examines the amount of wear that occurs on different materials. To do this, four materials were fed into a wear testing machine. The response is \(wear\) or the measured loss of weight in .1mm over the testing period. Random effects are defined by the \(run\) and \(position\) variables. This is a crossed random effect structure.
Here is a boxplot of \(wear\) by the different types of material:
boxplot(wear~material,data=abrasion)
Examining this, it appears that material A exhibits the most amount of wear and is different than the other three types.
Here is a model that finds the mean wear for each type and includes the crossed random effects:
library(lme4)
## Warning: package 'lme4' was built under R version 3.6.3
## Loading required package: Matrix
gmean5<-lmer(wear~0+material+(1|position)+(1|run),data = abrasion)
summary(gmean5)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wear ~ 0 + material + (1 | position) + (1 | run)
## Data: abrasion
##
## REML criterion at convergence: 100.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.08973 -0.30231 0.02697 0.42254 1.21052
##
## Random effects:
## Groups Name Variance Std.Dev.
## position (Intercept) 107.06 10.347
## run (Intercept) 66.90 8.179
## Residual 61.25 7.826
## Number of obs: 16, groups: position, 4; run, 4
##
## Fixed effects:
## Estimate Std. Error t value
## materialA 265.750 7.668 34.66
## materialB 220.000 7.668 28.69
## materialC 241.750 7.668 31.53
## materialD 230.500 7.668 30.06
##
## Correlation of Fixed Effects:
## matrlA matrlB matrlC
## materialB 0.740
## materialC 0.740 0.740
## materialD 0.740 0.740 0.740
Here is how to test whether the mean wear is significantly different between the four types:
comparemod<-lmer(wear~material+(1|position)+(1|run),data = abrasion)
nullmod<-lmer(wear~1+(1|position)+(1|run),data = abrasion)
## boundary (singular) fit: see ?isSingular
lrt<-2*(as.numeric(logLik(comparemod))-as.numeric(logLik(nullmod)))
pchisq(lrt,df=3,lower.tail = FALSE)
## [1] 2.574975e-08
This is comparing two models. The first is a grand mean model with the two crossed random effects for \(position\) and \(run\). The second model adds 3 additional coefficients for the different types of material. So a significant p-value would suggest that at least one type of material has a different amount of wear. By performing the Likelihood Ratio Test, it appears that the mean wear is in fact significantly different between the four types.