I will use bond data from Chapter2 of Littell et al. (2006) as a example.
#Load Packages
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(optimization)## Loading required package: Rcppdf.bonddf.bond <- data.frame(ingot = rep(1:7,each = 3),
metal = rep(c("n","i","c"),7),
pres = c(67,71.9,72.2,
67.5,68.8,66.4,
76,82.6,74.5,
72.7,78.1,67.3,
73.1,74.2,73.2,
65.8,70.8,68.7,
75.6,84.9,69))
df.bond$ingot <- factor(df.bond$ingot)
df.bond## ingot metal pres
## 1 1 n 67.0
## 2 1 i 71.9
## 3 1 c 72.2
## 4 2 n 67.5
## 5 2 i 68.8
## 6 2 c 66.4
## 7 3 n 76.0
## 8 3 i 82.6
## 9 3 c 74.5
## 10 4 n 72.7
## 11 4 i 78.1
## 12 4 c 67.3
## 13 5 n 73.1
## 14 5 i 74.2
## 15 5 c 73.2
## 16 6 n 65.8
## 17 6 i 70.8
## 18 6 c 68.7
## 19 7 n 75.6
## 20 7 i 84.9
## 21 7 c 69.0First we build linear mixed model by lme4 package with MLE method. Then we compare the parameter estimations by using equations from Some important proof of Mixed Model Equations and REML via optimization package.
lmm.bond <- lmer( pres ~ metal + (1|ingot),REML = FALSE,data = df.bond)
summary(lmm.bond)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: pres ~ metal + (1 | ingot)
## Data: df.bond
##
## AIC BIC logLik deviance df.resid
## 125.7 130.9 -57.9 115.7 16
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.45505 -0.81901 0.08048 0.52228 1.96114
##
## Random effects:
## Groups Name Variance Std.Dev.
## ingot (Intercept) 9.812 3.132
## Residual 8.890 2.982
## Number of obs: 21, groups: ingot, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 70.1857 1.6346 42.939
## metali 5.7143 1.5937 3.585
## metaln 0.9143 1.5937 0.574
##
## Correlation of Fixed Effects:
## (Intr) metali
## metali -0.488
## metaln -0.488 0.500anova(lmm.bond)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## metal 2 131.9 65.95 7.4186For linear mixed model
\[\begin{equation} \boldsymbol{Y} = \boldsymbol{X\beta} + \boldsymbol{Zu} + \boldsymbol{e} \end{equation}\]
\[\begin{equation} \begin{pmatrix} \boldsymbol{u}\\ \boldsymbol{e} \end{pmatrix} \sim \mathcal{N} \begin{pmatrix} \begin{pmatrix} \boldsymbol{0}\\ \boldsymbol{0} \end{pmatrix}, \begin{pmatrix} \boldsymbol{G} & \boldsymbol{0}\\ \boldsymbol{0} & \boldsymbol{R} \end{pmatrix} \end{pmatrix} \end{equation}\]
Design matrix \(\boldsymbol{X}\) is
matrix.x <- model.matrix(lmm.bond)
matrix.x## (Intercept) metali metaln
## 1 1 0 1
## 2 1 1 0
## 3 1 0 0
## 4 1 0 1
## 5 1 1 0
## 6 1 0 0
## 7 1 0 1
## 8 1 1 0
## 9 1 0 0
## 10 1 0 1
## 11 1 1 0
## 12 1 0 0
## 13 1 0 1
## 14 1 1 0
## 15 1 0 0
## 16 1 0 1
## 17 1 1 0
## 18 1 0 0
## 19 1 0 1
## 20 1 1 0
## 21 1 0 0
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$metal
## [1] "contr.treatment"
##
## attr(,"msgScaleX")
## character(0)Design matrix \(\boldsymbol{Z}\) is
df.z <- read.csv("files/Z.csv", header = FALSE)
matrix.z <- as.matrix(df.z)
colnames(matrix.z) <- paste0("ingot:",1:7)
matrix.z## ingot:1 ingot:2 ingot:3 ingot:4 ingot:5 ingot:6 ingot:7
## [1,] 1 0 0 0 0 0 0
## [2,] 1 0 0 0 0 0 0
## [3,] 1 0 0 0 0 0 0
## [4,] 0 1 0 0 0 0 0
## [5,] 0 1 0 0 0 0 0
## [6,] 0 1 0 0 0 0 0
## [7,] 0 0 1 0 0 0 0
## [8,] 0 0 1 0 0 0 0
## [9,] 0 0 1 0 0 0 0
## [10,] 0 0 0 1 0 0 0
## [11,] 0 0 0 1 0 0 0
## [12,] 0 0 0 1 0 0 0
## [13,] 0 0 0 0 1 0 0
## [14,] 0 0 0 0 1 0 0
## [15,] 0 0 0 0 1 0 0
## [16,] 0 0 0 0 0 1 0
## [17,] 0 0 0 0 0 1 0
## [18,] 0 0 0 0 0 1 0
## [19,] 0 0 0 0 0 0 1
## [20,] 0 0 0 0 0 0 1
## [21,] 0 0 0 0 0 0 1Define matrix \(\boldsymbol{G}\) as
\[\boldsymbol{G}_{7 \times 7} = \boldsymbol{I} \sigma^2_{\text{ingot}}\]
Define matrix \(\boldsymbol{R}\) as
\[\boldsymbol{R}_{21 \times 21} = \boldsymbol{I} \sigma^2_{\text{residual}}\]
Use equation(4.9) and equation (6.1) from Some important proof of Mixed Model Equations and REML
\[\begin{equation} \tilde{\boldsymbol{\beta}} = (\boldsymbol{X}^T \boldsymbol{V}^{-1} \boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{V}^{-1} \boldsymbol{Y} \tag{4.9} \end{equation}\]
\[\begin{equation} -2\ell(\boldsymbol{\theta};\boldsymbol{Y}) = \text{log}|\tilde{\boldsymbol{V}}(\boldsymbol{\theta})| + (\boldsymbol{Y}-\boldsymbol{X}\tilde{\boldsymbol{\beta}}(\boldsymbol{\theta}))^T \tilde{\boldsymbol{V}}(\boldsymbol{\theta})^{-1} (\boldsymbol{Y}-\boldsymbol{X}\tilde{\boldsymbol{\beta}}(\boldsymbol{\theta})) \tag{6.1} \end{equation}\]
Estimation of covariance parameters
loglikef <- function(x) {
vector.Y <- as.vector(df.bond$pres)
matrix.G <- x[1] * diag(1,nrow = 7)
matrix.R <- x[2] * diag(1,nrow = 21)
matrix.V <- matrix.z %*% matrix.G %*% t(matrix.z) + matrix.R
vector.Beta <- solve(t(matrix.x) %*% solve(matrix.V) %*% matrix.x) %*% t(matrix.x) %*% solve(matrix.V) %*% vector.Y
loglike <- log(det(matrix.V)) + t(vector.Y - matrix.x %*% vector.Beta) %*% solve(matrix.V) %*% (vector.Y - matrix.x %*% vector.Beta)
return(c(loglike))
}
optim_nm(fun = loglikef, k = 2, start = c(1,1),maximum= FALSE, tol = 0.0000000000001)$par##
## 9.812383 8.889932hat.cov.para <- optim_nm(fun = loglikef, k = 2, start = c(1,1),maximum= FALSE, tol = 0.0000000000001)$parWe have that
\[\sigma^2_{\text{ingot}} = 9.812383 \text{ and } \sigma^2_{\text{residual}} = 8.889932\]
This is identical to estimations from lme4.
Then estimate fix effect \(\boldsymbol{\beta}\)
hat.beta <- function(x) {
vector.Y <- as.vector(df.bond$pres)
matrix.G <- x[1] * diag(1,nrow = 7)
matrix.R <- x[2] * diag(1,nrow = 21)
matrix.V <- matrix.z %*% matrix.G %*% t(matrix.z) + matrix.R
vector.Beta <- solve(t(matrix.x) %*% solve(matrix.V) %*% matrix.x) %*% t(matrix.x) %*% solve(matrix.V) %*% vector.Y
return(vector.Beta)
}
round(hat.beta(x = hat.cov.para),4)## [,1]
## (Intercept) 70.1857
## metali 5.7143
## metaln 0.9143We have that \[\boldsymbol{\beta} = \begin{bmatrix} \beta_0\\ \beta_1\\ \beta_2 \end{bmatrix} = \begin{bmatrix} 70.1857\\ 5.7143\\ 0.9143 \end{bmatrix}\]
Again, this is identical to estimations from lme4.
Same as section 3, firstly we build linear mixed model by lme4 package with REML option. Then we compare โhandโ calculated estimations with estimations from lme4.
relmm.bond <- lmer( pres ~ metal + (1|ingot),REML = TRUE,data = df.bond)
summary(relmm.bond)## Linear mixed model fit by REML ['lmerMod']
## Formula: pres ~ metal + (1 | ingot)
## Data: df.bond
##
## REML criterion at convergence: 107.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.34712 -0.75825 0.07451 0.48354 1.81566
##
## Random effects:
## Groups Name Variance Std.Dev.
## ingot (Intercept) 11.45 3.383
## Residual 10.37 3.220
## Number of obs: 21, groups: ingot, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 70.1857 1.7655 39.754
## metali 5.7143 1.7214 3.320
## metaln 0.9143 1.7214 0.531
##
## Correlation of Fixed Effects:
## (Intr) metali
## metali -0.488
## metaln -0.488 0.500Using equation(7.3) from Some important proof of Mixed Model Equations and REML
\[\begin{equation} -2\ell_{REML}(\boldsymbol{\theta};\boldsymbol{Y}) = \text{log} |\tilde{\boldsymbol{V}}(\boldsymbol{\theta})| + \text{log} |\boldsymbol{X}^T \tilde{\boldsymbol{V}}^{-1} \boldsymbol{X}| + (\boldsymbol{Y} - \boldsymbol{X}\tilde{\boldsymbol{\beta}}(\boldsymbol{\theta}))^T \tilde{\boldsymbol{V}}^{-1} (\boldsymbol{Y} - \boldsymbol{X}\tilde{\boldsymbol{\beta}}(\boldsymbol{\theta})) \tag{7.3} \end{equation}\]
reloglikef <- function(x) {
vector.Y <- as.vector(df.bond$pres)
matrix.G <- x[1] * diag(1,nrow = 7)
matrix.R <- x[2] * diag(1,nrow = 21)
matrix.V <- matrix.z %*% matrix.G %*% t(matrix.z) + matrix.R
vector.Beta <- solve(t(matrix.x) %*% solve(matrix.V) %*% matrix.x) %*% t(matrix.x) %*% solve(matrix.V) %*% vector.Y
loglike <- log(det(matrix.V)) + log(det(t(matrix.x) %*% solve(matrix.V) %*% matrix.x)) + t(vector.Y - matrix.x %*% vector.Beta) %*% solve(matrix.V) %*% (vector.Y - matrix.x %*% vector.Beta)
return(c(loglike))
}
optim_nm(fun = reloglikef, k = 2, start = c(5,6),maximum= FALSE, tol = 0.0000000000001)$par##
## 11.44778 10.37159Again, parameter estimations are close enough.
Littell et al. (2006)