Installation

The package can be installed from GitHub using devtools.

# devtools::install_github("koekvall/lmmstest")
library(lmmstest)

Setting

Suppose a vector of responses \(Y\in \mathbb{R}^n\), matrix of predictors \(X\in \mathbb{R}^{n\times p}\), and design matrix \(Z \in \mathbb{R}^{n\times q}\) satisfy, for some \(\beta \in \mathbb{R}^p\), \[ Y = X\beta + ZU + E, \] where \(U = [U_1, \dots, U_q]^T\) is a multivariate normal vector of random effects with mean zero and the error vector \(E = [E_1, \dots, E_n]^T\) is multivariate normal with mean zero and covariance matrix \(\sigma^2 I_n\), independent of \(U\). Assume further that the elements of \(U\) are independent and that, for every \(j \in \{1, \dots, q\}\), there is a \(k \in \{1, \dots, d\}\) such that \(\mathrm{var}(U_j) = \lambda_k^2\). That is, the elements of \(\lambda = [\lambda_1, \dots, \lambda_d]^T \in [0, \infty)^d\) are the standard deviations of the random effects. Several random effects can have the same variance, in which case \(d < q\).

Background

Many popular tests perform poorly when one or more of the scale parameters \(\sigma, \lambda_1, \dots, \lambda_d\) are near or equal to zero. One way to address this problem is to notice that the Fisher information is singular when a scale parameter is zero. In particular, the element of the score function corresponding to such a scale parameter is equal to zero almost surely. The modified score test addresses this by replacing such elements, which are first order partial derivatives of the log-likelihood, by second order partial derivatives. When no scale parameter is equal to zero the modified score test statistic is the same as the usual score test statistic standardized by expected Fisher information; when one or more scale parameters are equal to zero the modified statistic can be defined as a limit of the usual one. In the linear mixed model, this limit has an analytical expression which facilitates our implementation.

Example

Suppose that for observational unit \(i = 1, \dots, N\) and time \(t = 1, \dots, N_T\), \[ Y_{i, t} = \beta_1 + \beta_2 X_{i, t} + U_{1, i} + U_{2, i} X_{i, t} + E_{i, t}. \] This model can be written in stacked form by letting \(Y = [Y_{1, 1}, \dots, Y_{N, N_T}]^T\) and defining \(X\), \(Z\), \(U\), and \(E\) accordingly. In particular, \(Z \in \mathbb{R}^{n \times 2N}\) with \(n = NN_T\). The random effects \(U_{1, i} \sim \mathcal{N}(0, \lambda_1)\) and \(U_{2, i} \sim \mathcal{N}(0, \lambda_2)\) are, respectively, a random intercept and a random slope, shared by observations from the same unit. In stacked notation, \[ U = \begin{bmatrix} U_{1, 1} \\ U_{1, 2} \\ \vdots \\ U_{2, N}\end{bmatrix} \sim \mathcal{N}\left(0, \mathrm{diag}(\lambda_1^2, \lambda_2^2) \otimes I_N\right), \] where \(\otimes\) is the Kronecker product.

Consider testing the null hypothesis \((\sigma, \lambda_1, \lambda_2) = (1, 0.5, 0)\). When \(\lambda_2 = 0\), the Fisher information is singular and common test statistics are unreliable. The following code generates one realization from the model with true parameters consistent with the null hypothesis and tests the null hypothesis using the proposed procedure. The predictors are drawn independently from a uniform distribution on \([-1, 2]\).

# For replicability
set.seed(2)

# Sample size
N <- 30
N_T <- 10
n <- N * N_T

# True parameters
sigma <- 1
lambda <- c(0.5, 0)
Beta <- c(1, 1)

# Generate matrix of data where rows index units and columns index time
XX <- matrix(runif(N * N_T, -1, 2) , nrow = N, ncol = N_T)
U1 <- rnorm(N, sd = lambda[1])
U2 <- rnorm(N, sd = lambda[2])
E <- matrix(rnorm(N * N_T, sd = sigma), nrow = N, ncol = N_T)
Y <-  Beta[1] + Beta[2] * XX + E
for(ii in 1:N){
  Y[ii, ] <- Y[ii, ] + U1[ii] + U2[ii] * XX[ii, ]
}

# Put data in tall format and add unit and time factors
ex_data <- data.frame(c(t(Y)), c(t(XX)))
ex_data$time <- rep(1:N_T, N)
ex_data$unit <- as.factor(rep(1:N, each = N_T))
names(ex_data)[1:2] <- c("y", "x")

# Inspect data
head(ex_data, 3)
##            y          x time unit
## 1 0.01218843 -0.4453532    1    1
## 2 1.54543509 -0.9687564    2    1
## 3 2.71624885  1.3100836    3    1
tail(ex_data, 3)
##             y         x time unit
## 298 0.9009858 0.4221224    8   30
## 299 1.3055184 0.6116278    9   30
## 300 0.1871057 0.5122688   10   30
# Create model matrices
X <- model.matrix(~x, data = ex_data)
Z <- model.matrix(~0 + unit + unit:x, data = ex_data)

# The first N columns of Z correspond to random intercepts and
# the next 30 to random slopes.
dim(Z)
## [1] 300  60
lam_idx <- rep(1:2, each = N)
L <- diag(lambda[lam_idx], 2 * N)

# Get estimate of coefficient vector under the null hypothesis.
# Since Sigma is known under the null hypothesis, this estimate is the
# generalized least squares estimate.
Sigma_null <- diag(sigma^2, n) + Z %*% L^2 %*% t(Z)
Beta_null <- c(qr.solve(crossprod(X, qr.solve(Sigma_null, X)),
                      crossprod(X, qr.solve(Sigma_null, ex_data$y))))
Beta_null
## [1] 1.059725 1.057824
# Compute modified score test statistic for elements 3:5 of theta =
# c(Beta, sigma, lambda[1], lambda_2), i.e. for (sigma, lambda[1], lambda_2).
# The supplied arguments for  (sigma, lambda[1], lambda_2) are the true values,
# so the null hypothesis is true.
lmmstest::score_test(y = ex_data$y,
                     X = X,
                     Z = Z,
                     Beta = Beta_null,
                     sigma = sigma,
                     lambda = lambda,
                     lam_idx = lam_idx,
                     test_idx = 3:5)
##   chi_sq       df    p_val 
## 1.075005 3.000000 0.783111

Comparison to common test statistics

The following figure shows quantile-quantile plots for the modified score (circles), likelihood ratio (triangles), and Wald (plus signs) test statistics from a simulation. The theoretical quantiles are from a chi-square distribution with 3 degrees of freedom. Data were generated according to the model from the previous section assuming the null hypothesis \((\sigma, \lambda_1, \lambda_2) = (1, 0, 0)\) (first plot), \((1, 0.01, 0.01)\) (second plot), or \((1, 1, 1)\) (third plot).

When no scale parameter is zero, classical theory says all three test statistics are asymptotically chi-square distributed with 3 degrees of freedom. The third plot indicates this asymptotic distribution is a good approximation when scale parameters are substantially different from zero.

The second plot shows the asymptotic distribution is a good approximation for the proposed procedure but not for the other two when scale parameters are near zero.

When some scale parameters are zero, the asymptotic distributions of the likelihood ratio and Wald test statistics are no longer chi-square distributions, and accordingly the first plot shows poor agreement between the quantiles. By contrast, the modified score test has an asymptotic chi-square distribution also in this setting, and the first plot shows the asymptotic distribution provides a good approximation.