
デスクワークをしながら、久しぶりにRを開きました。
私は最初、「言語さえ理解できれば応用が利く」と思って、統計にはRを利用していたのですが、周りの同僚達は皆、ラボが購入している統計ソフトを使っているので、最近はRから少し遠ざかっていました。
上司から、臨床よりのプロジェクトを一つ任されていたので、自宅で細々解析していたのですが、ラボの統計ソフトでは解析できないものも多く、ちょっと使いづらい……。
とゆー経緯で、やっぱりRに戻ってしまったのでした。
しかしながら、久しぶりにRを開いてみると、R言語をイロイロ忘れている事に気づき、愕然としました。
これからは忘れないように、ここにメモっておこうと思います。
今回の備忘録は、ロジスティック回帰分析。
いつもは青木先生の「Rによる統計処理」にかなりお世話になっているのですが、
95%CIを出すときなどにちょっと不都合があったので(私の読解力不足だと思います。青木先生、ごめんなさい。)、今回は他のサイトをイロイロ参考にしながら自分で出してみました。
以下、備忘録。
Rでロジスティック回帰分析
まず、ポジティブの症例数が因子の数x10以上である事を確認。
Excelでデータをコピー&ペースト。
Rに読み込む。(1行目はタイトルの場合、header=TRUE, データの場合はデフォルトでFALSEになっている)
入力例: あるデータ(データ名:data1, n=629, pathology (+=1, -=0), Sex(F=1, M=0), Age)を読み込む場合。
data1<-read.table(“clipboard”,header=TRUE)
ちゃんと読み込めているか、最初の6行を表示する。
>head(data1)
出力例:
Sex Age pathology
1 0 90 1
2 0 87 1
3 0 102 1
4 0 70 1
5 1 90 1
6 1 80 1
Rに用意されているglm関数を用いて多重ロジスティック回帰分析を行う。
>res<- glm(data1$目的変数 (pathology)~., data=data1, family=binomial)
res<- glm(data1$pathology~., data=data1, family=binomial)
もしくは
>res<- glm(data1$目的変数 (pathology)~説明変数1 (Sex)+説明変数2 (Age), data=data1, family=binomial)
res<- glm(data1$pathology~Sex+Age, data=data1, family=binomial)
そしてサマリー出力。
>summary(res)
出力例:
Call:
glm(formula = data1$pathology ~ Sex + Age, family = binomial, data = data1)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5701 -0.9785 -0.7176 1.2256 1.8862
Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.911966 0.720208 -6.820 9.09e-12 *** Sex 0.480557 0.173798 2.765 0.00569 ** Age 0.052668 0.009036 5.829 5.59e-09 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 |
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 830.31 on 628 degrees of freedom
Residual deviance: 781.38 on 626 degrees of freedom
AIC: 787.38
Number of Fisher Scoring iterations: 4
※ 解析結果の読み方は,基本的には線型回帰分析の場合と同じであり,「Coefficients」(赤枠)に注目すると,切片と係数の推定値(Estimate)とその標準誤差(Std. Error),P値(Pr( > |z|))を出すためのZ値(z value)が表示されている。
>coefficients(res)
出力例:
(Intercept) Sex Age
-4.91196580 0.48055663 0.05266811
>residuals(res)
出力例:
1 2 3 4 5 6 7 8 9
1.2512006 1.3203070 0.9892391 1.7220360 1.0494316 1.2712623 1.3874283 1.4817758 1.3640056
.....
ロジスティック回帰分析の結果の解釈
※ 線型回帰直線の場合と異なり,ロジスティック回帰分析では,得られた切片と係数の推定値に「一手間」かけないと,結果の解釈ができない。得られた推定値はexp(推定値) と変換することによって,オッズ比になる。
例えば、上記のデータの場合、Sexの推定値を入れると、Sexのオッズ比が出力される。
> exp(0.480557)
出力例:
[1] 1.616975
Ageの推定値を入れると、Ageのオッズ比が出力される。
> exp(0.052668)
出力例:
[1] 1.05408
Alternative: 間違いが起こらないように、解析結果resから係数をとる方が良さげ。
まず、解析結果resには何がどこに入っているか見るため、namesを入力。
> names(res)
出力例:
[1] "coefficients" "residuals" "fitted.values" "effects"
[5] "R" "rank" "qr" "family"
[9] "linear.predictors" "deviance" "aic" "null.deviance"
[13] "iter" "weights" "prior.weights" "df.residual"
[17] "df.null" "y" "converged" "boundary"
[21] "model" "call" "formula" "terms"
[25] "data1" "offset" "control" "method"
[29] "contrasts" "xlevels"
1番目に係数が格納されている事がわかったので、[1]から切片と係数の推定値を出力。
> res[[1]]
出力例:
(Intercept) Sex Age
-4.91196580 0.48055663 0.05266811
> exp(res[[1]][2])
出力例:
Sex
1.616974
さらに、この中の1番目の切片と、2番目のSex、3番目のAgeの係数を取り出すには、それぞれ、
> res[[1]][1]
出力例:
(Intercept)
-4.911966
> res[[1]][2]
出力例:
Sex
0.4805566
> res[[1]][3]
出力例:
Age
0.05266811
となる。
SexとAgeの係数をそれぞれオッズ比に変換するには、
> exp(res[[1]][2])
出力例:
Sex
1.616974
> exp(res[[1]][3])
出力例:
Age
1.05408
いっぺんに出力する場合は、
> exp(res[[1]])
出力例:
(Intercept) Sex Age
0.00735801 1.61697421 1.05407975
以上から、補正後のオッズ比は、Sex: 1.61, Age: 1.05となる。
95% CIまで一度に出力する方法
resのサマリーをres1として出力する。
> res1<-summary(res)
> res1
出力例:
Call:
glm(formula = data1$Pathology ~ Sex + Age, family = binomial, data = data1)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5701 -0.9785 -0.7176 1.2256 1.8862
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.911966 0.720208 -6.820 9.09e-12 ***
Sex 0.480557 0.173798 2.765 0.00569 **
Age 0.052668 0.009036 5.829 5.59e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 830.31 on 628 degrees of freedom
Residual deviance: 781.38 on 626 degrees of freedom
AIC: 787.38
Number of Fisher Scoring iterations: 4
res1の内容の中から、coefficientを取り出す。
> coe<-res1$coefficient
> coe
出力例:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.91196580 0.720208376 -6.820201 9.091330e-12
Sex 0.48055663 0.173797974 2.765030 5.691754e-03
Age 0.05266811 0.009035981 5.828710 5.585755e-09
オッズ比とオッズ比の95%信頼区間を出力する。
- オッズ比(OR)
- オッズ比の95%信頼区間(max: ORhigh, min: ORlow)
> OR<-exp(coe[,1])
> ORlow<-exp(coe[,1]-1.96*coe[,2])
> ORhigh<-exp(coe[,1]+1.96*coe[,2])
後から確認しやすいように、3つをまとめておく。
> res2 <- cbind(coe,OR,ORlow,ORhigh)
> res2
出力例:
Estimate | Std. Error | z value | Pr(>|z|) | OR | ORlow | ORhigh | |
(Intercept) | -4.91196580 | 0.720208376 | -6.820201 | 9.091330e-12 | 0.00735801 | 0.001793522 | 0.03018659 |
Sex | 0.48055663 | 0.173797974 | 2.765030 | 5.691754e-03 | 1.61697421 | 1.150173272 | 2.27322758 |
Age | 0.05266811 | 0.009035981 | 5.828710 | 5.585755e-09 | 1.05407975 | 1.035575787 | 1.07291435 |
以上から、補正後のオッズ比と95%CIは、Sex (COR, 1.62; 95%CI, 1.15-2.27), Age (COR, 1.05; 95%CI, 1.04-1.07)となる。
自分用プログラム
上記の内容をまとめて、時短の為にプログラムを作成。
- 目的変数:Y
- 説明変数:X1, X2
として入力。
data1<-read.table(“clipboard”,header=TRUE)
head(data1)
としてデータを取り込んだ後で、
res<-glm(data1$Y~.,data=data1,family=binomial)
res1<-summary(res)
res1
coe<-res1$coefficient
coe
OR<-exp(coe[,1])
ORlow<-exp(coe[,1]-1.96*coe[,2])
ORhigh<-exp(coe[,1]+1.96*coe[,2])
res2<-cbind(coe,OR,ORlow,ORhigh)
res2
出力値の解釈:
- P値:Pr(>|z|)
- オッズ比:OR
- 95%CI:ORlowーORhigh
今回、参考にしたサイトはこちら ▼
関連本 ▼