提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档
文章目录
- 前言
- 一、广义线性回归模型
- 二、数据集简介
- 三、使用步骤
- 1.引入库和读入数据
- 2.简单分析
- 3.详细分析
- 3.回归分析
- 4.计算混淆矩阵和整体预测率
前言
提示：glm被用于拟合广义线性模型，特别是通过给出对线性预测子的符号描述以及对误差分布的描述
一、广义线性回归模型
glm(formula, family = gaussian, data, weights, subset,na.action, start = NULL, etastart, mustart, offset,control = list(...), model = TRUE, method = "glm.fit",x = FALSE, y = TRUE, contrasts = NULL, ...)formula 数据关系/模型的描述，举例：group~性别+年龄+体重+身高+吸烟定性+饮酒定性，Y：group，X：性别，年龄，体重，身高，吸烟定性，饮酒定性
family 响应变量（因变量/Y）的分布类型，需要自己根据响应变量分布情况在参数里定义
family 种类(默认连接函数)binomial(link = "logit")#响应变量服从二项分布，连接函数为logit，即logistic回归
binomal(link=’probit’) #响应变量服从二项分布，连接函数为probit
gaussian(link = "identity")
Gamma(link = "inverse")
inverse.gaussian(link = "1/mu^2")
poisson(link = "log")#响应变量服从泊松分布，即泊松回归
quasi(link = "identity", variance = "constant")
quasibinomial(link = "logit")
quasipoisson(link = "log")data 可选择数据库，列表等
weights 是一个可选择向量，在模型拟合过程中可以“优化权重”，应该设置成NULL或者一列数字向量
subset 是一个可选择的向量，具体化在拟合过程的观察亚组
na.action 指示缺失值处理
start 线性预测子参数的起始值
etastart 线性预测子的起始值
mustart 均数向量的起始值
offset 可以被用来在拟合过程中具体化要被纳入线性预测子的优先已知部分，应该设置为NULL或者长度和案例一样的数字向量。
control 控制拟合过程参数的列表
model 指示model frame模型框是否应该被包含返回值的部分的逻辑变量
a logical value indicating whether model frame should be included as a component of the returned value.
method 用于拟合模型的方法，默认方法是"glm.fit"： 迭代重新加权最小二重法iteratively reweighted least squares (IWLS)，选择"model.frame"返回模型框并且不会拟合。
x, y 逻辑值指示被用于拟合过程的反应向量和模型矩阵是否应该被作为返回值的部分
contrasts 一个可选择的列表
intercept 逻辑，截距是否应该被包括在空模型里
object 可以从 "glm"类型里继承的对象
type 特征和偏好匹配被允许，从被拟合模型对象里提取出的权重类型，可以缩写
提示：以下是本篇文章正文内容，下面案例可供参考
二、数据集简介
示例数据来自同上文随机森林一样来源于网上比较流行的股票交易回报数据。下面是该数据集的说明：
The variables are:Year: The year that the observation was recorded
Lag1: Percentage return for previous week
Lag2: Percentage return for 2 weeks previous
Lag3: Percentage return for 3 weeks previous
Lag4: Percentage return for 4 weeks previous
Lag5: Percentage return for 5 weeks previous
Volume: Volume of shares traded (average number of daily shares traded in billions)
Today: Percentage return for this week
Direction: A factor with levels Down and Up indicating whether the market had a positive or negative return on a given week
三、使用步骤
1.引入库和读入数据
代码如下（示例）：
library(ISLR)
library(broom)
library(tidyverse)
library(ggplot2)
library(MASS)
library(class)
library(caret)
library(e1071)
data <- read.csv('weekly.csv')
2.简单分析
代码如下（示例）：
#描述性统计分析
summary(data)Year Lag1 Lag2 Lag3 Min. :1990 Min. :-18.1950 Min. :-18.1950 Min. :-18.1950 1st Qu.:1995 1st Qu.: -1.1540 1st Qu.: -1.1540 1st Qu.: -1.1580 Median :2000 Median : 0.2410 Median : 0.2410 Median : 0.2410 Mean :2000 Mean : 0.1506 Mean : 0.1511 Mean : 0.1472 3rd Qu.:2005 3rd Qu.: 1.4050 3rd Qu.: 1.4090 3rd Qu.: 1.4090 Max. :2010 Max. : 12.0260 Max. : 12.0260 Max. : 12.0260 Lag4 Lag5 Volume Today Min. :-18.1950 Min. :-18.1950 Min. :0.08747 Min. :-18.1950 1st Qu.: -1.1580 1st Qu.: -1.1660 1st Qu.:0.33202 1st Qu.: -1.1540 Median : 0.2380 Median : 0.2340 Median :1.00268 Median : 0.2410 Mean : 0.1458 Mean : 0.1399 Mean :1.57462 Mean : 0.1499 3rd Qu.: 1.4090 3rd Qu.: 1.4050 3rd Qu.:2.05373 3rd Qu.: 1.4050 Max. : 12.0260 Max. : 12.0260 Max. :9.32821 Max. : 12.0260 Direction Length:1089 Class :character Mode :character
#各个变量之间的相关性分析
cor(data[,-9])Year Lag1 Lag2 Lag3 Lag4
Year 1.00000000 -0.032289274 -0.03339001 -0.03000649 -0.031127923
Lag1 -0.03228927 1.000000000 -0.07485305 0.05863568 -0.071273876
Lag2 -0.03339001 -0.074853051 1.00000000 -0.07572091 0.058381535
Lag3 -0.03000649 0.058635682 -0.07572091 1.00000000 -0.075395865
Lag4 -0.03112792 -0.071273876 0.05838153 -0.07539587 1.000000000
Lag5 -0.03051910 -0.008183096 -0.07249948 0.06065717 -0.075675027
Volume 0.84194162 -0.064951313 -0.08551314 -0.06928771 -0.061074617
Today -0.03245989 -0.075031842 0.05916672 -0.07124364 -0.007825873Lag5 Volume Today
Year -0.030519101 0.84194162 -0.032459894
Lag1 -0.008183096 -0.06495131 -0.075031842
Lag2 -0.072499482 -0.08551314 0.059166717
Lag3 0.060657175 -0.06928771 -0.071243639
Lag4 -0.075675027 -0.06107462 -0.007825873
Lag5 1.000000000 -0.05851741 0.011012698
Volume -0.058517414 1.00000000 -0.033077783
Today 0.011012698 -0.03307778 1.000000000#绘制散点矩阵图
pairs(data[,-9])
该数据集共有1089行，9列。
从相关系数矩阵和散点矩阵图可以看出：滞后时间变量Lag1~Lag2之间没有显著性关系，但交易量Volume随时间不断有明显的增加。从交易量来看，自 90 年代以来，交易量显着增加。 这似乎在 2009 年左右达到顶峰，在 2010 年开始下降。
3.详细分析
#判断data中Lag1列往下移一行的数据与TOday列是否对应相等，从而判断数据是否按周增加
data %>% filter(lead(Lag1) != Today) %>% nrow()#按年分类并找出每年第一周的周序号
data$Week<-1:nrow(data)
Year_breaks<-data%>%group_by(Year)%>%summarise(Week=min(Week))#按周绘制交易量随时间的变化折线图
ggplot(data,aes(x=Week,y=Volume))+geom_line()+ #绘制折线图geom_smooth()+ #添加平滑趋势曲线theme_light() + #设置主题scale_x_continuous(breaks = Year_breaks$Week,minor_breaks = NULL,labels = Year_breaks$Year)+#如何按自己的意愿设置X轴的标签labs(title = "股票的每日交易量", x = "时间", y = "交易量")
#绘制堆积直方图
ggplot(data,aes(x=Year,fill=Direction))+ geom_bar(position = "fill")+geom_hline(yintercept = 0.5,col="black")+ #绘制y=0.5的水平参考线scale_x_continuous(breaks =seq(1990,2010),minor_breaks = NULL,labels = Year_breaks$Year )+scale_y_continuous(labels = scales::percent_format())+ #把y轴数值设为百分比制theme(axis.title.y =element_blank(),legend.position = "bottom")+ #取消y轴的标题ggtitle("每周的正回报与负回报情况")
这是成交回报随着时间变化的方向，似乎只有 4 年 >= 50% 的周没有看到正回报（分别是2000年、2001年、2002年、2008年）。
#分别计算出现Down和Up的概率
Week.probs<-prop.table(table(data$Direction))
Week.probs#绘制随时间变化的周波动
ggplot(data, aes(x = Week, y = Today/100 )) + #Today/100进行百分比化处理geom_line()+scale_x_continuous(breaks = Year_breaks$Week,minor_breaks = NULL,labels = Year_breaks$Year)+scale_y_continuous(labels = scales::percent_format(),breaks = seq(-0.2,0.2,0.05))+geom_hline(yintercept = 0,col="grey55")+ #绘制基准线theme_light()+labs(title = "每周百分比回报",x="时间",y="百分比回报")
可以看到，市场似乎经历了更高的变动/不稳定时期。 例如 2008 年 9 月在这里很突出。
3.回归分析
使用完整数据集执行逻辑回归，其中Direction作为响应变量，五个滞后变量加上Volume作为预测变量。
#进行logistic回归拟合
#把Direction变量由字符型转换成因子型
data$Direction <- as.factor(data$Direction)
glm.fit=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data = data,family = binomial)
summary(glm.fit)
Call:
glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, family = binomial, data = data)Deviance Residuals: Min 1Q Median 3Q Max
-1.6949 -1.2565 0.9913 1.0849 1.4579 Coefficients:Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.26686 0.08593 3.106 0.0019 **
Lag1 -0.04127 0.02641 -1.563 0.1181
Lag2 0.05844 0.02686 2.175 0.0296 *
Lag3 -0.01606 0.02666 -0.602 0.5469
Lag4 -0.02779 0.02646 -1.050 0.2937
Lag5 -0.01447 0.02638 -0.549 0.5833
Volume -0.02274 0.03690 -0.616 0.5377
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Dispersion parameter for binomial family taken to be 1)Null deviance: 1496.2 on 1088 degrees of freedom
Residual deviance: 1486.4 on 1082 degrees of freedom
AIC: 1500.4Number of Fisher Scoring iterations: 4tidy(glm.fit) #将模型输出结果转化为数据框
结果输出的z-statistic（z-value）的计算方法跟线性回归的T检验一样，z值的绝对值较大表明拒绝原假设H0：βj=0
由各个预测变量的P值可以看出，Lag2是显著性变量
4.计算混淆矩阵和整体预测率
#预测函数predict
#参数response告诉R只用输出概率P(Y=1|X)
#如果不给predict提供预测预测数据集，它会自动拟合logistic回归的训练数据的概率
glm.probs=predict(glm.fit,type = "response")
glm.probs[1:5] #这些值对应市场是上涨而不是下跌的概率1 2 3 4 5
0.6086249 0.6010314 0.5875699 0.4816416 0.6169013 #contrasts函数创建了一个哑变量
contrasts(data$Direction) UpDown 0Up 1
#将用logistic回归的训练数据的预测结果转化为变化方向
glm.pred=rep("Down",1089)
glm.pred[glm.probs>.5]='Up'
#计算预测结果与原来结果的混淆矩阵，从而计算预测一致的概率
attach(data)
table(glm.pred,Direction)Direction
glm.pred Down UpDown 54 48Up 430 557mean(glm.pred==Direction)
## [1] 0.5610652这略高于简单分类器（每次都预测 Up）所达到的基线准确度（55.56%）。
#用caret::confusionMatrix计算混淆矩阵
attach(data)
The following objects are masked from data:Direction, Lag1, Lag2, Lag3, Lag4, Lag5, Today, Volume, YearPredicted<-factor(ifelse(predict(glm.fit,type = "response")<.5,"Down","Up"))
confusionMatrix(Predicted,Direction,positive = "Up")
Confusion Matrix and StatisticsReference
Prediction Down UpDown 54 48Up 430 557Accuracy : 0.5611 95% CI : (0.531, 0.5908)No Information Rate : 0.5556 P-Value [Acc > NIR] : 0.369 Kappa : 0.035 Mcnemar's Test P-Value : <2e-16 Sensitivity : 0.9207 Specificity : 0.1116 Pos Pred Value : 0.5643 Neg Pred Value : 0.5294 Prevalence : 0.5556 Detection Rate : 0.5115 Detection Prevalence : 0.9063 Balanced Accuracy : 0.5161 'Positive' Class : Up prop.table(table(Predicted))
PredictedDown Up
0.09366391 0.90633609
用2009年之前的训练数据拟合logistic回归模型，其中只把Lag2作为预测变量，计算混淆矩阵的和测试集（2009和2010年）中总体预测准确率
#用1990-2008年的训练数据来拟合logistic回归模型，只把lag2作为预测变量，计算2009-2010的预测准确率
attach(data)
train=(Year<2009) #生成一个对应的布尔向量
#布尔向量可用于获取某个矩阵的行或子列
data.2009=data[!train,] #测试集数据
dim(data.2009) #查看该数据的维度
[1] 104 10glm.fit1<-glm(Direction~Lag2,data = data,family = binomial,subset = train)
summary(glm.fit1)
Call:
glm(formula = Direction ~ Lag2, family = binomial, data = data, subset = train)Deviance Residuals: Min 1Q Median 3Q Max
-1.536 -1.264 1.021 1.091 1.368 Coefficients:Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.20326 0.06428 3.162 0.00157 **
Lag2 0.05810 0.02870 2.024 0.04298 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1(Dispersion parameter for binomial family taken to be 1)Null deviance: 1354.7 on 984 degrees of freedom
Residual deviance: 1350.5 on 983 degrees of freedom
AIC: 1354.5Number of Fisher Scoring iterations: 4glm.probs1<-predict(glm.fit1,data.2009,type = "response")
glm.pred1<-rep("Down",104)
glm.pred1[glm.probs1>.5]="Up"
attach(data.2009)
table(glm.pred1,Direction)Direction
glm.pred1 Down UpDown 9 5Up 34 56
mean(glm.pred1==Direction)
[1] 0.625用caret::confusionMatrix计算混淆矩阵
> predicted1<-factor(ifelse(predict(glm.fit1,data.2009,type = "response")<.5,"Down","Up"))
> confusionMatrix(data = predicted1,reference=data.2009$Direction,positive = "Up")
Confusion Matrix and StatisticsReference
Prediction Down UpDown 9 5Up 34 56Accuracy : 0.625 95% CI : (0.5247, 0.718)No Information Rate : 0.5865 P-Value [Acc > NIR] : 0.2439 Kappa : 0.1414 Mcnemar's Test P-Value : 7.34e-06 Sensitivity : 0.9180 Specificity : 0.2093 Pos Pred Value : 0.6222 Neg Pred Value : 0.6429 Prevalence : 0.5865 Detection Rate : 0.5385 Detection Prevalence : 0.8654 Balanced Accuracy : 0.5637 'Positive' Class : Up
No Information Rate：0.5865，即表明测试数据的58.65%为最大分类（Up），因此这是分类器的一个基准线
Accuracy：0.625>0.5865
我们为单面测试提供了p值，以查看准确性是否优于“无信息率”。 P值[Acc> NIR]：0.2439> 0.05⟹没有明显的证据表明我们的分类器优于基线策略。
