R语言boot()输出解读：bootstrap重复测量ANOVA分析疑问

Hey there! Let's work through your bootstrap ANOVA problem step by step—you're absolutely right to turn to bootstrapping for non-normal reaction time data, but we need to adjust your approach to fit your within-subjects design, and clarify how to interpret the results.

1. First: Fix Your Bootstrap Code for Within-Subjects Data #

Your current code randomly samples rows from your dataframe, which breaks the structure of your within-subjects design. Each participant has multiple observations (one for every combination of test situation and cognitive load), so we need to resample participants, not individual rows, to preserve the repeated-measures structure.

First, make sure your dataframe has a unique SubjectID column (to identify each participant). Then use this revised code:

# Load required packages
library(boot)

# Bootstrap function tailored for within-subjects ANOVA
boot_within_anova <- function(data, indices) {
  # Resample participants with replacement
  sampled_subjects <- unique(data$SubjectID)[indices]
  sampled_data <- data[data$SubjectID %in% sampled_subjects, ]
  
  # Fit a within-subjects ANOVA (Error term accounts for participant variability)
  anova_fit <- aov(RT ~ Situation * Code + Error(SubjectID/(Situation * Code)), data = sampled_data)
  
  # Extract F-statistics for each effect
  f_values <- summary(anova_fit)[[1]][["F value"]]
  
  # Return F-stats for: Situation main effect, Cognitive Load main effect, Interaction
  return(c(Situation = f_values[1], Cognitive_Load = f_values[2], Interaction = f_values[3]))
}

# Run the bootstrap (resample participants, not rows)
boot_results <- boot(
  data = dataframe,
  statistic = boot_within_anova,
  R = 2000,
  strata = dataframe$SubjectID # Ensures we sample at the participant level
)

Why this works: #

• The Error(SubjectID/(Situation * Code)) term in aov properly models the variance from repeated measures on the same participant.
• Resampling participants (instead of rows) keeps all of a participant's observations together, which is critical for within-subjects designs.

2. Interpreting the Bootstrap Output #

Your original output was from a linear regression (lm), which is why you saw regression coefficients instead of ANOVA effect statistics. The revised output will look like this (example):

ORDINARY NONPARAMETRIC BOOTSTRAP

Call:
boot(data = dataframe, statistic = boot_within_anova, R = 2000, strata = dataframe$SubjectID)

Bootstrap Statistics :
          original        bias    std. error
t1*       4.231201  -0.1245678   1.890123
t2*      12.567890   0.0891234   3.123456
t3*       2.109876   0.0345678   1.234567

Here's what each column means:

• original: The F-statistic for each effect from your original dataset. t1* = Situation main effect, t2* = Cognitive Load main effect, t3* = Interaction effect.
• bias: The average difference between the bootstrap sample F-stats and your original F-stat (small bias is good!).
• std. error: The standard deviation of the bootstrap F-stats, which estimates the variability of the F-statistic.

3. Calculating p-Values for Effects #

Your intuition about counting extreme bootstrap statistics was on the mark! For each effect, we calculate how many bootstrap samples produced an F-statistic greater than or equal to your original F-stat—this gives us a p-value.

Use this code to compute p-values:

# Calculate p-values for each effect
p_situation <- mean(boot_results$t[,1] >= boot_results$t0[1])
p_cognitive_load <- mean(boot_results$t[,2] >= boot_results$t0[2])
p_interaction <- mean(boot_results$t[,3] >= boot_results$t0[3])

# Print results
cat("Situation Main Effect p-value:", round(p_situation, 3), "\n")
cat("Cognitive Load Main Effect p-value:", round(p_cognitive_load, 3), "\n")
cat("Interaction Effect p-value:", round(p_interaction, 3), "\n")

• A small p-value (e.g., < 0.05) means your original F-stat is extreme compared to the bootstrap distribution, indicating a significant effect.

4. Do You Need Multiple Comparisons Correction? #

Yes! You're testing 3 effects (2 main effects + 1 interaction), so you should adjust your significance threshold to avoid Type I errors. Two common options:

• Bonferroni correction: Divide your alpha level (e.g., 0.05) by the number of tests: 0.05 / 3 ≈ 0.0167. Only consider effects with p-values < 0.0167 as significant.
• Holm correction: A more powerful alternative to Bonferroni (less conservative). You can apply it using the p.adjust() function:

  adjusted_p <- p.adjust(c(p_situation, p_cognitive_load, p_interaction), method = "holm")
  names(adjusted_p) <- c("Situation", "Cognitive_Load", "Interaction")
  print(adjusted_p)
  

5. Why Your Original Code Didn't Work #

• You used lm instead of a within-subjects ANOVA model, which ignores the correlation between repeated measurements from the same participant.
• Resampling rows instead of participants breaks the within-subjects structure, leading to invalid bootstrap samples.

Don't worry—this is a common pitfall when learning bootstrap for repeated measures, and you're already ahead by recognizing the need for non-parametric methods with non-normal data!

内容的提问来源于stack exchange，提问作者Merle