This vignette demonstrates how to evaluate model performance using spatial cross-validation and compare different validation methods to understand the impact of spatial dependence.
spatialcvR provides standard model evaluation
metrics:
library(spatialcvR)
# Example predictions vs observations
observed <- c(10, 15, 20, 25, 30, 35, 40, 45, 50, 55)
predicted <- c(11, 14, 21, 24, 31, 34, 41, 44, 51, 54)
# Calculate metrics
metrics <- spatial_metrics(observed, predicted)
print(metrics)## Model Evaluation Metrics
## =======================
## Number of observations: 10
## RMSE: 1
## MAE: 1
## R2: 0.9952
## MAPE: 4.04 %
# Perfect prediction
observed_perfect <- c(10, 20, 30, 40, 50)
predicted_perfect <- c(10, 20, 30, 40, 50)
metrics_perfect <- spatial_metrics(observed_perfect, predicted_perfect)
print(metrics_perfect)## Model Evaluation Metrics
## =======================
## Number of observations: 5
## RMSE: 0
## MAE: 0
## R2: 1
## MAPE: 0 %
# Poor prediction
observed_poor <- c(10, 20, 30, 40, 50)
predicted_poor <- c(50, 40, 30, 20, 10)
metrics_poor <- spatial_metrics(observed_poor, predicted_poor)
print(metrics_poor)## Model Evaluation Metrics
## =======================
## Number of observations: 5
## RMSE: 28.2843
## MAE: 24
## R2: -3
## MAPE: 126 %
# Load sample data
data(sample_spatial_data)
# Create spatial folds
folds_block <- spatial_folds(
data = sample_spatial_data,
x = "longitude",
y = "latitude",
k = 5,
method = "block",
seed = 123
)
# Create random folds for comparison
folds_random <- spatial_folds(
data = sample_spatial_data,
x = "longitude",
y = "latitude",
k = 5,
method = "random",
seed = 123
)For this example, we’ll simulate model training and evaluation:
# Simulate model training and evaluation for each fold
evaluate_fold <- function(fold, data, target_var) {
train_idx <- fold$train
test_idx <- fold$test
# Simple linear model for demonstration
train_data <- data[train_idx, ]
test_data <- data[test_idx, ]
# Train model (using variable1 and variable2 as predictors)
model <- lm(target ~ variable1 + variable2, data = train_data)
# Make predictions
predictions <- predict(model, newdata = test_data)
observed <- test_data[[target_var]]
# Calculate metrics
metrics <- spatial_metrics(observed, predictions, na.rm = TRUE)
return(metrics)
}
# Evaluate spatial block CV
block_results <- lapply(folds_block$folds, function(fold) {
evaluate_fold(fold, sample_spatial_data, "target")
})
# Evaluate random CV
random_results <- lapply(folds_random$folds, function(fold) {
evaluate_fold(fold, sample_spatial_data, "target")
})# Aggregate metrics across folds
aggregate_metrics <- function(results_list) {
rmse_values <- sapply(results_list, function(x) x$RMSE)
mae_values <- sapply(results_list, function(x) x$MAE)
r2_values <- sapply(results_list, function(x) x$R2)
list(
RMSE_mean = mean(rmse_values),
RMSE_sd = sd(rmse_values),
MAE_mean = mean(mae_values),
MAE_sd = sd(mae_values),
R2_mean = mean(r2_values),
R2_sd = sd(r2_values)
)
}
block_summary <- aggregate_metrics(block_results)
random_summary <- aggregate_metrics(random_results)
cat("Spatial Block CV Summary:\n")## Spatial Block CV Summary:
## $RMSE_mean
## [1] 6.712582
##
## $RMSE_sd
## [1] 1.223729
##
## $MAE_mean
## [1] 5.477149
##
## $MAE_sd
## [1] 0.9768991
##
## $R2_mean
## [1] 0.7782187
##
## $R2_sd
## [1] 0.1693295
##
## Random CV Summary:
## $RMSE_mean
## [1] 6.83486
##
## $RMSE_sd
## [1] 0.5242892
##
## $MAE_mean
## [1] 5.522551
##
## $MAE_sd
## [1] 0.2974677
##
## $R2_mean
## [1] 0.8758788
##
## $R2_sd
## [1] 0.01977151
# Create summary metrics for comparison
block_summary_metrics <- list(
RMSE = block_summary$RMSE_mean,
MAE = block_summary$MAE_mean,
R2 = block_summary$R2_mean
)
random_summary_metrics <- list(
RMSE = random_summary$RMSE_mean,
MAE = random_summary$MAE_mean,
R2 = random_summary$R2_mean
)
# Compare methods
comparison <- compare_cv(
results_list = list(
spatial_block = block_summary_metrics,
random = random_summary_metrics
)
)
print(comparison)## Cross-Validation Method Comparison
## ===================================
## Number of methods compared: 2
## Best method (lowest RMSE): spatial_block
##
## Summary Table:
## method RMSE MAE R2 RMSE_diff RMSE_diff_pct
## spatial_block 6.712582 5.477149 0.7782187 0.0000000 0.000000
## random 6.834860 5.522551 0.8758788 0.1222781 1.821625
##
## Performance Comparison:
## Baseline method: spatial_block (RMSE = 6.7126)
## Best RMSE: 6.7126
## RMSE range: [6.7126, 6.8349]
## MAE range: [5.4771, 5.5226]
## R2 range: [0.7782, 0.8759]
# Get predictions from one fold for residual analysis
fold_idx <- 1
train_idx <- folds_block$folds[[fold_idx]]$train
test_idx <- folds_block$folds[[fold_idx]]$test
train_data <- sample_spatial_data[train_idx, ]
test_data <- sample_spatial_data[test_idx, ]
# Train model
model <- lm(target ~ variable1 + variable2, data = train_data)
# Make predictions
predictions <- predict(model, newdata = test_data)
observed <- test_data$target
# Remove NA values
valid_idx <- !is.na(predictions) & !is.na(observed)
predictions <- predictions[valid_idx]
observed <- observed[valid_idx]
coords <- cbind(
x = test_data$longitude[valid_idx],
y = test_data$latitude[valid_idx]
)
# Calculate spatial residuals
residuals <- spatial_residuals(observed, predictions, coords)
print(residuals)# Spatial scatter plot of residuals
plot_spatial_residuals(residuals, type = "scatter",
main = "Spatial Residuals - Scatter Plot")
# Histogram of residuals
# (Chunk disabled due to dependency on previous disabled chunk)
# plot_spatial_residuals(residuals, type = "histogram",
# main = "Spatial Residuals - Histogram")
# Q-Q plot for normality
# (Chunk disabled due to dependency on previous disabled chunk)
# plot_spatial_residuals(residuals, type = "qq",
# main = "Spatial Residuals - Q-Q Plot")# Simulate agricultural data
set.seed(456)
n <- 100
agri_data <- data.frame(
longitude = runif(n, 0, 1000),
latitude = runif(n, 0, 1000),
soil_quality = rnorm(n, 50, 10),
rainfall = rnorm(n, 800, 100),
temperature = rnorm(n, 20, 3),
yield = 50 + 0.3*soil_quality + 0.02*rainfall + 0.5*temperature +
rnorm(n, 0, 5)
)
# Add spatial structure to yield
agri_data$yield <- agri_data$yield +
0.01*agri_data$longitude + 0.02*agri_data$latitude
# Create spatial folds
agri_folds <- spatial_folds(
data = agri_data,
x = "longitude",
y = "latitude",
k = 5,
method = "block",
seed = 456
)
# Check for spatial leakage
agri_leakage <- detect_spatial_leakage(
data = agri_data,
folds = agri_folds,
x = "longitude",
y = "latitude"
)
print(agri_leakage)
# Evaluate model
agri_results <- lapply(agri_folds$folds, function(fold) {
evaluate_fold(fold, agri_data, "yield")
})
agri_summary <- aggregate_metrics(agri_results)
print(agri_summary)# Compare multiple spatial CV methods
methods <- c("block", "cluster", "random")
method_results <- list()
for (method in methods) {
folds <- spatial_folds(
data = sample_spatial_data,
x = "longitude",
y = "latitude",
k = 5,
method = method,
seed = 123
)
results <- lapply(folds$folds, function(fold) {
evaluate_fold(fold, sample_spatial_data, "target")
})
summary <- aggregate_metrics(results)
method_results[[method]] <- list(
RMSE = summary$RMSE_mean,
MAE = summary$MAE_mean,
R2 = summary$R2_mean
)
}
# Compare all methods
multi_comparison <- compare_cv(method_results)
print(multi_comparison)# Test different levels of spatial separation
k_values <- c(3, 5, 7, 10)
separation_results <- data.frame(
k = integer(),
mean_distance = numeric(),
RMSE = numeric()
# (Chunks disabled for CRAN build - they work interactively but not during build)
# For demonstration, please see the other vignettes or run the examples interactively
)Problem: Using random CV on spatial data without checking for spatial dependence.
Solution: Always compare spatial vs. random CV to assess the impact.
Problem: Using too much spatial information in predictors.
Solution: Keep spatial coordinates out of predictor variables unless explicitly modeling spatial effects.
Problem: Calculating distances on geographic coordinates.
Solution: Use projected CRS for distance-based methods, or acknowledge approximation.
Problem: Too few observations relative to spatial complexity.
Solution: Reduce k, use larger blocks, or collect more data.