-
Notifications
You must be signed in to change notification settings - Fork 82
/
one_factor_RE.R
105 lines (83 loc) · 2.79 KB
/
one_factor_RE.R
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
###################################################################################
### An approach for one 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; 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. The data regards crop yield from 10 randomly ###
### selected locations and three collections at each. See one_factor_RE.m and ###
### one_factor_RE.jl for the related Matlab and Julia files, and the respective ###
### twofactorRE.* for the associated two factor random effects examples. ###
###################################################################################
#####################
### Main function ###
#####################
one_factor_re_loglike = function(mu, sigma2_mu, sigma2){
# Args are 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)
}
###################
### Data set up ###
###################
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=T)
################################
### Starting values and test ###
################################
starts = list(
mu = mean(y),
sigma2_mu = var(rowMeans(y)),
sigma2 = mean(apply(y, 1, var))
)
### test
one_factor_re_loglike(starts[[1]], starts[[2]], starts[[3]])
#######################
### Run and compare ###
#######################
### bbmle has mle2 function for maximum likelihood estimation based on
### underlying R functions like optim. LBFGS-B is used to place lower bounds on
### the variance estimates
library(bbmle)
mlout = mle2(
one_factor_re_loglike,
start = starts,
method = 'L-BFGS-B',
lower = c(
mu = -Inf,
sigma2_mu = 0,
sigma2 = 0
),
trace = T
)
### 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())
lme = lmer(value ~ 1 | group, data = d, REML = F)
summary(mlout)
summary(lme)
-2 * logLik(lme)