ロジット変換
比率pを変換した値と説明変数xの間に線形関係を仮定する一般化線形モデルにおいて、比率pに対する変換として有名なものにロジット変換とプロビット変換がある。いずれもp→0でf(p)→-∞、p→1でf(p)→∞となるような変換であるが、ロジット変換の方がよく用いられる。
例えば、薬の投与量を、効果があらわれた人の割合をとしたとき、との間に
という関係を仮定する。
ロジスティック回帰分析
Rで上記例のパラメータとは以下のように推定する。
> dose <- 1:5 > effect <- c(2, 3, 9, 15, 18) > effect.prop <- effect / 20 > (res <- glm(cbind(effect, 20 - effect) ~ dose, family = binomial)) Call: glm(formula = cbind(effect, 20 - effect) ~ dose, family = binomial) Coefficients: (Intercept) dose -3.812 1.205 Degrees of Freedom: 4 Total (i.e. Null); 3 Residual Null Deviance: 45.34 Residual Deviance: 0.6074 AIC: 19.12
推定されたパラメータに基づく曲線をプロットに重ね合わせるには、次のようにする。
pred <- predict(res, newdata = data.frame(dose = seq(0, 5, 0.1))) fit <- res$family$linkinv(pred) plot(effect.prop ~ dose) lines(fit ~ seq(0, 5, 0.1))
説明変数が連続的に変化する場合のロジスティック回帰分析
説明変数が連続的に変化し、説明変数が特定の値であるときのサンプル数がほとんど1であるような場合にもロジスティック回帰分析は実行できる。
> library(vcd) > (res <- glm(Fail ~ Temperature, data = SpaceShuttle, family = binomial)) Call: glm(formula = Fail ~ Temperature, family = binomial, data = SpaceShuttle) Coefficients: (Intercept) Temperature 15.0429 -0.2322 Degrees of Freedom: 22 Total (i.e. Null); 21 Residual (1 observation deleted due to missingness) Null Deviance: 28.27 Residual Deviance: 20.32 AIC: 24.32
この例では目的変数は2つの水準を持つ因子型ベクトルであるが、分析時には1、0に自動的に対応付けられる。
因子型ベクトルをプロットに用いた時に変換される数値は必ずしも1、0ではないので、データポイントと推定された曲線を重ねてプロットする際には上手く調整してやる必要がある。
pred <- predict(res, newdata = data.frame(Temperature = seq(40, 100, 0.1))) fit <- res$family$linkinv(pred) SpaceShuttle$Failure <- as.numeric(SpaceShuttle$Fail) - 1 with(SpaceShuttle, plot(Temperature, Failure, ylab = "Failure(1:Yes, 0:No)")) lines(fit ~ seq(40, 100, 0.1))
集計データが与えられた場合
集計データは頻度をweight引数に与えることで解析できる。
> library(MASS) > head(housing) Sat Infl Type Cont Freq 1 Low Low Tower Low 21 2 Medium Low Tower Low 21 3 High Low Tower Low 28 4 Low Medium Tower Low 34 5 Medium Medium Tower Low 22 6 High Medium Tower Low 36 > result <- glm(Sat ~ Cont, data = housing, weight = Freq, family = binomial) > result Call: glm(formula = Sat ~ Cont, family = binomial, data = housing, weights = Freq) Coefficients: (Intercept) ContHigh 0.5431 0.2333 Degrees of Freedom: 71 Total (i.e. Null); 70 Residual Null Deviance: 2149 Residual Deviance: 2144 AIC: 2148
ちなみにこのデータは目的変数Satの水準が3つ(Low, Medium, High)なので、本来多項ロジスティック分析で解析すべきものであって係数の意味がよく分からなかったが、どうやらContがLowとHighそれぞれの値の時のSatのMediumとHighの合計の比率が推定されているらしい。最初の水準が0、他の水準が1として解析されるということだろうか。
> unique(result$family$linkinv(predict(result))) [1] 0.6325386 0.6849174 > > Ctable <- xtabs(Freq ~ Cont + Sat, data = housing) > ptable <- prop.table(Ctable, margin = 1) > ptable[, -1] Sat Cont Medium High Low 0.2496494 0.3828892 High 0.2768595 0.4080579 > rowSums(ptable[, -1]) Low High 0.6325386 0.6849174
多重ロジスティック回帰分析
> result <- glm(Fail ~ Temperature + Pressure, data = SpaceShuttle, + family = binomial) > summary(result) Call: glm(formula = Fail ~ Temperature + Pressure, family = binomial, data = SpaceShuttle) Deviance Residuals: Min 1Q Median 3Q Max -1.2130 -0.6089 -0.3870 0.3472 2.0928 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 14.359775 7.442899 1.929 0.0537 . Temperature -0.241540 0.109722 -2.201 0.0277 * Pressure 0.009534 0.008703 1.095 0.2733 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 28.267 on 22 degrees of freedom Residual deviance: 18.972 on 20 degrees of freedom (1 observation deleted due to missingness) AIC: 24.972 Number of Fisher Scoring iterations: 5
ステップワイズ法による変数選択がstepで実行できるが、このとき引数kに有効データ数のlogを与えると情報量基準にBICを使用できる。
> step(result) Start: AIC=24.97 Fail ~ Temperature + Pressure Df Deviance AIC - Pressure 1 20.315 24.315 <none> 18.972 24.972 - Temperature 1 27.503 31.503 Step: AIC=24.32 Fail ~ Temperature Df Deviance AIC <none> 20.315 24.315 - Temperature 1 28.267 30.267 Call: glm(formula = Fail ~ Temperature, family = binomial, data = SpaceShuttle) Coefficients: (Intercept) Temperature 15.0429 -0.2322 Degrees of Freedom: 22 Total (i.e. Null); 21 Residual (1 observation deleted due to missingness) Null Deviance: 28.27 Residual Deviance: 20.32 AIC: 24.32 > step(result, k = log(23)) Start: AIC=28.38 Fail ~ Temperature + Pressure Df Deviance AIC - Pressure 1 20.315 26.586 <none> 18.972 28.378 - Temperature 1 27.503 33.774 Step: AIC=26.59 Fail ~ Temperature Df Deviance AIC <none> 20.315 26.586 - Temperature 1 28.267 31.403 Call: glm(formula = Fail ~ Temperature, family = binomial, data = SpaceShuttle) Coefficients: (Intercept) Temperature 15.0429 -0.2322 Degrees of Freedom: 22 Total (i.e. Null); 21 Residual (1 observation deleted due to missingness) Null Deviance: 28.27 Residual Deviance: 20.32 AIC: 24.32
多項ロジスティック回帰分析、条件付きロジスティック回帰分析
理解が怪しいので保留。