by Resampling under the null hypothesis instead of case resampling will preserve the correlations and distributional characteristics of the data and allow you to perform hypothesis testing. For more information, please refer to Westfall and Young (1993).
Resampling under the null hypothesis
Selecting a null model for resampling
The very first step of this approach is to base the bootstrap on the null distribution, which is the population distribution that represents the null hypothesis.
In the context of ANOVA or ANCOVA designs, the following null distributions may be applicable:
- One-way ANOVA: intercept-only model
- One-way ANCOVA: covariates-only model
- In two-way factorial designs, we are testing three null hypotheses simultaneously based on two main effects and the interaction effect.
- ANOVA: intercept-only model
- ANCOVA: covariates-only model
Steps
- Fit the full model.
- Obtain the observed F-values from the full model.
- Fit the null model (i.e., intercept-only model for ANOVA or covariate-only model for ANCOVA).
- Optionally centre and rescale the residuals of the null model based on the leverage values of the observations.
- Our new dependent variable for the bootstrapping is the sum of the fitted values of the null model and the resampled residuals from the null model, but we use the predictors from the full model.
- Obtain the resampled F-values.
- Repeat steps 5-6, B times (e.g., 10,000).
- The bootstrapped p-value is the proportion of resampled F-values that are greater then the observed F-values.
Function
Required packages
carto fit the ANOVA and get the F-valuesbootto run the bootstrapping and calculate the BCa confidence intervals
Options
null.model: Null model of class lm full.model: Full model of class lm B: Number of bootstraps (default = 1000) scaled: Rescale residuals (default = TRUE) seed: Seed for replication (default = 1234) ci: Calculate confidence intervals (default = TRUE) ci.type: Confidence interval type (options = c(“bca”, “perc”)). cent: Percentile for confidence intervals (default = .95) dec: Number of decimal places (default = 5)
Code
| # Load required packages -------------------------------------------------- | |
| if (!require(pacman)) { | |
| install.packages("pacman") | |
| } | |
| pacman::p_load(car, boot) | |
| # Function ---------------------------------------------------------------- | |
| nullboot.Anova <- function(null.model, full.model, | |
| B = 1000, scaled = TRUE, seed = 1234, | |
| ci = TRUE, cent = .95, dec = 5, | |
| ci.type = c("bca", "perc")) { | |
| # Returns error if models are not 'lm' class | |
| stopifnot(class(null.model) == "lm" | class(full.model) == "lm") | |
| # Set seed | |
| set.seed(seed) | |
| # Extract data from the two models | |
| data_full <- model.frame(full.model) | |
| data_null <- model.frame(null.model) | |
| # Get residuals from null model and recenter and rescale based on | |
| # Davidson and MacKinnon (2004) if requested. | |
| if (isTRUE(scaled)) { | |
| # Number of observations | |
| N <- nrow(data_null) | |
| # Number of predictors (from full model) | |
| K <- length(labels(terms(full.model))) | |
| # Extract and scale residuals (if requested) from null model | |
| Er <- residuals(null.model) / sqrt(1 - hatvalues(null.model)) | |
| } else { | |
| Er <- residuals(null.model) | |
| } | |
| # Predicted values from null model | |
| Yhat <- fitted(null.model) | |
| # Get ANOVA statistics for full model | |
| mod_ANOVA <- Anova(full.model, type = 3) | |
| # Get observed F-values (excluding Intercept and NA [Residuals]) | |
| Fobs <- as.vector(na.omit(mod_ANOVA[-1, 3])) | |
| # Updated formula for bootstapping with Ystar as DV | |
| formula_resample <- update.formula(formula(full.model), Ystar ~ .) | |
| # Bootstrap function | |
| bs.Anova <- function(data, i) { | |
| data_full$Ystar <- Yhat + Er[i] | |
| # Fit linear model | |
| model_resample <- lm(formula_resample, data = data_full) | |
| # Get F-values from ANOVA | |
| Fstar <- as.vector(na.omit(Anova(model_resample, type = 3)[-1, 3])) | |
| return(Fstar) | |
| } | |
| # Run bootstrapping | |
| bootAnova <- boot( | |
| data_full, statistic = bs.Anova, R = B, | |
| parallel = "snow" | |
| ) | |
| # Bootstrapped F-values | |
| Fstar <- bootAnova$t | |
| # The proportion of resampled F-values greater than or equal to | |
| # the observed F-values | |
| pBoot <- vector() | |
| for (i in 1:length(Fobs)) { | |
| pBoot[i] <- (sum(Fstar[, i] >= Fobs[i]) + 1) / (length(Fstar[, i]) + 1) | |
| } | |
| # Create data.frame for outputting | |
| varNames <- labels(terms(full.model)) | |
| pRaw <- as.vector(na.omit(mod_ANOVA[-1, 4])) | |
| # Calculate confidence intervals if requested | |
| if (isTRUE(ci)) { | |
| ci.type <- match.arg(ci.type) | |
| FstarCI <- list() | |
| ciLB <- vector() | |
| ciUB <- vector() | |
| bcaLB <- vector() | |
| bcaUB <- vector() | |
| for (i in 1:length(Fobs)) { | |
| suppressWarnings(FstarCI[[i]] <- boot.ci( | |
| bootAnova, | |
| type = ci.type, | |
| index = i, | |
| conf = cent | |
| )) | |
| # Percentile confidence intervals | |
| if (ci.type == "perc") { | |
| ci_desc <- "Percentile" | |
| ciLB[i] <- FstarCI[[i]]$perc[, 4] | |
| ciUB[i] <- FstarCI[[i]]$perc[, 5] | |
| } | |
| # BCa confidence intervals | |
| if (ci.type == "bca") { | |
| ci_desc <- "BCa" | |
| ciLB[i] <- FstarCI[[i]]$bca[, 4] | |
| ciUB[i] <- FstarCI[[i]]$bca[, 5] | |
| } | |
| } | |
| bootOut <- data.frame( | |
| "F.value" = Fobs, | |
| "p.value" = pRaw, | |
| "p.boot" = pBoot, | |
| "CI.LB" = ciLB, | |
| "CI.UB" = ciUB, | |
| row.names = varNames | |
| ) | |
| } else { | |
| ci_desc <- "Not requested" | |
| bootOut <- data.frame( | |
| "F.value" = Fobs, | |
| "p.value" = pRaw, | |
| "p.boot" = pBoot, | |
| row.names = varNames | |
| ) | |
| } | |
| cat( | |
| "\n", "Bootstrapped ANOVA with Type III tests", | |
| "\n", | |
| "\n", "Resampling type: ", | |
| ifelse(isTRUE(scaled), | |
| "rescaled residuals", | |
| "residuals" | |
| ), | |
| "\n", "Number of bootstrap resamples: ", B, | |
| "\n", "Bootstrapped confidence interval type: ", ci_desc, | |
| "\n", "Confidence interval: ", 100 * cent, "%", | |
| "\n", | |
| "\n", sep = "" | |
| ) | |
| # Turn off scientific | |
| options(scipen = 999) | |
| print(round(bootOut, dec)) | |
| # Turn on scientific | |
| options(scipen = 0) | |
| } |
Example
To prevent issues with missing data I recommend fitting the full model first. Alternatively, you could just omit cases with missing values from your dataset.
full_model <- lm(dv ~ sex*group, data = df)
Then, fit the null model using the update() function, which is the intercept-only model for this example and using the model.frame() of the full_model.
null_model <- update(full_model, . ~ + 1, data = model.frame(full_model))
Finally, run the function. Note that I have not requested confidence intervals.
nullboot.Anova(null.model = null_model, full.model = full_model, ci = FALSE)
Output
The function will also output the BCa or percentile bootstrapped confidence intervals if requested.
Resampling type: rescaled residuals Number of bootstrap resamples: 1000 Bootstrapped confidence interval type: Not requested Confidence interval: 95% F.value p.value p.boot sex 0.42978 0.51425 0.51449 group 0.94369 0.33468 0.32867 group:sex 0.00101 0.97476 0.97203
Testing the null hypothesis
As we are resampling the residuals from the null model, the output below is the as above but with the addition of confidence intervals.
nullboot.Anova(null.model = null_model, full.model = full_model, ci = TRUE)
Output
Bootstrapped ANOVA with Type III tests
Resampling type: rescaled residuals
Number of bootstrap resamples: 1000
Bootstrapped confidence interval type: BCa
Confidence interval: 95%
F.value p.value p.boot CI.LB CI.UB
sex 0.42978 0.51425 0.51449 0.00073 5.23939
group 0.94369 0.33468 0.32867 0.03607 9.48676
group:sex 0.00101 0.97476 0.97203 0.00001 0.00111
Model-based resampling
If you are not interested in hypothesis testing, but would like to obtain confidence intervals for the F-statistic, then you should resample the residuals of the full model. This is also known as fixed-x resampling or the parametric bootstrap.
In which case, you would run the function with full_model for both the null.model and full.model options:
nullboot.Anova(null.model = full_model, full.model = full_model, ci = TRUE)
Parametric bootstrapping
Below is the parametric bootstrap where we’re resampling the residuals of the full model.
Bootstrapped ANOVA with Type III tests
Resampling type: rescaled residuals
Number of bootstrap resamples: 1000
Bootstrapped confidence interval type: BCa
Confidence interval: 95%
F.value p.value p.boot CI.LB CI.UB
sex 0.42978 0.51425 0.59540 0.00003 4.28447
group 0.94369 0.33468 0.50649 0.00095 8.67616
sex:group 0.00101 0.97476 0.97602 0.00000 0.00090
Notice that the p-values and confidence intervals differ depending on whether you’re resampling under the null hypothesis or are using parametric bootstrapping.
Recommendations
- If you are interesting in hypothesis testing (i.e., obtaining a bootstrapped p-value), then you should resample under the null hypothesis.
- If you are interested in obtaining confidence intervals for the F-values, then you should use parametric bootstrapping.
Reference
Westfall, P. H., & Young, S. S. (1993). Resampling-based multiple testing: Examples and methods for p-value adjustment (Vol. 279). John Wiley & Sons.
