For循环实现Bootstrap t检验分析蚂蚁物种数据的问题求助

Let's break down the problems in your code and walk through the fixes step by step. The core issues are incorrect index usage in your loop, improper grouping logic, and a function that doesn't return values properly—all leading to those frustrating NA p-values.

Key Problems in Your Original Code #

1. No index usage in the loop: You're reassigning ants = data.frame(mainbroca[,9:13]) every iteration and never using i to target a single ant species column. So you're testing all ants at once instead of one by one.
2. Incorrect grouping: ants == 0 creates a logical matrix for all ant columns, not a vector for the specific ant species you're testing. This leads to invalid subsetting of your response variables.
3. Result assignment bug: You're assigning all results to ac_pvals[1] instead of ac_pvals[i], so you overwrite the first position every time instead of filling the vector per species.
4. Function side effects: Your bootstrap_ttest assigns values to the global environment instead of returning them, so test1/test2 are NULL—hence the NA assignments.
5. Mismatched simulated data: Your fake data has ant columns starting at column 5, but your analysis targets columns 9-13, which breaks testing.

Fixed Bootstrap t-test Function #

First, let's rewrite the function to follow R best practices: return results explicitly, handle missing values, and avoid polluting the global environment.

bootstrap_ttest <- function(data1, data2, resamples = 1000){ 
  # Calculate observed difference in means (handle NAs)
  delta_real <- mean(data1, na.rm = TRUE) - mean(data2, na.rm = TRUE) 
  # Pool both datasets for null hypothesis resampling
  pooled_data <- c(data1, data2) 
  # Pre-allocate vector for null differences (faster than growing)
  null_differences <- numeric(resamples) 
  
  for(x in 1:resamples){ 
    # Sample with replacement to match original group sizes
    data1_null <- sample(pooled_data, size = length(data1), replace = TRUE) 
    data2_null <- sample(pooled_data, size = length(data2), replace = TRUE) 
    delta_null <- mean(data1_null, na.rm = TRUE) - mean(data2_null, na.rm = TRUE) 
    null_differences[x] <- delta_null 
  }
  
  # Two-tailed p-value calculation
  pvalue <- sum(abs(null_differences) > abs(delta_real)) / resamples 
  
  # Return all useful results in a list
  return(list(
    pvalue = pvalue,
    observed_delta = delta_real,
    null_distribution = null_differences
  ))
}

Corrected Loop for Batch Testing #

Now let's fix the loop to properly iterate over each ant species and response variable, with correct grouping and result storage.

# Define your target columns
ant_columns <- 9:13  # Columns for ant species in mainbroca
response_vars <- c("Attack_count", "Attack_percent", "Survival_count", "Survival_percent")

# Initialize a matrix to store p-values (rows = ant species, cols = responses)
pvalue_matrix <- matrix(
  nrow = length(ant_columns),
  ncol = length(response_vars),
  dimnames = list(colnames(mainbroca)[ant_columns], response_vars)
)

# Loop over each ant species
for(i in seq_along(ant_columns)){
  # Get the current ant species column
  current_ant <- mainbroca[, ant_columns[i]]
  # Create group filters: ant count == 0 vs ant count > 0
  group_zero <- current_ant == 0
  group_positive <- current_ant > 0
  
  # Loop over each response variable
  for(j in seq_along(response_vars)){
    # Extract response data for both groups
    response_data <- mainbroca[, response_vars[j]]
    data_zero <- response_data[group_zero]
    data_positive <- response_data[group_positive]
    
    # Skip if either group has no data (avoids errors/NA from empty vectors)
    if(length(data_zero) == 0 || length(data_positive) == 0){
      pvalue_matrix[i, j] <- NA
      next
    }
    
    # Run the bootstrap test and store the p-value
    test_result <- bootstrap_ttest(data_zero, data_positive)
    pvalue_matrix[i, j] <- test_result$pvalue
  }
}

# View your final p-value table
print(pvalue_matrix)

Reproducible Simulated Data #

Let's fix the fake data to match your column structure (9-13 are ant species) so you can test the code properly:

set.seed(123)  # Ensure reproducibility
mainbroca <- data.frame(
  # First 8 non-ant columns (fill with dummy data)
  matrix(rnorm(8*100), nrow = 100),
  # Ant species columns (9-13)
  ant1 = sample(c(0, 1, 2), 100, replace = TRUE),
  ant2 = sample(c(0, 3, 5), 100, replace = TRUE),
  ant3 = sample(c(0, 1), 100, replace = TRUE),
  ant4 = sample(c(0, 4), 100, replace = TRUE),
  ant5 = sample(c(0, 2, 3), 100, replace = TRUE),
  # Response variables (add meaningful differences for testing)
  Attack_count = rpois(100, lambda = ifelse(mainbroca$ant1 > 0, 5, 2)),
  Attack_percent = runif(100, 0, 100),
  Survival_count = rpois(100, lambda = ifelse(mainbroca$ant1 > 0, 10, 15)),
  Survival_percent = runif(100, 0, 100)
)

# Rename first 8 columns for clarity
colnames(mainbroca)[1:8] <- paste0("misc_col_", 1:8)

What This Fixes #

• Proper indexing: The loop now uses i to target one ant species at a time.
• Valid grouping: We subset response variables based on the current ant species' 0/non-0 status, not all ants at once.
• Clean result storage: The matrix organizes p-values by ant species and response variable, making it easy to interpret.
• Robust function: The revised function returns results safely and handles missing values.
• Testable data: The simulated data matches your column structure, so you can run the code immediately to verify it works.

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