Package 'mixedbiastest'

mixedbiastest: Bias Diagnostics for Linear Mixed Models

Description

The 'mixedbiastest' package provides a function to perform bias diagnostics on linear mixed models fitted with 'lmer' from the 'lme4' package. It implements permutation tests for assessing the bias of fixed effects, as described in Karl and Zimmerman (2021).

Details

The methods in this package are designed for Gaussian linear mixed models with diagonal random-effects covariance matrices and homoskedastic residual errors with covariance $\sigma^2 I_n$. Models with more general random-effects structures or residual covariance patterns are currently not supported.

While the bias diagnostic of Karl and Zimmerman (2021) is derived for general linear mixed models with arbitrary random-effects and residual covariance matrices $G$ and $R$, the current mixedbiastest implementation focuses on the practically important case described above. Extending the package to handle correlated random effects or more general residual covariance structures would require additional work on both the underlying linear algebra and the permutation scheme, and is therefore left for future research.

Functions

mixedbiastest

Performs the bias diagnostic test.

print.mixedbiastest

Prints the results of the bias diagnostic.

plot.mixedbiastest

Plots the permutation distributions and observed test statistics for each fixed effect.

list_fixed_effects

List Fixed Effects from an merMod Object.

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).

Author(s)

Maintainer: Andrew T. Karl [email protected] (ORCID)

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. 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. doi:10.17632/tmynggddfm.1

---

Example Dataset for the mixedbiastest Package

Description

A dataset containing 97 observations of three variables: y, x, and group.

Usage

example_data

Format

A data frame with 97 rows and 3 variables:

y

Numeric response variable.

x

Numeric predictor variable.

group

Integer indicating group membership.

---

List Fixed Effects from a merMod Object

Description

This function lists the fixed effects coefficients from an 'lmerMod'/'merMod' object, providing the index and name of each coefficient. This can help users construct contrast vectors ('k_list') for use with the 'mixedbiastest' function.

Usage

list_fixed_effects(model)

Arguments

model

An object that inherits from class merMod, such as an lmerMod or glmerMod object fitted using lmer or glmer from the lme4 package.

Value

A data frame with two columns:

Index

The index of each fixed effect coefficient.

Coefficient

The name of each fixed effect coefficient.

Acknowledgments

Development of this package was assisted by GPT o1-preview, 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).

Examples

if (requireNamespace("plm", quietly = TRUE) && requireNamespace("lme4", quietly = TRUE)) {
  library(lme4)
  data("Gasoline", package = "plm")
  # Fit a random effects model using lme4
  mixed_model <- lmer(lgaspcar ~ lincomep + lrpmg + lcarpcap + (1 | country),
                      data = Gasoline, REML = FALSE)

  # List fixed effects
  fixed_effects <- list_fixed_effects(mixed_model)
  print(fixed_effects)
}

---

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. 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. 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)
}

---

Plot Method for Bias Diagnostic Results

Description

Plots the permutation distributions and observed test statistics for each fixed effect.

Usage

## S3 method for class 'mixedbiastest'
plot(x, bins = 30, ...)

Arguments

x	An object of class "mixedbiastest".
bins	Integer, number of bins for the histograms (default 30).
...	Additional arguments (currently not used).

Value

A ggplot object (returned invisibly) showing permutation distributions for all fixed effects.

---

Print Method for Bias Diagnostic Results

Description

Prints the results of the bias diagnostic in a formatted table.

Usage

## S3 method for class 'mixedbiastest'
print(x, ...)

Arguments

x	An object of class '"mixedbiastest"'.
...	Additional arguments (currently not used).

Value

The input object, returned invisibly.

---