mixed-model-one-factor.Rmd

## 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