This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0/).
The R markdown source for this document is https://www.stat.umn.edu/geyer/5421/notes/random.Rmd.
Another name for random effects models is mixed models.
CRAN has a task view https://cran.r-project.org/view=MixedModels listing many packages doing some sort of random effects models.
We use four of them.
The version of R used to make this document is 4.6.1.
The version of the bookdown package used to make this document is
0.47.
The version of the rmarkdown package used to make this document is
2.31.
The version of the CatDataAnalysis package used to make this document is
0.1.5.
The version of the MASS package used to make this document is
7.3.66.
The version of the lme4 package used to make this document is
2.0.6.
The version of the glmm package used to make this document is
1.4.5.
The version of the parallel package used to make this document is
4.6.1.
The version of the aster package used to make this document is
1.5.0.
The version of the numDeriv package used to make this document is
2016.8.1.1.
library(CatDataAnalysis)
library(MASS)
library(lme4)## Loading required package: Matrixlibrary(glmm)## Loading required package: trust## Loading required package: mvtnorm## Loading required package: parallel## Loading required package: doParallel## Loading required package: foreach## Loading required package: iteratorslibrary(parallel)
library(aster)
# need bug fix
stopifnot(compareVersion(as.character(packageVersion("aster")), "1.3-8") >= 0)
library(numDeriv)Set random number generator seed for reproducibility.
set.seed(42)A statistical model becomes a random effects model when you start with a parametric statistical model and the reconsider some of the parameters turning them into random variables. To distinguish
the random variables that used to be parameters are now called random effects and
the parameters that are unchanged are now called fixed effects.
In order to complete the specification of a statistical model, we need a distribution for the random effects. Usually, they are considered independent and normally distributed but sometimes also dependent and normally distributed.
When we have both fixed effects and random effects, this is also called a mixed model. And when the model we started with was a linear model (LM) we say we have a linear mixed model (LMM).
LMM are not really in what this course is about, since they have continuous rather than categorical response. But they are much simpler, so we start there to get some intuition. We are following Oehlert (2010), Section 11.1 and 11.2.
Here is a simple example. Suppose we have data on patients and each patient is observed multiple times. Typically there will be some variability of the observations on one patient, and this may be much less than the between patient variability. So we have two levels of variability.
Here is the model in equations. Denote the components of the response vector as a matrix \(y_{i j}\) where \(i\) indexes patients and \(j\) indexes different measurements on each patient. Then we assume \[\begin{equation} y_{i j} = \mu_i + b_i + e_{i j} \tag{4.1} \end{equation}\] where
the \(\mu_i\) are the fixed effects,
the \(b_i\) are the random effects, assumed to be independent normal with mean zero and variance \(\sigma^2\), and
the \(e_{i j}\) are the errors, assumed to be independent normal with mean zero and variance \(\tau^2\).
The components of the response vector are still multivariate normal but are no longer independent. Rather we have \[\begin{alignat*}{2} \mathop{\rm var}(y_{i j}) & = \sigma^2 + \tau^2 \\ \mathop{\rm cov}(y_{i j}, y_{i k}) & = \sigma^2, & \qquad j \neq k \\ \mathop{\rm cov}(y_{i j}, y_{k m}) & = 0, & \qquad i \neq j \end{alignat*}\]
Estimating this is more or less like conventional ANOVA with minor differences. Oehlert (2010), Display 11.1, gives the following ANOVA table
| Source | DF | EMS |
| Patients | \(m - 1\) | \(\tau^2 + m \sigma^2\) |
| Error | \(N - m\) | \(\tau^2\) |
Here \(N\) is the total number of observations and \(m\) is the number of observations per individual.
This means the error sum of squares \[ \tfrac{1}{N - m} \sum\nolimits_{i j} (y_{i j} - \hat{\mu}_i)^2 \] is an unbiased estimator of \(\tau^2\), where \(\hat{\mu}_i\) is the average (sample mean) of the observations for patient \(i\). And the patient sum of squares \[ \tfrac{1}{m - 1} \sum\nolimits_i ( \hat{\mu}_i - \hat{\hat{\mu}})^2 \] where \(\hat{\hat{\mu}}\) is the average (sample mean) of all the observations, is an unbiased estimator of \(\tau^2 + m \sigma^2\).
But, of course, we don’t want to estimate that, so we subtract \[\begin{equation} \frac{1}{m} \left[ \left( \tfrac{1}{m - 1} \sum\nolimits_i ( \hat{\mu}_i - \hat{\hat{\mu}})^2 \right) - \left( \tfrac{1}{N - m} \sum\nolimits_{i j} (y_{i j} - \hat{\mu}_i)^2 \right) \right] \tag{4.2} \end{equation}\] is an unbiased estimate of \(\sigma^2\).
We already see a complication. Because of the subtraction in (4.2), it can give a negative estimate of variance, which is ridiculous.
So even in this simplest case, things are trickier than they seem at first sight. We should write down the log likelihood for the model and maximize it subject to the constraints \(\sigma^2 \ge 0\) and \(\tau^2 \ge 0\). We can always write down the log likelihood because the joint distribution of \(y\) and \(u\) is multivariate normal, hence so is the marginal distribution of \(y\) (because every marginal of a multivariate normal is again multivariate normal). And we can write as an R function the log likelihood for any multivariate normal distribution. It involves matrix inverses and determinants, but R knows how to do them. \[ l(\theta) = - \tfrac{1}{2} \log(\mathop{\rm det}(M(\theta))) - \tfrac{1}{2} (y - \mu(\theta))^T M(\theta)^{-1} (y - \mu(\theta)) \] where the distribution of \(y\) is multivariate normal with mean vector \(\mu(\theta)\) and variance matrix \(M(\theta)\).
The next step is generalized linear mixed models (GLMM), which were introduced by Stiratelli et al. (1984), although earlier authors had steps in that direction. We get GLMM from GLM in the same way we get LMM from LM: start with a parametric statistical model and the reconsider some of the parameters turning them into random variables.
Start with the model equation for a canonical affine submodel, which we copy here \[\begin{equation} \theta = a + M \beta \tag{4.3} \end{equation}\] where \(\theta\) is the saturated model canonical parameter vector, \(a\) the offset vector, \(M\) the model matrix, and \(\beta\) the submodel canonical parameter vector.
Now we want to turn some components into random effects. We denote this by writing what was \(\beta\) as a partitioned vector \((\beta, b)\) where now \(\beta\) is the fixed effects vector and \(b\) is the random effects vector. This means we have to also write the model matrix as a partitioned matrix \((M, Z)\) where now \(M\) is the model matrix for fixed effects and \(Z\) is the model matrix for random effects. \[\begin{equation} \theta = a + M \beta + Z b \tag{4.4} \end{equation}\] We could divide up the random effects even farther, writing \[\begin{equation} \theta = a + M \beta + \sigma_1 Z_1 u_1 + \cdots + \sigma_k Z_k u_k \tag{4.5} \end{equation}\] where now the random effects \(u_1\), \(\ldots,\) \(u_k\) all have independent standard normal components, and we see the standard deviation parameters explicitly.
In (4.4) the variance parameters are not explicit. They are implicit in \[ D = \mathop{\rm var}(b) \] Since the components of \(b\) are independent, \(D\) is a diagonal matrix. Let \(\nu_i = \sigma_i^2\). These parameters are called variance components. Thus the parameters of a random effects model are (to a frequentist)
the fixed effects (the components of \(\beta\)) and
the variance components (the \(\nu_i\)) or their square roots (the \(\sigma_i)\)
(of course, Bayesians also consider the random effects parameters). To keep straight which variance components go with which random effects, we define \[ E_i = \frac{\partial D}{\partial \nu_i} \] so \[ D = \sum_{i = 1}^k \nu_i E_i \] and we see that each diagonal component of \(D\) is some variance component and each \(E_i\) is a diagonal matrix with zero-or-one-valued diagonal components and \(E_i E_j = 0\) if \(i \neq j\) and \(\sum_i E_i\) is the identity matrix.
Note well that although we are assuming \(b\) has independent components (so \(D\) is diagonal), this does not mean that \(Z b\) has independent components. Its variance matrix is \[ \mathop{\rm var}(Z b) = E (Z b b^T Z^T) = Z \mathop{\rm var}(b b^T) Z^T = Z^T D Z \] and that is not (in general) a diagonal matrix. Hence the components of \(Z b\) are not independent.
The other parts of a GLMM are the same as a GLM (the link function, the family, and so forth).
All of the preceding section except the last sentence extends to arbitrary exponential family models Geyer et al. (2013), with one proviso. The full canonical parameter space of the model need not be a whole vector space (the negative binomial distribution provides an example). But the random effects (if assumed normal) take values that can be any real number. Thus the proviso is that the function (4.4) map into the full canonical parameter space for any real-valued vector \(b\).
A random effects model is a statistical in which some of the random variables are observed and other random variables are not. The latter are called random effects. But the same kind of model occurs in many contexts and these unobserved random variables have been given many different names in different contexts: missing data, latent variables, hidden variables, parameters (in so-called empirical Bayes)
It does not matter what they are called, the data comes in two parts \(x\) and \(y\), where \(y\) is the observed data and \(x\) is unobserved, whether unobservable (when we say random effects, latent variables, or hidden variables) or whether could have been observed but wasn’t (when we say missing data). And we assume it is easy to write down the joint distribution.
But likelihood theory says the likelihood is the probability of the observed data thought of as a function of the parameter vector. So we must eliminate \(x\), that is, we need the marginal distribution of \(y\). When \(x\) is a continuous random vector, we need to do an integral. \[\begin{equation} f_\theta(y) = \int f_\theta(x, y) \, d x \tag{5.1} \end{equation}\] where, if \(x\) is a vector, this is a multiple integral. Then taking logs gives the log likelihood \[\begin{equation} l(\theta) = \log \left( \int f_\theta(x, y) \, d x \right) \tag{5.2} \end{equation}\] There is no question that this is the Right Thing with a capital R and a capital T. Likelihood theory (our course notes here and here and Geyer (2013) and references cited therein) says that maximizers of (5.2) are consistent and asymptotically normal with asymptotic variance matrix inverse Fisher information matrix, and no estimators can be (asymptotically) better.
But there is a problem (except for LMM, where we can do the integral in (5.2) because every marginal of a multivariate normal distribution is again multivariate normal): we can almost never do the integral in (5.2); it is not found in any calculus book.
We can do the integral by Monte Carlo, either Markov chain Monte Carlo (MCMC) (Thompson and Guo (1991); Geyer (1994); Booth and Hobert (1999)) or ordinary Monte Carlo (OMC) (Sung and Geyer (2007); Knudson et al. (2021)).
We can do the integral by so-called numerical integration (https://cran.r-project.org/view=NumericalMathematics) but when the integral is high-dimensional, this is too slow. There is one special case where we can do numerical integration: when there is only one random effect involved in each observation. Then the multiple integral breaks up into a product of univarate integrals, and these can be done efficiently by numerical integration. There is even a special method for integrating out normal random effects called Gauss-Hermite quadrature (CRAN package fastGHQuad to do these univariate integrals. But this method is no help when there is more than one random effect involved in each observation (so-called crossed random effects).
We can do the integral by so-called Laplace approximation Wikipedia article following Breslow and Clayton (1993) and Geyer et al. (2013). This works as follows. Let \(l(\,\cdot\,)\) be the saturated model log likelihood and \(W(\,\cdot\,)\) be its second derivative. Then the integral (5.1) becomes \[\begin{equation*} \int \exp\left( l(a + M \beta + Z b) - \tfrac{1}{2} b^T D^{-1} b - \tfrac{1}{2} \log \mathop{\rm det}(D) \right) \, d b \end{equation*}\] The method of Laplace approximation procedes in two steps. First we find the point \(b\) where the integrand is highest; call that \(b^*\). Note that \(b^*\) is a function of \(\beta\) and \(\nu\) even though the notation does not indicate that. Then replace the log of the integrand with its best quadratic approximation (the Taylor series up to quadratic terms). This gives the integral \[\begin{equation} \int \exp\left( l(a + M \beta + Z b^*) + \tfrac{1}{2} (b - b^*)^T \bigl[ Z^T W(a + M \beta + Z b^*) Z - D^{-1} \bigr] (b - b^*) - \tfrac{1}{2} \log \mathop{\rm det}(D) \right) \, d b \tag{5.3} \end{equation}\] The point of this is that because now the integrand in (5.3) is now the exponential of a function quadratic in \(b\). Hence it is the kind of integral that occurs in calculating the marginal distribution of a multivariate normal distribution. Hence it can be done in closed form giving \[\begin{equation} l_\text{approx}(\beta, \nu) = l(a + M \beta + Z b^*) - \tfrac{1}{2} (b^*)^T D^{-1} b^* - \log \mathop{\rm det} \bigl[ Z^T W(a + M \beta + Z b^*) Z D + \text{Id} \bigr] \tag{5.4} \end{equation}\] where \(\text{Id}\) denotes the identity matrix of the appropriate dimension. But this is still too hard to maximize because we need derivatives of \(l_\text{approx}(\,\cdot\,)\) to find solutions and Fisher information, and don’t know derivatives of \(W(\,\cdot\,)\) (they do exist but are horribly messy). Hence we make, following Breslow and Clayton (1993), yet another approximation, treating \(W(\,\cdot\,)\) as a constant function having the value \(W^*\). This gives us \[\begin{equation} l_\text{approx}(\beta, \nu) = l(a + M \beta + Z b^*) - \tfrac{1}{2} (b^*)^T D^{-1} b^* - \log \mathop{\rm det} \bigl[ Z^T W^* Z D + \text{Id} \bigr] \tag{5.5} \end{equation}\] as our approximation to the log likelihood for the random effects model.
So we maximize (5.5) to find (approximate) maximum likelihood parameter estimates \(\hat{\beta}\) and \(\hat{\nu}\). But then we note that we used the wrong \(W^*\). So we set \[ W^* = W(a + M \hat{\beta} + Z b^*(\hat{\beta}, \hat{\nu})) \] and then re-optimize (5.5) and keep repeating this until \(W^*\) does not change (much).
This gives us estimates (even standard errors, Geyer et al. (2013), Section 2) but they are based on approximations of unknown validity. It is true that log likelihoods get approximately quadratic as sample size goes to infinity (so says the theory of maximum likelihood, Geyer (2013), Appendix C) but if \(b\) is a high-dimensional vector the sample size may need to be exceedingly large for Laplace approximation to be good. The second approximation, treating \(W(\,\cdot\,)\) as a constant function follows from the first: a quadratic function does have a constant second derivative function. But again sample size may need to be exceedingly large for this approximation to be good.
As always, if you distrust these large sample approximations, bootstrap! That what Geyer et al. (2013), Sections 5 and 8, recommend. Unfortunately, fitting these models is time-consuming, and bootstrap is even more so. Thus most users don’t bootstrap.
The main point of this section is not the mathematical details. It is that this is an approximation of unknown validity. We have changing from doing the actual random effects problem to an easier problem that is (hopefully) related.
R function glmmPQL in R package MASS uses Laplace approximation.
R function glmer in R package lme4 uses Laplace approximation if
there are crossed random effects and Gauss-Hermite quadrature otherwise.
R function glmm in R package glmm uses ordinary Monte Carlo.
R function reaster in R package aster uses Laplace approximation.
We use the following data from Agresti.
data(table_13.2)
head(table_13.2)## state pi n x proportion
## 1 AK 0.379 5 3 0.6000000
## 2 AL 0.387 29 9 0.3103448
## 3 AR 0.389 17 2 0.1176471
## 4 AZ 0.449 35 13 0.3714286
## 5 CA 0.609 207 129 0.6231884
## 6 CO 0.537 37 16 0.4324324resp <- with(table_13.2, cbind(x, n - x))Clearly what is \(T\) in the book is \(n\) in R.
What we have here is overdispersed binomial data.
Each state has a different proportion of voters voting for Obama.
So this is a competitor for analyses using R function glm with
family quasibinomial. The difference is that here we have a
real statistical model (the random effects model), and the quasi-likelihood
approach does not.
So we could have used any of the analyses below for problem 4.13 in Agresti (which was our homework problem 4-4).
gout <- glmmPQL(resp ~ 1, ~ 1 | state, family = binomial, data = table_13.2)## iteration 1## iteration 2## iteration 3summary(gout)## Linear mixed-effects model fit by maximum likelihood
## Data: table_13.2
## AIC BIC logLik
## NA NA NA
##
## Random effects:
## Formula: ~1 | state
## (Intercept) Residual
## StdDev: 0.1974425 1.429759
##
## Variance function:
## Structure: fixed weights
## Formula: ~invwt
## Fixed effects: resp ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 0.06679477 0.07459273 51 0.8954596 0.3747
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.1297601 -0.6454223 0.1034602 0.5145442 1.5567157
##
## Number of Observations: 51
## Number of Groups: 51That is clear as mud, but (from the examples in the help) the way this function is supposed to be used.
gout.too <- glmer(resp ~ 1 + 1 | state, family = binomial, data = table_13.2)
sout <- summary(gout.too)
sout## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: resp ~ 1 + 1 | state
## Data: table_13.2
##
## AIC BIC logLik -2*log(L) df.resid
## 288.9 292.8 -142.5 284.9 49
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9780 -0.3909 0.1121 0.3454 1.7044
##
## Random effects:
## Groups Name Variance Std.Dev.
## state (Intercept) 0.1721 0.4149
## Number of obs: 51, groups: state, 51
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.03027 0.08182 0.37 0.711This makes a lot more sense. The fixed effect estimate
sout$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.03027032 0.08181548 0.3699828 0.7113953corresponds to the mean value parameter
invlogit <- function(theta) 1 / (1 + exp(- theta))
sout$coefficients[ , "Estimate"] |> invlogit()## [1] 0.507567And we can make a Wald confidence interval.
conf.level <- 0.95
crit <- qnorm((1 + conf.level) / 2)
crit## [1] 1.959964(sout$coefficients[ , "Estimate"] + c(-1, 1) * crit *
sout$coefficients[ , "Std. Error"]) |> invlogit()## [1] 0.4675245 0.5475126using the trick of mapping the Wald interval from the canonical to the mean value parameter that we saw in Chapter Zero.
Note that this confidence interval is very different from the intro stats interval
pi.hat <- with(table_13.2, sum(x) / sum(n))
pi.hat## [1] 0.5245n <- with(table_13.2, sum(n))
(pi.hat + c(-1, 1) * crit * sqrt(pi.hat * (1 - pi.hat) / n))## [1] 0.5026133 0.5463867Both the point estimate and the confidence interval are different.
Our interval from glmer uses the heterogeneity of the states.
The intro stats interval assumes the states are homogenous
(all have exactly the same proportion of vote for this candidate,
which is an obviously false assumption).
gout.over <- glm(resp ~ 1, family = quasibinomial, data = table_13.2)
pout.over <- predict(gout.over, se.fit = TRUE, type = "response",
newdata = data.frame(state = "AK"))
# point estimate
pout.over$fit## 1
## 0.5245# 95% confidence interval
pout.over$fit + c(-1, 1) * crit * pout.over$se.fit## [1] 0.4900696 0.5589304We see, as explained in the notes about quasi-likelihood, that quasi-likelihood has the same point estimates as likelihood, even though this is wrong for any actual model of overdispersion. We have now looked at two models of overdispersion
negative binomial regression (our homework problem 2-4) and
random effects (this section)
and both change the point estimate, not just its standard error.
gout.too.too <- glmm(resp ~ 1, ~ 0 + state, data = table_13.2,
family = binomial.glmm, m = 1e6, varcomps.names = "state")
summary(gout.too.too)##
## Call:
## glmm(fixed = resp ~ 1, random = ~0 + state, varcomps.names = "state",
## data = table_13.2, family.glmm = binomial.glmm, m = 1e+06)
##
##
## Link is: "logit (log odds)"
##
## Fixed Effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0 0 0.359 0.72
##
##
## Variance Components for Random Effects (P-values are one-tailed):
## Estimate Std. Error z value Pr(>|z|)/2
## state 0 0 2.482 0.00653 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This looks like nonsense, but the estimates are not actually zero. A look inside the code reveals that it is trying to round the printout to the number of significant figures it accurately calculates, and here that is no significant figures. The estimates themselves are not actually zero.
gout.too.too$beta## (Intercept)
## 0.03045109gout.too.too$nu## state
## 0.1767598Nor are the estimates of standard errors
vcov(gout.too.too) |> diag() |> sqrt()## (Intercept) state
## 0.08490472 0.07121765Nor are the estimates of Monte Carlo standard errors
mcse(gout.too.too)## (Intercept) state
## NaN NaNHmmmm. Something wrong there. Maybe that is the source of the bug.
R package aster is not really designed for this, but can do it.
dat <- with(table_13.2,
data.frame(root = n, resp = x, state = state,
varb = factor(rep("only", length(state))),
id = seq_along(state)))
head(dat)## root resp state varb id
## 1 5 3 AK only 1
## 2 29 9 AL only 2
## 3 17 2 AR only 3
## 4 35 13 AZ only 4
## 5 207 129 CA only 5
## 6 37 16 CO only 6rout <- reaster(resp ~ 1, list(state = ~ 0 + state), pred = 0, fam = 1,
varb, id, root, data = dat)
summary(rout, standard.deviation = FALSE)##
## Call:
## reaster.formula(fixed = resp ~ 1, random = list(state = ~0 +
## state), pred = 0, fam = 1, varvar = varb, idvar = id, root = root,
## data = dat)
##
##
## Fixed Effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.03035 0.08047 0.377 0.706
##
## Variance Components (P-values are one-tailed):
## Estimate Std. Error z value Pr(>|z|)/2
## state 0.16512 0.07049 2.342 0.00958 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We see that the estimates and standard errors more or less agree between
R function glmer and R function reaster which is perhaps no suprise
since they may be using the same method (unless glmer is using
Gauss-Hermite quadrature, contrary to what the summary says).
But we were not able to compare with a function that does exact results for sufficient Monte Carlo sample size.
It is very unclear what “prediction” (in scare quotes) even means for
a random effects model. What should R generic function predict even
do for these model fits? Do methods even exist?
try(getS3method("predict", class(gout)[1]), silent = TRUE) |> is.function()## [1] TRUEtry(getS3method("predict", class(gout.too)[1]), silent = TRUE) |> is.function()## [1] FALSEtry(getS3method("predict", class(gout.too.too)[1]), silent = TRUE) |> is.function()## [1] FALSEtry(getS3method("predict", class(rout)[1]), silent = TRUE) |> is.function()## [1] FALSEOnly one has methods of predict for its results.
Prediction could mean
the unconditional expectation of the quantity being predicted,
the conditional expectation of the quantity being predicted given the observed data for that case,
a conditional mode if this is easier to do, or
many other ideas.
Here we will pick the last option to get per state predictions. R function
reaster output contains a component b that is the conditional modes of
the conditional distribution of the random effects given the observed data.
rout$b## stateAK stateAL stateAR stateAZ stateCA stateCO
## 0.063256937 -0.432795309 -0.651131912 -0.323091751 0.420667615 -0.181851750
## stateCT stateDC stateDE stateFL stateGA stateHI
## 0.106583691 0.279502729 0.126385522 0.211962776 -0.164241860 0.314222041
## stateIA stateID stateIL stateIN stateKS stateKY
## -0.226845271 -0.478529236 0.087427893 -0.079657141 -0.152219599 -0.323763255
## stateLA stateMA stateMD stateME stateMI stateMN
## -0.312895098 -0.216947551 0.551520963 0.389876483 0.136920711 -0.019575461
## stateMO stateMS stateMT stateNC stateND stateNE
## 0.124929007 0.076767385 -0.070860589 -0.470647614 -0.073633402 0.434353605
## stateNH stateNJ stateNM stateNV stateNY stateOH
## -0.293850221 0.098766325 0.043139713 -0.062597834 0.528552618 0.138405961
## stateOK stateOR statePA stateRI stateSC stateSD
## -0.628378814 -0.016273650 0.217890702 0.314222041 -0.432795309 0.126385522
## stateTN stateTX stateUT stateVA stateVT stateWA
## -0.206029866 -0.038943198 0.243990717 -0.045923802 0.200296849 0.441406524
## stateWI stateWV stateWY
## 0.183033773 0.047312733 -0.004301847adding the fixed effect to that gives
with(rout, alpha + b)## stateAK stateAL stateAR stateAZ stateCA stateCO
## 0.093608709 -0.402443537 -0.620780140 -0.292739979 0.451019387 -0.151499978
## stateCT stateDC stateDE stateFL stateGA stateHI
## 0.136935462 0.309854501 0.156737293 0.242314548 -0.133890088 0.344573813
## stateIA stateID stateIL stateIN stateKS stateKY
## -0.196493499 -0.448177464 0.117779665 -0.049305369 -0.121867827 -0.293411483
## stateLA stateMA stateMD stateME stateMI stateMN
## -0.282543326 -0.186595779 0.581872735 0.420228255 0.167272483 0.010776311
## stateMO stateMS stateMT stateNC stateND stateNE
## 0.155280779 0.107119157 -0.040508817 -0.440295842 -0.043281630 0.464705377
## stateNH stateNJ stateNM stateNV stateNY stateOH
## -0.263498449 0.129118097 0.073491485 -0.032246063 0.558904390 0.168757733
## stateOK stateOR statePA stateRI stateSC stateSD
## -0.598027042 0.014078122 0.248242474 0.344573813 -0.402443537 0.156737293
## stateTN stateTX stateUT stateVA stateVT stateWA
## -0.175678095 -0.008591426 0.274342489 -0.015572031 0.230648621 0.471758296
## stateWI stateWV stateWY
## 0.213385545 0.077664505 0.026049925and then mapping to the probability scale gives
e <- with(rout, alpha + b) |> invlogit() |> print()## stateAK stateAL stateAR stateAZ stateCA stateCO stateCT stateDC
## 0.5233851 0.4007254 0.3496040 0.4273332 0.6108816 0.4621973 0.5341805 0.5768497
## stateDE stateFL stateGA stateHI stateIA stateID stateIL stateIN
## 0.5391043 0.5602840 0.4665774 0.5853011 0.4510341 0.3897942 0.5294109 0.4876762
## stateKS stateKY stateLA stateMA stateMD stateME stateMI stateMN
## 0.4695707 0.4271689 0.4298304 0.4534859 0.6414982 0.6035379 0.5417209 0.5026941
## stateMO stateMS stateMT stateNC stateND stateNE stateNH stateNJ
## 0.5387424 0.5267542 0.4898742 0.3916705 0.4891813 0.6141298 0.4345039 0.5322348
## stateNM stateNV stateNY stateOH stateOK stateOR statePA stateRI
## 0.5183646 0.4919392 0.6361990 0.5420896 0.3547952 0.5035195 0.5617439 0.5853011
## stateSC stateSD stateTN stateTX stateUT stateVA stateVT stateWA
## 0.4007254 0.5391043 0.4561931 0.4978522 0.5681587 0.4961071 0.5574079 0.6157998
## stateWI stateWV stateWY
## 0.5531449 0.5194064 0.5065121And we can get standard errors
vout <- vcov(rout, re.too = TRUE, standard.deviation = FALSE)Make a function to post-process these estimates.
foo.factory <- function(v) {
is.alpha <- attr(v, "is.alpha")
is.bee <- attr(v, "is.b")
is.nu <- attr(v, "is.nu")
function(theta) {
stopifnot(is.numeric(theta))
stopifnot(is.finite(theta))
stopifnot(length(theta) == length(is.alpha))
alpha <- theta[is.alpha]
bee <- theta[is.bee]
invlogit(alpha + bee)
}
}
foo <- foo.factory(vout)Try it.
alphabeenu <- with(rout, c(alpha, b, nu))
all.equal(e, foo(alphabeenu))## [1] TRUEThe delta method is easily applied in this problem
using R function jacobian in CRAN package numDeriv to do the calculus.
jack <- jacobian(foo, alphabeenu)
vcov.e <- jack %*% vout %*% t(jack)Now assemble the results.
data.frame(Estimate = e, `Std. Error` = sqrt(diag(vcov.e))) |> round(3)## Estimate Std..Error
## stateAK 0.523 0.017
## stateAL 0.401 0.023
## stateAR 0.350 0.040
## stateAZ 0.427 0.016
## stateCA 0.611 0.005
## stateCO 0.462 0.011
## stateCT 0.534 0.011
## stateDC 0.577 0.030
## stateDE 0.539 0.019
## stateFL 0.560 0.005
## stateGA 0.467 0.008
## stateHI 0.585 0.029
## stateIA 0.451 0.016
## stateID 0.390 0.038
## stateIL 0.529 0.005
## stateIN 0.488 0.008
## stateKS 0.470 0.015
## stateKY 0.427 0.019
## stateLA 0.430 0.017
## stateMA 0.453 0.011
## stateMD 0.641 0.023
## stateME 0.604 0.031
## stateMI 0.542 0.006
## stateMN 0.503 0.007
## stateMO 0.539 0.008
## stateMS 0.527 0.012
## stateMT 0.490 0.017
## stateNC 0.392 0.014
## stateND 0.489 0.018
## stateNE 0.614 0.032
## stateNH 0.435 0.026
## stateNJ 0.532 0.007
## stateNM 0.518 0.013
## stateNV 0.492 0.013
## stateNY 0.636 0.010
## stateOH 0.542 0.005
## stateOK 0.355 0.035
## stateOR 0.504 0.009
## statePA 0.562 0.006
## stateRI 0.585 0.029
## stateSC 0.401 0.023
## stateSD 0.539 0.019
## stateTN 0.456 0.011
## stateTX 0.498 0.003
## stateUT 0.568 0.020
## stateVA 0.496 0.006
## stateVT 0.557 0.024
## stateWA 0.616 0.017
## stateWI 0.553 0.010
## stateWV 0.519 0.014
## stateWY 0.507 0.017For more about applying the delta method, with more examples, see the
vignette about the delta method in CRAN package aster.
As an example with crossed random effects we look at the salamander data
taken from the authoritative book on generalized linear models
McCullagh and Nelder (1983) (and ultimately coming from a biology paper cited therein)
and put in R package glmm
data(salamander)
class(salamander)## [1] "data.frame"sapply(salamander, class)## Mate Cross Female Male
## "numeric" "factor" "factor" "factor"unique(salamander$Mate)## [1] 1 0unique(salamander$Cross)## [1] R/R R/W W/R W/W
## Levels: R/R R/W W/R W/WThese data are about a salamander mating experiment. The response variable
Mate is Bernoulli (zero-or-one-valued). It is whether two salamanders
placed together by the experimenters mated. The salamanders came from two
locations called “Rough Butt” and “White Side” (names of locations, not
descriptions of salamanders). The Cross variable says which was mated with
which, female first, that is, type "R/W" indicates Rough Butt female with
White Side male. The variables Female and Male indicate the individual
salamanders in the experiment. Each salamander occurs in more than one line
of the data (was involved in more than one mating trial).
The model we fit has Cross as fixed effects and Female and Male
as crossed random effects.
lout <- glmer(Mate ~ 0 + Cross + (1 | Female) + (1 | Male),
data = salamander, family = binomial)
summary(lout)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Mate ~ 0 + Cross + (1 | Female) + (1 | Male)
## Data: salamander
##
## AIC BIC logLik -2*log(L) df.resid
## 430.6 453.9 -209.3 418.6 354
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0508 -0.6156 0.2714 0.5973 2.5506
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female (Intercept) 1.174 1.084
## Male (Intercept) 1.041 1.020
## Number of obs: 360, groups: Female, 60; Male, 60
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## CrossR/R 1.0082 0.3937 2.561 0.0105 *
## CrossR/W 0.3062 0.3747 0.817 0.4138
## CrossW/R -1.8959 0.4460 -4.251 2.13e-05 ***
## CrossW/W 0.9904 0.3912 2.532 0.0113 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## CrsR/R CrsR/W CrsW/R
## CrossR/W 0.280
## CrossW/R 0.112 -0.019
## CrossW/W 0.042 0.249 0.145time <- system.time(
gout <- glmm(Mate ~ 0 + Cross, list(~ 0 + Female, ~ 0 + Male),
data = salamander, varcomps.names = c("Female", "Male"),
family = binomial.glmm, m = 1e6,
cluster = makePSOCKcluster(detectCores()))
)
time## user system elapsed
## 38.352 12.339 1665.858summary(gout)##
## Call:
## glmm(fixed = Mate ~ 0 + Cross, random = list(~0 + Female, ~0 +
## Male), varcomps.names = c("Female", "Male"), data = salamander,
## family.glmm = binomial.glmm, m = 1e+06, cluster = makePSOCKcluster(detectCores()))
##
##
## Link is: "logit (log odds)"
##
## Fixed Effects:
## Estimate Std. Error z value Pr(>|z|)
## CrossR/R 1.03 0.42 2.435 0.0149 *
## CrossR/W 0.33 0.40 0.836 0.4030
## CrossW/R -1.91 0.45 -4.274 1.92e-05 ***
## CrossW/W 1.02 0.44 2.308 0.0210 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Variance Components for Random Effects (P-values are one-tailed):
## Estimate Std. Error z value Pr(>|z|)/2
## Female 1.40 0.59 2.371 0.00886 **
## Male 1.20 0.49 2.445 0.00724 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1So this took 27 minutes and 45.9 seconds of clock time on my computer.
dat <- transform(salamander, root = 1, varb = "only",
id = 1:nrow(salamander)) |> transform(varb = as.factor(varb))
head(dat)## Mate Cross Female Male root varb id
## 1 1 R/R 10 10 1 only 1
## 2 1 R/R 11 14 1 only 2
## 3 1 R/R 12 11 1 only 3
## 4 1 R/R 13 13 1 only 4
## 5 1 R/R 14 12 1 only 5
## 6 1 R/W 15 28 1 only 6rout <- reaster(fixed = Mate ~ 0 + Cross,
random = list(Female = ~ 0 + Female, Male = ~ 0 + Male),
pred = 0, fam = 1, varb, id, root, data = dat)
summary(rout, standard.deviation = FALSE)##
## Call:
## reaster.formula(fixed = Mate ~ 0 + Cross, random = list(Female = ~0 +
## Female, Male = ~0 + Male), pred = 0, fam = 1, varvar = varb,
## idvar = id, root = root, data = dat)
##
##
## Fixed Effects:
## Estimate Std. Error z value Pr(>|z|)
## CrossR/R 0.7782 0.3144 2.476 0.0133 *
## CrossR/W 0.2448 0.3033 0.807 0.4195
## CrossW/R -1.4822 0.3494 -4.242 2.22e-05 ***
## CrossW/W 0.7690 0.3135 2.453 0.0142 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Variance Components (P-values are one-tailed):
## Estimate Std. Error z value Pr(>|z|)/2
## Female 0.6544 0.3072 2.131 0.0166 *
## Male 0.5764 0.2874 2.005 0.0225 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
6.3.3 Comments
There does not seem to be any programmatic access to the estimate of the variance component and its standard error (which the summary function does print but does not put in the returned object). This is bad design.
It says it uses Laplace approximation rather than Gauss-Hermite quadrature, even though the function help says the opposite. What does it do? Are these accurate estimates or not?