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

No.05383　ロジスティック解析について　　【resident】　2008/01/15(Tue) 18:06

はじめまして。突然で申し訳ありませんが，質問させてください。
市中病院で研修医をしておりまして，ある疾患に対する治療薬に反応するかしないかを予測するために多重ロジスティック解析をしたいと考えているものです。説明変数としては，年齢，性別，発病から薬剤投与までの日数，各血液検査データです。

質問
ロジットの式を作ることによって，薬に反応するかしないかの確率は出せると思うのですが，日常の臨床で毎回あの式に代入するのは大変なので，もうすこし簡易化して，スコア化したいと考えています。たとえば，10歳以上なら1点，日数が5日以内なら1点，白血球数が20000以上なら1点というように設定し，その和の大小で結果を予想するという様にです。そのカットオフの値の決め方がわからないのです。カットオフ値をその値に設定したときのオッズ比が最大になるときということは感覚として理解できるのですが，そのような事もロジスティック解析できるのでしょうか？若しくは，別の手法により求めることができるのでしょうか？

素人なので用語などもつたなく，意味がよくわからないかもしれません。
もうしわけありませんが，よろしくおねがいします。

No.05384　Re: ロジスティック解析について　　【青木繁伸】　2008/01/15(Tue) 19:08

元データの独立変数をあるカットオフ値で区切って 0/1 データにして，それを用いて多重ロジスティックモデルを適用するとき，AIC が最小になるようなカットオフ値を求めるというようなことをやればよいわけで，R でやると，超簡単です。

http://aoki2.si.gunma-u.ac.jp/R/lr.data
のデータを例として，ちょっと以下のような関数を書いてみる。
func <- function(par) # 独立変数を二値データに変換してglmでaicを求める
{
	df$z1 <- df$x1 > par[1]
	df$z2 <- df$x2 > par[2]
	return(summary(glm(y~z1+z2, df, family=binomial))$aic)
}
par <- numeric(2)
par[1] <- mean(df$x1) # 初期値は平均値にでもしておく
par[2] <- mean(df$x2)
optim(par, func) # 最適解を求める
結果は
$par
[1] 148.0209 226.0251 # 2変数のカットオフ値

$value
[1] 74.90795 # AIC

$counts
function gradient 
      29       NA 

$convergence
[1] 0

$message
NULL
このあとどうしましょうか。
> df$z1 <- df$x1 > par[1]
> df$z2 <- df$x2 > par[2]
> summary(glm(y~z1+z2, df, family=binomial))

Call:
glm(formula = y ~ z1 + z2, family = binomial, data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7085  -0.6754  -0.3552  -0.3371   2.4066  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.7317     0.6374  -4.286 1.82e-05 ***
z1TRUE       -0.1074     0.6172  -0.174   0.8618    
z2TRUE        1.4774     0.6981   2.116   0.0343 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 76.714  on 97  degrees of freedom
Residual deviance: 71.443  on 95  degrees of freedom
AIC: 77.443

Number of Fisher Scoring iterations: 5
z1=1, z2=1 でも，合計点は-2.7317-0.1074+1.4774=-1.3617 でそれを確率に直すには 1/(1+exp(1.3617))=0.2039641みたいに計算しなくてはなりません。
小数点以下一桁くらいの足し算にしたり，代表的なλについて前もって1/(1+exp(λ))を計算しておくなどと，簡略化することはできるでしょうけど。

でも，それくらいの簡略化なら，元のデータから計算した係数をそのまま使ってエクセルか何かを使って計算するようにしておけば*正*確*に*値が求められますし，情報のロスもないわけですし。。。。

労多くして益なしではないでしょうか。

No.05419　Re: ロジスティック解析について　　【ありがとうございます！】　2008/01/17(Thu) 13:30

　青木先生，お答えいただきましてありがとうございます。
　先生のおしゃったように，以下のようにやってみましたが，うまくいきません。
　Rを初めて使うので，あまりに初歩的なミスがあるのでしょうが，
　それすらわかりません。日本語のマニュアルを読み漁りましたが，
　理解できませんでした。あと何が足りないのでしょうか。
　お忙しいところ恐縮ですが，お時間のあるときで結構ですので，お願いいたします。

> x <- read.table(stdin(), header=TRUE)
0: x1x2y
1: 2.913.90
2: 0.82.21
3: 0.43.91
4: 0.513.61
5: 0.718.51
6: 0.519.21
7: 0.311.31
8: 1.212.81
9: 1.93.81
10: 0.52.11
11: 0.41.71
12: 0.714.31
13: 0.614.31
14: 0.48.91
15: 1.29.11
16: 3.320.81
17: 0.581
18: 0.47.41
19: 0.419.81
20: 1.728.10
21: 0.37.61
22: 0.43.91
23: 0.94.51
24: 0.411.11
25: 0.44.11
26: 0.924.50
27: 0.913.31
28: 0.45.81
29: 0.95.31
30: 0.716.81
31: 0.78.31
32: 0.52.51
33: 0.37.51
34: 0.711.21
35: 1.47.30
36: 0.22.31
37: 0.53.41
38: 0.213.41
39: 0.59.71
40: 2.84.60
41: 0.6250
42: 2.64.60
43: 0.411.51
44: 0.36.51
45: 0.32.11
46: 0.43.61
47: 0.43.61
48: 0.58.91
49: 0.56.91
50: 0.815.21
51: > x
x1x2y
1 2.913.90
2 0.82.21
3 0.43.91
4 0.513.61
5 0.718.51
6 0.519.21
7 0.311.31
8 1.212.81
9 1.93.81
10 0.52.11
11 0.41.71
12 0.714.31
13 0.614.31
14 0.48.91
15 1.29.11
16 3.320.81
17 0.581
18 0.47.41
19 0.419.81
20 1.728.10
21 0.37.61
22 0.43.91
23 0.94.51
24 0.411.11
25 0.44.11
26 0.924.50
27 0.913.31
28 0.45.81
29 0.95.31
30 0.716.81
31 0.78.31
32 0.52.51
33 0.37.51
34 0.711.21
35 1.47.30
36 0.22.31
37 0.53.41
38 0.213.41
39 0.59.71
40 2.84.60
41 0.6250
42 2.64.60
43 0.411.51
44 0.36.51
45 0.32.11
46 0.43.61
47 0.43.61
48 0.58.91
49 0.56.91
50 0.815.21
51
> func <- function(par)
+ {
+ df$z1 <- df$x1 > par[1]
+ df$z2 <- df$x2 > par[2]
+ return(summary(glm(y~z1+z2, df, family=binomial))$aic)
+ }
> par <- numeric(2)
> par[1] <- mean(df$x1)
Warning message:
引数は数値でも論理値でもありません。NA 値を返します in: mean.default(df$x1)
> par[2] <- mean(df$x2)
Warning message:
引数は数値でも論理値でもありません。NA 値を返します in: mean.default(df$x2)
> optim(par, func)

　

No.05420　Re: ロジスティック解析について　　【青木繁伸】　2008/01/17(Thu) 14:19

入力行で，数値が区切られていないので，変数名が「 x1x2y」になり，数値も3つの数値がひとつながりになって「2.913.90」のようになっています（これでは数値と認識されませんね）。
また，データフレームは x としたようですね。私のプログラムでは df であると仮定しているので，データがちゃんと読めたとしても mean(df$x1) をしようとしたときに「そういうオブジェクトはない」という別のエラーメッセージになるでしょう。

まずは，データがちゃんと読めることが前提です。

No.05422　Re: ロジスティック解析について　　【resident】　2008/01/17(Thu) 19:32

　　恐縮です。こんな初歩的なことに答えていただき，本当にありがとうございます。
　そもそもデータフレームが3列の表になっていまいのですね。お恥ずかしい。さらには，ないデータフレームから計算しようとしていたとは・・・。引数を少し勉強して，結果をだせました。ありがとうございます。
　　簡易スコア化として，カットオフ値ごとの0，1のデータにして再度ロジスティック回帰分析をして，説明変数の個別のオッズ比に重み付けをして，ポイントを1点や2点とつけて和が何点以上なら感度・特異度はこれくらいと見積もろうかと思っています。

　本当にありがとうございました。