一般化线性模型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. 线形预测子系数的含义