ロジスティック回帰(2)

ロジスティック回帰 (1)の続きです。例題はロジスティック回帰 (1)と同じですが、データセットの構造が異なります。例題の内容や結果については色々とご意見があると思いますが、統計学の練習のためのサンプルですのでご了承ください。

参考記事 ↑↑↑
例題
リハビリを受けている入院中の患者に対してリハビリに関する満足度調査を実施した。1週間に割り当てられた単位数をもとにランダムにインタビューした。「満足」と回答する割合と1週間のリハビリ単位数には関係があるのでしょうか?

データの準備と要約

データセット log01log02 をダウンロードします(ファイルの読み込み方

dat1 <- read.csv("log01.csv", header=T, fileEncoding = "UTF-8")
dat2 <- read.csv("log02.csv", header=T, fileEncoding = "UTF-8")

dat1の確認

dat1

単位数=1週間のリハ単位数
患者数=各単位数に該当する患者数
満足=「満足」と回答した患者数

dat2の確認

dim(dat2)

100行あるので20行のみを表示してみます

#20行のみ表示
dat2[0:20,]

ロジスティック回帰

それぞれのセットでロジスティック回帰を実行してみます

目的変数の設定が異なるので注意してください

fit1: dat1 を使用したロジスティック回帰

fit1 <- glm(
    cbind(満足, 患者数-満足) ~ 単位数, 
    family=binomial, 
    data=dat1
)
summary(fit1)

fit2: dat2 を使用したロジスティック回帰

fit2 <- glm(
    cbind(満足, 1-満足) ~ 単位数, 
    family=binomial, 
    data=dat2
)
summary(fit2)

fit1とfit2で異なる部分はdevianeceの部分です

逸脱度(deviance) = -2 × 対数尤度

以下はRの出力のまとめです(Rは尤度比検定に必要なNull devianceResidual devianceを出力します)

Null deviance:切片モデルのdevianceから尤も当てはまりのよいモデル(Full model)のdevianceを引いた値

Residual deviance:回帰モデルのdevianceから尤も当てはまりのよいモデル(Full model)のdevianceを引いた値

fit1fit2
Null deviance43.6268138.59
Null devianceの自由度899
Residual deviance6.3246101.29
Residual devianceの自由度798
AIC30.394105.29

この違いから、fit1とfit2では対数尤度が異なることが理解できます

deviance

fit1のdeviance

$y_i=$満足の数$、n_i=$患者数$、i=1,2,3, \dotsm ,9$ (ロジスティック回帰 (1)より)

対数尤度

$=L_{\theta1}=log(\prod_{i=1}^{n} (P_i^{y_i}(1-P_i)^{n_i-y_i}))$

$\quad=log(\prod_{i=1}^{9} (P_i^{y_i}(1-P_i)^{n_i-y_i}))$

$\quad=\sum(y_i(\beta_0+\beta_1x_i)-n_ilog(1+e^{\beta_0+\beta_1x_i}))+Const.$

fit1の最大逸脱度

fit1_null <- glm(
    cbind(満足, 患者数-満足) ~ 1, 
    family=binomial, 
    data=dat1
)

Dmax1 <- -2*logLik(fit1_null)
print(Dmax1)

fit1モデルの逸脱度

Dreg1 <- -2*logLik(fit1)
print(Dreg1)

fit1の最小逸脱度

二項分布より全ての確率から尤度を算出して対数尤度を求めます

Dmin1 <- -2*sum(log(dbinom(dat1$満足, dat1$患者数, dat1$満足/dat1$患者数)))
print(Dmin1)

最大逸脱度と最小逸脱度との差がNull devianceになります

Rでは Null devianceResidual deviance が算出されます

Null deviance(最大逸脱度$-$最小逸脱度)

Dmax1 - Dmin1

Residual deviance(回帰分析モデルの逸脱度$-$最小逸脱度)

Dreg1 - Dmin1

fit2のdeviance

$y_i=$満足の数$(0, 1)、n_i=1、i=1,2,3, \dotsm ,100$

fit2_null <- glm(
    cbind(満足, 1-満足) ~ 1, 
    family=binomial, 
    data=dat2
)
summary(fit2_null)

fit2の最大逸脱度

Dmax2 <- -2*logLik(fit2_null)
print(Dmax2)

fit2モデルの逸脱度

Dreg2 <- -2*logLik(fit2)
print(Dreg2)

fit2の最小逸脱度

dat2の場合はグループデータではなく全てのデータが対象となっており対数尤度=0となります

したがって最小逸脱度=0

#回帰分析モデルの逸脱度
-2*logLik(fit2)

#-2*対数尤度からも同じ答えが算出されます
b0 <- summary(fit2)$coef[1, 1] 
b1 <- summary(fit2)$coef[2, 1]
print(c(b0, b1))
-2*sum(dat2$満足*(b0+b1*dat2$単位数) - log(1 + exp(b0+b1*dat2$単位数)))

AIC(赤池情報量規準: Akaike’s Information Criterion)

$AIC=-2*($最大対数尤度$-$最尤推定したパラメタ数$)$

例)fit1のAIC

-2*(logLik(fit1) - 2)

参考文献:久保拓弥. (2012). データ解析のための統計モデリング入門 (Vol. 267). 岩波書店.

リハビリテーション研究に必要な統計学について、R(Windows, ChatGPT)を使って紹介してます。サンプルは全て架空のデータで作成しています。したがって解析結果は事実とは異なりますのでご了承ください。間違いなどのご指摘はコメント欄にご記入いただければ助かります。

統計学備忘録をフォローする

ダメ出し 間違い、分かりにくい部分などのご意見をお待ちします

タイトルとURLをコピーしました