mixed-model-two-factor.Rmd

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