简介这套R语言实战案例聚焦大气污染数据分析面向有一定R基础的学员和环境数据爱好者适合希望从数据导入、清洗、统计到可视化走通全流程的读者。压缩包共10个文件包含Rmd分析脚本、HTML结果预览、TXT说明文档以及Rdata样本数据整体大小仅1.24MB。其中3个Rmd脚本分别演示变异系数热图、浓度柱状图与浓度热图的制作方法HTML文件对应输出结果TXT文件则交代了数据来源、变量含义及分析目标。已有9748人学习下载便于初学者快速上手。案例覆盖缺失值与异常值处理、dplyr数据整理、ggplot2绘图等核心操作也涉及利用forecast、caret等包开展趋势预测与建模的思路。Rdata数据文件可直接载入运行有助于读者边学边练理解污染物时空分布的分析方法。 很多人下载了公开的空气质量数据之后第一反应就是plot()、ggplot()上去猛画几张图然后盯着图发呆——图是真的好看但下一步该干什么数据里到底藏着什么信息却说不上来。这个大气污染数据分析案例就是冲着这个痛点去的。它不是一个“点一下出结果”的傻瓜式模板而是一套从数据体检、清洗、可视化到建模的完整R语言分析流程用的是一份真实结构的国控站点监测数据含PM2.5、PM10、SO2、NO2、CO、O3六项污染物和气象要素。无论你是刚学R的学生、做环境相关研究的科研人员还是想转行数据分析的职场人都能从这套案例里拿到可以直接改改就跑的代码以及比代码更重要的——分析思路。1. 拿到污染物数据后先别急着画图我在处理这份数据时踩过最大的坑就是“想一口气把图全画完”。后来发现急着画图的代价是要花三倍时间回头改数据格式。所以这套案例的第一步刻意安排成了数据体检而且是用最朴素的方式。1.1 数据结构和字段含义案例用的数据是CSV格式包含日期、站点、六项污染物浓度、四项气象要素。字段大致长这样# 用data.table读取速度更快 library(data.table) air_data - fread(air_data.csv) str(air_data) # 输出示例 # $ date : chr 2023-01-01 2023-01-02 ... # $ station : chr 站点A 站点A ... # $ pm25 : num 45.2 78.3 ... # $ pm10 : num 80.1 112.5 ... # $ so2 : num 8.2 10.1 ... # $ no2 : num 32.4 45.2 ... # $ co : num 0.8 1.2 ... # $ o3 : num 62.3 55.8 ... # $ temperature : num 3.2 1.5 ... # $ humidity : num 55 68 ... # $ windspeed : num 2.1 1.4 ... # $ wind_direction: num 180 220 ... # $ pressure : num 1023 1021 ...先看结构不是走形式。这一步能暴露大多数问题日期是不是字符型、有没有把站点读成因子、污染物浓度是不是数值型。如果哪列被读成了字符后面所有计算都会报错或者得出荒谬结果。我曾经遇到过一次co列因为有ND未检出标记被整体读成字符型后面算相关系数时直接全部NA排查了很久。1.2 缺失值和异常值检查大气污染监测数据有个特点不是每天每个站点都能出数。设备校准、通讯故障、停电都会导致缺失。更麻烦的是有些缺失会以NA形式存在有些却会以-9999、0甚至空格存在。# 缺失值概览 colSums(is.na(air_data)) # 更严格一点查看每列的“异常”情况 summary(air_data[, .(pm25, pm10, so2, no2, co, o3)])我的习惯是重点看summary()里的Min.和Max.。PM2.5出现0值要警惕因为很多仪器在浓度极低时会报0但更常见的情况是缺测被填成0SO2出现负值也可能是仪器零点漂移。对于后续要算均值、做相关分析的任务这些异常值会实实在在影响结果。这个案例里我的处理原则是先boxplot画出分布再用IQR规则标记极端值但在删之前一定先去看看原始记录不能机械地“非黑即白”。数据体检不是炫技是在帮后面的所有分析排雷。2. 数据清洗日期解析、缺失值和合并气象数据数据清洗是整个分析流程里最不性感但最重要的一步。大气污染数据的清洗有几个固定的坑这套案例里全部处理了一遍。2.1 日期解析别让时间轴变成乱麻现实里的日期格式五花八门2023/1/1、2023-01-01、01-01-2023以及Excel导出的45000这种序列值。如果时间轴没有正确解析时间序列图上的横坐标会乱成一团。案例里用lubridate包统一解析library(lubridate) air_data[, date : ymd(date)] # 或者按年月日拆开方便后面做月季分析 air_data[, :( year year(date), month month(date), season case_when( month %in% c(3,4,5) ~ 春季, month %in% c(6,7,8) ~ 夏季, month %in% c(9,10,11) ~ 秋季, TRUE ~ 冬季 ) )]这里有个细节ymd()会尝试自动识别格式但如果你的数据里有2023/1/1这种格式建议先统一成干净格式再解析。跨年数据更要小心我曾经因为日期列里混入了一个2023-02-30导致整个时间序列的排序错位画出来的线图在2月底突然断崖下跌花了很久才发现是脏数据。2.2 缺失值的处理策略区分“真缺”和“假缺”对于污染物时间序列缺失值处理我很少用粗暴的均值填充。原因很简单污染物浓度有明显的日变化和季节变化1月凌晨的PM2.5和7月下午的PM2.5差了十倍不止用一个全域均值填进去等于硬造虚假规律。案例里采用了两步法先删除关键污染物的缺失行如果缺失比例很低比如小于5%然后用线性插值填补剩余缺口。# 计算各列的缺失比例 missing_ratio - colMeans(is.na(air_data)) %% round(3) print(missing_ratio) # 对pm25用na.approx做线性插值 library(zoo) air_data[order(date), pm25_filled : na.approx(pm25, na.rm FALSE)]注意na.approx处理中间缺口没问题但首尾的NA它填不了需要单独用na.locf或者最近邻补齐。实务中如果首尾长期缺测我最倾向直接删掉那段数据避免人为制造“假数值”。2.3 气象数据合并把污染物放到环境下看大气污染不是孤立存在的气象条件直接决定了污染物的积累和扩散。风速大时PM2.5浓度通常低湿度高时往往利于二次颗粒物生成。案例把污染物数据和气象数据按日期合并成一张宽表为后面的相关性分析和建模做准备met_data - fread(met_data.csv) met_data[, date : ymd(date)] combined - merge(air_data, met_data, by c(date, station), all.x TRUE)我一直建议把“合并”当成一次验证的机会merge之后立刻检查行数有没有变多那是产生了笛卡尔积或变少那是丢失了匹配项。案例数据里两个文件日期是严格对齐的但你自己处理新数据时这一步最容易出问题。3. 时空分布可视化一张图讲清污染格局可视化不是用来“好看”的是用来快速定位问题的。大气污染数据最常见也最有价值的三个可视化方向这套案例里都写了代码。3.1 浓度时间序列先看整体趋势再看异常峰值折线图是第一张必须画的图。看什么看趋势、看季节、看极值。案例用ggplot2画了一张全年逐日变化曲线并用geom_smooth()叠加了平滑趋势线一眼就能看出秋冬季浓度明显高于夏季library(ggplot2) ggplot(combined, aes(x date, y pm25)) geom_line(colour grey60, linewidth 0.4, alpha 0.6) geom_smooth(method loess, colour #D73027, se FALSE, linewidth 1) labs(title PM2.5 日平均浓度变化趋势, x 日期, y 浓度 (μg/m³)) theme_classic(base_size 14)画完这张图我往往会顺手标出最高峰值对应的日期——那大概率对应某次重污染过程值得回去翻翻当时的天气和边界层条件。这种“图→现象→原因”的链条才是数据分析的价值所在。3.2 站点对比与箱线图不同站点的污染画像如果是多站点数据箱线图比折线图更能体现站点之间的差异。案例里按站点分组画了PM2.5的箱线图能清楚看到哪个站点的中位数高、哪个站点的离群值多ggplot(combined, aes(x station, y pm25, fill station)) geom_boxplot(outlier.colour black, outlier.size 0.8, alpha 0.7) labs(title 各站点 PM2.5 浓度分布对比, x 站点, y 浓度 (μg/m³)) theme_classic(base_size 14) theme(legend.position none)很多分析到此就停了其实还可以继续往下走把站点按区域类型工业区、居民区、交通点着上不同颜色往往能看出更有意思的规律。这个案例数据带了站点分类字段代码直接在原基础上加了一行fill station_type。3.3 污染物相关矩阵先猜后证不要盲目上模型做相关分析之前先画一张相关矩阵图能帮自己建立“数据直觉”。比如PM2.5和PM10因为同源相关性必然高CO和NO2都来自燃烧排放也高度相关。案例里用corrplot画了六项污染物加气象因子的相关矩阵library(corrplot) cor_matrix - combined[, .(pm25, pm10, so2, no2, co, o3, temperature, humidity, windspeed)] | cor(use complete.obs) corrplot(cor_matrix, method color, type upper, tl.col black, addCoef.col black, number.cex 0.7)这张图的价值在于后续建模选特征时你会很清楚哪些变量存在共线性哪些变量是相对独立的。比如PM2.5和PM10相关系数如果超过0.9那建模时就不该把它们同时放进去做普通线性回归否则会带来多重共线性的麻烦。4. 影响PM2.5浓度的气象因子量化分析大气污染数据分析如果只停在“画图”层面说服力还不够。气象因子到底怎么影响污染物浓度这个案例用回归和相关分析做了量化拆解。4.1 单因子相关分析哪个气象因素影响最大先做单因子检验算出PM2.5与温度、湿度、风速、气压的相关系数和p值with(combined, cor.test(pm25, temperature, method spearman)) with(combined, cor.test(pm25, humidity, method spearman)) with(combined, cor.test(pm25, windspeed, method spearman)) with(combined, cor.test(pm25, pressure, method spearman))这里我特意用了spearman而不是pearson因为PM2.5浓度和气象要素的关系通常不是线性的比如风速对PM2.5的“清除效应”在低风速段明显到高风速段趋于平缓。用秩相关更稳健。实测下来风速和PM2.5往往呈显著负相关湿度则常常呈正相关高湿促进二次颗粒物生成温度的影响在冬季和夏季方向是相反的所以做全样本相关时会被抵消掉这也是为什么后面需要再分层看。4.2 回归模型用一个方程描述多重影响单因子相关只能看“两两关系”实际大气过程是多因子共同作用的。案例用多元线性回归量化各气象因子的贡献model_lm - lm(pm25 ~ temperature humidity windspeed pressure, data combined) summary(model_lm)重点看三样东西Estimate的正负方向是否符合物理常识、Pr(|t|)显著的变量有哪些、Adjusted R-squared整体解释力如何。这个案例运行下来windspeed和pressure的系数通常显著为负humidity系数显著为正说明静稳、高湿、高压的气象条件确实利于PM2.5累积。回归不是为了得到一个完美的预测公式而是为了把“天气影响污染”这个模糊感觉变成一个可量化的系数。4.3 分层分析分季节看效应全样本拟合完之后我强烈建议做一层分季节的亚组分析。因为夏季PM2.5主要受O3和光化学反应影响冬季则受燃煤排放和静稳天气影响气象因子的作用方向和强度完全不同。案例里用split按季节拆开再分别跑回归season_models - lapply(split(combined, combined$season), function(df) { lm(pm25 ~ temperature humidity windspeed pressure, data df) }) lapply(season_models, summary)我见过最典型的情况是风速在全年回归里负效应明显但到夏季变得不显著因为夏季本身的扩散条件整体较好风速的边际效应被削弱了。这种分层结果写进报告里比只报一个全样本模型要专业得多。5. 预测模型的构建从线性回归到随机森林分析的最后一步是建模。模型不是越复杂越好而是越匹配问题越好。案例里我特意安排了三个层次的模型方便不同需求的读者对号入座。5.1 线性回归基准模型跑通一个可解释的底稿用污染物和气象变量预测PM2.5浓度线性回归作为基准模型有不可替代的价值可解释性强、计算快、稳定性好。代码上需要把类别变量转成哑变量并检查VIF方差膨胀因子确认没有严重的共线性library(car) model_lm_full - lm(pm25 ~ pm10 so2 no2 co temperature humidity windspeed pressure, data combined) vif(model_lm_full) # 如果VIF 10说明存在严重共线性需要剔除变量VIF的判断很重要。我在这个案例里遇到PM2.5和PM10的VIF超过20的情况后来把PM10从特征里去掉模型的稳定性立刻好了起来。这是很多初学者最容易忽视的环节。5.2 随机森林抓非线性关系随机森林的优势是不用预设函数形式能自动捕捉非线性关系和变量交互。案例用randomForest包训练了一个包含500棵树的模型library(randomForest) # 去除含NA的行并选择建模需要的列 model_data - combined[complete.cases(combined), .(pm25, so2, no2, co, o3, temperature, humidity, windspeed, pressure)] set.seed(42) rf_model - randomForest(pm25 ~ ., data model_data, ntree 500, importance TRUE) print(rf_model)随机森林跑起来不需要太多调参ntree500基本上够用关键是要看变量重要性图varImpPlot()。这个案例里排在前面的通常是pressure、humidity和so2说明静稳天气条件和燃煤源对PM2.5浓度的影响很突出。5.3 模型评估不看R²只看RMSE和泛化能力训练集上的R²再高也不代表真实预测能力。案例里用caret的createDataPartition做了七三拆分在测试集上比较不同模型的RMSE和R²library(caret) set.seed(123) train_index - createDataPartition(model_data$pm25, p 0.7, list FALSE) train_data - model_data[train_index, ] test_data - model_data[-train_index, ] # 线性模型在测试集上的评估 lm_fit - lm(pm25 ~ ., data train_data) lm_pred - predict(lm_fit, newdata test_data) lm_rmse - sqrt(mean((test_data$pm25 - lm_pred)^2)) lm_r2 - cor(test_data$pm25, lm_pred)^2 # 随机森林在测试集上的评估 rf_fit - randomForest(pm25 ~ ., data train_data, ntree 500) rf_pred - predict(rf_fit, newdata test_data) rf_rmse - sqrt(mean((test_data$pm25 - rf_pred)^2)) rf_r2 - cor(test_data$pm25, rf_pred)^2 cat(线性回归 RMSE:, lm_rmse, R²:, lm_r2, \n) cat(随机森林 RMSE:, rf_rmse, R²:, rf_r2, \n)这个案例实测下来随机森林比线性回归的RMSE低10%-20%但优势没有想象中夸张。这说明PM2.5浓度里很大一部分方差是被排放源决定的气象因子只能解释其中一部分这也符合大气科学的常识。另外如果数据里含有站点固定效应建议加一个station变量进来或者用混合效应模型能进一步提升预测精度。提示建模过程中set.seed()一定不要省。随机森林和训练集划分都涉及随机数不设种子的话你每次跑的结果都会不一样复现性就没法保证了。一点个人体会这套案例做下来我最大的收获不是哪个模型分数高而是整个分析链条的顺畅感先体检数据再清洗整合然后可视化找模式和异常接着用相关和回归量化关系最后建模验证。每一步都在为下一步铺路没有一步是白做的。你在跑这份代码的时候如果遇到报错优先检查数据格式和缺失值八成的问题都出在数据上。做数据分析这个事情慢就是快前期把数据弄干净后面所有步骤都会很顺。本文还有配套的精品资源点击获取 
长河君 · 2026/9/8 22:12:18 · 阅读 1.2k R语言大气污染数据分析：从数据清洗到随机森林建模实战
简介这套R语言实战案例聚焦大气污染数据分析面向有一定R基础的学员和环境数据爱好者适合希望从数据导入、清洗、统计到可视化走通全流程的读者。压缩包共10个文件包含Rmd分析脚本、HTML结果预览、TXT说明文档以及Rdata样本数据整体大小仅1.24MB。其中3个Rmd脚本分别演示变异系数热图、浓度柱状图与浓度热图的制作方法HTML文件对应输出结果TXT文件则交代了数据来源、变量含义及分析目标。已有9748人学习下载便于初学者快速上手。案例覆盖缺失值与异常值处理、dplyr数据整理、ggplot2绘图等核心操作也涉及利用forecast、caret等包开展趋势预测与建模的思路。Rdata数据文件可直接载入运行有助于读者边学边练理解污染物时空分布的分析方法。 很多人下载了公开的空气质量数据之后第一反应就是plot()、ggplot()上去猛画几张图然后盯着图发呆——图是真的好看但下一步该干什么数据里到底藏着什么信息却说不上来。这个大气污染数据分析案例就是冲着这个痛点去的。它不是一个“点一下出结果”的傻瓜式模板而是一套从数据体检、清洗、可视化到建模的完整R语言分析流程用的是一份真实结构的国控站点监测数据含PM2.5、PM10、SO2、NO2、CO、O3六项污染物和气象要素。无论你是刚学R的学生、做环境相关研究的科研人员还是想转行数据分析的职场人都能从这套案例里拿到可以直接改改就跑的代码以及比代码更重要的——分析思路。1. 拿到污染物数据后先别急着画图我在处理这份数据时踩过最大的坑就是“想一口气把图全画完”。后来发现急着画图的代价是要花三倍时间回头改数据格式。所以这套案例的第一步刻意安排成了数据体检而且是用最朴素的方式。1.1 数据结构和字段含义案例用的数据是CSV格式包含日期、站点、六项污染物浓度、四项气象要素。字段大致长这样# 用data.table读取速度更快 library(data.table) air_data - fread(air_data.csv) str(air_data) # 输出示例 # $ date : chr 2023-01-01 2023-01-02 ... # $ station : chr 站点A 站点A ... # $ pm25 : num 45.2 78.3 ... # $ pm10 : num 80.1 112.5 ... # $ so2 : num 8.2 10.1 ... # $ no2 : num 32.4 45.2 ... # $ co : num 0.8 1.2 ... # $ o3 : num 62.3 55.8 ... # $ temperature : num 3.2 1.5 ... # $ humidity : num 55 68 ... # $ windspeed : num 2.1 1.4 ... # $ wind_direction: num 180 220 ... # $ pressure : num 1023 1021 ...先看结构不是走形式。这一步能暴露大多数问题日期是不是字符型、有没有把站点读成因子、污染物浓度是不是数值型。如果哪列被读成了字符后面所有计算都会报错或者得出荒谬结果。我曾经遇到过一次co列因为有ND未检出标记被整体读成字符型后面算相关系数时直接全部NA排查了很久。1.2 缺失值和异常值检查大气污染监测数据有个特点不是每天每个站点都能出数。设备校准、通讯故障、停电都会导致缺失。更麻烦的是有些缺失会以NA形式存在有些却会以-9999、0甚至空格存在。# 缺失值概览 colSums(is.na(air_data)) # 更严格一点查看每列的“异常”情况 summary(air_data[, .(pm25, pm10, so2, no2, co, o3)])我的习惯是重点看summary()里的Min.和Max.。PM2.5出现0值要警惕因为很多仪器在浓度极低时会报0但更常见的情况是缺测被填成0SO2出现负值也可能是仪器零点漂移。对于后续要算均值、做相关分析的任务这些异常值会实实在在影响结果。这个案例里我的处理原则是先boxplot画出分布再用IQR规则标记极端值但在删之前一定先去看看原始记录不能机械地“非黑即白”。数据体检不是炫技是在帮后面的所有分析排雷。2. 数据清洗日期解析、缺失值和合并气象数据数据清洗是整个分析流程里最不性感但最重要的一步。大气污染数据的清洗有几个固定的坑这套案例里全部处理了一遍。2.1 日期解析别让时间轴变成乱麻现实里的日期格式五花八门2023/1/1、2023-01-01、01-01-2023以及Excel导出的45000这种序列值。如果时间轴没有正确解析时间序列图上的横坐标会乱成一团。案例里用lubridate包统一解析library(lubridate) air_data[, date : ymd(date)] # 或者按年月日拆开方便后面做月季分析 air_data[, :( year year(date), month month(date), season case_when( month %in% c(3,4,5) ~ 春季, month %in% c(6,7,8) ~ 夏季, month %in% c(9,10,11) ~ 秋季, TRUE ~ 冬季 ) )]这里有个细节ymd()会尝试自动识别格式但如果你的数据里有2023/1/1这种格式建议先统一成干净格式再解析。跨年数据更要小心我曾经因为日期列里混入了一个2023-02-30导致整个时间序列的排序错位画出来的线图在2月底突然断崖下跌花了很久才发现是脏数据。2.2 缺失值的处理策略区分“真缺”和“假缺”对于污染物时间序列缺失值处理我很少用粗暴的均值填充。原因很简单污染物浓度有明显的日变化和季节变化1月凌晨的PM2.5和7月下午的PM2.5差了十倍不止用一个全域均值填进去等于硬造虚假规律。案例里采用了两步法先删除关键污染物的缺失行如果缺失比例很低比如小于5%然后用线性插值填补剩余缺口。# 计算各列的缺失比例 missing_ratio - colMeans(is.na(air_data)) %% round(3) print(missing_ratio) # 对pm25用na.approx做线性插值 library(zoo) air_data[order(date), pm25_filled : na.approx(pm25, na.rm FALSE)]注意na.approx处理中间缺口没问题但首尾的NA它填不了需要单独用na.locf或者最近邻补齐。实务中如果首尾长期缺测我最倾向直接删掉那段数据避免人为制造“假数值”。2.3 气象数据合并把污染物放到环境下看大气污染不是孤立存在的气象条件直接决定了污染物的积累和扩散。风速大时PM2.5浓度通常低湿度高时往往利于二次颗粒物生成。案例把污染物数据和气象数据按日期合并成一张宽表为后面的相关性分析和建模做准备met_data - fread(met_data.csv) met_data[, date : ymd(date)] combined - merge(air_data, met_data, by c(date, station), all.x TRUE)我一直建议把“合并”当成一次验证的机会merge之后立刻检查行数有没有变多那是产生了笛卡尔积或变少那是丢失了匹配项。案例数据里两个文件日期是严格对齐的但你自己处理新数据时这一步最容易出问题。3. 时空分布可视化一张图讲清污染格局可视化不是用来“好看”的是用来快速定位问题的。大气污染数据最常见也最有价值的三个可视化方向这套案例里都写了代码。3.1 浓度时间序列先看整体趋势再看异常峰值折线图是第一张必须画的图。看什么看趋势、看季节、看极值。案例用ggplot2画了一张全年逐日变化曲线并用geom_smooth()叠加了平滑趋势线一眼就能看出秋冬季浓度明显高于夏季library(ggplot2) ggplot(combined, aes(x date, y pm25)) geom_line(colour grey60, linewidth 0.4, alpha 0.6) geom_smooth(method loess, colour #D73027, se FALSE, linewidth 1) labs(title PM2.5 日平均浓度变化趋势, x 日期, y 浓度 (μg/m³)) theme_classic(base_size 14)画完这张图我往往会顺手标出最高峰值对应的日期——那大概率对应某次重污染过程值得回去翻翻当时的天气和边界层条件。这种“图→现象→原因”的链条才是数据分析的价值所在。3.2 站点对比与箱线图不同站点的污染画像如果是多站点数据箱线图比折线图更能体现站点之间的差异。案例里按站点分组画了PM2.5的箱线图能清楚看到哪个站点的中位数高、哪个站点的离群值多ggplot(combined, aes(x station, y pm25, fill station)) geom_boxplot(outlier.colour black, outlier.size 0.8, alpha 0.7) labs(title 各站点 PM2.5 浓度分布对比, x 站点, y 浓度 (μg/m³)) theme_classic(base_size 14) theme(legend.position none)很多分析到此就停了其实还可以继续往下走把站点按区域类型工业区、居民区、交通点着上不同颜色往往能看出更有意思的规律。这个案例数据带了站点分类字段代码直接在原基础上加了一行fill station_type。3.3 污染物相关矩阵先猜后证不要盲目上模型做相关分析之前先画一张相关矩阵图能帮自己建立“数据直觉”。比如PM2.5和PM10因为同源相关性必然高CO和NO2都来自燃烧排放也高度相关。案例里用corrplot画了六项污染物加气象因子的相关矩阵library(corrplot) cor_matrix - combined[, .(pm25, pm10, so2, no2, co, o3, temperature, humidity, windspeed)] | cor(use complete.obs) corrplot(cor_matrix, method color, type upper, tl.col black, addCoef.col black, number.cex 0.7)这张图的价值在于后续建模选特征时你会很清楚哪些变量存在共线性哪些变量是相对独立的。比如PM2.5和PM10相关系数如果超过0.9那建模时就不该把它们同时放进去做普通线性回归否则会带来多重共线性的麻烦。4. 影响PM2.5浓度的气象因子量化分析大气污染数据分析如果只停在“画图”层面说服力还不够。气象因子到底怎么影响污染物浓度这个案例用回归和相关分析做了量化拆解。4.1 单因子相关分析哪个气象因素影响最大先做单因子检验算出PM2.5与温度、湿度、风速、气压的相关系数和p值with(combined, cor.test(pm25, temperature, method spearman)) with(combined, cor.test(pm25, humidity, method spearman)) with(combined, cor.test(pm25, windspeed, method spearman)) with(combined, cor.test(pm25, pressure, method spearman))这里我特意用了spearman而不是pearson因为PM2.5浓度和气象要素的关系通常不是线性的比如风速对PM2.5的“清除效应”在低风速段明显到高风速段趋于平缓。用秩相关更稳健。实测下来风速和PM2.5往往呈显著负相关湿度则常常呈正相关高湿促进二次颗粒物生成温度的影响在冬季和夏季方向是相反的所以做全样本相关时会被抵消掉这也是为什么后面需要再分层看。4.2 回归模型用一个方程描述多重影响单因子相关只能看“两两关系”实际大气过程是多因子共同作用的。案例用多元线性回归量化各气象因子的贡献model_lm - lm(pm25 ~ temperature humidity windspeed pressure, data combined) summary(model_lm)重点看三样东西Estimate的正负方向是否符合物理常识、Pr(|t|)显著的变量有哪些、Adjusted R-squared整体解释力如何。这个案例运行下来windspeed和pressure的系数通常显著为负humidity系数显著为正说明静稳、高湿、高压的气象条件确实利于PM2.5累积。回归不是为了得到一个完美的预测公式而是为了把“天气影响污染”这个模糊感觉变成一个可量化的系数。4.3 分层分析分季节看效应全样本拟合完之后我强烈建议做一层分季节的亚组分析。因为夏季PM2.5主要受O3和光化学反应影响冬季则受燃煤排放和静稳天气影响气象因子的作用方向和强度完全不同。案例里用split按季节拆开再分别跑回归season_models - lapply(split(combined, combined$season), function(df) { lm(pm25 ~ temperature humidity windspeed pressure, data df) }) lapply(season_models, summary)我见过最典型的情况是风速在全年回归里负效应明显但到夏季变得不显著因为夏季本身的扩散条件整体较好风速的边际效应被削弱了。这种分层结果写进报告里比只报一个全样本模型要专业得多。5. 预测模型的构建从线性回归到随机森林分析的最后一步是建模。模型不是越复杂越好而是越匹配问题越好。案例里我特意安排了三个层次的模型方便不同需求的读者对号入座。5.1 线性回归基准模型跑通一个可解释的底稿用污染物和气象变量预测PM2.5浓度线性回归作为基准模型有不可替代的价值可解释性强、计算快、稳定性好。代码上需要把类别变量转成哑变量并检查VIF方差膨胀因子确认没有严重的共线性library(car) model_lm_full - lm(pm25 ~ pm10 so2 no2 co temperature humidity windspeed pressure, data combined) vif(model_lm_full) # 如果VIF 10说明存在严重共线性需要剔除变量VIF的判断很重要。我在这个案例里遇到PM2.5和PM10的VIF超过20的情况后来把PM10从特征里去掉模型的稳定性立刻好了起来。这是很多初学者最容易忽视的环节。5.2 随机森林抓非线性关系随机森林的优势是不用预设函数形式能自动捕捉非线性关系和变量交互。案例用randomForest包训练了一个包含500棵树的模型library(randomForest) # 去除含NA的行并选择建模需要的列 model_data - combined[complete.cases(combined), .(pm25, so2, no2, co, o3, temperature, humidity, windspeed, pressure)] set.seed(42) rf_model - randomForest(pm25 ~ ., data model_data, ntree 500, importance TRUE) print(rf_model)随机森林跑起来不需要太多调参ntree500基本上够用关键是要看变量重要性图varImpPlot()。这个案例里排在前面的通常是pressure、humidity和so2说明静稳天气条件和燃煤源对PM2.5浓度的影响很突出。5.3 模型评估不看R²只看RMSE和泛化能力训练集上的R²再高也不代表真实预测能力。案例里用caret的createDataPartition做了七三拆分在测试集上比较不同模型的RMSE和R²library(caret) set.seed(123) train_index - createDataPartition(model_data$pm25, p 0.7, list FALSE) train_data - model_data[train_index, ] test_data - model_data[-train_index, ] # 线性模型在测试集上的评估 lm_fit - lm(pm25 ~ ., data train_data) lm_pred - predict(lm_fit, newdata test_data) lm_rmse - sqrt(mean((test_data$pm25 - lm_pred)^2)) lm_r2 - cor(test_data$pm25, lm_pred)^2 # 随机森林在测试集上的评估 rf_fit - randomForest(pm25 ~ ., data train_data, ntree 500) rf_pred - predict(rf_fit, newdata test_data) rf_rmse - sqrt(mean((test_data$pm25 - rf_pred)^2)) rf_r2 - cor(test_data$pm25, rf_pred)^2 cat(线性回归 RMSE:, lm_rmse, R²:, lm_r2, \n) cat(随机森林 RMSE:, rf_rmse, R²:, rf_r2, \n)这个案例实测下来随机森林比线性回归的RMSE低10%-20%但优势没有想象中夸张。这说明PM2.5浓度里很大一部分方差是被排放源决定的气象因子只能解释其中一部分这也符合大气科学的常识。另外如果数据里含有站点固定效应建议加一个station变量进来或者用混合效应模型能进一步提升预测精度。提示建模过程中set.seed()一定不要省。随机森林和训练集划分都涉及随机数不设种子的话你每次跑的结果都会不一样复现性就没法保证了。一点个人体会这套案例做下来我最大的收获不是哪个模型分数高而是整个分析链条的顺畅感先体检数据再清洗整合然后可视化找模式和异常接着用相关和回归量化关系最后建模验证。每一步都在为下一步铺路没有一步是白做的。你在跑这份代码的时候如果遇到报错优先检查数据格式和缺失值八成的问题都出在数据上。做数据分析这个事情慢就是快前期把数据弄干净后面所有步骤都会很顺。本文还有配套的精品资源点击获取