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

一、线性回归

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

代码语言：javascript复制# 加载数据

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")

二、逻辑回归

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

代码语言：javascript复制# 加载数据

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")三、多项式回归

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

代码语言：javascript复制# 创建多项式特征

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")四、泊松回归

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

代码语言：javascript复制# 加载数据

data(warpbreaks)

# 建立泊松回归模型

model <- glm(breaks ~ wool + tension, data = warpbreaks, family = poisson)

# 查看模型摘要

summary(model)

五、负二项回归

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

代码语言：javascript复制library(MASS)

# 建立负二项回归模型

model <- glm.nb(breaks ~ wool + tension, data = warpbreaks)

# 查看模型摘要

summary(model)

六、贝叶斯回归

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

代码语言：javascript复制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")

七、生存分析

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

代码语言：javascript复制# 安装并加载必要的包

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)

八、决策树

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

代码语言：javascript复制# 安装并加载必要的包

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）。

代码语言：javascript复制# 加载示例数据集

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（汽车的每加仑油行驶的英里数） 。

代码语言：javascript复制# 加载示例数据集

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是一种强大的分类算法，特别适合处理高维数据。

代码语言：javascript复制# 加载示例数据集

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), "%"))

十一、神经网络

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

代码语言：javascript复制# 加载数据

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是线性模型的扩展，可以处理非正态分布的响应变量。

代码语言：javascript复制# 建立GLM模型 (以逻辑回归为例)

model <- glm(am ~ wt + hp, data = mtcars, family = binomial)

# 查看模型摘要

summary(model)

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

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

代码语言：javascript复制# 加载必要包

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用于降维和探索性数据分析。

代码语言：javascript复制# 执行PCA

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

# 查看结果

summary(pca_result)

# 可视化

biplot(pca_result)

十五、因子分析

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

代码语言：javascript复制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个组。

代码语言：javascript复制# 执行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)十七、层次聚类

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

代码语言：javascript复制# 计算距离矩阵

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模型用于分析和预测时间序列数据。

代码语言：javascript复制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 包进行时间序列交叉验证。

代码语言：javascript复制# 安装和加载必要的包

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包进行时间序列交叉验证

代码语言：javascript复制# 安装 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）用于分析具有分层结构的数据。

代码语言：javascript复制# 多层次模型（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)