R语言回归分析实战:预测首尔自行车共享需求

发布时间:2026/9/25 21:15:22
R语言回归分析实战:预测首尔自行车共享需求 简介面向需要在R环境中完成回归建模与需求预测的数据分析学习者这是一份首尔自行车共享需求预测完整项目资源。资源围绕天气、时间、假期、季节等多种因素对每小时租车量的影响展开提供从数据探索、变量重要性分析到CUBIST、随机森林、CART、KNN、条件推断树等多种模型构建与评估的R代码并配套报告文档便于理解规则集成模型如何实现约95%的R²解释力。包体共4个文件含R脚本、CSV数据集、Markdown说明及Word分析报告压缩包大小1.71MB结构精简适合直接运行复现。目前已有256人学习。通过该资源可获取完整R代码、原始CSV数据集及分析报告适合课程设计、案例复现或作为回归预测入门参考能帮助快速掌握模型对比、交叉验证与变量重要性评估的实操流程。1. 首尔自行车共享需求回归从“每小时要备多少辆车”说起首尔自行车共享需求数据是拿 R 做回归分析入门到实战之间最好的中间站——它有公开数据集、有时间结构、有天气和日历混杂出的真实噪声但又不至于像工业场景那样脏到没法收拾。标题里的活儿其实就是一件事用 R 读入按小时记录的数据把气温、湿度、风速、降雨、节假日这些特征组织起来训练一个能预测下一个小时需要多少辆自行车的回归模型。适合两类人一类刚学完线性模型想找个完整项目练手另一类在做共享出行调度、用户消费预测这类计数回归场景想看看经典回归在真实业务数据上能承受到什么程度。先说一个反直觉的结论这类按小时记录的数据如果像普通回归那样随机划分训练集和验证集预测误差会非常好看但一上线马上翻车。原因是同一天的不同小时会被拆到两边模型等于提前“见过了当天的天气”。要测真实水平后面必须按时间拆。2. 数据长什么样理解字段与建模前处理这一步决定回归天花板2.1 字段清单与目标变量的“代理”性质先看字段。这个数据集的记录跨度约一年每行是一条小时级记录目标列是当小时实际租出的自行车数。严格说真正的“需求”不可观测租出量只是需求的代理下雨天可能有人想骑但没人出门租出量就会低估需求。业务解读时要留这个心眼建模时则把租出量当目标即可。字段含义类型建模注意Date日期字符转成 Date 类型拆出周几、月份Rented Bike Count小时租出量数值目标变量右偏带大量 0Hour小时 0–23数值核心时间特征对需求量影响最大Temperature气温 ℃数值与需求呈倒 U 型考虑平方项Humidity湿度 %数值高湿常伴降雨但信息不重复Wind speed风速 m/s数值强风明显压制骑行需求Visibility能见度数值单位是 10m不是米建模前换算Dew point temperature露点温度数值与湿度相关性高注意共线性Solar Radiation太阳辐射 MJ/m²数值白天信号强与 Hour 有交互Rainfall降雨量 mm数值下雨需求骤降0 值有业务含义Snowfall降雪量 cm数值冬季气候和温度强相关Seasons季节类别转成因子不要保留为文本Holiday是否节假日类别转成 0/1Functioning Day是否运营日类别非运营日租出量恒为 0要单独处理这里最容易忽略的是 Visibility 的单位。原始数据里能见度出现几百甚至上千的数值直接当米解释会以为数据异常实际上单位是 10m1000 代表 10000 米。另一个是 Functioning Day它表示当天系统是否运营等于 No 的日子租出量基本是 0这部分样本混进回归会把均值往下拉后面专门讲怎么处理。2.2 用 dplyr 完成清洗与时间特征构造先读数据再看实际列名。不同渠道下载的首尔自行车数据列名略有差异有的带空格、有的带括号和单位读进来之后 R 会自动做 check.names 处理。我一般先打印一遍名字再统一改成短列名省得后面写 formula 时到处加反引号。library(dplyr) bike_raw - read.csv(SeoulBikeData.csv, fileEncoding UTF-8) cat(names(bike_raw), sep \n) # 观察实际列名按自己那份数据来对 # 统一转小写把空格、括号、斜杠替换成下划线便于后续引用 names(bike_raw) - names(bike_raw) %% tolower() %% gsub([^a-z0-9], _, .) %% gsub(_, _, .) %% sub(_$, , .) bike - bike_raw %% rename( rented rented_bike_count, temp temperature, hum humidity, wind wind_speed, vis visibility_10m, dew dew_point_temperature, solar solar_radiation_mj_m2, rainfall rainfall_mm, snowfall snowfall_cm )fileEncodingUTF-8是 Windows 上常见的坑不指定的话某些镜像下载的 CSV 会用系统本地编码读取中文标签直接乱码。正则清洗列名时注意gsub([^a-z0-9], _, .)会把所有非字母数字字符统一压成下划线双下划线再合并一次最后去掉末尾的下划线。rename 里的新列名是后面所有代码的基础如果某份数据的列名对不上以cat(names())打出来的结果为准改映射。接着构造时间特征和业务哑变量。日期格式在源数据里是日/月/年as.Date的 format 必须写%d/%m/%Y顺手用 lubridate 的wday判断周末避免依赖系统 locale 的weekdays()。library(lubridate) bike - bike %% mutate( date as.Date(date, format %d/%m/%Y), hour as.numeric(hour), month as.numeric(format(date, %m)), is_weekend ifelse(wday(date, week_start 1) %in% c(6, 7), 1, 0), raining ifelse(rainfall 0, 1, 0), snowing ifelse(snowfall 0, 1, 0), holiday ifelse(holiday Yes, 1, 0), functioning ifelse(functioning_day Yes, 1, 0) )wday(date, week_start 1)返回 1–71 是周一6、7 是周六日这样判断与系统语言无关。holiday和functioning从字符转成 0/1是因为后面放进lm()和glmnet时数值哑变量比字符因子更直观也避免因子水平顺序带来的模型解释困惑。2.3 检查缺失值与分布顺带拆一次训练/验证清洗完先看 summary重点不是均值而是有没有不可能出现的值、缺失值有没有被 R 读成 NA、目标变量右偏到什么程度。小时租出量的典型分布是凌晨大量 0早晚高峰两个峰整体右偏。回归模型对右偏目标往往拟合不好第 3 章的对数变换就是为这个做准备。summary(bike) hist(bike$rented, breaks 50, main 小时租出量分布) # 非运营日的租出量到底长什么样 bike %% filter(functioning 0) %% summarise(n n(), mean_rented mean(rented), max_rented max(rented))非运营日样本的租出量均值接近 0这符合业务定义。处理策略我一般看建模目的如果目标是预测运营状态下的需求直接把非运营日过滤掉模型更干净如果目标是端到端预测“系统实际要承担多少租出量”则保留 functioning 作为特征预测时输入运营状态即可。两种做法没有绝对对错但要在同一份代码里保持一致不要建模时又混进又剔除。拆训练/验证这一步强烈建议在清洗之后就做不要等模型调完再拆。关键点是按日期排序后切分而不是随机抽样dates - sort(unique(bike$date)) cut_idx - round(length(dates) * 0.8) train - bike %% filter(date dates[cut_idx]) test - bike %% filter(date dates[cut_idx])这里 0.8 是常见比例但对共享单车这类强周期业务我更推荐按“最后 30 天做测试”来切因为业务上线时预测的永远是未来不是过去。切分完成后后面所有特征工程和模型调参都只允许看 traintest 要一直留到最后一章做外推验证。3. 建立基线回归用 lm() 跑通第一版靠残差决定下一步3.1 特征怎么进模型数值型、哑变量、交互项线性回归的第一步不是把字段全塞进 formula而是想清楚每个特征的形态。温度对骑行需求是典型的倒 U 型太冷不出门太热也不骑所以加一个I(temp^2)平方项。Hour 本身是 0–23 的整数直接放进线性模型会被当成“越晚需求线性越低”实际上是早晚双峰这里先用线性近似后面换树模型时自动处理非线性。降雨和降雪已经转成 0/1 哑变量保留原始降雨量数值会引入大量 0 值干扰。交互项我习惯先加两个业务上确定的temp:hour表示不同时段对温度的敏感度不同raining:is_weekend表示雨天在周末的影响可能比工作日更弱。交互项不要一上来加一堆否则 VIF 立刻爆表自己都分不清是共线性还是真信号。library(car) fit_lm - lm( rented ~ temp I(temp^2) hum wind vis solar raining snowing hour is_weekend holiday functioning temp:hour raining:is_weekend, data train ) summary(fit_lm) vif(fit_lm)vif()来自 car 包输出大于 5 就要警惕大于 10 基本可以断定这组变量在抢同一个信息。首尔数据里最容易爆的是露点温度和湿度它俩物理含义高度重叠如果 VIF 高我一般先去掉 dew保留 hum因为湿度对骑行决策的解释更直接。3.2 第一版线性模型与残差诊断summary()里先看 F 统计量的 p 值再看每个变量的显著性。第一次跑这个模型时你会发现温度、湿度、小时、降雨系数全显著但 Adjusted R² 可能只在 0.3–0.4 徘徊。别急着调参先画残差图。par(mfrow c(2, 2)) plot(fit_lm)四张图按顺序读Residuals vs Fitted 看有没有喇叭口Normal Q-Q 看残差是否正态Scale-Location 看方差的稳定性Residuals vs Leverage 找影响点。这个数据集上最常见的形态是拟合值增大时残差散开成漏斗形QQ 图两端明显偏离直线。这是计数数据的典型异方差——租出量均值越高波动越大线性模型强行用同一个方差去拟合导致小预测值被高估、大预测值被低估。残差图是黑匣子但它告诉你下一步该怎么走。两个方向一是对目标变量做对数变换让方差更稳二是换用泊松回归或树模型。泊松回归在概念上更贴合计数数据但首尔这份数据存在过度离散直接glm(..., family poisson)会让标准误被低估。所以先把对数变换跑通再上正则化和树模型。3.3 对数变换后重跑处理极端值log1p是log(x 1)的简写专门处理含 0 的计数数据。直接log(0)会得到负无穷log1p把 0 映射成 0语义上说得通深夜无人租车对数需求就是 0。train$log_rented - log1p(train$rented) fit_lm_log - lm( log_rented ~ temp I(temp^2) hum wind vis solar raining snowing hour is_weekend holiday functioning temp:hour raining:is_weekend, data train ) summary(fit_lm_log) hist(resid(fit_lm_log), breaks 50)对比两次 summary你会看到对数模型的 R² 明显上升残差直方图也更接近正态。但这不意味着模型真的更好只是目标尺度变了评估指标也跟着变。后续比较模型时一定要把预测值还原到原始租出量尺度再算 RMSE 或 MAE否则两个模型根本不在一个比较维度上。还原时的常见错误是直接exp(predict(...)) - 1。对数变换让模型拟合的是对数空间的中位数还原后整体会系统性偏低。这个偏差在业务上很要命第 5 章专门讲怎么用 smearing 修正。4. 把预测误差压下去glmnet 弹性网络与随机森林怎么选4.1 什么时候该从 lm 切到正则化模型线性模型跑通之后下一步不是立刻上深度学习而是先用正则化回归和树模型把误差压一压。首尔数据里有几个特征天生高度相关露点温度与湿度、太阳辐射与气温、风速与能见度。lm 对这组相关特征非常敏感系数估计方差大换个训练集系数就漂移。glmnet 的弹性网络elasticnet同时混入 L1 和 L2 惩罚L1 帮忙做变量选择L2 帮忙稳定相关变量的系数正好治这个病。在 R 里 glmnet 的alpha参数控制 L1/L2 比例alpha 1是 lassoalpha 0是 ridge中间值就是弹性网络。默认alpha 1很多人就这么用但对相关特征较多的数据纯 lasso 会随机挑一个代表变量丢掉另一个稳定性差。我习惯把alpha从 0 到 1 扫一遍让cv.glmnet自己选。4.2 用 cv.glmnet 扫 alpha 并选 lambda 的完整命令glmnet的接口和lm最大的区别是它不接受 formula必须自己用model.matrix构造特征矩阵。这一步有两个坑一是因子列会自动展开成哑变量二是展开结果自带截距列喂给 glmnet 前要删掉否则会跟 glmnet 内部加的截距重复。library(glmnet) # 用 train 数据构造特征矩阵-1 去掉截距列 x_train - model.matrix( ~ temp I(temp^2) hum wind vis solar raining snowing hour is_weekend holiday functioning temp:hour raining:is_weekend, data train )[, -1] y_train - train$log_rented # 扫 alpha对每个 alpha 做 5 折交叉验证选 lambda alphas - seq(0, 1, by 0.1) cv_list - lapply(alphas, function(a) { cv.glmnet(x_train, y_train, alpha a, nfolds 5) }) # 取交叉验证误差最小的那一组 best_idx - which.min(sapply(cv_list, function(m) min(m$cvm))) alpha_best - alphas[best_idx] fit_en - cv_list[[best_idx]] plot(fit_en) coef(fit_en, s fit_en$lambda.min)alphas按 0.1 步长扫 11 个值对这份数据量完全够用不用更细。cv.glmnet的nfolds 5是常见默认数据量小可以改 10但这数据有近一年的小时记录5 折足够稳定。lambda.min是交叉验证误差最小的 lambdalambda.1se是误差在一个标准误之内的最大 lambda系数更稀疏。我一般两个都看一眼上线求稳用lambda.1se竞赛刷指标用lambda.min。plot(fit_en)画的是不同 lambda 下系数收缩路径能直观看到哪些变量最后被压成 0。coef(fit_en, s fit_en$lambda.min)输出选定 lambda 下的系数如果某个特征系数是 0说明它对预测没有边际贡献。4.3 随机森林补位做三模型横评正则化回归解决了共线性但解决不了特征的非线性和交互。Hour 的早晚双峰、温度倒 U 型这些 lm 要用平方项和交互项手工拼而树模型天生能切分这种关系。R 里我优先用 ranger 而不是 randomForest因为数据量稍大时 randomForest 慢到让人怀疑人生ranger 在同等精度下快一个量级。library(ranger) fit_rf - ranger( log_rented ~ temp I(temp^2) hum wind vis solar raining snowing hour is_weekend holiday functioning temp:hour raining:is_weekend, data train, num.trees 500, mtry 5, importance impurity, seed 42 ) print(fit_rf) sort(fit_rf$variable.importance, decreasing TRUE)num.trees 500对这份数据足够再高耗时增加明显但精度提升有限。mtry 5是每个分裂节点随机抽 5 个特征默认是特征数的平方根这里特征展开后约 20 个默认 4–5 都合理。importance impurity输出的是 Gini 重要性计算快适合看特征排序permutation更准但慢最后定稿时可以用。三模型横评的骨架是这样# 测试集同样构造特征矩阵 x_test - model.matrix( ~ temp I(temp^2) hum wind vis solar raining snowing hour is_weekend holiday functioning temp:hour raining:is_weekend, data test )[, -1] y_test - test$log_rented pred_lm - predict(fit_lm_log, newdata test) pred_en - predict(fit_en, s fit_en$lambda.min, newx x_test) pred_rf - predict(fit_rf, data test)$predictions calc_rmse - function(actual, pred) sqrt(mean((actual - pred)^2)) calc_mae - function(actual, pred) mean(abs(actual - pred)) # 在 log 尺度上横评还原到原始尺度的偏差修正见第 5 章 sapply(list(lm pred_lm, enet pred_en, rf pred_rf), function(p) { c(rmse calc_rmse(y_test, p), mae calc_mae(y_test, p)) })经验上在这类数据里随机森林的 RMSE 通常明显低于线性模型弹性网络介于两者之间但弹性网络的优势是可解释性和稳定性——特征只有几十个系数还能打印出来讲给业务听。树模型一旦跑起来就是纯黑匣子业务方问“为什么这个小时预测这么高”你很难一句话讲清楚。5. 避坑首尔自行车需求回归里最常见的 5 个翻车点5.1 随机抽样导致预测虚高同一天的天气被模型提前看到现象训练/验证用sample()随机切分验证集 RMSE 低得喜人发布上线后每天预测都偏业务方直接质疑模型没学过当天数据。原因小时级数据里同一天的 24 小时共享同一组天气和节假日状态。随机切分会把同一天的一部分小时分进训练集另一部分分进验证集模型等于先看了当天的天气答案再去预测当天剩下的小时信息泄露非常隐蔽。解决切分必须按日期整体切。用第 2 章的sort(unique(date))方法保证训练集日期全部早于验证集日期。头条原则任何涉及时间序列的数据测试集日期必须全部在训练集之后这是回归预测类项目的第一纪律。5.2 非运营日混进回归早上 8 点的预测值被系统性拉低现象模型整体 R² 还行但分时段看误差时工作日早高峰的预测值普遍偏低而且偏低幅度稳定。原因Functioning Day 等于 No 的日子租出量为 0这些样本占比不小把回归均值往下拉。模型学到的是“存在一整天完全没人骑车的日子”但它不知道是哪天于是把这种不确定性摊到所有样本上高峰期预测被拖低。解决按业务目标二选一。只做运营日需求预测就在清洗时filter(functioning 1)做端到端租出量预测就保留 functioning 特征预测时显式传入运营状态。最忌讳的是建模时不处理、解释时又说“非运营日本来就是 0”两头不靠。5.3 用 exp() 还原对数预测总量预测出现系统性低估现象把 log 模型的预测值exp()还原后分小时误差看着还行但按天聚合的总量预测比实际少一到两成。原因log1p建模后残差在对数空间是零均值但exp()是非线性变换还原后残差的均值不再为 0。exp(predict)得到的是对数空间的中位数不是原始尺度的均值。业务要的是总量中位数一定偏小。解决加一个 smearing 修正项用训练集残差计算mean(exp(resid))乘到预测值上再减 1smear - mean(exp(resid(fit_lm_log))) pred_original - exp(predict(fit_lm_log, newdata test)) * smear - 1这个修正项的意义是训练集上模型平均低估的比例就是未来预测要补偿的比例。代价是整体方差被放大一点但总量预测的偏差能明显收回来。5.4 星期几特征依赖系统语言换台电脑结果就变现象同样的代码在 A 机器跑is_weekend判断正常在 B 机器跑出来的周末标记全反了模型结果对不上。原因基础 R 的weekdays()返回星期几的文本这个文本由系统 locale 决定。中文系统返回“星期六”英文系统返回“Saturday”如果代码里写weekdays(date) %in% c(Saturday, Sunday)在中文机器上永远是 FALSE。解决不要用weekdays()做判断改用lubridate::wday(date, week_start 1)它返回数字 1–7与语言无关。或者用format(date, %u)也能拿到 ISO 周几数字。这类与环境相关的问题最坑人因为它不在你的代码逻辑里而在运行环境的配置里。5.5 把能见度原始数值当米解释异常值其实是单位问题现象做异常值清洗时发现 Visibility 有大量超过 1000 的“异常值”准备一刀切删除结果模型训练集和验证集分布对不上。原因这个数据集里 Visibility 的单位是 10m原始值 1000 表示 10000 米。把它当米处理会以为晴天的能见度上万是数据错误其实不仅正常而且是有用的天气信号。雨雪天的能见度原始值会掉到几百甚至以下这种差异对预测很有区分度。解决读数据后先确认字段单位能见度转成公里再存vis_km vis / 100。删除异常值前先看单位定义别让单位问题毁掉一个真正有效的特征。这类坑在公开数据集里尤其常见下载页的说明文档和数据字典值得先读一遍再动手。6. 从模型到结论外推验证与特征效应解读6.1 按时间切分的验证与滚动窗口随机切分在第 5 章已经否掉了但一次性按时间切也还不够。首尔这个城市有明显的季节周期春秋骑行量大、冬夏差很多只切一刀可能训练集里没有见过足够多的冬天样本验证集又刚好落在冬天误差会被高估。更稳的做法是滚动窗口模拟真实上线时“每周重训一次”的节奏。dates - sort(unique(bike$date)) # 以 30 天为窗口每天滚动一天 for (i in 1:(length(dates) - 30)) { trn - bike %% filter(date dates[i], date dates[i 30]) tst - bike %% filter(date dates[i 30]) # 这里放你的模型训练代码记录当天预测误差到结果表 # err_day[i] - calc_rmse(tst$rented, pred) }这个循环跑起来会重训很多次模型计算量大但得到的是模型在一年各个季节里真实的外推表现。跑完后按月份聚合误差如果 1 月和 7 月的误差明显高于春秋说明模型还没学到温度极值对骑行需求的非线性压制这就是下一步要加的交互项或新特征。6.2 用特征效应而不是变量重要性解释模型树模型的变量重要性只告诉你哪些特征参与分裂多不告诉你怎么影响预测。Hour 重要性最高是因为它把早晚高峰切开了但你从重要性数字里看不出“早上 8 点需求爬升、下午 6 点另一个峰”。要把结论讲给业务听得看偏依赖图。library(pdp) pd_hour - partial(fit_rf, pred.var hour, train train) plot(pd_hour, type b)partial计算的是固定其他特征为均值时Hour 变化带来的预测均值变化。画出来能看到清晰的早晚双峰。这个图就是调度业务最需要的那个结论早高峰备车、午间低谷、晚高峰再备一次。数据现象模型表现调度动作工作日早 7–9 点需求峰Hour 偏依赖图第一个峰早高峰前 1 小时从仓库调车到地铁口降雨量 0 时需求骤降raining 系数为负且与周末交互显著雨天减少路侧调度转入车辆检修气温超过 30℃ 后需求回落temp 平方项负系数高温天减少露天站点车辆避免暴晒损耗最后说一个我自己的习惯模型定稿前强制自己看一眼按天聚合的误差图。如果误差和气温、节假日出现清晰模式说明还有信息没进模型如果误差平稳地绕零波动才敢往上线走。这套流程在首尔自行车数据上成立换到用户消费预测、销量预测这些同样带时间结构的场景也一样。希望帮到你。本文还有配套的精品资源点击获取