| Title: | Bias Diagnostic for Linear Mixed Models |
|---|---|
| Description: | Provides a function to perform bias diagnostics on linear mixed models fitted with lmer() from the 'lme4' package. Implements permutation tests for assessing the bias of fixed effects, as described in Karl and Zimmerman (2021) <doi:10.1016/j.jspi.2020.06.004>. Karl and Zimmerman (2020) <doi:10.17632/tmynggddfm.1> provide R code for implementing the test using 'mvglmmRank' output. Development of this package was assisted by 'GPT o1-preview' for code structure and documentation. |
| Authors: | Andrew T. Karl [cre, aut]
|
| Maintainer: | Andrew T. Karl <[email protected]> |
| License: | GPL-3 |
| Version: | 0.3.0 |
| Built: | 2024-12-13 10:54:38 UTC |
| Source: | CRAN |
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).
mixedbiastestPerforms the bias diagnostic test.
print.mixedbiastestPrints the results of the bias diagnostic.
plot.mixedbiastestPlots the permutation distributions and observed test statistics for each fixed effect.
list_fixed_effectsList Fixed Effects from an lmerMod Object.
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).
Maintainer: Andrew T. Karl [email protected] (ORCID)
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
A dataset containing 97 observations of three variables: y, x, and group.
example_dataexample_data
A data frame with 97 rows and 3 variables:
Numeric response variable.
Numeric predictor variable.
Integer indicating group membership.
This function lists the fixed effects coefficients from an 'lmerMod' object, providing the index and name of each coefficient. This can help users construct contrast vectors ('k_list') for use with the 'mixedbiastest' function.
list_fixed_effects(model)list_fixed_effects(model)
model |
An object of class 'lmerMod' fitted using 'lmer' from the 'lme4' package. |
A data frame with two columns:
IndexThe index of each fixed effect coefficient.
CoefficientThe name of each fixed effect coefficient.
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).
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) # Suppose we want to test the contrast: lincomep - lcarpcap p <- nrow(fixed_effects) k <- rep(0, p) k[fixed_effects$Index[fixed_effects$Coefficient == "lincomep"]] <- 1 k[fixed_effects$Index[fixed_effects$Coefficient == "lcarpcap"]] <- -1 # Run the bias test with the custom contrast result <- mixedbiastest(mixed_model, k_list = list(k)) print(result) plot(result) } else { message("Please install 'plm' and 'lme4' packages to run this example.") }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) # Suppose we want to test the contrast: lincomep - lcarpcap p <- nrow(fixed_effects) k <- rep(0, p) k[fixed_effects$Index[fixed_effects$Coefficient == "lincomep"]] <- 1 k[fixed_effects$Index[fixed_effects$Coefficient == "lcarpcap"]] <- -1 # Run the bias test with the custom contrast result <- mixedbiastest(mixed_model, k_list = list(k)) print(result) plot(result) } else { message("Please install 'plm' and 'lme4' packages to run this example.") }
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"'.
mixedbiastest(model, n_permutations = 10000, k_list = NULL, verbose = FALSE)mixedbiastest(model, n_permutations = 10000, k_list = NULL, verbose = FALSE)
model |
An object of class 'lmerMod' fitted using 'lmer' from the 'lme4' package. |
n_permutations |
Integer. Number of permutations to perform (default is 10000). |
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. |
**Note:** This function currently supports only models with diagonal random effects covariance matrices (i.e., the G matrix is diagonal). The methodology for non-diagonal G matrices is described in Karl and Zimmerman (2021), but is not implemented in this version of the package.
See the list_fixed_effects function if you would like to estimate the bias of a contrast of fixed effects.
An object of class "mixedbiastest" containing:
results_tableA data frame with the test results for each fixed effect or contrast, including bias estimates and p-values.
permutation_valuesA list of numeric vectors containing permutation values for each fixed effect or contrast.
modelThe original lmerMod model object provided as input.
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).
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
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) 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) }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) 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) }
Plots the permutation distributions and observed test statistics for each fixed effect.
## S3 method for class 'mixedbiastest' plot(x, ...)## S3 method for class 'mixedbiastest' plot(x, ...)
x |
An object of class '"mixedbiastest"'. |
... |
Additional arguments (currently not used). |
A ggplot object showing permutation distributions for all fixed effects.
Prints the results of the bias diagnostic in a formatted table.
## S3 method for class 'mixedbiastest' print(x, ...)## S3 method for class 'mixedbiastest' print(x, ...)
x |
An object of class '"mixedbiastest"'. |
... |
Additional arguments (currently not used). |
No return value. This function is called for its side effects (printing the results).