Bootstrap加权中位数及置信区间：R代码实现咨询

Got it, let's walk through how to implement a weighted median Bootstrap confidence interval using R's boot package. I'll build out your partial code into a complete, working solution, and explain the key parts so you understand what's going on.

Step 1: Load Packages and Prepare Data #

First, let's get our data and weights set up (I've kept your original values intact):

library(boot)

# Your original data and weight vector
data <- c(6, 11, 7, 8, 3, 9, 4, 1, 1, 8, 2, 2, 5, 3, 1)
weight <- c(0.839432605459112, 0.774215027235327, 0.709256693551626, 
            0.809376516981207, 0.809698716683444, 0.880849581474519, 
            0.829263837448813, 1.80390621483409, 1.12749447791778, 
            0.93389158146594, 1.07832286911631, 0.79541512406283, 
            1.06708509325217, 0.946752658104578, 0.968003233015867)

Step 2: Define a Weighted Median Function #

Base R's median() doesn't support weights, so we'll write a simple function to calculate the weighted median. This works by sorting the data, then finding the first value where the cumulative weight proportion reaches 0.5:

weighted_median <- function(x, weights) {
  # Sort the data and corresponding weights
  sorted_idx <- order(x)
  x_sorted <- x[sorted_idx]
  w_sorted <- weights[sorted_idx]
  
  # Calculate cumulative weight proportions
  cum_weight_pct <- cumsum(w_sorted) / sum(w_sorted)
  
  # Grab the first value where cumulative weight hits 50%
  x_sorted[which(cum_weight_pct >= 0.5)[1]]
}

Step 3: Build the Bootstrap Statistic Function #

We need a function that the boot package can use to repeatedly calculate the weighted median on resampled data. We'll combine our data and weights into a data frame so we can resample rows together:

# Combine data and weights into a single data frame
df <- data.frame(value = data, weight = weight)

# Bootstrap statistic function: resample rows, then compute weighted median
boot_stat <- function(df, resample_idx) {
  # Extract the resampled subset
  resampled_data <- df[resample_idx, ]
  # Calculate weighted median for this resample
  weighted_median(resampled_data$value, resampled_data$weight)
}

Step 4: Run the Bootstrap and Calculate Confidence Intervals #

Now we'll run 1000 resamples (you can adjust R for more precision), then compute both percentile and BCa (bias-corrected accelerated) confidence intervals:

# Set seed for reproducibility
set.seed(123)

# Run the Bootstrap
Mboot <- boot(data = df, statistic = boot_stat, R = 1000)

# View the basic Bootstrap results
print(Mboot)

# Calculate 95% confidence intervals
boot.ci(Mboot, type = c("perc", "bca"))

Alternative: Weighted Resampling Approach #

If your weights represent the relative probability of each observation being selected (e.g., survey sampling weights), you can instead resample observations proportional to their weights, then compute the regular median. Here's how that looks:

# Bootstrap function with weighted resampling
boot_stat_weighted <- function(df, dummy_idx) {
  # Resample observations proportional to their weights
  resampled_idx <- sample(nrow(df), size = nrow(df), replace = TRUE, prob = df$weight)
  # Compute regular median on the weighted resample
  median(df$value[resampled_idx])
}

# Run weighted resampling Bootstrap
set.seed(123)
Mboot_weighted <- boot(data = df, statistic = boot_stat_weighted, R = 1000)

# Get confidence intervals
boot.ci(Mboot_weighted, type = c("perc", "bca"))

Key Notes: #

• The first method (resample rows, compute weighted median) is ideal when weights represent the importance of each observation.
• The second method (weighted resampling, compute regular median) works best when weights represent how many times an observation would appear in an unweighted dataset.

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