企业官网建设流程全解析

简介：这套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，但更常见的情况是缺测被填成0；SO2出现负值也可能是仪器零点漂移。对于后续要算均值、做相关分析的任务，这些异常值会实实在在影响结果。这个案例里，我的处理原则是：先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)

随机森林跑起来不需要太多调参，ntree=500基本上够用，关键是要看变量重要性图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()一定不要省。随机森林和训练集划分都涉及随机数，不设种子的话，你每次跑的结果都会不一样，复现性就没法保证了。

一点个人体会

这套案例做下来，我最大的收获不是哪个模型分数高，而是整个分析链条的顺畅感：先体检数据，再清洗整合，然后可视化找模式和异常，接着用相关和回归量化关系，最后建模验证。每一步都在为下一步铺路，没有一步是白做的。你在跑这份代码的时候如果遇到报错，优先检查数据格式和缺失值，八成的问题都出在数据上。做数据分析这个事情，慢就是快，前期把数据弄干净，后面所有步骤都会很顺。

本文还有配套的精品资源，点击获取