6

現在、HardinandHilbeの著書「GeneralizedLinearModelsandExtension」（第2版、2007年）を使用しています。著者は、OLSモデルの代わりに、「ログリンクは通常、連続スケールで正の値のみを取る応答データに使用される」と提案しています。もちろん、アイデンティティリンクを使用した「通常の」線形モデルを引き続き使用できるかどうかを確認するための残差プロットも提案しています。

私は彼らがSTATAの本で何をしているのかをRで再現しようとしています。確かに、私はログリンクでSTATAに問題はありません。ただし、Rのglm-functionを使用して同じモデルを呼び出すが、指定するfamily=gaussian(link="log")と、開始値を指定するように求められます。それらをすべてゼロに設定すると、アルゴリズムが収束しなかったというメッセージが常に表示されます。他の値を選択すると、メッセージが同じになる場合がありますが、多くの場合、次のようになります。

Error in glm.fit(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  :
     NA/NaN/Inf in 'x'

私が言ったように、STATAでは、開始値を設定したり、エラーを発生させたりすることなく、これらのモデルを実行できます。多くの異なるモデルと異なるデータセットを試しましたが、問題は常に同じです（単一の独立変数のみを含めない限り）。なぜこれが当てはまるのか、私が間違っているのか、または本から提案されたモデルが適切でない可能性があるのか​​、誰かに教えてもらえますか？助けていただければ幸いです、ありがとう！

編集：エラーを再現する例として、ここからダウンロードできるデータセットを考えてみましょう。このデータセットを読み込んだ状態で、次のモデルを実行します。

mod <- glm(betaplasma ~ age + vituse, family=gaussian(link="log"), data=data2, start=c(0,0,0))

これにより、アルゴリズムが収束しなかったという警告メッセージが生成されます。

Edit2：そのモデルのSTATA出力も提供するように依頼されました。ここにあります：

. glm betaplasma age vituse, link(log)

Iteration 0:   log likelihood = -2162.1385  
Iteration 1:   log likelihood = -2096.4765  
Iteration 2:   log likelihood = -2076.2465  
Iteration 3:   log likelihood = -2076.2244  
Iteration 4:   log likelihood = -2076.2244  

Generalized linear models                          No. of obs      =       315
Optimization     : ML                              Residual df     =       312
                                                   Scale parameter =  31384.51
Deviance         =  9791967.359                    (1/df) Deviance =  31384.51
Pearson          =  9791967.359                    (1/df) Pearson  =  31384.51

Variance function: V(u) = 1                        [Gaussian]
Link function    : g(u) = ln(u)                    [Log]

                                                   AIC             =  13.20142
Log likelihood   = -2076.224437                    BIC             =   9790173

------------------------------------------------------------------------------
             |                 OIM
  betaplasma |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         age |   .0056809   .0032737     1.74   0.083    -.0007354    .0120972
      vituse |   -.273027   .0650773    -4.20   0.000    -.4005762   -.1454779
       _cons |   5.467577   .2131874    25.65   0.000     5.049738    5.885417
------------------------------------------------------------------------------
4

2 に答える 2

13

コメントで述べたように、StataがRよりも（統計的な意味ではなく数値的に）GLMフィッティングを備えていることはおそらく真実です。とはいえ、この特定のデータセットのフィッティングはそれほど難しくはないようです。

データの読み取り：

data2 <- read.table("http://lib.stat.cmu.edu/datasets/Plasma_Retinol",
         skip=30,nrows=315)
dnames <- c("age","sex","smokstat","quetelet","vituse","calories","fat","fiber",
           "alcohol","cholesterol","betadiet","retdiet","betaplasma","retplasma")
names(data2) <- dnames

データをプロットします。

par(mfrow=c(1,2),las=1,bty="l")
with(data2,plot(betaplasma~age))
with(data2,boxplot(betaplasma~vituse))

ここに画像の説明を入力してください

切片パラメーターの開始値を妥当な値（つまり、対数スケールのデータの平均に近い値：これらのいずれかが機能するもの）に設定することで、これらを適合させるのはかなり簡単です。

mod <- glm(betaplasma ~ age + vituse, family=gaussian(link="log"), data=data2,
           start=c(10,0,0))
mod <- glm(betaplasma ~ age + vituse, family=gaussian(link="log"), data=data2,
           start=c(log(mean(data2$betaplasma)),0,0))

後者の場合は、おそらくログリンクフィットを開始するための合理的なデフォルト戦略です。結果（少し省略）は、Stataの結果と非常によく一致しています。

summary(mod)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.467575   0.218360  25.039  < 2e-16 ***
## age          0.005681   0.003377   1.682   0.0935 .  
## vituse      -0.273027   0.065552  -4.165 4.03e-05 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
## 
## (Dispersion parameter for gaussian family taken to be 31385.26)
## 
##     Null deviance: 10515638  on 314  degrees of freedom
## Residual deviance:  9791967  on 312  degrees of freedom
## AIC: 4160.4
## 
## Number of Fisher Scoring iterations: 9

confint(mod)
##                     2.5 %      97.5 %
## (Intercept)  5.0364648709  5.87600710
## age         -0.0007913795  0.01211007
## vituse      -0.4075213916 -0.14995759

（Rはp値と（？）信頼区間にZ統計ではなくtを使用しています）

ただし、このモデルをこれらのデータに適合させない理由はいくつかあります。特に、（ガウスモデルに関連する）一定の分散の仮定はあまり合理的ではありません。これらのデータは、対数正規モデル（または、標準のガウスモデルを使用した対数変換と分析）に適しているようです。

スケールでのプロットlog(1+x)（データにはゼロエントリがあります）：

with(data2,plot(log(1+betaplasma)~age))
with(data2,boxplot(log(1+betaplasma)~vituse))

ここに画像の説明を入力してください

プロット（これは、加法モデルをフィッティングするのではなく、ggplotの値ごとに別々の線にフィッティングします）vituse

library(ggplot)
theme_set(theme_bw())
(g1 <- qplot(age,1+betaplasma,colour=factor(vituse),data=data2)+
    geom_smooth(method="lm")+
    scale_y_log10())

ここに画像の説明を入力してください

'外れ値'なしで表示：

g1 %+% subset(data2,betaplasma>0)

ここに画像の説明を入力してください

他の2つのポイント：（1）このデータセットに値0の応答があるのは少し奇妙です。不可能ではありませんが、奇妙です。vituse（2）数値ではなく因子として扱われるべきであるように見えます（「1 =はい、かなり頻繁に、2 =はい、あまり頻繁ではない、3 =いいえ」）-おそらく序数。

于 2012-11-27T14:16:56.150 に答える
5

正常ではない可能性のあるエラーである可能性があることをお勧めします。同意する場合（またはデータが同意する場合）は、次の構造を検討してください。

?family
?glm
?binomial
lfit <- glm( dep <- indep1 + indep2, data=dat, family=binomial(link="probit")

これにより、アイデンティティインクモデルの周りに二項誤差が生じるはずです。これの利点は、変数の元のスケールで推定を解釈しやすくなることです。family = poissonプロビットリンクで使用する以前の誤った提案についてお詫びします。データや分布の説明さえ提示したことがないことを忘れないでください。明らかに、二項誤差は@BenBolkerが提供するデータセットには適切ではありません。

対数正規分布のエラーを伴う非整数値がある場合は、準ポアソンモデルを検討する必要があります。Ben Bolkerが提示したデータでこのモデルを実行し、gaussian（link = "log）モデルを比較すると、ほとんど区別がつかず、開始値は必要ありません。

> mod2 <- glm(betaplasma ~ age + vituse, family=quasipoisson, data=data2         )
> mod2

Call:  glm(formula = betaplasma ~ age + vituse, family = quasipoisson, 
    data = data2)

Coefficients:
(Intercept)          age       vituse  
   5.452014     0.006096    -0.276679  

Degrees of Freedom: 314 Total (i.e. Null);  312 Residual
Null Deviance:      37270 
Residual Deviance: 33420    AIC: NA 

> glm(betaplasma ~ age + vituse, family=gaussian(link="log"), data=data2,
+            start=c(10,0,0))

Call:  glm(formula = betaplasma ~ age + vituse, family = gaussian(link = "log"), 
    data = data2, start = c(10, 0, 0))

Coefficients:
(Intercept)          age       vituse  
   5.467575     0.005681    -0.273027  

Degrees of Freedom: 314 Total (i.e. Null);  312 Residual
Null Deviance:      10520000 
Residual Deviance: 9792000  AIC: 4160

vituseは明らかに3つのレベルの要素であるため、おそらくもう少し複雑なモデルを使用する必要があります。

> mod2 <- glm(betaplasma ~ age + factor(vituse), family=quasipoisson, data=data2         )
> mod2

Call:  glm(formula = betaplasma ~ age + factor(vituse), family = quasipoisson, 
    data = data2)

Coefficients:
    (Intercept)              age  factor(vituse)2  factor(vituse)3  
       5.151076         0.006359        -0.224107        -0.562727  

Degrees of Freedom: 314 Total (i.e. Null);  311 Residual
Null Deviance:      37270 
Residual Deviance: 33380    AIC: NA
于 2012-11-26T16:36:52.033 に答える