R 语言统计建模大全：20 个经典模型实战解析，速收藏！

天意生信云

发布于 2025-02-27 14:28:48

1.9K0

统计建模是数据科学中至关重要的一部分，帮助分析和预测数据中的趋势与模式。在数据科学中，常用的统计模型有回归分析、时间序列分析、分类模型、聚类模型等，每种模型有其独特的应用场景。在R语言中，我们可以通过丰富的统计包，如lm()进行线性回归分析，glm()用于广义线性模型，arima()进行时间序列建模等。这些模型能够帮助我们从数据中提取信息并做出科学决策，成为数据分析中的强大工具。

一、线性回归

线性回归是最基本也是最常用的统计模型之一，用于分析因变量与一个或多个自变量之间的线性关系。

# 加载数据
data(mtcars)

# 建立线性回归模型
model <- lm(mpg ~ wt + hp, data = mtcars)

# 查看模型摘要
summary(model)

# 可视化
plot(mtcars$wt, mtcars$mpg, main = "Weight vs MPG", xlab = "Weight", ylab = "MPG")
abline(lm(mpg ~ wt, data = mtcars), col = "red")

二、逻辑回归

逻辑回归用于预测二分类结果，常用于分类问题。

# 加载数据
data(mtcars)

# 创建二分类变量
mtcars$am_factor <- factor(mtcars$am, levels = c(0, 1), labels = c("Automatic", "Manual"))

# 建立逻辑回归模型
model <- glm(am_factor ~ wt + hp, data = mtcars, family = "binomial")

# 查看模型摘要
summary(model)

# 预测
predictions <- predict(model, type = "response")

三、多项式回归

多项式回归用于建模非线性关系。

# 创建多项式特征
mtcars$wt_squared <- mtcars$wt^2

# 建立多项式回归模型
model <- lm(mpg ~ wt + wt_squared, data = mtcars)

# 查看模型摘要
summary(model)

# 可视化
plot(mtcars$wt, mtcars$mpg, main = "Weight vs MPG (Polynomial)", xlab = "Weight", ylab = "MPG")
curve(coef(model)[1] + coef(model)[2]*x + coef(model)[3]*x^2, add = TRUE, col = "red")

四、泊松回归

泊松回归用于建模计数数据。

# 加载数据
data(warpbreaks)

# 建立泊松回归模型
model <- glm(breaks ~ wool + tension, data = warpbreaks, family = poisson)

# 查看模型摘要
summary(model)

五、负二项回归

负二项回归用于处理过度离散的计数数据。

library(MASS)

# 建立负二项回归模型
model <- glm.nb(breaks ~ wool + tension, data = warpbreaks)

# 查看模型摘要
summary(model)

六、贝叶斯回归

贝叶斯线性回归（Bayesian Linear Regression）使用贝叶斯方法对模型的参数进行推断，与传统的线性回归方法不同，贝叶斯方法为参数提供了概率分布而不是点估计。这意味着贝叶斯方法能够处理模型参数的不确定性，提供更为丰富的信息。

library(rstanarm)

# 加载数据
data(mtcars)
head(mtcars)

# 使用 rstanarm 进行贝叶斯线性回归
# 我们使用 mpg 作为目标变量，disp 和 wt 作为自变量
bayes_model <- stan_glm(mpg ~ disp + wt, data = mtcars, 
                        family = gaussian(), 
                        prior = normal(0, 10), # 设置先验分布
                        prior_intercept = normal(0, 10), 
                        chains = 4, iter = 2000)

# 查看模型摘要
summary(bayes_model)

# 预测
predictions <- predict(bayes_model, newdata = mtcars)

# 可视化预测结果
plot(mtcars$mpg, predictions, main = "贝叶斯线性回归预测 vs 真实值", 
     xlab = "真实值", ylab = "预测值")
abline(0, 1, col = "red")

七、生存分析

生存分析用于分析时间到事件数据。

# 安装并加载必要的包
install.packages("survival")
install.packages("survminer")

library(survival)
library(survminer)

# 加载示例数据集
data(lung)
head(lung)

# 创建生存对象
surv_obj <- Surv(time = lung$time, event = lung$status)
print(surv_obj)

# 拟合 Kaplan-Meier 生存曲线（整体数据）
fit <- survfit(surv_obj ~ 1, data = lung)
summary(fit)

# 绘制 Kaplan-Meier 生存曲线
ggsurvplot(fit, data = lung, conf.int = TRUE, pval = TRUE,
           risk.table = TRUE, surv.median.line = "hv",
           ggtheme = theme_minimal())

# 按变量（性别 sex）分组进行生存分析
fit_sex <- survfit(Surv(time, status) ~ sex, data = lung)
ggsurvplot(fit_sex, data = lung, pval = TRUE, conf.int = TRUE,
           risk.table = TRUE, ggtheme = theme_minimal())

# Cox 比例风险回归模型
cox_model <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)
summary(cox_model)

# 绘制 Cox 回归森林图
ggforest(cox_model, data = lung)

八、决策树

决策树是一种直观的机器学习模型，常用于分类和回归问题。

# 安装并加载必要的包
install.packages("rpart")
install.packages("rpart.plot")

library(rpart)
library(rpart.plot)

# 加载示例数据集 iris
data(iris)
head(iris)

# 构建决策树模型（分类）
tree_model <- rpart(Species ~ ., data = iris, method = "class")

# 可视化决策树
rpart.plot(tree_model, type = 3, extra = 101, fallen.leaves = TRUE)

# 进行预测（使用训练数据集）
predictions <- predict(tree_model, iris, type = "class")

# 计算模型准确率
accuracy <- sum(predictions == iris$Species) / nrow(iris)
print(paste("模型准确率：", round(accuracy * 100, 2), "%"))

# 构建决策树回归模型
set.seed(123)
data(mtcars)
tree_reg <- rpart(mpg ~ ., data = mtcars, method = "anova")

# 可视化回归决策树
rpart.plot(tree_reg, type = 3, extra = 101, fallen.leaves = TRUE)

# 进行预测
pred_mpg <- predict(tree_reg, mtcars)

# 输出预测值和实际值的对比
comparison <- data.frame(Actual = mtcars$mpg, Predicted = pred_mpg)
print(head(comparison))

九、随机森林

（1）使用随机森林进行分类

我们使用iris数据集进行分类任务，预测鸢尾花的种类（Species）。

# 加载示例数据集
data(iris)
head(iris)

# 构建随机森林分类模型
rf_model_class <- randomForest(Species ~ ., data = iris, ntree = 100)

# 查看模型摘要
print(rf_model_class)

# 绘制误差率图
plot(rf_model_class)

# 进行预测
pred_class <- predict(rf_model_class, iris)

# 计算准确率
accuracy_class <- sum(pred_class == iris$Species) / nrow(iris)
print(paste("分类模型准确率：", round(accuracy_class * 100, 2), "%"))

（2）使用随机森林进行回归

我们使用mtcars数据集进行回归任务，预测mpg（汽车的每加仑油行驶的英里数）。

# 加载示例数据集
data(mtcars)
head(mtcars)

# 构建随机森林回归模型
rf_model_reg <- randomForest(mpg ~ ., data = mtcars, ntree = 100)

# 查看模型摘要
print(rf_model_reg)

# 绘制误差率图
plot(rf_model_reg)

# 进行预测
pred_reg <- predict(rf_model_reg, mtcars)

# 输出预测值和实际值的对比
comparison_reg <- data.frame(Actual = mtcars$mpg, Predicted = pred_reg)
print(head(comparison_reg))

# 计算回归模型的R平方值
rsq_reg <- cor(mtcars$mpg, pred_reg)^2
print(paste("回归模型R平方值：", round(rsq_reg, 2)))

十、支持向量机

SVM是一种强大的分类算法，特别适合处理高维数据。

# 加载示例数据集
data(iris)
head(iris)

# 构建支持向量机分类模型
svm_model <- svm(Species ~ ., data = iris, kernel = "radial", cost = 1, scale = TRUE)

# 查看模型摘要
summary(svm_model)

# 进行预测
pred_svm <- predict(svm_model, iris)

# 计算准确率
accuracy_svm <- sum(pred_svm == iris$Species) / nrow(iris)
print(paste("SVM分类模型准确率：", round(accuracy_svm * 100, 2), "%"))

十一、神经网络

神经网络是深度学习的基础，能够处理复杂的模式识别任务。

# 加载数据
data(iris)
head(iris)

# 将因子变量转换为二进制指示变量（one-hot encoding）
iris$Species <- as.numeric(iris$Species) 

# 划分训练集和测试集
set.seed(123)
index <- sample(1:nrow(iris), 0.7 * nrow(iris))
train_data <- iris[index, ]
test_data <- iris[-index, ]

# 构建神经网络模型
nn_model <- neuralnet(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, 
                      data = train_data, hidden = c(5, 3), linear.output = FALSE)

# 绘制神经网络结构
plot(nn_model)

# 进行预测
nn_pred <- compute(nn_model, test_data[, -5])$net.result
nn_pred_class <- round(nn_pred)  # 取整获得分类结果

# 计算准确率
accuracy_nn <- sum(nn_pred_class == test_data$Species) / nrow(test_data)
print(paste("神经网络分类模型准确率：", round(accuracy_nn * 100, 2), "%"))

十二、广义线性模型 (GLM)

GLM是线性模型的扩展，可以处理非正态分布的响应变量。

# 建立GLM模型 (以逻辑回归为例)
model <- glm(am ~ wt + hp, data = mtcars, family = binomial)

# 查看模型摘要
summary(model)

十三、广义加性模型 (GAM)

GAM允许预测变量与响应变量之间的非线性关系。

# 加载必要包
library(mgcv)

# 构造数据集：这里我们生成 32 个样本
set.seed(123)  # 保证结果可重复
n <- 32
mpg <- rnorm(n, mean = 20, sd = 5)   # mpg：均值20，标准差5的正态分布
wt  <- rnorm(n, mean = 3,  sd = 0.5)  # wt：均值3，标准差0.5的正态分布
hp  <- rnorm(n, mean = 150, sd = 30)  # hp：均值150，标准差30的正态分布

# 创建数据框，命名为 mtcars（与内置数据集同名）
mtcars <- data.frame(mpg = mpg, wt = wt, hp = hp)

# 建立 GAM 模型
model <- gam(mpg ~ s(wt) + s(hp), data = mtcars)

# 查看模型摘要
summary(model)

# 可视化模型中平滑函数的估计
plot(model)

十四、主成分分析 (PCA)

PCA用于降维和探索性数据分析。

# 执行PCA
pca_result <- prcomp(mtcars[, c("mpg", "disp", "hp", "drat", "wt")], scale. = TRUE)

# 查看结果
summary(pca_result)

# 可视化
biplot(pca_result)

十五、因子分析

因子分析用于探索变量之间的潜在结构。

library(psych)

# 执行因子分析
fa_result <- fa(mtcars[, c("mpg", "disp", "hp", "drat", "wt")], nfactors = 2, rotate = "varimax")

# 查看结果
print(fa_result)

# 可视化
fa.diagram(fa_result)

十六、聚类分析 (K-means)

K-means是一种常用的聚类算法，用于将数据分成K个组。

# 执行K-means聚类
kmeans_result <- kmeans(mtcars[, c("mpg", "wt")], centers = 3)

# 可视化结果
plot(mtcars$wt, mtcars$mpg, col = kmeans_result$cluster, pch = 19,
     main = "K-means Clustering", xlab = "Weight", ylab = "MPG")
points(kmeans_result$centers[, c("wt", "mpg")], col = 1:3, pch = 8, cex = 2)

十七、层次聚类

层次聚类创建一个树状结构来表示数据的聚类。

# 计算距离矩阵
dist_matrix <- dist(mtcars[, c("mpg", "wt")])

# 执行层次聚类
hc_result <- hclust(dist_matrix, method = "ward.D2")

# 可视化结果
plot(hc_result, main = "Hierarchical Clustering Dendrogram", xlab = "", sub = "")

十八、时间序列分析 (ARIMA)

ARIMA模型用于分析和预测时间序列数据。

library(forecast)

# 创建时间序列对象
ts_data <- ts(AirPassengers, frequency = 12)

# 拟合ARIMA模型
model <- auto.arima(ts_data)

# 查看模型摘要
summary(model)

# 预测
forecast_result <- forecast(model, h = 12)
plot(forecast_result)

十九、时间序列交叉验证

使用 caret 包进行时间序列交叉验证。

# 安装和加载必要的包
install.packages("caret")
library(caret)
library(tseries)

# 示例数据
data(AirPassengers)
ts_data <- AirPassengers

# 定义时间序列的长度
train_size <- length(ts_data) * 0.8
train_data <- ts_data[1:train_size]
test_data <- ts_data[(train_size + 1):length(ts_data)]

# 设置交叉验证的参数
train_control <- trainControl(method = "timeslice", 
                              initialWindow = 36,  # 初始训练窗口大小（例如36个月）
                              horizon = 12,        # 每次验证期为12个月
                              fixedWindow = TRUE)  # 固定窗口（滚动窗口）

# 使用 ARIMA 模型进行训练和验证
model <- train(train_data ~ 1, 
               method = "auto", 
               trControl = train_control)

# 查看模型结果
print(model)

使用tscv包进行时间序列交叉验证

# 安装 tsccv 包
install.packages("tscv")
library(tscv)

# 示例数据
data(AirPassengers)
ts_data <- AirPassengers

# 定义时间序列交叉验证设置
cv_results <- tscv(ts_data, 
                  k = 5,                # 将数据集分为5个时间段
                  window_type = "rolling", # 使用滚动窗口
                  horizon = 12,           # 每次验证期为12个月
                  initial_window = 36)    # 初始窗口为36个月

# 查看交叉验证结果
print(cv_results)

二十、多层次模型

多层次模型（也称为混合效应模型，Mixed Effects Model）用于分析具有分层结构的数据。

# 多层次模型（Multilevel Models, MLM）在 R 中的应用

# 加载必要的包
install.packages("lme4")
library(lme4)
install.packages("nlme")
library(nlme)

# 生成示例数据
set.seed(123)
school_data <- data.frame(
  score = rnorm(1000, mean = 75, sd = 10),  # 学生成绩
  teacher = factor(rep(1:50, each = 20)),   # 50 名教师，每位教师教授 20 名学生
  school = factor(rep(1:10, each = 100)),    # 10 所学校，每所学校有 100 名学生
  study_hours = rnorm(1000, mean = 5, sd = 1) # 学习时间
)

# 1. 仅包含随机截距的模型
model1 <- lmer(score ~ (1 | school) + (1 | teacher), data = school_data)
summary(model1)

# 2. 添加固定效应变量 study_hours
model2 <- lmer(score ~ study_hours + (1 | school) + (1 | teacher), data = school_data)
summary(model2)

# 3. 允许斜率在学校层面随机变化
model3 <- lmer(score ~ study_hours + (study_hours | school) + (1 | teacher), data = school_data)
summary(model3)

# 4. 使用 nlme 包进行建模
model_nlme <- lme(score ~ study_hours, random = ~1 | school/teacher, data = school_data)
summary(model_nlme)

# 5. 比较模型的 AIC 值
AIC(model1, model2, model3)

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2025-02-27，如有侵权请联系 cloudcommunity@tencent.com 删除

模型