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

No.12367　ロジスティック回帰分析のエラー　　【雪】　2010/04/02(Fri) 13:33

20匹の動物の障害の有無と，あるタンパク質x1の量について，ロジスティック回帰分析を行っています。データは以下のような感じです。
左列…データフレームの行番号＝動物個体の番号
真ん中の列…障害の有無(1/0)
右列…x1の測定値　　　　　　です。

   Injury   x1
1   0   1.0122456
2   0   2.0412369
3   0   0.4611751
4   0   0.7499913
5   0   0.7353509
6   1   0.9793186
7   1   2.9590673
8   1   2.0849206
9   1   1.1778935
10  1   1.3341162
11  1   7.4915782
12  1   45.8988990
13  1   13.7614024
14  1   7.0526272
15  1   21.3706421
16  0   1.2741729
17  0   2.9629299
18  0   0.5090419
19  0   0.8374493
20  0   1.6199026

このデータについて
result <- glm(Injury ~ x1, data, family = binomial)
とすると，
In glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart, : 数値的に 0 か 1 である確率が生じました。
というエラーが帰ってきます。

Complete Separationに該当しているのかと思いましたが，6番目の個体を見ると，変数x1だけで完全に説明できていないと思います。
私の勉強不足なだけで，これもComplete SeparationやQuasi-Complete Separationの一例なのでしょうか？
12番目の個体のx1の値が大きすぎるせいのように見えて，これを小さい値に修正すると分析できるようになるので，Complete Separationに近い問題かとは思っているのですが…

よろしくお願いいたします。

No.12369　Re: ロジスティック回帰分析のエラー　　【波音】　2010/04/02(Fri) 16:36

理屈は分かりませんが，少なくともそのようなエラーが出る原因は12番目のデータにあるようです。試しに　12 1 45.8988990　という部分のデータのみを取り除いて解析してみれば，上手く行きます。
> inj
 [1] 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0
Levels: 0 1
> x1
 [1]  1.0122456  2.0412369  0.4611751  0.7499913  0.7353509  0.9793186
 [7]  2.9590673  2.0849206  1.1778935  1.3341162  7.4915782 13.7614024
[13]  7.0526272 21.3706421  1.2741729  2.9629299  0.5090419  0.8374493
[19]  1.6199026
> result <- glm(inj ~ x1, binomial)
> summary(result)

Call:
glm(formula = inj ~ x1, family = binomial)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4590  -0.7443  -0.6207   0.5459   1.6636  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -1.9527     1.1071  -1.764   0.0778 .
x1            0.8755     0.6170   1.419   0.1559

No.12370　Re: ロジスティック回帰分析のエラー　　【雪】　2010/04/02(Fri) 17:46

お返事ありがとうございます。

そうなんです。
12番目のデータを除くか，x1の値を36以下ぐらいにすると，解析できるようになります。しかし，その原因がよくわかりません。

今，怪しいと思っているのは，summary(result)の一番下の行です。
Number of Fisher Scoring iterations: 7

「Rによる医療統計学」P.183によると，"interationの数が10を超える場合には，glmでは当てはめを停止する。この閾値を自分で設定することも可能である"と書かれていますが，具体的な方法は記されていません。
?glmを見れば一発で分かりそうなのですが，ヘルプを見ても閾値の変え方が分からないのです。

もしご存知の方がおられれば，教えて頂けると嬉しいです。

No.12371　Re: ロジスティック回帰分析のエラー　　【青木繁伸】　2010/04/02(Fri) 17:52

エラーではなく，警告（worning）であるところに注意。

この警告は，glm が呼ぶ glm.fit の最後の方
        if (family$family == "binomial") {
            if (any(mu > 1 - eps) || any(mu < eps)) 
                warning("fitted probabilities numerically 0 or 1 occurred")
        }
で出しているのですね。eps は 2.22044604925031e-15，mu は予測値（ $fitted.values）
    list(coefficients = coef, residuals = residuals, fitted.values = mu, 以下略
で，何を意味しているかは文字通り予測値が1や0に非常に近いものがありますよ。ということです。「numerically 0 or 1」というのと「数値的に 0 か 1」というのは若干のニュアンスの違いがあるような（英文は，『数値としてほとんど 0 か 1 と解釈される』というような感じ。日本語は断定。そう感じるのは私だかかな？）。

で，扱いをどうするかですが，そのような意味なので無視してもよいでしょう。警告の意味は，こんなに1や0に近い値になるのは，独立変数のどれかが外れ値かも知れませんよということでしょう。実際にそのような値が観察されていることもあるだろうし，もしかしたらデータ入力に間違いがあるのかもしれないので，データを見直したらいかが？ということでしょう。今回はそう言うことはないのでしょうけど，もしかしたら 45.898899 ではなく，5.898899 ならこの警告は出ないでしょうからね。

よくわからないことがあるときには，ソースを見れば解決できることがありますね。そこが「自由なソフトウェア」のよいところです。

No.12372　Re: ロジスティック回帰分析のエラー　　【青木繁伸】　2010/04/02(Fri) 18:01

> ?glmを見れば一発で分かりそうなのですが，ヘルプを見ても閾値の変え方が分からないのです。

glm の引数に，control = glm.control(...) がありますよ。
説明を読むと，a list of parameters for controlling the fitting process. See the documentation for glm.control for details. ということで，glm.control という関数で設定できるとありますね。

そして，glm.control を見ると，maxit を設定することができるのが分かります。

> iteration の数が 10 を超える場合には，glm では当てはめを停止する

というのはうそ（？）で，デフォルトは 25 だということも分かります。

また，glm が返す結果も ? glm の下の方（Value という項）に書いてありますよ。そのなかに，収束したかどうかは converged でかえされる（logical. Was the IWLS algorithm judged to have converged?）ということもわかります。

オンラインヘルプを見るのが一番です。

No.12374　Re: ロジスティック回帰分析のエラー　　【雪】　2010/04/02(Fri) 19:36

青木先生

解説をありがとうございました。よく分かりました。
恥ずかしながら，options(warn=-1)で警告を切れることも，今まで知りませんでした。
とても助かりました。

No.12377　Re: ロジスティック回帰分析のエラー　　【芳賀】　2010/04/04(Sun) 12:49

　データをグラフ化して眺めると，そのままロジスティック曲線をあてはめるのは不自然です．x1を対数変換してロジスティック曲線を当てはめるときれいな曲線を描くことができます．
　確率が50%になるxの値は1.80，傾きの係数は1.909となります．
　なお，xの対数に対してロジスティック曲線があてはまるときは，通常のロジスティック曲線の式（1/(1+exp(-(a+bx))) ではなく，x^b/(x^b+x50^b) という式を用いるとp=0.5となるxが直接求められます．この式は，医薬品分野ではEmaxモデルなどと呼ばれています．

No.12382　Re: ロジスティック回帰分析のエラー　　【青木繁伸】　2010/04/04(Sun) 21:32

まあ，実質的な独立変数をどのように設定するかという，固有科学的な問題なのでしょうかね。