読者です 読者をやめる 読者になる 読者になる

Rでロジスティック回帰分析

R

Rでロジスティック回帰分析をしたので、おおまかな流れをメモ。
無駄にオブジェクトや命令文が増えているが、自分にとっての分かりやすさを重視で。

分析

分析自体はデフォルトの関数(glm)を用いるが、追加でpsclパッケージを用いた。psclに含まれるpR2()で、対数尤度やMacFaddenの疑似R二乗値を計算してくれる。

他にもDesignパッケージのlrm関数というのもモデルの検定なんかに便利らしいので、今後また調べるつもり。

library(pscl)			#今回必要なパッケージ

#関数本体。dfはデータフレーム、yは0か1の従属変数、x1〜xnが独立変数。
ans <- glm(y ~ x1+x2+xn, family=binomial(link="logit"), df)

#切片のみモデル
ans0 <- glm(y ~ 1, family=binomial(link="logit"), df)

summary(ans)			#結果を表示
exp(ans$coef)			#オッズ比
pR2(ans)			#尤度比統計(詳細はヘルプファイル)
anova(ans0, ans, test="Chisq")	#尤度比検定(-2対数尤度でカイ二乗検定)

データを整えてファイル出力

ついでに結果をまとめるためのデータ操作も書いておく。

library(Hmisc)				#ここで使うパッケージ

##データの整形
(temp <- summary(ans)$coe)		#結果のスコアをデータフレームに
(usb <- temp[,1:2])			#偏回帰係数と標準誤差
(or <- data.frame(exp(ans$coef)))	#オッズ比
(table <- data.frame(usb, or))		#表の統合
(table <- round(table, 3))		#表を小数点以下3桁で揃える
(logL <- round(pR2(ans), 3))		#尤度比統計も桁を揃える
summary(ans)			#有意確率など、スコアの確認用に出力

ここまでやれば、あとはExcelなりLaTeXなりで使えるように出力して、表を整える。サンプルサイズや有意確率も忘れずに。