-
Notifications
You must be signed in to change notification settings - Fork 17
/
mixed-model-two-factor.Rmd
120 lines (74 loc) · 2.98 KB
/
mixed-model-two-factor.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
## Two-factor Mixed Model
An approach for two factor random effects model via maximum likelihood in R Matlab and Julia. It's based on Statistical Modeling and Computation (2014) Chapter 10, example 10.10. See the [previous chapter][One-factor Mixed Model] for a one factor model, and the [Supplemental Section][Two-factor] for the Matlab and Julia versions of this example. Note that the text has a typo on the `sigma2` variance estimate (value should be .0023 not .023).
### Data Setup
The data regards the breeding value of a set of five sires in raising pigs. Each
sire is mated to a random group of dams, with the response being the average
daily weight gain in pounds of two piglets in each litter.
```{r two-factor-re-setup}
library(tidyverse)
y = c(1.39,1.29,1.12,1.16,1.52,1.62,1.88,1.87,1.24,1.18,
.95,.96,.82,.92,1.18,1.20,1.47,1.41,1.57,1.65)
# for use in lme4, but also a more conceptual representation of the data
d = expand.grid(sire = rep(1:5, 2), dam = 1:2)
d = data.frame(d[order(d$sire), ], y)
```
### Function
The function takes the log variances `eta` as input to keep positive.
```{r two-factor-re}
two_factor_re <- function(mu, eta_alpha, eta_gamma, eta) {
# Args
# mu: intercept
# eta_alpha: random effect one
# eta_gamma: random effect two
# eta: residual variance of y
sigma2_alpha = exp(eta_alpha)
sigma2_gamma = exp(eta_gamma)
sigma2 = exp(eta)
n = length(y)
# covariance matrix of observations
Sigma = sigma2 * diag(n) + sigma2_alpha * tcrossprod(Xalpha) +
sigma2_gamma * tcrossprod(Xgamma)
# log likelihood
ll = -n / 2 * log(2 * pi) - sum(log(diag(chol(Sigma)))) -
.5 * t(y - mu) %*% chol2inv(chol(Sigma)) %*% (y - mu)
return(-ll)
}
```
### Estimation
Starting values and test.
```{r two-factor-re-starts}
starts = list(
mu = mean(y),
eta_alpha = var(tapply(y, d$sire, mean)),
eta_gamma = var(y) / 3,
eta = var(y) / 3
)
Xalpha = diag(5) %x% rep(1, 4)
Xgamma = diag(10) %x% rep(1, 2)
```
Estimate at starting values.
```{r two-factor-re-est}
two_factor_re(
mu = starts[[1]],
eta_alpha = starts[[2]],
eta_gamma = starts[[3]],
eta = starts[[4]]
)
```
Package <span class="pack" style = "">bbmle</span> has an <span class="func" style = "">mle2</span> function for maximum likelihood estimation based on underlying R functions like <span class="func" style = "">optim</span>, and produces a nice summary table. LBFGS-B is used to place lower bounds on the variance estimates.
```{r two-factor-re-bbmle}
library(bbmle)
fit_mle = mle2(two_factor_re, start = starts, method = 'BFGS')
```
### Comparison
We can compare to the <span class="pack" style = "">lme4</span> model result.
```{r two-factor-re-compare}
### lme4 comparison
library(lme4)
fit_mer = lmer(y ~ (1 | sire) + (1 | dam:sire), d, REML = FALSE)
summary(fit_mle)
exp(coef(fit_mle)[-1])
summary(fit_mer)
```
### Source
Original code available at https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/Mixed%20Models/two_factor_RE.R