美文网首页
一般化线性模型GLM学习笔记

一般化线性模型GLM学习笔记

作者: Jason数据分析生信教室 | 来源:发表于2021-06-26 23:20 被阅读0次

Chapter 2 一般化线性模型GLM及其要素

説明変数 x (objective variable) , 目的変数 y(response)

  • 数据分析根据说明变数x是量的变数还是质的变数的不同,用法而不同。
  • 根据目的变数分散的不同,用法也不同。比如,方差分析和直线回归就要求是等分散的正规分布。如果非等分散,就得用到非参数分析,比如Mann-Whitney

-- GLM可以看破说明变数对目的变数的影响的分析方法的结构。不管说明变数,目的变数是质是量都无所谓。

-- GLM覆盖了传统的回归分析,应用范围更广。

  • リンク関数(目的変数の予測値)可以让数据保持自然,比方说让计数值永远为正,让百分比值永远在0到1之间。
    -- log(目的変数の予測値)
    -- log{(目的変数の予測値/(1-目的変数の予測値)}

1.1 GLM的3要素

  • リンク関数(目的変数の予測値)=線形予測値
  • 目的変数の予測値=リンク関数の逆関数(線形予測値)
    -- log(目的変数の予測値)=線形予測子
    -- 目的変数の予測値=exp(線形予測子)
  • 三要素
    -- 線形予測子
    -- リンク関数
    -- 誤差構造

1.1.2 リンク函数

  • 目的変数の予測値=線形予測子
    -- identity(恒等)リンク

1.1.3 誤差構造

  • 计数数据: 平均值越大,分散越大
  • 百分比数据 : 0,1两端的分散比较小,0.5附近的分散比较大。
  • GLM根据数据分布特征选择误差构造。

1.2 GLM和以往方法的关系

- 误差构造 关联函数
直线回归 正规分布 identity
重回归 正规分布 identity
分散分析 正规分布 identity
共分散分析 正规分布 identity
ロジスティック回帰 二項分布 ロジット
ブロビット分析 二項分布 ブロビット
ボアソン回帰 ボアソン分布 対数
対数線形モデル ボアソン分布 対数

1.3 概率函数和GLM

  • i.i.d. independent and identically distributed
  • 说明变数固定不变,目的变数的概率模型。

Chapter 2 最尤法和一般线性模型

2.1 概率分布和一般化线形模型

2.1.1 经常使用的例子

0820

  • 二項分布
    --有/无,生存/死亡,雄/雌

  • ポアソン分布
    -- 单位时间里事件发生的概率相同。某单位时间长度里发生y次事件的概率分布
    -- 平均=分散

  • 指数分布(负的指数分布)
    -- 从某个时间点发生的事件,单位时间里终止的概率一定。
    -- 事件持续时间的概率密度
    -- 平均=1/γ, 分散=1/γ2,変動係数=1

  • 正規分布

  • 指数型分布
    --\gamma 分布、逆\gamma分布、\beta分布

2.1.2 GLM的概率模型

0821

  • 目的变数的期待值就是说明变数的函数

例1. 一个说明变数,poisson分布,对数link
目的变数的期待值= e^ {\beta_0+\beta_1x_1}代入公式


y: 事件发生的次数,\gamma : 平均值
柏松分布特征:
1.y=0,1,2,3,4,5....的概率总和为1
2.概率由且仅有期待值决定。

--例3. 两个说明变数,正规分布,identity link

2.2 尤度和最尤法

2.2.1 尤度

概率分布的参数决定以后,就可以推测出某个数据值的概率。
反过来,给出一堆数据,参数一定的情况下,能得到给出数据的概率,就是尤度(likelihood)

  • 尤度


    柏松尤度

2.2.2 尤度的最大化和最尤推定

  • 对数尤度


    对数尤度
  • 对数尤度偏微分(Socre) Score=0的偏微分求解


    对数尤度偏微分
  • 最尤推定


    最尤推定

总结一下就是对概率密度函数进行偏微分,求解socre为0时的方程就是期待值?
最尤推定。

2.2.3 最尤推定值的性质

  • 样本尺寸越大,最尤推定值的分布就越接近正规分布。

利用这个性质,不仅可以推定得到最尤值,还可以计算信赖空间。

2.2.4 最尤法和检定

  • Wald test: z value= Estimate/ Std.error,评价Estimate值是不是有意义
  • Score test : 评价尤度有没有到最大
  • 尤度比检定: 与归无模型的尤度进行比较。评价模型与数据的吻合度
  • Deviance: 2x{(无制约模型对数尤度)-(有制约模型的对数尤度)} 模型与数据的吻合度
    --对数尤度差的分布和卡方分布类似。(自由度就是被制约的参数量)

其他: Residual Deviance是个啥

  • Deviance: D = −2logL∗(最大化对数尤度)
  • Residual Deviance: D − (理论最小 D)

理论最小D:Full model的D

2.2.5 最尤法和推定

  • Wald test: 最尤推定值的点推和区间推。平均0,分散1

2.2.6 预测和AIC

  • AIC =-2 x Deviance + 2 × (参数量)
  • AIC 越小,预测精度越高

2.4 R和glm函数

2.4.1 glm函数

x1<-c(6.5,3.8,3.4,2.4,3.0,5.5,2.4,6.6)
x2<-c(3.7,4.9,1.0,1.8,4.6,4.8,3.8,2.7)
y1<-c(8,5,2,0,1,11,4,9)
# x1, x2 说明变数, y1应答变数
> glm(formula=y1~x1+x2,family=poisson(link="log"))

Call:  glm(formula = y1 ~ x1 + x2, family = poisson(link = "log"))

Coefficients:
(Intercept)           x1           x2  
    -1.3277       0.4039       0.2770  

Degrees of Freedom: 7 Total (i.e. Null);  5 Residual
Null Deviance:      26.78 
Residual Deviance: 7.085    AIC: 36.69

link(目的变数的期待值)= -1.3277+0.4039x1+0.2770x2
log(目的变数的期待值)= -1.3277+0.4039x1+0.2770x2
目的变数的期待值= exp(-1.3277+0.4039x1+0.2770x2)

  • dataframe
> data_gp01<- as.data.frame(cbind(x1,x2,y1))
> resgp01_12<- glm(formula = y1~x1+x2,family=poisson(link="log"),data=data_gp01)
> resgp01_12

Call:  glm(formula = y1 ~ x1 + x2, family = poisson(link = "log"), data = data_gp01)

Coefficients:
(Intercept)           x1           x2  
    -1.3277       0.4039       0.2770  

Degrees of Freedom: 7 Total (i.e. Null);  5 Residual
Null Deviance:      26.78 
Residual Deviance: 7.085    AIC: 36.69
> summary(resgp01_12)

Call:
glm(formula = y1 ~ x1 + x2, family = poisson(link = "log"), data = data_gp01)

Deviance Residuals: 
       1         2         3         4         5         6         7         8  
-0.71757   0.09939   0.49348  -1.51703  -1.43278   0.56161   1.24101   0.32697  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.3277     0.9081  -1.462 0.143716    
x1            0.4039     0.1084   3.726 0.000195 ***
x2            0.2770     0.1561   1.775 0.075949 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 26.777  on 7  degrees of freedom
Residual deviance:  7.085  on 5  degrees of freedom
AIC: 36.688

Number of Fisher Scoring iterations: 5

模型的预测值

> fitted(resgp01_12)
        1         2         3         4         5         6         7         8 
10.204733  4.781036  1.380816  1.150689  3.184833  9.240961  2.002539  8.054392 
> predict(resgp01_12,type="response")
        1         2         3         4         5         6         7         8 
10.204733  4.781036  1.380816  1.150689  3.184833  9.240961  2.002539  8.054392 
> predict(resgp01_12,type="link")
        1         2         3         4         5         6         7         8 
2.3228516 1.5646572 0.3226747 0.1403607 1.1583999 2.2236459 0.6944161 2.0862176 
> x1<- c(1,2)
> x2<- c(3,4)
> new_data1<- as.data.frame(cbind(x1,x2))
> predict(resgp01_12,newdata=new_data1,type="link")
          1           2 
-0.09271676  0.58824713
  • 残差
> # residuals
> residuals(resgp01_12,type="response")
         1          2          3          4          5          6          7          8 
-2.2047332  0.2189643  0.6191839 -1.1506887 -2.1848333  1.7590385  1.9974606  0.9456077
  • Residuals Deviance
    -- Deviance: D = −2logL∗(最大化对数尤度)
    -- Residual Deviance: D − (理论最小 D)
> residuals(resgp01_12)
          1           2           3           4           5           6           7           8 
-0.71756855  0.09939095  0.49347962 -1.51702916 -1.43278283  0.56161229  1.24100973  0.32697208 
  • Pearson残差:(目的变数的实测值-预测值)/分散的平方根
> residuals(resgp01_12,type="pearson")
         1          2          3          4          5          6          7          8 
-0.6901686  0.1001411  0.5269286 -1.0727016 -1.2242636  0.5786511  1.4115221  0.3331920 
  • 默认link

不指定link函数

> glm(y1~x1+x2,family=poisson,data=data_gp01)

Call:  glm(formula = y1 ~ x1 + x2, family = poisson, data = data_gp01)

Coefficients:
(Intercept)           x1           x2  
    -1.3277       0.4039       0.2770  

Degrees of Freedom: 7 Total (i.e. Null);  5 Residual
Null Deviance:      26.78 
Residual Deviance: 7.085    AIC: 36.69

统计检定

-- Wald test

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.3277     0.9081  -1.462 0.143716    
x1            0.4039     0.1084   3.726 0.000195 ***
x2            0.2770     0.1561   1.775 0.075949 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

--尤度比检定

> anova(resgp01_12,test="LRT")
Analysis of Deviance Table

Model: poisson, link: log

Response: y1

Terms added sequentially (first to last)


     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                     7     26.777              
x1    1  16.2114         6     10.566 5.665e-05 ***
x2    1   3.4807         5      7.085   0.06209 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

但是这个检定只检定了x1,x1+x2的情况。写formula的时候反一下就好。

> resgp01_21<- glm(y1~x2+x1, family = poisson(link="log"),data=data_gp01)
> resgp01_21

Call:  glm(formula = y1 ~ x2 + x1, family = poisson(link = "log"), data = data_gp01)

Coefficients:
(Intercept)           x2           x1  
    -1.3277       0.2770       0.4039  

Degrees of Freedom: 7 Total (i.e. Null);  5 Residual
Null Deviance:      26.78 
Residual Deviance: 7.085    AIC: 36.69
> summary(resgp01_21)

Call:
glm(formula = y1 ~ x2 + x1, family = poisson(link = "log"), data = data_gp01)

Deviance Residuals: 
       1         2         3         4         5         6         7         8  
-0.71757   0.09939   0.49348  -1.51703  -1.43278   0.56161   1.24101   0.32697  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.3277     0.9081  -1.462 0.143716    
x2            0.2770     0.1561   1.775 0.075949 .  
x1            0.4039     0.1084   3.726 0.000195 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 26.777  on 7  degrees of freedom
Residual deviance:  7.085  on 5  degrees of freedom
AIC: 36.688

Number of Fisher Scoring iterations: 5
  • drop1()
    -- 指定模型 vs 删除一个说明变量的模型比较

0822

> drop1(resgp01_12,test="Chisq")
Single term deletions

Model:
y1 ~ x1 + x2
       Df Deviance    AIC     LRT  Pr(>Chi)    
<none>       7.085 36.688                      
x1      1   22.802 50.405 15.7171 7.355e-05 *** #(x1+x2) vs x2 =  删了x1
x2      1   10.566 38.168  3.4807   0.06209 .  #(x1+x2) vs x1 = 删了x2
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

不指定算法的话

> drop1(resgp01_12)
Single term deletions

Model:
y1 ~ x1 + x2
       Df Deviance    AIC
<none>       7.085 36.688
x1      1   22.802 50.405
x2      1   10.566 38.168
  • add1()
    与drop1相反,在指定模型的基础上添加变量
> add1(resgp01_0,~x1,test="Chisq")
Single term additions

Model:
y1 ~ 1
       Df Deviance    AIC    LRT  Pr(>Chi)    
<none>      26.777 52.380                     
x1      1   10.566 38.168 16.211 5.665e-05 ***  #常数项 vs x1
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> add1(resgp01_0,~x1+x2,test="Chisq")
Single term additions

Model:
y1 ~ 1
       Df Deviance    AIC    LRT  Pr(>Chi)    
<none>      26.777 52.380                     
x1      1   10.566 38.168 16.211 5.665e-05 ***   #常数项 vs x1
x2      1   22.802 50.405  3.975   0.04618 *      #常数项 vs x2
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

对数尤度和Deviance

> logLik(resgp01_12)
'log Lik.' -15.34385 (df=3)    #对数尤度
> deviance(resgp01_12)         #Deviance
[1] 7.084974

Null Deviance: 2 x (仅含有切片的模型对数尤度- Full Model对数尤度)
Residual deviance: 2 x (候选模型对数尤度- Full Model对数尤度)

AIC

> AIC(resgp01_12)
[1] 36.68771

> stepAIC(resgp01_12)
Start:  AIC=36.69
y1 ~ x1 + x2

       Df Deviance    AIC
<none>       7.085 36.688
- x2    1   10.566 38.168
- x1    1   22.802 50.405

Call:  glm(formula = y1 ~ x1 + x2, family = poisson(link = "log"), data = data_gp01)

Coefficients:
(Intercept)           x1           x2  
    -1.3277       0.4039       0.2770  

Degrees of Freedom: 7 Total (i.e. Null);  5 Residual
Null Deviance:      26.78 
Residual Deviance: 7.085    AIC: 36.69

指定线形预测子

> glm(y1~1,family = poisson(link="log"),data = data_gp01) ##仅含有常数项

Call:  glm(formula = y1 ~ 1, family = poisson(link = "log"), data = data_gp01)

Coefficients:
(Intercept)  
      1.609  

Degrees of Freedom: 7 Total (i.e. Null);  7 Residual
Null Deviance:      26.78 
Residual Deviance: 26.78    AIC: 52.38
> glm(y1~x1-1,family = poisson(link="log"),data = data_gp01) ##无常数项

Call:  glm(formula = y1 ~ x1 - 1, family = poisson(link = "log"), data = data_gp01)

Coefficients:
   x1  
0.348  

Degrees of Freedom: 8 Total (i.e. Null);  7 Residual
Null Deviance:      91.53 
Residual Deviance: 10.75    AIC: 36.35

2.4.2 线形预测子

A. formula可以是函数
0826

> glm(exp(Temp/Wind)~Day, family=gaussian(link="identity"),data=airquality)

Call:  glm(formula = exp(Temp/Wind) ~ Day, family = gaussian(link = "identity"), 
    data = airquality)

Coefficients:
(Intercept)          Day  
 -4.941e+16    1.412e+16  

Degrees of Freedom: 152 Total (i.e. Null);  151 Residual
Null Deviance:      6.734e+38 
Residual Deviance: 6.71e+38     AIC: 13350

上面的例子里应答变数就是exp(Temp/Wind),说明变数是Day

B. offset
offset的系数是1

C. 名义尺度的说明变数
0或者1的ダミー变数

2.4.3 说明变量不能线性相关

2.4.5 Link Function

Chapter 3 目的变数为离散数据的GLM

3.1 以比例为目的变数的分析

  • oddsオッズ
    level=2的时候,就是p 与 1-p。 odds:p/(1-p)
  • 对数odds
    log{p/(1-p)}=logp - log(1-p)
    p=0.5的话,对数odds就是0。
    odds越大,对数odds越大。
    odds永远为正值。对数odds负无穷到正无穷,没有限制。
  • odds rattio odds比
    两个条件的odds比,概率p1,p2
odds odds ratio
p1/(1-p1) p1(1-p2)/(1-p1)p2
p2/(1-p2) p2(1-p1)/(1-p2)p1

如果两个条件的odds没有差的话,odds比就是1,对数odds比是0.

3.1.1 罗辑回归 Logistic Regression

  • 误差构造: 二项分布
  • link函数: Logit
  • y是比例值说明变数,y的期待值的logit= 线形预测子

log(E[y]/(1-E[y]))=ax1+b
A. 线形预测子系数的含义

相关文章

网友评论

      本文标题:一般化线性模型GLM学习笔记

      本文链接:https://www.haomeiwen.com/subject/dbjtultx.html