R语言：如何用apply系列函数向量化代码，替代循环提升效率

Hey there! I totally feel your pain—running a for loop over 23k rows is a surefire way to lock up your 16GB RAM and waste hours waiting. Let's ditch the slow loop and use vectorized approaches from the apply family (plus some spatial data tricks) to get this done in a fraction of the time.

The key here is to leverage R's optimized vectorized functions instead of iterating row-by-row. We'll handle spatial distance calculations, date filtering, and state-based grouping in bulk, which is way more efficient.

Step 1: Preprocess Your Data First #

First, make sure your date and spatial data are in the right format—this avoids headaches later:

# Ensure date is properly formatted as Date type
event$date <- as.Date(event$date)

# Use the sf package for fast spatial calculations (way better than manual distance formulas)
library(sf)
# Convert your data to a spatial object (WGS84 is standard for lat/lon)
event_sf <- st_as_sf(event, coords = c("LON", "LAT"), crs = 4326)
# Transform to a planar CRS (like UTM) to calculate distances in meters (adjust CRS to match your data's region)
event_sf_utm <- st_transform(event_sf, crs = 32610) # Example: UTM Zone 10N for western US

Step 2: Calculate Spatial + Temporal Event Counts (No More For Loops) #

We'll use vapply (from base R's apply family) to process each row, but the heavy lifting is done by vectorized spatial/date operations under the hood.

Option 1: Base R with vapply #

# Define a function to count matching events for a single row
count_matching_events <- function(row_idx, data, max_dist_km, max_days_future) {
  # Grab current row's key values
  current_state <- data$STATEid[row_idx]
  current_date <- data$date[row_idx]
  current_point <- data$geometry[row_idx]
  
  # Filter to only same-state records first (cuts down on calculations)
  same_state_subset <- data[data$STATEid == current_state, ]
  
  # Calculate distances from current point to all same-state points
  distances <- st_distance(current_point, same_state_subset$geometry)
  # Keep only records within the distance threshold (convert km to meters)
  dist_pass <- distances <= max_dist_km * 1000
  
  # Keep only records in the future time window (exclude current row if needed)
  date_pass <- same_state_subset$date > current_date & same_state_subset$date <= current_date + max_days_future
  
  # Count how many records meet both criteria
  sum(dist_pass & date_pass)
}

# Calculate counts for 50km + 30 days
event$count_50km_30d <- vapply(
  1:nrow(event_sf_utm),
  count_matching_events,
  FUN.VALUE = integer(1),
  data = event_sf_utm,
  max_dist_km = 50,
  max_days_future = 30
)

# Repeat for other combinations (100km/30d, 50km/60d, 100km/60d)
event$count_100km_30d <- vapply(
  1:nrow(event_sf_utm),
  count_matching_events,
  FUN.VALUE = integer(1),
  data = event_sf_utm,
  max_dist_km = 100,
  max_days_future = 30
)

Option 2: Tidyverse with purrr (Alternative to base apply) #

If you prefer the tidyverse syntax, purrr::map_int works just as well and is more readable for some:

library(purrr)

event$count_50km_30d <- map_int(1:nrow(event_sf_utm), ~{
  current_row <- event_sf_utm[.x, ]
  same_state <- event_sf_utm[event_sf_utm$STATEid == current_row$STATEid, ]
  dist_ok <- st_distance(current_row$geometry, same_state$geometry) <= 50 * 1000
  date_ok <- same_state$date > current_row$date & same_state$date <= current_row$date + 30
  sum(dist_ok & date_ok)
})

Step 3: Add Case/Viold Counts for Before/After Time Windows #

For your supplementary request, we'll adjust the date filter to include both past and future days, then sum the CASE and viold fields. Again, vectorized operations are key:

# Function to count CASE and viold for a row's time window
count_case_viold_window <- function(row_idx, data, day_window) {
  current_state <- data$STATEid[row_idx]
  current_date <- data$date[row_idx]
  
  # Filter same-state records within the time window
  same_state_window <- data[
    data$STATEid == current_state &
    data$date >= current_date - day_window &
    data$date <= current_date + day_window,
  ]
  
  # Sum the two fields
  c(
    CASE_total = sum(same_state_window$CASE),
    viold_total = sum(same_state_window$viold)
  )
}

# Calculate 30-day window counts
case_viold_30 <- t(vapply(
  1:nrow(event),
  count_case_viold_window,
  FUN.VALUE = integer(2),
  data = event,
  day_window = 30
))
# Merge results back to original data
event$CASE_30d <- case_viold_30[, "CASE_total"]
event$viold_30d <- case_viold_30[, "viold_total"]

# Repeat for 60-day window
case_viold_60 <- t(vapply(
  1:nrow(event),
  count_case_viold_window,
  FUN.VALUE = integer(2),
  data = event,
  day_window = 60
))
event$CASE_60d <- case_viold_60[, "CASE_total"]
event$viold_60d <- case_viold_60[, "viold_total"]

Pro Tip: Even Faster Performance with Grouped Precomputation #

If you still need more speed, precompute distance and date matrices for each state group. This avoids repeating the same calculations over and over:

# Split data into groups by STATEid
state_groups <- split(event_sf_utm, event_sf_utm$STATEid)

# Function to process a single state group
process_state_group <- function(group_data, max_dist_km, max_days_future) {
  n_rows <- nrow(group_data)
  counts <- integer(n_rows)
  
  # Precompute all pairwise distances in the group
  dist_matrix <- st_distance(group_data$geometry)
  # Precompute all pairwise date checks (future window)
  date_matrix <- outer(
    group_data$date,
    group_data$date,
    FUN = function(x, y) y > x & y <= x + max_days_future
  )
  
  # Use precomputed matrices to count matches quickly
  for(i in 1:n_rows) {
    counts[i] <- sum(dist_matrix[i, ] <= max_dist_km * 1000 & date_matrix[i, ])
  }
  
  group_data$count <- counts
  return(group_data)
}

# Process all groups and combine back into one dataset
optimized_results <- lapply(
  state_groups,
  process_state_group,
  max_dist_km = 50,
  max_days_future = 30
)
event_final <- do.call(rbind, optimized_results)

This method cuts down on redundant work because we only calculate distances and date checks once per state, not once per row.

Why This Works Better Than Your Original Loop #

• Vectorized functions: st_distance, outer, and sum are all optimized in low-level languages (C/Fortran) instead of R's interpreted loop.
• Reduced redundancy: Grouping by state means we don't re-filter same-state records for every single row.
• Efficient spatial handling: The sf package is built for fast spatial operations, way better than manual Haversine distance calculations.

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