广义线性回归模型

提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档


前言

提示: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.007825873
               Lag5      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.4

Number of Fisher Scoring iterations: 4

tidy(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) 
      Up
 Down  0
 Up    1
#将用logistic回归的训练数据的预测结果转化为变化方向
glm.pred=rep("Down",1089)
glm.pred[glm.probs>.5]='Up'
#计算预测结果与原来结果的混淆矩阵,从而计算预测一致的概率
attach(data)
table(glm.pred,Direction)
        Direction
glm.pred Down  Up
    Down   54  48
    Up    430 557
    
mean(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, Year

Predicted<-factor(ifelse(predict(glm.fit,type = "response")<.5,"Down","Up"))
confusionMatrix(Predicted,Direction,positive = "Up")
Confusion Matrix and Statistics

          Reference
Prediction Down  Up
      Down   54  48
      Up    430 557
                                         
               Accuracy : 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))
Predicted
      Down         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  10

glm.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.5

Number of Fisher Scoring iterations: 4

glm.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 Up
     Down    9  5
     Up     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 Statistics

          Reference
Prediction Down Up
      Down    9  5
      Up     34 56
                                         
               Accuracy : 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⟹没有明显的证据表明我们的分类器优于基线策略。

  • 0
    点赞
  • 11
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值