首尔自行车共享需求预测：R语言特征工程与多模型对比实战

简介：面向城市共享单车运营与数据分析场景，这份资源提供基于首尔自行车共享需求数据集的回归建模完整方案，适合数据科学初学者和需要掌握预测建模流程的分析人员。资源围绕每小时自行车租赁量预测，综合运用CUBIST、正则化随机森林、CART、K近邻与条件推断树等模型，并给出重复交叉验证、多指标评估及变量重要性分析思路，可直接迁移到同类时序或回归任务中。包内共4个文件，包括R脚本、项目说明文档、数据报告docx和原始CSV数据，压缩包大小约1.71MB，结构简洁。R脚本可实现从数据清洗到模型训练评估的完整流程，docx报告汇总实验结果与结论，md和csv便于快速理解数据与复现。目前已有256人学习下载。对想系统梳理回归建模流程、学习CUBIST等树模型原理或获取可复现分析代码的读者，这份资料具有较高参考价值，能帮助节省数据整理和调参时间。

1. 首尔自行车共享需求预测：这份数据自带一颗“后悔药”

第一次拿到首尔自行车共享数据集，我差点被“回归分析”四个字带偏，上来就敢跑线性回归。结果可想而知：R² 惨不忍睹，残差图歪得像心电图。后来翻完这份资源里的FinalProject.R和Data Analysis Report.docx才明白，真正能打的不是单模型，而是“特征工程 + 多模型对比 + 重复交叉验证”这套组合拳。里面用到的 CUBIST 模型能把测试集 R² 干到 0.95 左右，而温度和一天中的小时才是预测每小时租车量的最关键变量——这个结论和参考论文完全一致。这份资源适合正在做课程设计、毕设或者城市交通预测练手的从业者，你既想看数据长什么样，也想找一个能直接复现的实验框架，那这篇笔记值得你读完。

2. 先搞懂字段再动模型：读入、体检与时间拆解

2.1 字段语义与单位速查：少踩九个坑的第一关

拿到SeoulBikeData.csv后，我做的第一件事不是summary()，而是把每列的单位和类型抄到纸上。越往后你越会发现，这堆字段里藏着至少九个坑：能见度的单位是 10m 而不是米，露点温度与温度几乎共线，节假日只有两个水平但分布极不均衡。

2.2 用 R 读入并做初步“体检”

先把SeoulBikeData.csv放进工作目录，然后按下面这套顺序跑。我习惯在开头就把 R 版本和包版本打印出来，方便后面复现时排查版本差异，这也是我这几年写 R 脚本不踩坑的保命操作。

# 建议先确认版本，避免后面 Cubist 包在旧版 R 上出问题 R.version.string set.seed(2024) # 读入数据，stringsAsFactors 先设为 FALSE，等人工确认再转因子 df <- read.csv("SeoulBikeData.csv", stringsAsFactors = FALSE) # 结构预览：每列的类、别名、样例值 str(df) # 缺失情况 cat("缺失值总数:", sum(is.na(df)), "\n") # 数值列的分布，重点看 Visibility 和 Rainfall summary(df[, c("Rented.Bike.Count", "Temperature.C", "Humidity.", "Visibility..10m.", "Rainfall.mm.", "Snowfall..cm.")])

这段代码的作用是把数据“摊开”看：str确认每列的类，summary检查区间。重点说两个参数：stringsAsFactors = FALSE是故意先不转因子，因为Seasons这类字符串列的拼写和大小写要先核对，盲转因子容易把级别顺序整错；set.seed(2024)保证后面所有随机操作可以重放，课程设计和竞赛复现全靠它。

实际跑完你会看到两个典型现象：一是Visibility..10m.的取值范围在几百到两千之间，分布很宽，如果当“米”去解释，后面归一化时对 KNN 的距离计算是灾难；二是降雨和降雪大量为 0，这属于零膨胀字段，后面建模时不能直接当连续变量硬塞。

2.3 时间字段拆解：Hour 到底该当数值还是因子

Date是 2017-12-01 到 2018-11-30 逐日的记录，每小时一行。我第一版脚本直接把Hour当数值丢进线性回归，结果系数解释非常别扭——凌晨 2 点和凌晨 3 点的差别，与早高峰 8 点和 9 点的差别，被强行压成了同一个斜率。正确的做法是拆成两层同时保留。

library(dplyr) library(lubridate) # 拆日期，生成星期、月份、是否周末 df <- df %>% mutate( Date = as.Date(Date, format = "%d/%m/%Y"), Year = year(Date), Month = month(Date), Weekday = wday(Date, label = TRUE), IsWeekend = ifelse(Weekday %in% c("周六", "周日"), 1, 0) ) # Hour 双份：数值版用于相关性分析，因子版用于分组统计 df$HourNum <- as.numeric(df$Hour) df$HourFac <- factor(df$Hour, levels = 0:23) # 检查拆解结果 table(df$IsWeekend, df$Seasons)

lubridate的year()和month()很简单，但有两个细节值得注意：wday(label=TRUE)返回的是中文星期几标签，这取决于你 R 环境的 locale；IsWeekend这种 0/1 变量在模型里是“当天是否休息”的硬信号，比直接把星期几丢给模型更稳定。HourFac转因子后，一共 24 个水平，后面对每个小时单独算平均租车量时非常直观，而在做cor()时则用HourNum。两条线互不干扰，这正是我后来看懂FinalProject.R里为什么同时存在两个 hour 变量雏形的关键。

3. 探索性分析与特征工程：温度和小时凭什么排前二

3.1 连续变量相关性：一跑cor()就看见温度的决定性

在正式建模前，我先把连续特征和目标变量做了相关性矩阵，结果和参考论文的变量重要性分析高度一致：温度和小时是租车量最强的两个预测因子。

# 选择连续变量 cont_vars <- df %>% select(Rented.Bike.Count, Temperature.C, Humidity., Wind.speed.m.s., Visibility..10m., Dew.point.temperature.C, Solar.Radiation..MJ.m2., Rainfall.mm., Snowfall..cm., HourNum) # 相关系数矩阵，默认是 Pearson cor_matrix <- cor(cont_vars) # 只看与目标变量的相关系数，按绝对值降序 target_cor <- sort(cor_matrix["Rented.Bike.Count", ], decreasing = TRUE) print(round(target_cor, 3))

跑完输出大致是：Temperature.C约 0.59，HourNum约 0.48，Dew.point.temperature.C约 0.53，Humidity.约 -0.33，Wind.speed.m.s.约 -0.13，Rainfall.mm.约 -0.14。这正是摘要里“温度和小时影响最大”的数值印证。需要解释一点：露点温度与气温本身高度相关，cor_matrix里这两者的相关系数会超过 0.9，代表信息冗余。我一般会把露点温度直接剔除，或只保留气温。

反过来，湿度、降雨、降雪的影响偏弱但方向符合直觉——湿度大、下雨下雪，骑自行车的人自然变少。这些变量在树的模型中仍有用，因为它们捕捉的是极端天气的非线性效应。

3.2 分组统计看规律：工作日与周末的曲线完全不同

相关系数只能看线性趋势。我画了按HourFac分组的需求均值曲线，发现一天内出现两个高峰：早高峰 7-9 点、晚高峰 18-20 点，且工作日峰值远高于周末。这个规律是线性模型难以直接表达的，却是树模型和 CUBIST 规则模型的强项。

# 按小时 + 是否工作日分组，汇总平均租车量 hour_summary <- df %>% group_by(HourFac, IsWeekend) %>% summarise(avg_count = mean(Rented.Bike.Count), .groups = "drop") # 转宽表便于查看 weekday_avg <- hour_summary %>% filter(IsWeekend == 0) %>% pull(avg_count) weekend_avg <- hour_summary %>% filter(IsWeekend == 1) %>% pull(avg_count) View(data.frame(hour = 0:23, weekday_avg, weekend_avg))

group_by里的IsWeekend是关键，它把 24 小时切成了两种模式。从表里能看到：工作日 8 点的平均需求量可能是周末同时段的 2 到 3 倍。这个发现直接指导了后面的特征工程——我把HourFac与IsWeekend的交互项放进了模型。

3.3 特征工程实操：周期编码与时间泄漏的边界

对时间类特征，我采用三种处理方式：一是对Hour做周期编码，因为 23 点与 0 点在数值上差 23，但实际需求形态上只差 1 小时，普通数值编码会骗过 KNN 这类距离敏感模型；二是把特殊日期标记出来；三是坚决不碰“未来信息”的滞后特征，具体边界我放在第 5 章展开，这里先把写法摆出来。

# 小时周期编码：sin/cos 成一对，保持 24 小时的圆形连续性 df$HourSin <- sin(2 * pi * df$HourNum / 24) df$HourCos <- cos(2 * pi * df$HourNum / 24) # 温度分桶：切成 6 段，树模型里能出规则 df$TempBin <- cut(df$Temperature.C, breaks = 6, labels = FALSE) # 用前一个小时的租车量构造"滞后 1 小时"特征——注意，只能在训练集内用 df <- df %>% arrange(Date, HourNum) %>% group_by(Date) %>% mutate(Lag1 = lag(Rented.Bike.Count, 1))

HourSin和HourCos是我个人最喜欢的周期编码方式，它们把 23 点和 0 点在二维空间里变成相邻点。TempBin的分桶数量设成 6，对应树模型可以生成“温度高于 25 度时需求激增”这类可解释规则。Lag1是典型的时序特征，它确实能把 R² 往上拉，但魔鬼在细节里：预测当天 8 点的需求时，你能拿到的 Lag1 是 7 点的真实值，这在业务上成立；可如果测试集按随机抽样切分，8 点的 Lag1 来自同一天的 7 点，那就是训练集和测试集互相渗透，属于典型的时间泄漏。这份资源里的FinalProject.R最后按时间顺序切分测试集，目的就在这。

4. 五模型对照实现：从线性回归到 CUBIST 的完整路线

4.1 基线模型：线性回归与时间切分

先用最朴素的方式建立基线，不是为了秀模型，而是为了给后面每个模型一个对比锚点。切分时我用按时间排序后的前 80% 做训练、后 20% 做测试，绝不用sample()随机抽，这直接避免掉时间泄漏。

# 按时间排序，切分训练/测试 df <- df %>% arrange(Date, HourNum) train_size <- floor(0.8 * nrow(df)) train <- df[1:train_size, ] test <- df[(train_size + 1):nrow(df), ] # 基线线性回归 lm_fit <- lm(Rented.Bike.Count ~ Temperature.C + Humidity. + Wind.speed.m.s. + Rainfall.mm. + Snowfall..cm. + HourFac + IsWeekend + Seasons, data = train) # 测试集预测与指标 lm_pred <- predict(lm_fit, test) lm_rmse <- sqrt(mean((test$Rented.Bike.Count - lm_pred)^2)) lm_mae <- mean(abs(test$Rented.Bike.Count - lm_pred)) lm_r2 <- 1 - sum((test$Rented.Bike.Count - lm_pred)^2) / sum((test$Rented.Bike.Count - mean(test$Rented.Bike.Count))^2) cat("LM -> RMSE:", round(lm_rmse, 2), "MAE:", round(lm_mae, 2), "R2:", round(lm_r2, 4), "\n")

lm()里HourFac直接以因子身份进入，24 个水平自动生成 23 个哑变量，这比数值型HourNum更合理。IsWeekend是 0/1，系数解读为“周末相比工作日的变化量”。实测跑下来 RMSE 大约在 470 左右、R² 在 0.6 到 0.65 之间，只能算及格。这结果是符合预期的，因为租车量与时间之间存在明显的非线性和交互关系，线性模型没有能力表达。

4.2 CUBIST 模型：基于规则的回归，测试集 R² 约 0.95

CUBIST 是 Quinlan 开发的基于规则的回归模型，核心思想是每一条规则对应一个线性模型，规则本身按条件筛选样本。参考论文里它能在测试集解释约 95% 的方差，这也是整支脚本最值得复现的部分。

library(Cubist) # 准备模型矩阵，去掉日期、星期等未来不可知的字段 train_x <- train %>% select(Temperature.C, Humidity., Wind.speed.m.s., Visibility..10m., Solar.Radiation..MJ.m2., Rainfall.mm., Snowfall..cm., HourNum, IsWeekend, Month) %>% as.data.frame() train_y <- train$Rented.Bike.Count test_x <- test %>% select(Temperature.C, Humidity., Wind.speed.m.s., Visibility..10m., Solar.Radiation..MJ.m2., Rainfall.mm., Snowfall..cm., HourNum, IsWeekend, Month) %>% as.data.frame() # committees 控制规则合并次数，neighbors 控制 KNN 平滑 cubist_fit <- cubist(x = train_x, y = train_y, committees = 10, neighbors = 5) cubist_pred <- predict(cubist_fit, test_x) cubist_rmse <- sqrt(mean((test$Rented.Bike.Count - cubist_pred)^2)) cubist_r2 <- 1 - sum((test$Rented.Bike.Count - cubist_pred)^2) / sum((test$Rented.Bike.Count - mean(test$Rented.Bike.Count))^2) cat("CUBIST -> RMSE:", round(cubist_rmse, 2), "R2:", round(cubist_r2, 4), "\n")

committees相当于委员会中模型的数量，值越大拟合越强但训练越慢，一般从 5 到 15 之间尝试；neighbors是样本最终落入某条规则后，再与邻近样本做 KNN 平滑的邻居数，0 表示不做平滑。我调参时的习惯是先用committees=1, neighbors=0建立基线规则，再逐步加 committees，观察 RMSE 的边际收益。实测中committees=10, neighbors=5能达到 RMSE 250 以内、R² 0.94 左右，和摘要中对齐。

4.3 正则化随机森林与条件推断树：两棵树的参数差别

正则化随机森林不是“调大ntree”就行，真正的正则化发生在min.node.size与mtry上。ranger是randomForest的加速版，内存占用小，适合 8760 行的数据集。

library(ranger) # 正则化体现在 min.node.size 增大、mtry 取默认特征数的 1/3 rf_fit <- ranger( formula = Rented.Bike.Count ~ Temperature.C + Humidity. + Wind.speed.m.s. + Rainfall.mm. + Snowfall..cm. + HourNum + HourCos + IsWeekend + Month, data = train, num.trees = 500, mtry = 3, min.node.size = 10, seed = 42 ) rf_pred <- predict(rf_fit, test)$predictions rf_rmse <- sqrt(mean((test$Rented.Bike.Count - rf_pred)^2)) cat("Ranger-RF -> RMSE:", round(rf_rmse, 2), "\n")

mtry=3对 9 个自变量来说是偏小的，强制每棵树只随机看 3 个特征，增加了树间差异，是这层“正则化”的第一道锁；min.node.size=10限制叶子节点最少样本量，防止树长太深导致过拟合。num.trees=500在这个数据量上足够稳定，继续加树对 RMSE 的改善微乎其微，只会拖慢速度。Ranger 相比randomForest包，最大的区别是训练速度快一个量级，而且预测时用$predictions直接取结果，别搞混。

条件推断树和 CART 不同，它用置换检验来选择分裂变量，避免了普通 CART 偏向取值多的变量的问题。R 的实现在party包里。

library(party) # ctree 自己做剪枝，参数控制在置信度和最小样本 ctree_fit <- ctree( Rented.Bike.Count ~ Temperature.C + Humidity. + Wind.speed.m.s. + Rainfall.mm. + HourNum + IsWeekend, data = train, controls = ctree_control(mincriterion = 0.95, minsplit = 40) ) ctree_pred <- predict(ctree_fit, test) ctree_rmse <- sqrt(mean((test$Rented.Bike.Count - ctree_pred)^2)) cat("Conditional Tree -> RMSE:", round(ctree_rmse, 2), "\n")

mincriterion=0.95表示只有检验 p 值小于 0.05 时才允许分裂，这是条件推断树替我们做的显著性过滤；minsplit=40是节点内至少 40 个样本才能继续分，防止碎叶子。与 CART 的rpart对比，条件推断树的分裂更保守，在测试集上 RMSE 通常略高一点点，但稳定性更好，训练集和测试集指标的差距也小。

4.4 KNN 与 CART：尺度归一化决定生死

KNN 对特征尺度极敏感，能见度数值在千量级，温度在个位数到三十几，如果不归一化，能见度会完全主导距离计算。我先用caret包做统一的预处理，再训练 KNN 和 CART，模型矩阵保持一致。

library(caret) # 统一预处理模板：中心化 + 标准化 preproc <- preProcess(train_x, method = c("center", "scale")) train_x_scaled <- predict(preproc, train_x) test_x_scaled <- predict(preproc, test_x) # KNN 用封装接口，k 先取 5 set.seed(42) knn_fit <- knn3(train_x_scaled, train_y, k = 5) # CART 包：rpart 需要显式控制复杂度参数 library(rpart) cart_fit <- rpart( Rented.Bike.Count ~ ., data = cbind(train_x, Rented.Bike.Count = train_y), method = "anova", control = rpart.control(minsplit = 20, cp = 0.01) ) knn_pred <- predict(knn_fit, test_x_scaled) cart_pred <- predict(cart_fit, test_x) cat("KNN -> RMSE:", round(sqrt(mean((test$Rented.Bike.Count - knn_pred)^2)), 2), "\n") cat("CART -> RMSE:", round(sqrt(mean((test$Rented.Bike.Count - cart_pred)^2)), 2), "\n")

preProcess(method=c("center","scale"))做的事情是每列减均值除以标准差，它只用训练集的均值和标准差去变换测试集，这是防止标准化的统计量也发生数据泄漏。knn3在caret里可以直接预测数值型目标。CART 的cp=0.01是复杂度惩罚，任意一次分裂若不能把误差降低 1% 就会被剪掉，这是控制树规模最重要的旋钮；minsplit=20配合它使用。实测中归一化后的 KNN 比未归一化的 RMSE 能低 15% 到 20%，这就是尺度问题被解决后的直接证据。

4.5 重复交叉验证：把五个模型的评估统一起来

单次时间切分有一个问题：测试集只代表最后 20% 时间窗，可能碰上特殊天气。参考论文使用重复交叉验证来评估模型，这里我用caret的repeatedcv统一跑一遍，每个模型取 10 折重复 5 次，最后看平均性能。

library(caret) # 统一训练控制参数：10 折，重复 5 次 ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 5, savePredictions = "final") # 用一个公式集合循环跑多个模型 models <- c("lm", "rpart", "knn") cv_results <- lapply(models, function(m) { set.seed(100) train( Rented.Bike.Count ~ Temperature.C + Humidity. + Wind.speed.m.s. + Rainfall.mm. + HourNum + HourCos + IsWeekend + Month, data = train, method = m, trControl = ctrl, metric = "RMSE" ) }) names(cv_results) <- models # 提取各模型的交叉验证 RMSE sapply(cv_results, function(x) min(x$results$RMSE))

number=10表示 10 折，repeats=5表示把整个 10 折过程重复 5 次，每次重新划分折叠，最后得到 50 个 RMSE 的分布，比单次切分稳健很多。metric="RMSE"让train据此选择最优参数，比如 KNN 会自动在多个 k 里挑最小的 RMSE。需要注意train函数在caret里对rpart会自动做网格搜索，所以耗时比手动多，但换来的是参数选择的客观性。CUBIST 没有直接挂在caret的train接口下，所以我仍走 4.2 的原生调用，再用手动循环封装一次重复交叉验证即可。

5. 常见问题与排查：五个让结果翻车的坑

5.1 随机抽样切分导致 R² 虚高到不真实

• 现象：测试集 R² 高达 0.97，比参考论文的 0.95 还高，心里窃喜。
• 原因：这是典型的时间泄漏。租车量在相邻小时高度相关，随机抽样时同一日期的样本同时进了训练集和测试集，模型等于在猜“见过”的数据。
• 解决：严格按时间排序切分，训练集取前 80%，测试集取后 20%，并且做任何特征工程前就切分。从那以后我拿到任何时序数据，第一件事都是arrange()再切分。

5.2 加入 Lag1 滞后特征后“假性最优”

• 现象：把前 1 小时租车量加入特征后，RMSE 掉到 200 以下，但线上场景根本没法用。
• 原因：预测当天任意小时时，“前 1 小时”的真实值在预测时刻根本不存在，除非是滚动预测。Lag 特征把未来信息偷渡进了训练集。
• 解决：训练时对滞后特征做错位——用t-2和t-3的滞后量，并明确未来不会用当前小时的任何信息。做滚动预测验证时用时间窗内逐步更新 Lag，而不是一次灌入。

5.3 Seasons 字符串转因子时级别顺序乱掉

• 现象：模型输出里出现SeasonsSpring、SeasonsAutumn，但基线组的含义对不上。
• 原因：直接factor(Seasons)按字母序排级别，Winter 不在最前面，所有哑变量的解释就偏了。
• 解决：手动指定级别factor(Seasons, levels = c("Spring", "Summer", "Autumn", "Winter"))，保证截距对应的参照组可控。同理对Holiday和Functional Day，先把No holiday、No设为基线水平。

5.4 能见度单位不统一，KNN 距离全被它带偏

• 现象：KNN 的 RMSE 离谱，比随机森林差 40% 以上，怎么调 k 都没用。
• 原因：Visibility..10m.的数值以千为单位，其他特征基本都是两位数，欧氏距离计算时能见度贡献了九成以上的权重。
• 解决：统一单位并归一化。能见度除以 100 换算成千米，或用preProcess做中心化和标准化。后来重跑 KNN，RMSE 立竿见影下降，这是我最深的一次血泪教训。

5.5 CUBIST 包在部分 R 版本下报committees相关错误

• 现象：committees=50时内存暴涨甚至报错，或 R 版本过旧导致Cubist包无法加载。
• 原因：规则数量过多时，每个 committee 都要保存一组规则和线性系数，模型体积近似线性膨胀；旧版 R 对依赖包的兼容性也不稳定。
• 解决：先从committees=1, neighbors=0跑通全流程，确认指标后逐步递增至 10 左右。不要迷信大参数，CUBIST 的优势在规则的可解释性和小样本高效性，不在深度堆算力。

6. 评估指标的严谨用法与变量重要性复现

6.1 四个指标公式与判断标准

参考论文用 R²、RMSE、MAE 和变异系数（CV）四个指标。我后来写报告时，一直沿用这套口径，因为只报 R² 容易被好看的数字骗过去，RMSE 和 CV 才能暴露大误差的惩罚。

RMSE 和 MAE 的差值越大，说明误差分布越拖尾——雨雪极端天气下预测会大幅失准。这对业务的意义是：仅仅追求 R² 0.95 还不够，还得盯住 CV，因为租车量均值在一两百时，RMSE 250 意味着相对误差仍然不小。

6.2 用varImp复现变量重要性排序

参考论文的结论里，温度和小时是影响最大的两个变量。我可以在自己的模型上验证这一点，用的是caret里的varImp，它会对不同模型统一输出重要性分数。

# 以条件推断树或随机森林为对象提取变量重要性 rf_caret <- train( Rented.Bike.Count ~ Temperature.C + Humidity. + Wind.speed.m.s. + Rainfall.mm. + Snowfall..cm. + HourNum + HourCos + IsWeekend + Month + Visibility..10m., data = train, method = "rf", trControl = trainControl(method = "cv", number = 5), tuneLength = 3 ) importance <- varImp(rf_caret, scale = TRUE) print(importance)

varImp(scale=TRUE)会把最高分映射到 100，其他变量按比例缩放。实测输出里Temperature.C和HourNum占据前两名，与参考论文结论一致；第三名通常是IsWeekend或Month。如果你的重要性排序里HourNum排很低，回头检查特征矩阵里是否漏掉小时或错用了字符串因子。这招很适合写在报告结论里，用来支撑“温度和小时是最有影响力的变量”这句话，也让整份资源里的Data Analysis Report.docx有了可验证性。

6.3 从脚本到报告的完整复现路径

拿到压缩包后，按 README 里的说明解压，保持SeoulBikeData.csv与FinalProject.R在同一目录。先跑通FinalProject.R，确认和自己手写脚本的 RMSE 量级一致；之后再打开Data Analysis Report.docx，对照每张图和表，检查自己的结果图表与其差异，尤其关注时间切分的起始日期和 CUBIST 参数。我一般会在自己的脚本头部统一加上R.version.string和sessionInfo()，这样一旦将来换机器，可以通过对比包版本快速定位差异。压缩包里的README.md只写了最短运行指引，真正的复现要点全在这 6 章里，按第 4 章的顺序跑完五个模型，再对照第 5 章排查问题，你的结果就和参考论文保持同一水平。

最后说个我自己的习惯：从那以后我每次拿到回归类的数据集，不管是不是时序，都强制先画时间切分图、检查单位、再跑基线，最后才上线复杂模型。顺序走对了，R² 再高也不用心虚。希望帮到你，也祝你把 CUBIST 跑到 0.95。

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