mixedbiastest: Bias Diagnostic for Linear Mixed Models

In mixedbiastest: Bias Diagnostic for Linear Mixed Models

Bias Diagnostic for Linear Mixed Models

Description

Performs a permutation test to assess the bias of fixed effects in a linear mixed model fitted with lmer. This function computes the test statistic and performs the permutation test, returning an object of class "mixedbiastest".

Usage

mixedbiastest(model, n_permutations = 10000, k_list = NULL, verbose = FALSE)

Arguments

model	An object of class lmerMod (or a subclass, such as lmerModLmerTest) fitted using lmer from the lme4 package. Models fitted with glmer or with non-identity residual covariance structures (for example, non-unit weights) are not currently supported.
n_permutations	Integer. Number of permutations to perform (default is 10000). Must be a positive integer.
k_list	Optional list of numeric vectors. Each vector specifies a linear combination of fixed effects to test. If NULL, each fixed effect is tested individually.
verbose	Logical. If TRUE, prints detailed messages during execution.

Details

The implementation follows Karl and Zimmerman (2021) and is currently restricted to:

• Gaussian linear mixed models fitted by lmer (no GLMMs).

• Diagonal random-effects covariance matrices (i.e., a block-diagonal G with scalar blocks for each random-effect coefficient).

• Homoskedastic residual errors with covariance \sigma^2 I_n (no observation weights or residual correlation structures).

In particular, models fitted with glmer or with non-identity residual covariance structures (for example, non-unit observation weights) are beyond the scope of the current implementation.

While the diagnostic of Karl and Zimmerman (2021) is formulated for general linear mixed models with arbitrary covariance matrices G and R, this function implements the special case of Gaussian lmer models with diagonal G and homoskedastic residual errors. Extending mixedbiastest() to correlated random effects or more general residual covariance structures would require substantial additional work on the underlying linear algebra and permutation scheme, and is left for future research.

See the list_fixed_effects function if you would like to construct contrasts of fixed effects to be used as k_list.

Value

An object of class "mixedbiastest" containing:

results_table

A data frame with the test results for each fixed effect or contrast, including bias estimates and permutation p-values.

permutation_values

A list of numeric vectors containing permutation values for each fixed effect or contrast.

model

The original lmerMod model object provided as input.

Acknowledgments

Development of this package was assisted by GPT o1-preview and GPT 5 Pro, which helped in constructing the structure of much of the code and the roxygen documentation. The code is based on the R code provided by Karl and Zimmerman (2020).

References

Karl, A. T., & Zimmerman, D. L. (2021). A diagnostic for bias in linear mixed model estimators induced by dependence between the random effects and the corresponding model matrix. Journal of Statistical Planning and Inference, 212, 70-80. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.jspi.2020.06.004")}

Karl, A., & Zimmerman, D. (2020). Data and Code Supplement for "A Diagnostic for Bias in Linear Mixed Model Estimators Induced by Dependence Between the Random Effects and the Corresponding Model Matrix". Mendeley Data, V1. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.17632/tmynggddfm.1")}

Examples

if (requireNamespace("plm", quietly = TRUE) && requireNamespace("lme4", quietly = TRUE)) {
library(lme4)
data("Gasoline", package = "plm")
mixed_model <- lmer(lgaspcar ~ lincomep + lrpmg + lcarpcap + (1 | country),
                    data = Gasoline)
result <- mixedbiastest(mixed_model)
print(result); plot(result)
}
if (requireNamespace("lme4", quietly = TRUE)) {
library(lme4)
example_model <- lmer(y ~ x + (1 | group), data = example_data)
result2 <- mixedbiastest(example_model)
print(result2); plot(result2)

# Simulate data
set.seed(123)
n_groups <- 30
n_obs_per_group <- 10
group <- rep(1:n_groups, each = n_obs_per_group)
x <- runif(n_groups * n_obs_per_group)
beta0 <- 2; beta1 <- 5
sigma_u <- 1; sigma_e <- 0.5
u <- rnorm(n_groups, 0, sigma_u)
e <- rnorm(n_groups * n_obs_per_group, 0, sigma_e)
y <- beta0 + beta1 * x + u[group] + e

data_sim <- data.frame(y = y, x = x, group = factor(group))
model3 <- lmer(y ~ x + (1 | group), data = data_sim)
result3 <- mixedbiastest(model3, verbose = TRUE)
plot(result3)
}

mixedbiastest documentation built on Dec. 1, 2025, 1:08 a.m.