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 -
正規分布
-
指数型分布
-- 分布、逆分布、分布
2.1.2 GLM的概率模型
0821
- 目的变数的期待值就是说明变数的函数
例1. 一个说明变数,poisson分布,对数link
目的变数的期待值= 代入公式
y: 事件发生的次数, : 平均值
柏松分布特征:
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. 线形预测子系数的含义
网友评论