simulate_wald_pvalues: Generate Wald test P-Values with Monte Carlo Simulations

Description

This function performs Monte Carlo simulations to generate p-values for model coefficients by refitting the model with new simulated responses and computing the Wald test statistics for each simulation. It's standard behavior verify if the type I error from Wald tests are under control, considering the provided model as "true", i.e., the model assumptions are valid. It supports univariate and joint Wald tests, using chi-squared distributions to calculate p-values.

Usage

simulate_wald_pvalues(
  model,
  nsim = 1000,
  responses = NULL,
  test_coefficients = NULL,
  refit_fn = NULL,
  coef_fn = get_fixef,
  vcov_fn = get_vcov,
  show_progress = TRUE,
  plot.it = TRUE,
  ...
)

Value

An object of class AD_pvalues, which contains the following components:

test_coefficients: Vector of coefficients being tested.
pvalues_matrix: Matrix of p-values where each column corresponds to a simulation and each row corresponds to a coefficient.
pvalues_joint: Vector containing the joint p-values obtained from each simulation.
simulation_fixef: List of fixed effect coefficient estimates from each simulation.
simulation_vcov: List of covariance matrices estimated from each simulation.
simulation_warning: Vector of boolean indicating if a simulation threw a warning.
converged: Logical vector indicating whether model refitting converged for each simulation.
responses: Simulated responses used for refitting the model.

Arguments

model: A fitted model object that will be used to simulate responses.
nsim: The number of simulations to perform.
responses: An optional list of values to be used as response variables to refit the model.
test_coefficients: Numeric vector. A vector with values to be used to compute the test statistic. It should be the coefficients that was used to compute the fitted values of the response. If NULL defaults to coef_fn(model).
refit_fn: Function to refit the model with new responses. If NULL, defaults to get_refit(model, y, ...).
coef_fn: Function that retrieves the coefficients of the model.
vcov_fn: Function that computes the variance-covariance matrix for the models adjusted in the simulations.
show_progress: Display a progress bar for the simulation iteration.
plot.it: Logical. Generate ecdf plot for joint Wald test.
...: Additional arguments to be passed to refit_fn.

Details

If responses is provided, the function refits the model with these new response vectors. Otherwise, it generates new responses with stats::simulate().

For each new response calls get_refit() to generate a new model with the new response. It gets the fixed effects and the variance and covariance matrix with get_fixef() and get_vcov().

Each simulated model is refitted using the specified refit_fn (or the default refit function) and the fixed effects coefficients and variance-covariance matrix are extracted using coef_fn and vcov_fn, respectively. The univariate Wald test is computed from the Wald statistic for each coefficient, while the joint Wald test uses the inverse variance-covariance matrix to compute a Wald statistic for the test_coefficients. P-values are calculated from a chi-squared distribution with appropriate degrees of freedom.

Examples

Run this code

# from help("glm")
counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)
model <- glm(counts ~ outcome + treatment, family = poisson())
new_responses <- replicate(100, MASS::rnegbin(fitted.values(model), theta = 4.5), simplify = FALSE)
simulate_wald_pvalues(model, responses = new_responses, nsim = 100)

## Using custom refit_fn
if (require("survival")) {
  fit <- survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian,
    dist = "exponential"
  )
  fitted_rate <- 1 / fitted(fit)
  new_responses <- replicate(100, rexp(length(fitted_rate), fitted_rate), simplify = FALSE)
  refit_surv_ovarian <- function(.y) {
    survreg(Surv(.y, fustat) ~ ecog.ps + rx, ovarian, dist = "exponential")
  }
  simulate_wald_pvalues(fit, responses = new_responses, refit_fn = refit_surv_ovarian)
}