| ポアソン回帰におけるパラメータ推定 |
|
| 投稿者：kamo 2017/07/09(Sun) 16:15 No. 22415 |
お世話になります。
Rのoptim関数を用いて、ポアソン回帰（説明変数2つ)の係数を推定したいのですが、うまくいきません。
> x1 <- as.factor(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1))
> x2 <- as.factor(c(50,50,50,60,60,60,60,40,50,50,50,60,60,60,60,60,60,60,60,50,50,50,50,60,50,40,40,40,40,50,50,50,50,50,50,50,60,50,50,50,50))
> y <- c(0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,3,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1)
>
>
> logL <- function(par){
+ beta0 <- par[1]
+ beta1 <- par[2]
+ beta2 <- par[3]
+ sum(dpois(y,exp(beta0+beta1*x1+beta2*x2),log = TRUE))
+ }
>
> opt <- optim(c(1,1,1),logL,control = list(fnscale=-1),method = "BFGS",hessian = TRUE)
Error in optim(c(1, 1, 1), logL, control = list(fnscale = -1), method = "BFGS", :
initial value in 'vmmin' is not finite
In addition: Warning messages:
1: In Ops.factor(beta1, x1) : * not meaningful for factors
2: In Ops.factor(beta2, x2) : * not meaningful for factors
>
ご教示頂ければ幸いです。よろしくお願い致します。 |
| Re: ポアソン回帰におけるパラメータ推定 |
| 投稿者：青木繁伸 2017/07/09(Sun) 18:15 No. 22416 |
> x1 <- as.factor(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1))
> x2 <- as.factor(c(50,50,50,60,60,60,60,40,50,50,50,60,60,60,60,60,60,60,60,50,50,50,50,60,50,40,40,40,40,50,50,50,50,50,50,50,60,50,50,50,50)) x1, x2 が factor なので，その数値演算は無理です。
任意の beta0, beta1, beta2 とすると，> beta0 = 0
> beta1 = -0.1
> beta2 = -0.11
> beta0+beta1*x1+beta2*x2
[1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
警告メッセージ:
1: Ops.factor(beta1, x1) で: ‘*’ は因子に対しては無意味です
2: Ops.factor(beta2, x2) で: ‘*’ は因子に対しては無意味です
3: Ops.factor(beta1, x1) で: ‘*’ は因子に対しては無意味です
4: Ops.factor(beta2, x2) で: ‘*’ は因子に対しては無意味です
> exp(beta0+beta1*x1+beta2*x2)
[1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
警告メッセージ:
1: Ops.factor(beta1, x1) で: ‘*’ は因子に対しては無意味です
2: Ops.factor(beta2, x2) で: ‘*’ は因子に対しては無意味です
> dpois(y,exp(beta0+beta1*x1+beta2*x2),log = TRUE)
[1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
警告メッセージ:
1: Ops.factor(beta1, x1) で: ‘*’ not meaningful for factors
2: Ops.factor(beta2, x2) で: ‘*’ not meaningful for factors
> sum(dpois(y,exp(beta0+beta1*x1+beta2*x2),log = TRUE))
[1] NA
警告メッセージ:
1: Ops.factor(beta1, x1) で: ‘*’ not meaningful for factors
2: Ops.factor(beta2, x2) で: ‘*’ not meaningful for factors となるので，optim では求まりません。
独立変数が factor のとき，lm や glm などでは自動的にダミー変数が作成されるのですが，今回のような場合は自分でダミー変数を作って対処しないとだめでしょう。> # x1 は既にダミー変数
> x1 = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
> x2 = c(50,50,50,60,60,60,60,40,50,50,50,60,60,60,60,60,60,60,60,50,50,50,50,60,50,40,40,40,40,50,50,50,50,50,50,50,60,50,50,50,50)
> x2.50 = (x2 == 50)+0 # かっこは必要, +0 で logical を numeric に
> x2.60 = (x2 == 60)+0
> # x2 == 40 の場合は，x2.50 = 0, x2.60 = 0
>
> logL <- function(par){
+ sum(dpois(y,exp(par[1]+par[2]*x1+par[3]*x2.50+par[4]*x2.60),log = TRUE)) # par を直接使って問題ない
+ }
> opt <- optim(rep(1, 4),logL,control = list(fnscale=-1),method = "BFGS",hessian = TRUE)
> opt
$par
[1] -0.6495380 -1.4754486 0.9360592 0.3696557
$value
[1] -38.44917
$counts
function gradient
32 19
$convergence
[1] 0
$message
NULL
$hessian
[,1] [,2] [,3] [,4]
[1,] -28.998906 -4.0006148 -1.799968e+01 -9.999199e+00
[2,] -4.000615 -4.0006148 -3.350026e+00 -1.728501e-01
[3,] -17.999681 -3.3500262 -1.799968e+01 1.776357e-09
[4,] -9.999199 -0.1728501 1.776357e-09 -9.999199e+00
> # 検証
> par = opt$par
> exp(par[1]+par[2]*x1+par[3]*x2.50+par[4]*x2.60)
[1] 1.3317864 1.3317864 1.3317864 0.7558728 0.7558728 0.7558728 0.7558728 0.5222870 1.3317864 1.3317864 1.3317864 0.7558728 0.7558728
[14] 0.7558728 0.7558728 0.7558728 0.7558728 0.7558728 0.7558728 1.3317864 1.3317864 1.3317864 1.3317864 0.7558728 1.3317864 0.1194346
[27] 0.1194346 0.1194346 0.1194346 0.3045477 0.3045477 0.3045477 0.3045477 0.3045477 0.3045477 0.3045477 0.1728500 0.3045477 0.3045477
[40] 0.3045477 0.3045477
> dpois(y,exp(par[1]+par[2]*x1+par[3]*x2.50+par[4]*x2.60),log = TRUE)
[1] -1.3317864 -1.3317864 -1.3317864 -0.7558728 -0.7558728 -0.7558728 -0.7558728 -1.1718250 -1.0452652 -1.0452652 -1.0452652 -1.0357550
[13] -1.0357550 -1.0357550 -1.0357550 -1.0357550 -1.0357550 -1.0357550 -1.0357550 -1.4518912 -1.4518912 -1.4518912 -1.4518912 -2.0087844
[25] -2.2639823 -0.1194346 -0.1194346 -0.1194346 -0.1194346 -0.3045477 -0.3045477 -0.3045477 -0.3045477 -0.3045477 -0.3045477 -0.3045477
[37] -0.1728500 -1.4934752 -1.4934752 -1.4934752 -1.4934752
> sum(dpois(y,exp(par[1]+par[2]*x1+par[3]*x2.50+par[4]*x2.60),log = TRUE))
[1] -38.44917
|
| Re: ポアソン回帰におけるパラメータ推定 |
| 投稿者：kamo 2017/07/09(Sun) 22:03 No. 22417 |
詳細なご回答、誠にありがとうございます。
ダミー変数・デザイン行列まで作成してから、関数に放り込むべきであったのですね。
とてもよく理解できました。
ありがとうございました。 |
| Re: ポアソン回帰におけるパラメータ推定 |
| 投稿者：kamo 2017/07/10(Mon) 15:06 No. 22418 |
度々失礼します。
上記質問の続きとしまして、
「残差逸脱度Residual Deviance」の算出方法を知りたいです。
ご回答の最終行
> sum(dpois(y,exp(par[1]+par[2]*x1+par[3]*x2.50+par[4]*x2.60),log = TRUE))
[1] -38.44917
がこのモデルにおける対数尤度と理解しておりますが、
ここから残差逸脱度というものが計算できますでしょうか？
ご教示頂ければ幸いです。 |
| Re: ポアソン回帰におけるパラメータ推定 |
| 投稿者：青木繁伸 2017/07/11(Tue) 15:04 No. 22419 |
| Re: ポアソン回帰におけるパラメータ推定 |
| 投稿者：kamo 2017/07/15(Sat) 15:17 No. 22420 |
書籍等で調べると、先生のおっしゃるように残差逸脱度＝対数尤度の-2倍との記述がございました。
しかし、
> x1 <- as.factor(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1))
> x2 <- as.factor(c(50,50,50,60,60,60,60,40,50,50,50,60,60,60,60,60,60,60,60,50,50,50,50,60,50,40,40,40,40,50,50,50,50,50,50,50,60,50,50,50,50))
> y <- c(0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,3,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1)
> fm <- glm(y~x1+x2,family="poisson"(link="log"))
> summary(fm)
Call:
glm(formula = y ~ x1 + x2, family = poisson(link = "log"))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.6321 -0.7804 -0.3009 0.5385 1.2393
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.6495 1.0356 -0.627 0.53058
x11 -1.4757 0.5638 -2.617 0.00886 **
x250 0.9360 1.0405 0.900 0.36832
x260 0.3697 1.0805 0.342 0.73225
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 40.539 on 40 degrees of freedom
Residual deviance: 28.838 on 37 degrees of freedom
AIC: 84.898
Number of Fisher Scoring iterations: 5
Residual deviance: 28.838と出力され、（これが残差逸脱度ですよね？）
-2*(-38.44917) = 76.89833
と一致しません。
どこかで誤解しているでしょうか？
ご教示いただければ幸いです。 |
| Re: ポアソン回帰におけるパラメータ推定 |
| 投稿者：青木繁伸 2017/07/15(Sat) 21:55 No. 22421 |
glm でできるなら，わざわざ変なことして optim でやる必要ないでしょう
full model の 対数尤度は
> sum(dpois(y, lambda=y, log=TRUE)
[1] -24.03019
したがって，その -2 倍は full model の deviance
> sum(dpois(y, lambda=y, log=TRUE))*-2
[1] 48.06037
当該モデルの residual deviance が 28.838 なので，deviance は
> sum(dpois(y, lambda=y, log=TRUE))*-2+28.838
[1] 76.89837
逆にたどれば，当該モデルの deviance が 76.89837 で full model の deviance が48.06037 なので，residual deviance は
76.89837 - 48.06037 = 28.838 |
| Re: ポアソン回帰におけるパラメータ推定 |
| 投稿者：kamo 2017/07/16(Sun) 19:07 No. 22422 |
青木先生
ご回答ありがとうございます。
glm関数におけるポアソン回帰の計算手順を理解することが、このプログラミングの目的でした。お手数おかけして、すみません。
とてもよく理解することができました。ありがとうございました。
|
| Re: ポアソン回帰におけるパラメータ推定 |
| 投稿者：青木繁伸 2017/07/16(Sun) 19:14 No. 22423 |
> glm関数におけるポアソン回帰の計算手順を理解すること
そのような場合には，ソースコードを読むことをお勧めします
一点の曖昧さもなく，完璧に記述されていますから |
| Re: ポアソン回帰におけるパラメータ推定 |
| 投稿者：kamo 2017/07/17(Mon) 15:02 No. 22424 |
ご返信ありがとうございます。
ソースコードを読解するのは中々ハードルが高そうで、敬遠していましたが...
やはり、計算手順の理解には、その方法が確実ということですよね。
チャレンジしてみます。
またつまづくようなことがあったときは、この掲示板にお世話になるかと思いますが、今後ともご指導よろしくお願い致します。 |
|
