一、线性回归分析
1、lm()函数
lm()函数是用于拟合线性模型(Linear Models)的主要函数。线性模型是一种统计方法,用于描述一个或多个自变量(预测变量、解释变量)与因变量(响应变量)之间的关系。它可以处理简单的线性回归、多元线性回归以及带有分类预测变量的回归(通过创建虚拟变量或指示变量)。
基本格式:
lm(formula, data, subset, weights, ...)
- formula:描述因变量与自变量间关系的符号表达式。
- data:包含公式中所有变量的数据框(data frame)或列表(list)。若未明确指定,R 将在全局环境中搜索变量。
- subset(子集):逻辑向量或表达式,用于从数据中筛选用于模型拟合的观测值。默认为NULL,即使用全部数据。
- weights(权重):可选参数,用于为各观测值分配权重。默认为
NULL
,即所有观测值权重相等。- ...(其他参数):lm函数还接受其他多个参数,这些参数通常与模型的拟合与优化相关。例如,na.action参数可用于定义缺失值(NA)的处理方式,method参数可用于指定拟合方法(尽管对于普通线性模型,此参数通常设为默认值 "qr" 即可)。
2、简单线性回归
用R语言内置的cars数据集做演示,此数据集记录了汽车的速度(speed)和停车距离(dist),一共50条记录。
head(cars, n=5)
# 简单线性模型拟合
fit <- lm(dist ~ speed, data=cars)
# 拟合结果的详细信息
summary(fit)
# 模型参数
coeffcients(fit)
# 回归系数置信区间
confint(fit)
# 模型预测值
fitted(fit)
# 模型的残差
residuals(fit)
从上面结果可知,拟合得到的模型参数的截距项为-17.5791,回归系数是3.9324,调整的多重R2(Adjusted R-squared)为0.6438,说明该模型能解释停车距离为64.38%的变异。方差分析结果也显示整个模型是显著的(p=1.49e-12 < 0.05)。因为简单线性回归只有一个自变量,所以模型的F检验和回归系数的t检验的结果是相同的。
plot(cars)
lines(x=cars$speed, y=fitted(fit), col="red")
3、多重线性回归
多重线性回归包含多个自变量。
下面使用R语言内置的数据集mtcars进行演示,此数据集包含了32种汽车的11种基本性能数据。通过汽车排量(disp),总功率(hp),后桥速比(drat)和车重(wt)四个变量来预测汽车油耗指数(mpg),mpg越大,油耗越低。
head(mtcars, n=5)
fit <- lm(mpg ~ disp + hp + drat + wt, data=mtcars)
summary(fit)
从以上结果可知:汽车排量和后桥速比与汽车油耗指数正相关,而汽车总功率和车重于汽车油耗指数负相关。在多重线性回归中,回归系数表示当1个自变量每增加1个单位,且其它自变量不变时,因变量所增加或减少的数量,例如,车重的回归系数为-3.479668,表示当排量、总功率和后桥速比不变时,车重每增加1个单位,汽车油耗指数将下降约3.48个单位。方差分析结果表明,整个回归模型是显著的(F=34.82,p=2.704e-10<0.01)。在截距项和回归系数显著性检验中,截距项(Intercept)、总功率(hp)和车重(wt)的回归系数显著(Pr<0.05) ,排量(disp)和后桥速比(drat)的回归系数不显著。整个模型能解释油耗指数81.36%的变异。
4、plot()函数
R语言中有一个实用的基础函数plot(),可以生成四种回归模型诊断图:残差图、正态QQ图、尺度-位置图和残差-杠杆图。
fit <- lm(mpg ~ disp+hp+drat+wt, data=mtcars)
# 将四种形态组合成一张图
par(mfrow=c(2,2))
plot(fit)
5、多重共线性
如果自变量之间为多重共线性,即自变量之间有较强的相关性,将使回归系数的估计产生非常严重的误差,以至于估计出来的回归系数没有任何意义。如果要判断回归模型是否存在严重的多重共线性,可以使用方差膨胀因子。
library(car)
fit <- lm(mpg ~ disp+hp+drat+wt, data=mtcars)
vif <- vif(fit)
vif
# 查看哪些变量膨胀因子大于10
vif > 10
# 查看哪些变量膨胀因子的开方大于2
sqrt(vif) > 2
从上面结果可知,如果以方差膨胀因子是否大于10来作为判断准则,那么该回归模型中不存在严重的多重共线性;如果以方差膨胀因子的开方大于2为判断准则,那么该回归模型中存在disp和wt两个变量时,存在严重的多重共线性。
二、判别分析
判别分析就是利用若干个特征来表征事物,通过对这些特征的定量分析,最终将事物判定为某一已知总体。
常见的判别分析有如下三种。
1、距离判别
距离判别(Distance-based Discriminant Analysis)对空间中的某个点进行类属判别,最容易想到的是使用该点与各已知总体的距离远近来进行判别。
对数据进行距离判别,有很多种选择:借助mahalanobis()函数得到马氏距离,接着自编函数进行距离判别;使用WMDB扩展包的wmd()函数,此函数可以进行加权或非加权的马氏距离判别;使用WeDiBaDis扩展包的WDBdisc()函数,此函数也可以进行加权或非加权的马氏距离判别。
以下是如何在R中实现基于距离的分类的基本步骤:
1.1 准备数据
确保你的数据集已经加载并准备好。数据集应该包含特征变量(用于计算距离)和目标变量(类别标签)。
1.2 计算类别中心
对于每个类别,计算其所有样本的均值(或其他代表点),这将作为该类别的中心。
1.3 计算距离
对于新的未知样本,计算它到每个类别中心的距离。可以使用欧氏距离、马氏距离等。
1.4 分类
将样本分类到距离最小的类别中。
1.5 评估模型
使用测试集评估模型的性能,通常通过混淆矩阵、准确率等指标。
1.6 示例
使用R语言中内置的iris数据集进行演示,此数据集包含了3类鸢尾花(setosa、versicolor和virginica)的4个特征,从150条记录。
# 先查看数据信息
head(iris)
str(iris)
library(iris)
describe(iris)
# 从iris数据集中随机抽取3种鸢尾花的数据各一条作为测试集,剩余的作为训练集
# 设定随机种子
set.seed(1234)
# 随机抽取测试集
data <- cbind(rownames = rownames(iris),iris) # 将行名添加为数据框的一列
library(dplyr)
test_data <- data %>% group_by(Species) %>% sample_n(1)
# 剩余数据作为训练集
train_data <- filter(data, !(rownames %in% test_data$rownames))
# 移除行名列以进行后续计算
test_data <- test_data[,-1] %>% ungroup()
test_data
# 移除行名列以进行后续计算
train_data <- train_data[,-1] %>% ungroup()
head(train_data,n=10)
使用WDBdisc()函数进行马氏距离判别:
4.4.2版本的R语言不支持安装WeDiBaDis扩展包。
# 将数据框转换为矩阵
library(dplyr)
test_data1 <- mutate(test_data, Species=as.numeric(Species)) %>%
as.matrix()
train_data1 <- mutate(train_data, Species=as.numeric(Species)) %>%
as.matrix()# 进行马氏距离判别
library(WeDiBaDis)
fit1 <- WDBdisc(data=train_data1, datatype="m", classcol=5, distance="Mahalanobis", method="DB")
summary(fit1)
如下使用欧氏距离进行基于距离的分类:
# 查看数据集
head(iris, n=5)
# 加载数据集
data(iris)# 拆分数据集为训练集和测试集
set.seed(12345)
index <- sample(1:nrow(iris), 0.7 * nrow(iris)) # 70%训练集
train_data <- iris[index, -5] # 训练集,去掉最后的类别标签用于计算中心
train_labels <- iris[index, 5]test_data <- iris[-index, -5] # 测试集
test_labels <- iris[-index, 5]# 计算类别中心
centers <- aggregate(train_data, by=list(Species=train_labels), FUN=mean)# 定义一个函数来计算欧氏距离
euclidean_distance <- function(x, y) {sqrt(sum((x - y)^2))
}# 对测试集中的每个样本进行分类
predictions <- apply(test_data, 1, function(row) {distances <- sapply(split(centers[, -1], centers$Species), function(center) {euclidean_distance(row, center)})# 返回距离最小的类别names(which.min(distances))
})# 评估模型性能
conf_matrix <- table(Predicted=predictions, Actual=test_labels)
accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix)
print(conf_matrix)
print(paste("Accuracy:", round(accuracy, 2)))
2、Fisher判别
Fisher判别分析(Fisher Discriminant Analysis, FDA),也被称为线性判别分析(Linear Discriminant Analysis, LDA)在统计模式识别领域有着广泛的应用。尽管“Fisher判别分析”和“线性判别分析”在术语上存在些许差异,但在大多数情况下,它们指的是同一种方法。FDA/LDA的目标是找到一个线性组合(或投影)方向,使得在这个方向上,不同类别之间的样本投影点尽可能分开,而同一类别内的样本投影点尽可能紧凑。
使用MASS扩展包的lda()函数做演示:
library(MASS) # 包含lda函数
library(ggplot2) # 可视化# 使用经典鸢尾花数据集
data(iris)
head(iris)# 数据预处理
set.seed(12) # 设置随机种子保证可重复性
train_index <- sample(1:nrow(iris), nrow(iris)*0.9) # 90%训练集
train_data <- iris[train_index, ]
test_data <- iris[-train_index, ]# 执行Fisher判别分析(LDA);Species ~ .表示使用所有特征预测品种
lda_model <- lda(Species ~ ., data = train_data)# 查看模型概要
print(lda_model)
lda()函数会输出各类别的先验概率(Prior probabilities)、分组均值(Group means)、判别函数系数(Coefficients of linear discriminants)和迹的比重(Proportion of trace)。其中,LD1能解释总变异的99.22%,LD2只能解释总变异的0.78%,故LD1就是所需要的线性函数。
# 模型预测
predictions <- predict(lda_model, newdata = test_data)# 生成混淆矩阵
confusion_matrix <- table(Predicted = predictions$class, Actual = test_data$Species)
print(confusion_matrix)
# 计算准确率
accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)
cat("\n测试集准确率:", round(accuracy*100, 2), "%\n")
从上面结果可知,总共15种预测,全都预测成功。
# 可视化判别结果
projected_data <- data.frame(LD1 = predictions$x[,1],LD2 = predictions$x[,2],Species = test_data$Species
)ggplot(projected_data, aes(x = LD1, y = LD2, color = Species)) +geom_point(size = 3) +stat_ellipse(level = 0.95) +labs(title = "Fisher判别投影结果",x = "第一判别函数",y = "第二判别函数") +theme_minimal()
3、Bayes判别
使用klaR扩展包中的NaiveBayes()函数。
library(klaR)
set.seed(12) # 设置随机种子保证可重复性
train_index <- sample(1:nrow(iris), nrow(iris)*0.9) # 90%训练集
train_data <- iris[train_index, ]
test_data <- iris[-train_index, ]
# 首先建立先验概率相等的Bayes判别模型
data1 <- NaiveBayes(Species ~ ., data=train_data)
# 建立先验概率分别为0.3,0.5,0.2的Bayes判别模型
data2 <- NaiveBayes(Species ~., data=train_data, prior=c(3/10, 5/10, 2/10))
# 查看data1和data2的结构
str(data1)
str(data2)
# 计算两个模型的混淆矩阵
x <- table(Actual = train_data$Species, predicted = predict(data1, train_data)$class)
y <- table(Actual = train_data$Species, predicted = predict(data2, train_data)$class)
x
y# 计算正确率
sum(diag(prop.table(x)))
sum(diag(prop.table(y)))
从上面结果可知,先验概率相等时,有6朵花判错;先验概率不等时,也有6朵花判错。但两者的概率相等,都是95.556%。