使用lm()函数建立多元线性回归模型，使用glm()建立广义线性回归模型。
1. 青少年赌博情况数据集
# 加载包
library(faraway)
library(ggcorrplot)
library(ggplot2)
# 加载数据
data(teengamb)
head(teengamb)
teengamb数据集
# 建立回归模型
## 一元线性回归模型
lm1 <- lm(gamble~income,data = teengamb)
summary(lm1)
ggplot(teengamb,aes(x = income,y = gamble)) +
theme_bw() +
geom_point()
一元线性回归结果
## 探索各个变量之间的相关性
gamble_cor <- cor(teengamb)
ggcorrplot(gamble_cor,method = "square",lab = TRUE)+ # square表示矩形热图，lab=TRUE表示添加相关系数
theme(axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10))
## 多元线性回归模型
lm1 <- lm(gamble~sex+income+verbal,data = teengamb)
summary(lm1)
## 剔除无关变量后的多元线性回归模型
lm2 <- lm(gamble~sex+income,data = teengamb)
summary(lm2)
## sex和income跟gamble相关性更强
变量之间的相关性热图
多元线性回归结果
# 利用所有变量进行逐步回归
## 分为训练集和测试集
set.seed(12)
index <- sample(nrow(teengamb),round(nrow(teengamb)*0.7))
train <- teengamb[index,]
test <- teengamb[-index,]
teengamblm <- lm(gamble~.,data = train)
summary(teengamblm)
## 预测回归模型
library(Metrics)
prelm <- predict(teengamblm,test)
sprintf("均方根误差为：%f",mse(test$gamble,prelm))
# '均方根误差为：322.534961'
## 使用step()函数进行逐步回归
gamblestep <- step(teengamblm,direction = "both")
summary(gamblestep)
prestep <- predict(gamblestep,test)
sprintf("均方根误差：%f",mse(test$gamble,prestep))
## 逐步回归后的均方根误差反而增大了，Adjusted R-squared值有稍许增加，
## 说明逐步回归的结果相比之前的模型结果并没有很好的提升
# Logistic回归模型
## 在训练集上使用所有变量进行逻辑回归
teengamblm <- glm(sex~.,data = teengamb,family = "binomial")
## 对逻辑回归模型进行逐步回归，来筛选变量
teengambstep <- step(teengamblm,direction = "both")
summary(teengambstep)
逻辑回归模型结果
2. Adult收入数据集
# 数据处理
## 读取数据
## 因为含有字符型变量，所以为了回归分析，需要把变量变为因子变量
adult <- read.csv("adult.csv",header = F,stringsAsFactors = T)
head(adult)
## 添加变量名
colname<-c("age","workclass","fnlwgt","education","education.num",
"marital.status","occupation","relationship",
"race","sex","capital.gain","capital.loss","hours.per.week",
"native.country","class")
colnames(adult)<-colname
## 查看各变量类型
str(adult)
数据集
table()函数--R语言
## 缺失值识别
sum(is.na(adult))
## 判断是否存在非NA型的缺失值情况
## table()函数统计每个因子的频数
table(adult$workclass,exclude = NULL)
缺失值识别
## 缺失值处理
## ?字符前面有一个空格
## 可以看到结果当中还有？这一项，再找找有没有更合适的处理方法
adult$occupation[adult$occupation==" ?"] <- NA
adult$native.country[adult$native.country==" ?"] <- NA
table(adult$workclass,exclude = NULL)
## 缺失值删除
adult1 <- na.omit(adult)
nrow(adult)
nrow(adult1)
head(adult1)
缺失值处理结果
# 可视化分析数据信息
## 相关系数热力图,只能选取数值型变量
library(dplyr)
adult_cor <- cor(adult1 %>% select(age,fnlwgt,education.num,
capital.gain,capital.loss,hours.per.week))
ggcorrplot(adult_cor,method = "square",lab = TRUE)+
theme(axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10))
相关系数热力图
## 不同教育程度的收入差异（直方图）
library(patchwork)
p1 <- ggplot(adult1,aes(x = education,fill=class)) +
theme_bw(base_family = "STKaiti") +
geom_histogram(stat = "count") +
coord_flip() +
labs(x="教育程度",y="人数",title = "不同教育程度的收入差异情况") +
theme(axis.text.x = element_text(vjust = 0.5),plot.title = element_text(hjust = 0.5))
## geom_histogram()函数默认position="stack",show.legend=T
## x是因子变量,这里应该是stat="count"，如果是数值变量，这里应该是stat="identity"
## 给了我启示，字符型变量的count不需要先计算生成一列num再绘图，
## 转换成因子变量可以直接绘图
p2 <- ggplot(adult1,aes(x = education,fill=class)) +
theme_bw(base_family = "STKaiti") +
geom_histogram(stat = "count",position = "fill") +
coord_flip() +
labs(x="教育程度",y="人数",title = "不同教育程度的收入差异情况") +
theme(axis.text.x = element_text(vjust = 0.5),plot.title = element_text(hjust = 0.5))
options(repr.plot.width=12, repr.plot.height=6)
p1+p2+plot_layout(guides = 'collect')
## 博士和硕士的人数虽然较少，但博士和硕士学历大多数都能拿到较高的工资
不同教育程度的收入差距
## 受教育年限不同的收入差异
p1 <- ggplot(adult1,aes(x = education.num,fill=class)) +
theme_bw(base_family = "STKaiti") +
geom_histogram(stat = "count") +
coord_flip() +
labs(x="教育年限",y="人数",title = "不同教育年限的收入差异情况") +
theme(axis.text.x = element_text(vjust = 0.5),plot.title = element_text(hjust = 0.5))
p2 <- ggplot(adult1,aes(x = education.num,fill=class)) +
theme_bw(base_family = "STKaiti") +
geom_histogram(stat = "count",position = "fill") +
coord_flip() +
labs(x="教育年限",y="人数",title = "不同教育年限的收入差异情况") +
theme(axis.text.x = element_text(vjust = 0.5),plot.title = element_text(hjust = 0.5))
p1+p2+plot_layout(guides = 'collect')
## 接受教育的年限越长，其收入较高的可能性越大
不同教育年限的收入差异
## 不同性别的收入差异
p1 <- ggplot(adult1,aes(x = sex,fill=class)) +
theme_bw(base_family = "STKaiti") +
geom_histogram(stat = "count") +
labs(x="性别",y="人数",title = "不同性别的收入差异情况") +
theme(axis.text.x = element_text(vjust = 0.5),plot.title = element_text(hjust = 0.5))
p2 <- ggplot(adult1,aes(x = sex,fill=class)) +
theme_bw(base_family = "STKaiti") +
geom_histogram(stat = "count",position = "fill") +
coord_flip() +
labs(x="性别",y="人数",title = "不同性别的收入差异情况") +
theme(axis.text.x = element_text(vjust = 0.5),plot.title = element_text(hjust = 0.5))
## 不同年龄的收入差异
p3 <- ggplot(adult1,aes(x = class,y = age,fill=class)) +
theme_bw(base_family = "STKaiti") +
geom_boxplot() +
labs(x="收入",y="年龄",title = "不同年龄的收入差异情况") +
theme(axis.text.x = element_text(vjust = 0.5),plot.title = element_text(hjust = 0.5))
## 不同工作时间的收入差异
p4 <- ggplot(adult1,aes(x = class,y = hours.per.week,fill=class)) +
theme_bw(base_family = "STKaiti") +
geom_boxplot() +
labs(x="收入",y="工作时间",title = "不同工作时间的收入差异情况") +
theme(axis.text.x = element_text(vjust = 0.5),plot.title = element_text(hjust = 0.5))
options(repr.plot.width=12, repr.plot.height=12)
(p1+p2)/(p3+p4)+plot_layout(guides = 'collect')
不同因素的收入差异
# Logistic回归模型对年收入进行预测
adult1$class <- factor(adult1$class,levels=c(" <=50K"," >50K"),labels=c(0,1))
## 数据集切分为训练集和测试集
index <- sample(nrow(adult1),round(nrow(adult1)*0.7))
train <- adult1[index,]
test <- adult1[-index,]
## 在训练集上使用变量进行逻辑回归
adult1train <- glm(class~.,data = train,family = "binomial")
## 逐步回归筛选变量
adult1step <- step(adult1train,direction = "both")
summary(adult1step)
3. Prostate前列腺病人的数据集
# 读取数据
data(prostate)
head(prostate)
# 多元线性回归
## 探索各个变量之间的相关性
prostate_cor <- cor(prostate)
options(repr.plot.width=6, repr.plot.height=6)
ggcorrplot(prostate_cor,method = "square",lab = TRUE)+ # square表示矩形热图，lab=TRUE表示添加相关系数
theme(axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10))
## 多元线性回归模型
lm <- lm(lpsa~.,data = prostate)
summary(lm)
数据集
变量之间的相关系数热力图
# 逐步回归
lmstep <- step(lm,direction = "both")
summary(lmstep)
## Adjusted R-squared值基本没变化，但去掉了三个变量lcp,gleason,pgg45
# Ridge回归
## 加载包
library(glmnet)
## 数据集拆分为训练集和测试集
set.seed(123)
index <- sample(nrow(prostate),round(nrow(prostate)*0.7))
train <- prostate[index,]
test <- prostate[-index,]
## 数据标准化
train_s <- scale(train,center = TRUE,scale = TRUE)
test_s <- scale(test,center = TRUE,scale = TRUE)
head(train_s)
head(test_s)
## 利用训练集和cv.glmnet()函数分析模型参数的影响
lambdas <- seq(0,5,length.out = 200)
x <- train_s[,1:8]
y <- train_s[,9]
set.seed(1234)
ridge_model <- cv.glmnet(x,y,alpha=0,lambda=lambdas,nfolds = 3)
summary(ridge_model)
p1 <- plot(ridge_model)
p2 <- plot(ridge_model$glmnet.fit,"lambda",label = T)
p1+p2
## 建立ridge模型
ridge_min <- ridge_model$lambda.min
ridge_min
ridge_best <- glmnet(x,y,alpha = 0,lambda = ridge_min)
coef(ridge_best)
## 验证预测结果
test_pre <- predict(ridge_best,as.matrix(test_s[,1:8]))
sprintf("标准化后平均绝对误差为：%f",mae(test_s[,9],test_pre))
# '标准化后平均绝对误差为：0.416944'
p1
p2
# Lasso回归
set.seed(1234)
lasso_model <- cv.glmnet(x,y,alpha=1,lambda = lambdas,nfolds = 3)
plot(lasso_model)
plot(lasso_model$glmnet.fit,"lambda",label = T)
## 建立lasso模型
lasso_min <- lasso_model$lambda.min
lasso_best <- glmnet(x,y,alpha = 1,lambda = lasso_min)
coef(lasso_best)
## 验证预测结果
test_pre <- predict(lasso_best,test_s[,1:8])
sprintf("标准化后平均绝对误差为：%f",mae(test_s[,9],test_pre))
# '标准化后平均绝对误差为：0.421402'
日常废话：这章内容好多啊，包括很多公式推导，还得好好消化一番~~~
