「統計学関連なんでもあり」の過去ログ--- 048

No.22415　ポアソン回帰におけるパラメータ推定　　【kamo】　2017/07/09(Sun) 16:15

お世話になります。
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

ご教示頂ければ幸いです。よろしくお願い致します。

No.22416　Re: ポアソン回帰におけるパラメータ推定　　【青木繁伸】　2017/07/09(Sun) 18:15

> 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

No.22417　Re: ポアソン回帰におけるパラメータ推定　　【kamo】　2017/07/09(Sun) 22:03

詳細なご回答，誠にありがとうございます。
ダミー変数・デザイン行列まで作成してから，関数に放り込むべきであったのですね。
とてもよく理解できました。
ありがとうございました。

No.22418　Re: ポアソン回帰におけるパラメータ推定　　【kamo】　2017/07/10(Mon) 15:06

度々失礼します。
上記質問の続きとしまして，
「残差逸脱度Residual Deviance」の算出方法を知りたいです。

ご回答の最終行
> sum(dpois(y,exp(par[1]+par[2]*x1+par[3]*x2.50+par[4]*x2.60),log = TRUE))
[1] -38.44917
がこのモデルにおける対数尤度と理解しておりますが，
ここから残差逸脱度というものが計算できますでしょうか？

ご教示頂ければ幸いです。

No.22419　Re: ポアソン回帰におけるパラメータ推定　　【青木繁伸】　2017/07/11(Tue) 15:04

「対数尤度の -2 倍がデビアンス」では？

No.22420　Re: ポアソン回帰におけるパラメータ推定　　【kamo】　2017/07/15(Sat) 15:17

書籍等で調べると，先生のおっしゃるように残差逸脱度＝対数尤度の-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と一致しません。
どこかで誤解しているでしょうか？
ご教示いただければ幸いです。

No.22421　Re: ポアソン回帰におけるパラメータ推定　　【青木繁伸】　2017/07/15(Sat) 21:55

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

No.22422　Re: ポアソン回帰におけるパラメータ推定　　【kamo】　2017/07/16(Sun) 19:07

青木先生
ご回答ありがとうございます。
glm関数におけるポアソン回帰の計算手順を理解することが，このプログラミングの目的でした。お手数おかけして，すみません。
とてもよく理解することができました。ありがとうございました。

No.22423　Re: ポアソン回帰におけるパラメータ推定　　【青木繁伸】　2017/07/16(Sun) 19:14

> glm関数におけるポアソン回帰の計算手順を理解すること

そのような場合には，ソースコードを読むことをお勧めします
一点の曖昧さもなく，完璧に記述されていますから

No.22424　Re: ポアソン回帰におけるパラメータ推定　　【kamo】　2017/07/17(Mon) 15:02

ご返信ありがとうございます。
ソースコードを読解するのは中々ハードルが高そうで，敬遠していましたが...
やはり，計算手順の理解には，その方法が確実ということですよね。
チャレンジしてみます。
またつまづくようなことがあったときは，この掲示板にお世話になるかと思いますが，今後ともご指導よろしくお願い致します。