-
Notifications
You must be signed in to change notification settings - Fork 17
/
mixed-model-one-factor.Rmd
132 lines (94 loc) · 2.99 KB
/
mixed-model-one-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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
## One-factor Mixed Model
The following is an approach for one factor random effects model via maximum
likelihood in R (and Matlab and Julia in the [Supplemental
Section][One-factor]). It's based on Statistical Modeling and Computation (2014)
Chapter 10, example 10.10. Unfortunately I did this before knowing they had both
Matlab and R code on their website, though the R code here is a little cleaner
and has comments.
### Data Setup
The data regards crop yield from 10 randomly selected
locations and three collections at each location.
```{r one-factor-re-setup}
library(tidyverse)
y = matrix(c(22.6, 20.5, 20.8,
22.6, 21.2, 20.5,
17.3, 16.2, 16.6,
21.4, 23.7, 23.2,
20.9, 22.2, 22.6,
14.5, 10.5, 12.3,
20.8, 19.1, 21.3,
17.4, 18.6, 18.6,
25.1, 24.8, 24.9,
14.9, 16.3, 16.6),
10, 3, byrow = TRUE)
```
### Function
The estimating function.
```{r one-factor-re}
one_factor_re <- function(mu, sigma2_mu, sigma2){
# Args
# mu: intercept
# sigma2_mu: variance of intercept
# sigma2: residual variance of y
# I follow their notation for consistency
d = nrow(y)
ni = ncol(y)
# covariance matrix of observations
Sigmai = sigma2 * diag(ni) + sigma2_mu * matrix(1, ni, ni)
# log likelihood
l = rep(NA, 10)
# iterate over the rows
for(i in 1:d){
l[i] = .5 * t(y[i, ] - mu) %*% chol2inv(chol(Sigmai)) %*% (y[i, ] - mu)
}
ll = -(ni*d) / 2*log(2*pi) - d / 2*log(det(Sigmai)) - sum(l)
return(-ll)
}
```
### Estimation
Starting values.
```{r one-factor-re-starts}
starts = list(
mu = mean(y),
sigma2_mu = var(rowMeans(y)),
sigma2 = mean(apply(y, 1, var))
)
```
Estimate at the starting values.
```{r one-factor-re-est}
one_factor_re(mu = starts[[1]],
sigma2_mu = starts[[2]],
sigma2 = starts[[3]])
```
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 one-factor-re-bbmle}
library(bbmle)
fit_mle = mle2(
one_factor_re ,
start = starts,
method = 'L-BFGS-B',
lower = c(
mu = -Inf,
sigma2_mu = 0,
sigma2 = 0
),
trace = TRUE
)
```
### Comparison
We can compare to the <span class="pack" style = "">lme4</span> model result.
```{r one-factor-re-compare}
library(lme4)
library(tidyverse)
d = data.frame(y) %>%
pivot_longer(everything(), names_to = 'x', values_to = 'value') %>%
arrange(x) %>%
group_by(x) %>%
mutate(group = 1:n())
fit_mer = lmer(value ~ 1 | group, data = d, REML = FALSE)
summary(fit_mle)
summary(fit_mer)
-2 * logLik(fit_mer)
```
### Source
Original code available at https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/Mixed%20Models/one_factor_RE.R