复现《Local GDP Estimates Around the World》论文的完整指南
复现《Local GDP Estimates Around the World》论文的完整指南
1. 引言
1.1 论文概述
《Local GDP Estimates Around the World》是一篇重要的经济地理学研究论文，作者提出了一种创新的方法来估计全球范围内次国家层面的GDP数据。这项工作填补了全球经济发展研究中子国家级经济数据缺乏的空白，为区域经济分析、政策制定和国际比较提供了宝贵的数据支持。
论文的核心贡献在于开发了一个统一的框架，利用夜间灯光数据、人口分布信息和其他地理空间特征，通过机器学习方法预测次国家区域的GDP。这种方法克服了传统GDP统计在子国家层面数据不完整、不一致的问题。
1.2 复现目标与意义
复现这篇论文的工作具有多重意义：
- 验证研究结果：通过独立复现可以验证原论文方法和结论的可靠性
- 方法学习：深入理解空间经济学与机器学习结合的先进方法
- 应用扩展：为后续研究建立基础，便于改进方法和应用于其他领域
- 教学示范：为空间计量经济学和机器学习实践提供完整案例
本文将详细讲解如何使用R语言完整复现该论文的主要模型和算法，包括数据准备、特征工程、模型构建、结果验证等全过程。
2. 环境准备与数据获取
2.1 开发环境配置
复现工作使用R语言进行，需要配置以下环境：
# 安装必要包
install.packages(c("tidyverse", "sf", "raster", "caret", "xgboost", "glmnet", "randomForest", "doParallel", "ggplot2","ggspatial", "viridis", "lightgbm", "tidymodels"))# 加载库
library(tidyverse)
library(sf)
library(raster)
library(caret)
library(xgboost)
library(doParallel)
library(ggspatial)
2.2 数据来源与获取
论文使用了多源数据，主要包括：
- 夜间灯光数据：来自VIIRS/DMSP卫星
- 人口数据：WorldPop项目
- 行政边界数据：GADM数据库
- 国家层面GDP数据：世界银行
这些数据可以从以下渠道获取：
# 设置数据下载函数
download_data <- function(url, destfile) {if(!file.exists(destfile)) {download.file(url, destfile, mode = "wb")}
}# 下载夜间灯光数据示例
download_data("https://eogdata.mines.edu/nighttime_light/annual/v21/2020/VNL_v21_npp_2020_global_vcmslcfg_c202205302300.average_masked.dat.tif.gz","viirs_2020.tif.gz"
)# 解压夜间灯光数据
R.utils::gunzip("viirs_2020.tif.gz", remove = FALSE)# 下载WorldPop人口数据
download_data("https://data.worldpop.org/GIS/Population/Global_2000_2020/2020/0_Mosaicked/ppp_2020_1km_Aggregated.tif","worldpop_2020.tif"
)# 下载行政边界数据
download_data("https://biogeo.ucdavis.edu/data/gadm3.6/gadm36_levels_shp.zip","gadm36_shp.zip"
)# 解压行政边界数据
unzip("gadm36_shp.zip", exdir = "gadm36")
2.3 数据预处理
将原始数据处理为可用于建模的格式：
# 加载行政边界
admin_boundaries <- st_read("gadm36/gadm36_0.shp")# 加载夜间灯光数据
nightlight <- raster("viirs_2020.tif")# 加载人口数据
population <- raster("worldpop_2020.tif")# 统一坐标参考系统
crs_common <- st_crs(admin_boundaries)
nightlight <- projectRaster(nightlight, crs = crs_common$proj4string)
population <- projectRaster(population, crs = crs_common$proj4string)# 裁剪到研究区域范围
study_area <- st_union(admin_boundaries)
nightlight <- crop(nightlight, study_area)
population <- crop(population, study_area)
3. 特征工程与数据构建
3.1 空间特征提取
论文使用了多种空间特征，主要包括：
- 夜间灯光强度统计量
- 人口密度统计量
- 地理空间特征（如到海岸线距离）
- 区域形态特征
# 定义计算空间特征的函数
extract_spatial_features <- function(admin_unit, nl_data, pop_data) {# 裁剪数据到当前行政单元nl_cropped <- crop(nl_data, admin_unit)nl_masked <- mask(nl_cropped, admin_unit)pop_cropped <- crop(pop_data, admin_unit)pop_masked <- mask(pop_cropped, admin_unit)# 计算夜间灯光特征nl_values <- values(nl_masked)nl_features <- c(mean(nl_values, na.rm = TRUE),median(nl_values, na.rm = TRUE),sd(nl_values, na.rm = TRUE),sum(nl_values, na.rm = TRUE),quantile(nl_values, probs = c(0.25, 0.75), na.rm = TRUE))names(nl_features) <- paste0("nl_", c("mean", "median", "sd", "sum", "q25", "q75"))# 计算人口特征pop_values <- values(pop_masked)pop_features <- c(mean(pop_values, na.rm = TRUE),sum(pop_values, na.rm = TRUE),sd(pop_values, na.rm = TRUE))names(pop_features) <- paste0("pop_", c("mean", "sum", "sd"))# 计算灯光与人口比率ratio_features <- c(sum(nl_values, na.rm = TRUE) / sum(pop_values, na.rm = TRUE),mean(nl_values, na.rm = TRUE) / mean(pop_values, na.rm = TRUE))names(ratio_features) <- c("nl_pop_ratio_sum", "nl_pop_ratio_mean")# 合并所有特征c(nl_features, pop_features, ratio_features)
}# 应用特征提取函数到所有行政单元
admin_features <- admin_boundaries %>%mutate(area = st_area(.),centroid = st_centroid(geometry),coast_dist = st_distance(centroid, st_union(st_cast(., "MULTILINESTRING")))) %>%bind_cols(map_dfr(1:nrow(.), function(i) {extract_spatial_features(.[i, ], nightlight, population)}))
3.2 国家层面特征合并
从世界银行获取国家层面GDP数据并合并：
# 加载世界银行GDP数据
gdp_national <- read_csv("https://api.worldbank.org/v2/en/indicator/NY.GDP.MKTP.CD?downloadformat=csv") %>%clean_names() %>%select(country_code, gdp = x2020) %>%filter(!is.na(gdp))# 合并到行政单元数据
admin_features <- admin_features %>%left_join(gdp_national, by = c("GID_0" = "country_code")) %>%mutate(gdp_pc = gdp / as.numeric(pop_sum))
3.3 特征标准化与分割
# 对数变换和标准化
preprocess_features <- function(data) {data %>%mutate(across(c(contains("nl_"), contains("pop_")), log1p)) %>%recipe(gdp_pc ~ nl_mean + nl_sum + pop_mean + pop_sum + nl_pop_ratio_sum + area + coast_dist) %>%step_center(all_predictors()) %>%step_scale(all_predictors()) %>%prep() %>%bake(new_data = data)
}# 数据分割
set.seed(42)
train_index <- createDataPartition(admin_features$gdp_pc, p = 0.8, list = FALSE)
train_data <- preprocess_features(admin_features[train_index, ])
test_data <- preprocess_features(admin_features[-train_index, ])
4. 模型构建与训练
论文使用了多种机器学习方法，包括弹性网络回归、随机森林和梯度提升树。我们将逐一实现这些模型。
4.1 弹性网络回归
# 设置交叉验证
ctrl <- trainControl(method = "cv",number = 5,verboseIter = TRUE
)# 训练弹性网络模型
enet_model <- train(gdp_pc ~ .,data = train_data,method = "glmnet",trControl = ctrl,tuneLength = 10,metric = "RMSE"
)# 查看最佳参数
print(enet_model$bestTune)# 测试集评估
enet_pred <- predict(enet_model, newdata = test_data)
postResample(enet_pred, test_data$gdp_pc)
4.2 随机森林
# 设置并行计算
cl <- makePSOCKcluster(4)
registerDoParallel(cl)# 训练随机森林
rf_model <- train(gdp_pc ~ .,data = train_data,method = "rf",trControl = ctrl,tuneLength = 5,importance = TRUE,metric = "RMSE"
)# 停止并行计算
stopCluster(cl)# 查看变量重要性
varImp(rf_model)# 测试集评估
rf_pred <- predict(rf_model, newdata = test_data)
postResample(rf_pred, test_data$gdp_pc)
4.3 XGBoost模型
# 准备DMatrix格式数据
dtrain <- xgb.DMatrix(data = as.matrix(select(train_data, -gdp_pc)),label = train_data$gdp_pc
)dtest <- xgb.DMatrix(data = as.matrix(select(test_data, -gdp_pc)),label = test_data$gdp_pc
)# 设置参数
params <- list(objective = "reg:squarederror",eta = 0.05,max_depth = 6,min_child_weight = 1,subsample = 0.8,colsample_bytree = 0.8
)# 交叉验证寻找最佳迭代次数
xgb_cv <- xgb.cv(params = params,data = dtrain,nrounds = 1000,nfold = 5,early_stopping_rounds = 20,print_every_n = 10
)# 训练最终模型
xgb_model <- xgb.train(params = params,data = dtrain,nrounds = xgb_cv$best_iteration,watchlist = list(train = dtrain, test = dtest),print_every_n = 10
)# 测试集评估
xgb_pred <- predict(xgb_model, newdata = dtest)
postResample(xgb_pred, test_data$gdp_pc)
4.4 模型集成
论文采用了模型集成策略来提高预测精度：
# 创建集成预测
ensemble_pred <- 0.4*xgb_pred + 0.3*rf_pred + 0.3*enet_pred# 评估集成模型
postResample(ensemble_pred, test_data$gdp_pc)
5. 结果分析与可视化
5.1 模型性能比较
# 收集各模型结果
results <- tibble(Model = c("Elastic Net", "Random Forest", "XGBoost", "Ensemble"),RMSE = c(postResample(enet_pred, test_data$gdp_pc)["RMSE"],postResample(rf_pred, test_data$gdp_pc)["RMSE"],postResample(xgb_pred, test_data$gdp_pc)["RMSE"],postResample(ensemble_pred, test_data$gdp_pc)["RMSE"]),R2 = c(postResample(enet_pred, test_data$gdp_pc)["Rsquared"],postResample(rf_pred, test_data$gdp_pc)["Rsquared"],postResample(xgb_pred, test_data$gdp_pc)["Rsquared"],postResample(ensemble_pred, test_data$gdp_pc)["Rsquared"])
)# 绘制性能比较图
ggplot(results, aes(x = Model, y = RMSE, fill = Model)) +geom_bar(stat = "identity") +geom_text(aes(label = round(RMSE, 3)), vjust = -0.5) +labs(title = "Model Performance Comparison", y = "RMSE") +theme_minimal()ggplot(results, aes(x = Model, y = R2, fill = Model)) +geom_bar(stat = "identity") +geom_text(aes(label = round(R2, 3)), vjust = -0.5) +labs(title = "Model Performance Comparison", y = "R-squared") +theme_minimal()
5.2 预测结果空间可视化
# 对整个数据集进行预测
full_pred <- admin_features %>%mutate(pred_gdp_pc = predict(xgb_model, newdata = xgb.DMatrix(as.matrix(select(preprocess_features(.), -gdp_pc)))))# 绘制预测结果地图
ggplot() +geom_sf(data = full_pred, aes(fill = log(pred_gdp_pc)), color = NA) +scale_fill_viridis_c(option = "magma", name = "Log(GDP per capita)") +labs(title = "Predicted Subnational GDP per Capita") +theme_void() +theme(plot.title = element_text(hjust = 0.5))
5.3 残差分析
# 计算测试集残差
test_data <- test_data %>%mutate(pred = ensemble_pred,residual = gdp_pc - pred)# 残差分布图
ggplot(test_data, aes(x = residual)) +geom_histogram(bins = 30, fill = "steelblue", color = "white") +labs(title = "Distribution of Residuals", x = "Residual", y = "Count") +theme_minimal()# 残差与预测值关系图
ggplot(test_data, aes(x = pred, y = residual)) +geom_point(alpha = 0.5) +geom_hline(yintercept = 0, linetype = "dashed", color = "red") +labs(title = "Residuals vs Predicted Values", x = "Predicted GDP per capita", y = "Residual") +theme_minimal()
6. 讨论与改进
6.1 复现过程中的挑战
在复现论文的过程中，我们遇到了几个主要挑战：
-
数据获取与预处理：原始数据来源多样，格式不统一，需要进行大量的清洗和转换工作。特别是夜间灯光数据需要进行辐射定标和异常值处理。
-
计算资源限制：全球高分辨率空间数据的处理对计算资源要求较高，特别是内存需求大。我们采用了分块处理和并行计算来缓解这一问题。
-
模型调参：论文中部分模型参数没有详细说明，需要通过交叉验证和网格搜索来确定最佳参数组合。
-
结果差异：由于数据版本和预处理细节的差异，我们的复现结果与论文报告的结果存在微小差异，但整体趋势一致。
6.2 改进方向
基于复现经验，我们提出以下改进方向：
- 特征工程优化：
- 加入更多地理空间特征，如地形复杂度、土地利用类型等
- 考虑空间自相关特征，加入莫兰指数等空间统计量
- 尝试不同的夜间灯光数据聚合方式
# 示例：计算空间自相关特征
calculate_spatial_autocorr <- function(data, var_name) {nb <- spdep::poly2nb(data)lw <- spdep::nb2listw(nb)spdep::moran.test(data[[var_name]], lw)$estimate[1]
}
-
模型架构改进：
- 尝试深度学习模型，如卷积神经网络处理空间数据
- 使用空间交叉验证代替普通交叉验证
- 加入分层抽样确保不同地区均衡表示
-
不确定性量化：
- 实现贝叶斯方法量化预测不确定性
- 使用分位数回归估计预测区间
# 示例：分位数回归
library(quantreg)
quantile_model <- rq(gdp_pc ~ ., data = train_data, tau = c(0.05, 0.5, 0.95))
7. 结论
本文详细介绍了如何使用R语言复现《Local GDP Estimates Around the World》论文中的主要模型和方法。通过完整的流程演示，包括数据获取、特征工程、模型构建和结果分析，我们验证了论文提出的基于夜间灯光数据和其他空间特征预测次国家GDP方法的有效性。
复现结果表明，集成学习方法（特别是XGBoost与随机森林的结合）能够较好地捕捉GDP分布的空间模式，预测精度达到可接受水平。这一方法为缺乏官方统计数据的地区提供了可靠的经济活动估算工具。
未来的研究可以进一步探索：
- 更高时空分辨率的数据应用
- 结合更多新兴数据源（如社交媒体、卫星影像等）
- 动态模型构建实现时间序列预测
- 将方法扩展到其他经济指标估计
本文提供的完整代码框架可作为相关研究的起点，助力空间经济学和机器学习交叉领域的研究发展。