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

カテゴリカルデータ解析読書メモ(第6章)

R

ロジット変換

比率pを変換した値と説明変数xの間に線形関係を仮定する一般化線形モデルにおいて、比率pに対する変換として有名なものにロジット変換とプロビット変換がある。いずれもp→0でf(p)→-∞、p→1でf(p)→∞となるような変換であるが、ロジット変換の方がよく用いられる。
例えば、薬の投与量をx、効果があらわれた人の割合をpとしたとき、xpの間に
\log \frac{p}{1-p} = a + bx
という関係を仮定する。

ロジスティック回帰分析

Rで上記例のパラメータabは以下のように推定する。

> 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))

f:id:Rion778:20160727225542p:plain

説明変数が連続的に変化する場合のロジスティック回帰分析

説明変数が連続的に変化し、説明変数が特定の値であるときのサンプル数がほとんど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))

f:id:Rion778:20160727234850p:plain

集計データが与えられた場合

集計データは頻度を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

BICを指定してもAICという記述になっているが、一番最後のAIC以外はBICになっているらしい。

多項ロジスティック回帰分析、条件付きロジスティック回帰分析

理解が怪しいので保留。