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

Rで逆推定(2)

はてなから1年前のブログとかいうメールが来てこんな記事を書いていたのを思い出した。
rion778.hatenablog.com
再び少し調べたところ、これはinvestrというパッケージで実行できるらしい。

データとモデルオブジェクトの作成

データは前回と同じものを使うが、データフレームにしておく。これは、モデルオブジェクト中の$dataがNULLだと逆推定を行う際にデータの指定が必要となるためだ。

temp <- rep(c(15, 20, 25), c(5, 5, 5))
day <- c(34, 33, 36, 37, 35,
         30, 28, 29, 25, 28,
         23, 25, 22, 24, 26)
d <- data.frame(temp, day)
result <- lm(day ~ temp, data = d)

逆推定

逆推定はinvest()関数で行う。さらに、plotFit()関数で推定区間と共にプロットできる。

library(investr)
(inv <- invest(result, y0 = 30, interval = "inversion", mean.response = TRUE))
# estimate    lower    upper 
# 19.09091 18.09070 19.99693 
plotFit(result, interval = "confidence", shade = TRUE)
abline(h = 30, v = c(inv$lower, inv$estimate, inv$upper), lty = 2)

f:id:Rion778:20160806220620p:plain
ちなみに、JMPで実行した場合の推定値は

  • 予測値: 19.09091
  • 下側95%:18.09069
  • 上側95%:19.99694

であるので、ほぼ一致している。

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

対応分析

カテゴリカルデータにおいて各カテゴリーの変数に適当な数字を割り当てると、2つのカテゴリ間で散布図を書いたり相関係数を計算することができる。このとき、相関係数が最大となるような数値をカテゴリー変数に割り当てることで、カテゴリー間の関係を表現、探索する方法が対応分析である。
まず、カテゴリー変数に適当に数字を割り当ててプロットしてみる。テキストでは円の大きさを度数の対数に対応させているが、分かりにくいので普通に度数に比例させた。

# カテゴリーに数値を割り当て、度数に円の大きさを対応させてプロットする
plot(c(0,6), c(0,4), type = "n")
for(j in 1:5){
  for(i in 1:3){
    points(j,i,cex = tabledata[i,j]/mean(tabledata))
  }
}

f:id:Rion778:20160731005044p:plain
この状態では行と列の間に相関はほとんど見られない。
相関係数の最大化は、MASSパッケージのcorresp()関数で実行できる。

# corresp()を使って相関係数を最大化する
result <- corresp(tabledata, nf = 1)
result
plot(-3:2, -3:2, type = "n", xlab = "Urban", ylab = "Status")
for(i in 1:3){
  for(j in 1:5){
    points(result$cscore[j], result$rscore[i],
           cex = (tabledata[i, j])/mean(tabledata))
  }
}

plotしてみると、度数の多いセルが対角線に近い位置になるよう並んでいることが分かる。
f:id:Rion778:20160731005149p:plain
同様の図は単にcorresp()の結果をplotするだけでも得られる。

plot(result)

f:id:Rion778:20160731005211p:plain
nf = 2とすることで、各カテゴリーを二次元データで表現できる。これをplotすると自動的にバイプロットが作成され、カテゴリ間の関係を視覚的に確認することができる。

result <- corresp(tabledata, nf = 2)
plot(result)

f:id:Rion778:20160731005518p:plain

多重対応分析

テキストのデータだと少し分かりにくいので、MASSパッケージのfarmsデータを例にする。

> summary(farms)
 Mois   Manag  Use    Manure
 M1:7   BF:3   U1:7   C0:6  
 M2:4   HF:5   U2:8   C1:3  
 M4:2   NM:6   U3:5   C2:4  
 M5:7   SF:6          C3:4  
                      C4:3  

farmsは牧草地の環境や管理、利用方法についてのデータで、4つのカテゴリカル変数を持つ。
テキストではMASSパッケージのmca()関数を使用しているが、ここではFactoMineRパッケージのMCA()関数を使用してみる。

library(FactoMineR)
result <- MCA(farms)

MCA()ではplot = FALSEと明示的にOFFにしなければ実行時にプロットが複数表示される。カテゴリーによるバイプロットの結果の部分を示そう。軸名の部分には固有値の割合も同時に示されている。
f:id:Rion778:20160731012703p:plain
この例では、第2軸までではデータの変動の40%程度までしか説明できていない。固有値と累積割合はresult$eigで確認できる。

> result$eig
         eigenvalue percentage of variance
dim 1  6.499174e-01           2.166391e+01
dim 2  5.551954e-01           1.850651e+01
dim 3  5.169428e-01           1.723143e+01
dim 4  3.819977e-01           1.273326e+01
dim 5  3.102940e-01           1.034313e+01
dim 6  2.208944e-01           7.363148e+00
dim 7  1.332712e-01           4.442372e+00
dim 8  8.908661e-02           2.969554e+00
dim 9  7.744688e-02           2.581563e+00
dim 10 4.752489e-02           1.584163e+00
dim 11 1.742866e-02           5.809553e-01
dim 12 2.442300e-32           8.140999e-31
       cumulative percentage of variance
dim 1                           21.66391
dim 2                           40.17043
dim 3                           57.40185
dim 4                           70.13511
dim 5                           80.47825
dim 6                           87.84139
dim 7                           92.28377
dim 8                           95.25332
dim 9                           97.83488
dim 10                          99.41904
dim 11                         100.00000
dim 12                         100.00000

また、各カテゴリーがどの程度軸に寄与しているかはresult$var$contribで確認できる。

> result$var$contrib
          Dim 1      Dim 2        Dim 3      Dim 4      Dim 5
M1  1.997032543  6.3987899 6.653719e-04  2.8555287 21.7714144
M2  1.332351513  6.1585438 9.578125e+00 11.0329611  0.1482187
M4  1.942796388  2.6231723 1.344942e+00  4.5402859 46.1550354
M5  9.185471009  2.3084618 8.911297e+00  0.1010620  0.5528613
BF  1.253592590  8.4488996 5.772589e+00 12.9978140 11.0988568
HF  0.474068031 12.2938426 3.762828e+00  3.7546643  2.2497386
NM 20.591151562  3.0315099 1.913005e+00  1.6407989  0.6060488
SF  9.718852258 12.3551134 1.719340e+00  0.2505076  3.1151790
U1  6.648880395  1.2997710 1.072406e+01  0.2151174  2.9499596
U2  9.699621774  6.0048359 2.352817e-02  4.1868505  0.9439940
U3  0.789421444  3.0648808 1.655492e+01  9.8408916  0.6452055
C0 20.591151562  3.0315099 1.913005e+00  1.6407989  0.6060488
C1  0.237784097 12.3464560 7.134416e+00 17.8537482  5.6832374
C2  5.424873542  4.5789897 4.094891e-01 24.3956217  2.2873335
C3  0.002754719  0.1617089 2.943082e+01  3.1784495  0.8359484
C4 10.110196574 15.8935143 8.069635e-01  1.5148998  0.3509199

これを見ると、第1軸にはNM、C0、C4などの影響が大きく、第2軸にはHF、SF、C1、C4などの影響が大きい。

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

対数線形モデル

対数線形モデルにおいては、分割表における各セルの期待度数を予測する。すなわち、複数のカテゴリカル変数の中から特定の目的変数を設定するのでなく、全ての変数について相互の関係を調べることを目的としている。

lonlin()とloglm()

対数線形モデルによる解析を行う関数としてloglin()とloglm()が用意されている。
loglin()は以下のように使用する。

> library(vcd)
> tabledata <- xtabs(Freq ~ Status + Urban, data = DanishWelfare)
> result <- loglin(tabledata, margin = list(1, 2), param = TRUE, fit = TRUE)
2 iterations: deviation 2.273737e-13 
> result$para
$`(Intercept)`
[1] 5.383959

$Status
     Widow    Married  Unmarried 
-0.7693390  1.0654469 -0.2961079 

$Urban
   Copenhagen SubCopenhagen     LargeCity          City       Country 
   -0.4836304    -0.3771835    -0.4102991     0.6787275     0.5923855 

対してloglm()は次のようにモデル式で入力する。

> result <- loglm(~ Status + Urban, data = tabledata)
> result
Call:
loglm(formula = ~Status + Urban, data = tabledata)

Statistics:
                      X^2 df P(> X^2)
Likelihood Ratio 153.9505  8        0
Pearson          158.1145  8        0
> coef(result)
$`(Intercept)`
[1] 5.383959

$Status
     Widow    Married  Unmarried 
-0.7693390  1.0654469 -0.2961079 

$Urban
   Copenhagen SubCopenhagen     LargeCity          City       Country 
   -0.4836304    -0.3771835    -0.4102991     0.6787275     0.5923855 

計算される値はどちらも同じだが、入出力の理解のしやすさはloglm()の方が上の様に思う。以降、loglm()を使用する。

期待度数と残差を求める

期待度数はfitted()、残差はresiduals()を用いて計算できる。

> round(fitted(result), 1)
Re-fitting to get fitted values
           Urban
Status      Copenhagen SubCopenhagen LargeCity   City Country
  Widow           62.2          69.2      67.0  199.0   182.5
  Married        389.9         433.6     419.5 1246.5  1143.4
  Unmarried       99.9         111.1     107.5  319.4   293.0
> round(residuals(result), 1)
Re-fitting to get frequencies and fitted values
           Urban
Status      Copenhagen SubCopenhagen LargeCity City Country
  Widow            6.1           2.5      -1.5  1.0    -6.3
  Married         -5.4          -0.2       0.0  0.0     3.1
  Unmarried        4.5          -1.8       1.2 -0.8    -1.7

三元表での対数線形モデル

変数が3つ以上になると作成できるモデルの数がかなり多くなる。
交互作用項について細かく検討するのであればupdate()を利用すると入力の手間を節約できる。

> Modelindependent <- loglm(~ Gender + Admit + Dept, data = UCBAdmissions)
> Modelindependent
Call:
loglm(formula = ~Gender + Admit + Dept, data = UCBAdmissions)

Statistics:
                      X^2 df P(> X^2)
Likelihood Ratio 2097.671 16        0
Pearson          2000.328 16        0
> # GenderとAdmitの交互作用を加える
> ModelGA <- update(Modelindependent, .~. + Gender:Admit)
> ModelGA
Call:
loglm(formula = . ~ Gender + Admit + Dept + Gender:Admit, data = UCBAdmissions)

Statistics:
                      X^2 df P(> X^2)
Likelihood Ratio 2004.222 15        0
Pearson          1748.160 15        0

テキストでは全てのモデルを作成して比較しているが、面倒なのでstep()を使用してモデル選択をしてみよう。
まずフルモデルを作成する。

> ModelFull <- loglm(~ Gender * Admit * Dept, data = UCBAdmissions)
> ModelFull
Call:
loglm(formula = ~Gender * Admit * Dept, data = UCBAdmissions)

Statistics:
                 X^2 df P(> X^2)
Likelihood Ratio   0  0        1
Pearson            0  0        1

フルモデルでは期待度数と観測度数が一致する。
では、AICによるモデル選択を試みる。

> step(ModelFull)
Start:  AIC=48
~Gender * Admit * Dept

                    Df    AIC
<none>                 48.000
- Gender:Admit:Dept  5 58.204
Call:
loglm(formula = ~Gender * Admit * Dept, data = UCBAdmissions)

Statistics:
                 X^2 df P(> X^2)
Likelihood Ratio   0  0        1
Pearson            0  0        1

AICではフルモデルが選択された。
次にBICによるモデル選択を試みる。

> step(ModelFull, k = log(sum(UCBAdmissions)))
Start:  AIC=202.02
~Gender * Admit * Dept

                    Df    AIC
- Gender:Admit:Dept  5 180.14
<none>                 202.02

Step:  AIC=180.14
~Gender + Admit + Dept + Gender:Admit + Gender:Dept + Admit:Dept

               Df     AIC
- Gender:Admit  1  173.25
<none>             180.14
- Admit:Dept    5  901.45
- Gender:Dept   5 1266.75

Step:  AIC=173.25
~Gender + Admit + Dept + Gender:Dept + Admit:Dept

              Df     AIC
<none>            173.25
- Admit:Dept   5  986.49
- Gender:Dept  5 1351.78
Call:
loglm(formula = ~Gender + Admit + Dept + Gender:Dept + Admit:Dept, 
    data = UCBAdmissions, evaluate = FALSE)

Statistics:
                      X^2 df    P(> X^2)
Likelihood Ratio 21.73551  6 0.001351993
Pearson          19.93841  6 0.002840164

BICを基準にした場合は、Admitとgenderの相互作用がモデルに含まれていない。これは第5章での解析結果とも一致する。
AICBICで結果が異なったが、データ数が多い場合にAICを基準とすると、フルモデルを選択しやすい傾向があり、BICをモデル選択に使用する場合が多い、とのことである。

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

ポアソン分布

適合度検定はvcdパッケージのgoodfit()関数で行う。

> t <- 0:5
> x <- c(8, 6, 8, 3, 4, 1)
> x <- matrix(c(x, t), nr = 6)
> library(vcd)
> result <- goodfit(x, type = "poisson", method = "MinChisq")
> result

Observed and fitted values for poisson distribution
with parameters estimated by `MinChisq' 

 count observed    fitted pearson residual
     0        8 5.1375937      1.262848928
     1        6 9.0658624     -1.018235178
     2        8 7.9988673      0.000400486
     3        3 4.7049673     -0.786027501
     4        4 2.0756111      1.335733618
     5        1 0.7325299     -0.016953894

結果をplot()に渡すことで、結果を視覚的に判断することもできる。
f:id:Rion778:20160730184840p:plain

ポアソン回帰分析

データセットInsectSpraysは殺虫剤ごとの死虫数*1をまとめたものである。このデータは平均値が大きいほど分散が大きく、ポアソン分布の特徴を持っているように見える。

ggplot(InsectSprays, aes(x = spray, y = count)) + geom_boxplot() 

f:id:Rion778:20160730185448p:plain
glm()による解析を行ってみる。

> result <- glm(count ~ spray, data = InsectSprays, family = "poisson")
> summary(result)

Call:
glm(formula = count ~ spray, family = "poisson", data = InsectSprays)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3852  -0.8876  -0.1482   0.6063   2.6922  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.67415    0.07581  35.274  < 2e-16 ***
sprayB       0.05588    0.10574   0.528    0.597    
sprayC      -1.94018    0.21389  -9.071  < 2e-16 ***
sprayD      -1.08152    0.15065  -7.179 7.03e-13 ***
sprayE      -1.42139    0.17192  -8.268  < 2e-16 ***
sprayF       0.13926    0.10367   1.343    0.179    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 409.041  on 71  degrees of freedom
Residual deviance:  98.329  on 66  degrees of freedom
AIC: 376.59

Number of Fisher Scoring iterations: 5

この結果から推定値を得るためには、例えばIntercept(=sprayA)であれば\exp(2.67415) = 14.5、sprayBであれば\exp(2.67415+0.05588) = 15.3というようにやるのだが、まとめて求めるにはpredict()関数を使用する。

> pred <- predict(result, newdata = data.frame(spray = LETTERS[1:6]),
+ type = "response")
> pred
        1         2         3         4         5         6 
14.500000 15.333333  2.083333  4.916667  3.500000 16.666667 

newdataには説明変数と同じ名前を持つ列に、予測したい説明変数の値をセットしたデータフレームを指定する。省略すると元のデータフレームの説明変数の値全てに対応する予測値が表示される。typeに"response"を指定することで、実際の予測値(上でexp()を使って計算した後の値)が表示される。指定しない場合は線形予測子の値(上でexp()に与えた中身)になる。
ついでに95%予測区間を求めてプロットしてみる。

> upr <- qpois(0.975, pred)
> lwr <- qpois(0.025, pred)
> ggplot(InsectSprays, aes(x = spray, y = count)) + geom_point() +
+ geom_pointrange(data = data.frame(spray = LETTERS[1:6], count = pred),
+ aes(ymin = lwr, ymax = upr), color = "blue")

f:id:Rion778:20160730190809p:plain
ちなみにこの予測区間はパラメータの標準誤差を考慮していないので、実際の予測区間よりたぶんちょっと狭い。そもそもこのような方法で予測区間を構成して良いものかどうかよく分からない。

オフセットによる調整法

まず、都道府県別の交通事故発生件数の一部を対象にポアソン回帰分析を行ってみる。

> accident <- c(45703, 8906, 7938, 12091, 7327, 9820, 11526, 6525)
> pref <- c("A福岡", "B佐賀", "C長崎", "D熊本", 
+           "E大分", "F宮崎", "G鹿児島", "H沖縄")
> result1 <- glm(accident ~ pref, family = "poisson")
> result1

Call:  glm(formula = accident ~ pref, family = "poisson")

Coefficients:
(Intercept)    prefB佐賀    prefC長崎    prefD熊本    prefE大分  
     10.730       -1.635       -1.751       -1.330       -1.831  
  prefF宮崎  prefG鹿児島    prefH沖縄  
     -1.538       -1.378       -1.947  

Degrees of Freedom: 7 Total (i.e. Null);  0 Residual
Null Deviance:	    60920 
Residual Deviance: 5.994e-12 	AIC: 105.1

この結果を図示して解釈するためにまたpredict()を使ってみるが、typeに"terms"を指定すると、係数をスケールとした値が得られる。
具体的には"constant"属性に1つの値を持つベクトルが返ってくるが、constantの値とベクトルの各値を足したものが線形予測子である(result1のCoefficientsの値と比較してみると良い)。

> pred <- predict(result1, type = "terms")
> pred
         pref
1  1.42600826
2 -0.20943048
3 -0.32449433
4  0.09630569
5 -0.40458953
6 -0.11173456
7  0.04844967
8 -0.52051473
attr(,"constant")
[1] 9.303911

これをプロットしてみる。

ggplot(data.frame(pref, pred), 
       aes(x = pref, y = pred)) +
  geom_bar(stat = "identity") +
  theme_bw(base_family = "HiraKakuProN-W3") #Macで日本語を表示するために指定

f:id:Rion778:20160730192121p:plain
圧倒的に福岡の係数が大きいが、これをもって福岡は交通事故が多いと判断して良いだろうか?
ここで、人口当たりの交通事故発生件数で解析することを試みる。
交通事故発生件数の平均を\lambda、人口をnとすると、次のようなモデルを解析することになる。
\log\lambda = \log n + a + b_1\mathrm{prefB} + b_2\mathrm{prefC} + ...
真数の商は対数の差であるので、右辺に\log nを含めて解析することで\frac{\lambda}{n}を解析することができる。\log{n}の様な項をオフセットと呼び、glmではoffset=に指定する。

> pop <- c(5056, 859, 1453, 1828, 1203, 1143, 1730, 1373)
> result2 <- glm(accident ~ pref, offset = log(pop), family = poisson)
> result2

Call:  glm(formula = accident ~ pref, family = poisson, offset = log(pop))

Coefficients:
(Intercept)    prefB佐賀    prefC長崎    prefD熊本    prefE大分  
    2.20159      0.13712     -0.50356     -0.31235     -0.39484  
  prefF宮崎  prefG鹿児島    prefH沖縄  
   -0.05082     -0.30510     -0.64295  

Degrees of Freedom: 7 Total (i.e. Null);  0 Residual
Null Deviance:	    5699 
Residual Deviance: 8.338e-12 	AIC: 105.1

先ほどと同様にプロットしてみると、結果は違ったものとなる。

> pred <- predict(result2, type = "terms")
> ggplot(data.frame(pref, pred),
+ aes(x = pref, y = pred)) +
+ geom_bar(stat = "identity") +
+ theme_bw(base_family = "HiraKakuProN-W3")

f:id:Rion778:20160730193421p:plain
この例ではもともとがフルモデルによる解析なのでoffsetを指定してもモデルが変化せず、従ってAICも変化しない*2

過分散である場合の解析方法

ポアソン分布におけるパラメータ\lambdaが、観測のたびにガンマ分布に従って変動すると仮定する。
このような仮定にもとづいて乱数を発生させると、平均よりも分散が大きくなる。

> lambda <- rgamma(n, shape = 5, scale = 3)
> y <- rpois(n, lambda)
> mean(y)
[1] 13.15
> var(y)
[1] 49.08158

このような状況を過分散と呼ぶ。
次に、過分散する用量反応データを人工的に作成する。

n = 20
v = 3
dose <- 0:4 * 5
x <- rep(dose, n)
a <- exp(3 + 0.05 * dose) / v
lambda <- rgamma(5 * n, shape = a, scale = v)
y <- rpois(5 * n, lambda)
dataNB <- data.frame(x, y)

xで層別してプロットするとともに平均・分散を求めてみる。

> ggplot(dataNB, aes(x = x, y = y)) + geom_boxplot(aes(group = x))
> library(dplyr)
> dataNB %>% group_by(x) %>% summarize(mean(y), var(y))
# A tibble: 5 x 3
      x mean(y)    var(y)
  <dbl>   <dbl>     <dbl>
1     0   18.60  43.09474
2     5   23.80  97.74737
3    10   30.70  56.22105
4    15   49.45 134.47105
5    20   52.70 212.64211

f:id:Rion778:20160730194443p:plain
そのように作成したので当然ではあるが、過分散が確認できる。
このようなデータの解析には負の二項分布モデルを使用することができ、RではMASSパッケージのglm.nb()関数を用いる。使用方法はglm()とほとんど同じである。

> # 負の二項分布モデルによる分析
> result <- glm.nb(y ~ x, data = dataNB)
> # ポアソン回帰分析
> result2 <- glm(y ~ x, data = dataNB, family = poisson)
> result

Call:  glm.nb(formula = y ~ x, data = dataNB, init.theta = 15.79541194, 
    link = log)

Coefficients:
(Intercept)            x  
    2.91633      0.05636  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    266.5 
Residual Deviance: 102.7 	AIC: 747.9
> result2

Call:  glm(formula = y ~ x, family = poisson, data = dataNB)

Coefficients:
(Intercept)            x  
     2.9285       0.0553  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    825.1 
Residual Deviance: 314.4 	AIC: 845.6

ポアソン回帰分析に比べ、残差逸脱度、AIC共に小さくなっていることが分かる。
これもまたプロットしておこう。

# プロット
ux <- seq(0, 20, by = 0.1)
# 平均(mu)の計算
pred <- predict(result, newdata = data.frame(x = ux),
                type = "response")
# 推定された形状パラメータ(theta)から予測区間を構成
# ※推定値の標準誤差を考慮していない
upr <- qnbinom(0.975, size = result$theta, mu = pred)
lwr <- qnbinom(0.025, size = result$theta, mu = pred)
# plot
ggplot(data.frame(ux, pred),aes(x = ux, y = pred)) +
  geom_line() +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.25) +
  geom_point(data = dataNB, aes(x = x, y = y))

f:id:Rion778:20160730201831p:plain
プロットの結果から見てもさほど問題は無いと思うが、予測区間はパラメータの標準誤差を考慮していないのできっと少し狭い。

*1:原典にあたってないが、helpにはInsect countとしか書いてないので生虫数かもしれない

*2:フルモデルで解析しても元の値が返ってくるだけなので、そもそもglmで解析する必要があるのかはよく分からない

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

ロジット変換

比率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になっているらしい。

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

理解が怪しいので保留。

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

シンプソンのパラドックス

層別をしないで解析すると関連が見られるが、層別をして解析すると関連が見られなくなるような現象。層別に用いるカテゴリカル変数が、他のカテゴリカル変数全てに影響を与えているような場合に発生する、見かけ上の相関。連続変数で言うところの擬似相関と同様のもの。2つの変数の関係を調べる際には両方の変数に影響する他の変数が無いかどうかに注意する必要がある。

層別2 x 2表の解析

カリフォルニア州立大学バークレー校の1973年の入試結果データであるUCBAdmissionsについて、各クラスの性別ごとの合否結果を合算して集計すると、男性の方が合格率が高いという結果が得られる。

> GAdata <- apply(UCBAdmissions, c(2, 1), sum)
> GAdata
        Admit
Gender   Admitted Rejected
  Male       1198     1493
  Female      557     1278
> chisq.test(GAdata)

	Pearson's Chi-squared test with Yates' continuity correction

data:  GAdata
X-squared = 91.61, df = 1, p-value < 2.2e-16

> mosaicplot(GAdata)

f:id:Rion778:20160726215732p:plain
しかし、層別にプロットすると、A学科を除いて合格者と不合格者の男女比はあまり変わらない。

old <- par(mfrow = c(2, 3), mar = c(3, 3, 3, 3), oma = c(2, 2, 2, 2))
for(i in 1:6){
  mosaicplot(~ Admit + Gender, data = UCBAdmissions[,,i],
             main = paste("Department", LETTERS[i]))
}
mtext("Student admissions at UC Barkeley", outer = TRUE)
par(old)

f:id:Rion778:20160726215922p:plain

条件付き独立

2つの変数間の関係について、3つめの変数で層別した場合に関係が見られなくなるような場合、最初の2つの変数は条件付き独立であるという。
3つめの変数で層別された2 x 2表について、条件付き独立であることを帰無仮説とする検定(マンテル・ヘンセル検定)を、Rではmantelhaen.test()で実行できる。

> mantelhaen.test(UCBAdmissions)

	Mantel-Haenszel chi-squared test with continuity correction

data:  UCBAdmissions
Mantel-Haenszel X-squared = 1.4269, df = 1, p-value = 0.2323
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 0.7719074 1.0603298
sample estimates:
common odds ratio 
        0.9046968 

ここでcommon odds ratioとして推定されているのは共通オッズ比と呼ばれるもので、層別した場合の各層におけるオッズ比は等しいと仮定した場合のオッズ比を指している。

オッズ比の均一性の検定

それぞれの層のオッズ比が等しいかどうかを検定するために、パッケージepiDisplayのmhor()関数が使用できる*1

> library(epiDisplay)
> mhor(mhtable = UCBAdmissions)

Stratified analysis by  Dept 
                OR lower lim. upper lim.  P value
Dept A       0.350      0.197      0.592 1.67e-05
Dept B       0.803      0.294      2.004 6.77e-01
Dept C       1.133      0.845      1.516 3.87e-01
Dept D       0.921      0.679      1.250 5.99e-01
Dept E       1.221      0.806      1.839 3.60e-01
Dept F       0.828      0.433      1.576 5.46e-01
M-H combined 0.905      0.772      1.060 2.17e-01

M-H Chi2(1) = 1.52 , P value = 0.217 
Homogeneity test, chi-squared 5 d.f. = 17.96 , P value = 0.003 

f:id:Rion778:20160726221638p:plain
また、テキストに記述のある「マンテル・ヘンセル統計量を用いた別の均一性の検定法(Breslow and Day, 1980)」は、現在はDescToolsパッケージのBreslowDayTest()関数で実行できる様子。

> library(DescTools)
> BreslowDayTest(UCBAdmissions)

	Breslow-Day test on Homogeneity of Odds Ratios

data:  UCBAdmissions
X-squared = 18.826, df = 5, p-value = 0.002071

層別したオッズ比を求める

層別した後のオッズ比はvcdパッケージのoddsratio()関数で計算でき、confint()による信頼区間の計算やplot()によるグラフ出力にも対応している。
当時よりバージョンが上がっているのか、テキストに記載のものに比べると出力が見やすくなっている印象。

> library(vcd)
> lor <- oddsratio(UCBAdmissions, log = FALSE)
> lor
 odds ratios for Admit and Gender by Dept 

        A         B         C         D         E         F 
0.3492120 0.8025007 1.1330596 0.9212838 1.2216312 0.8278727 
> confint(lor)
                                    2.5 %    97.5 %
Admitted:Rejected/Male:Female|A 0.2086756 0.5843954
Admitted:Rejected/Male:Female|B 0.3403815 1.8920166
Admitted:Rejected/Male:Female|C 0.8545328 1.5023696
Admitted:Rejected/Male:Female|D 0.6863345 1.2366620
Admitted:Rejected/Male:Female|E 0.8250748 1.8087848
Admitted:Rejected/Male:Female|F 0.4552059 1.5056335

f:id:Rion778:20160726222138p:plain

層別I x J表の解析

層別I x J表においても、条件付き独立性の検定にはmantelhaen.test()が、線形関連の検定にはcoinパッケージのlbl_test()関数が同様に使用できる。

> library(coin)
> mantelhaen.test(jobsatisfaction)

	Cochran-Mantel-Haenszel test

data:  jobsatisfaction
Cochran-Mantel-Haenszel M^2 = 10.2, df = 9, p-value = 0.3345

> lbl_test(jobsatisfaction)

	Asymptotic Linear-by-Linear Association Test

data:  Job.Satisfaction (ordered) by
	 Income (<5000 < 5000-15000 < 15000-25000 < >25000) 
	 stratified by Gender
Z = 2.5736, p-value = 0.01006
alternative hypothesis: two.sided

*1:テキストに記載のepicalcパッケージの新しいもの。epicalcパッケージに入っていた関数はそのまま使える様子。cf:r - Why was package 'epicalc' removed from CRAN? - Stack Overflow cf2:整形した結果を表示する - 盆栽日記

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

オッズとオッズ比

2つの割合p_1p_2に対して、オッズ\frac{p}{1-p}の比\frac{p_2(1-p_1)}{p_1(1-p_2)}を考えると、このオッズ比が1より大ならp_1\lt p_2、1より小ならp_1\gt p_2である。
オッズ比は、ロジスティック回帰において回帰係数がオッズ比の対数に一致することや、患者対照研究のようにオッズ比のみが推定できる研究方法があるためにしばしば用いられる。pが小さい場合には、オッズ比は割合の比に近似するので、リスクが何倍になる、といった文脈で使われる場合もある。

独立性のカイ二乗検定

2つの二項割合p_1p_2について、帰無仮説をp_1=p_2とする。このとき、割合の推定量の差を分散のルートで補正した量の2乗をピアソンのカイ二乗統計量\chi ^2と呼び、データ数が大きい時に\chi^2が自由度1のカイ二乗分布で近似されるため、この近似によって棄却域やp値が計算できる。
このとき、より近似を良くするためにイェーツの補正と呼ばれる補正を行う場合がある。Rのchisq.test()ではデフォルトである。

> Ctable <- matrix(c(512, 1385, 567, 2419), nr = 2)
> chisq.test(Ctable)

	Pearson's Chi-squared test with Yates' continuity correction

data:  Ctable
X-squared = 42.679, df = 1, p-value = 6.449e-11

この補正を行ってもデータ数が少ないと近似が不正確になり、警告メッセージが表示される場合がある。

> Ctable <- matrix(c(2, 10, 3, 15), nr = 2)
> chisq.test(Ctable)

	Pearson's Chi-squared test

data:  Ctable
X-squared = 0, df = 1, p-value = 1

 警告メッセージ: 
 chisq.test(Ctable):   カイ自乗近似は不正確かもしれません 

フィッシャーの直接確率法

紅茶が先かミルクが先かという話は、「ミルクが先」が美味しいということで決着が付いているらしい*1が、ともかく、サンプル数が少ない場合に超幾何分布モデルによって結果が偶然に起こる確率を直接計算する方法をフィッシャーの直接確率法と呼ぶ。
Rではfisher.test()で行う。なお、信頼区間や推定値はオッズ比として表示される。

> table <- matrix(c(3, 1, 1, 3), nr = 2)
> fisher.test(table)

	Fisher's Exact Test for Count Data

data:  table
p-value = 0.4857
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
   0.2117329 621.9337505
sample estimates:
odds ratio 
  6.408309 

かつてはデータ数が大きくなると計算量が増えるという点が問題であったが、2x2程度の表であれば多少データ数が大きくても実行に支障はない。

> Ctable <- matrix(c(512, 1385, 567, 2419), nr = 2)
> fisher.test(Ctable)

	Fisher's Exact Test for Count Data

data:  Ctable
p-value = 8.63e-11
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.372396 1.811955
sample estimates:
odds ratio 
  1.576932 

関連性の指標

2x2表が独立ではない場合にどの程度の関連性があるかを調べる指標をまとめて計算できる関数assocstats()がパッケージvcdに入っている。また、oddsratio()でオッズ比も計算できる。

> library(MASS)
> library(vcd)
> Ctable <- xtabs(~ Sex + W.Hnd, data = survey)
> assocstats(Ctable)
                     X^2 df P(> X^2)
Likelihood Ratio 0.54629  1  0.45984
Pearson          0.54351  1  0.46098

Phi-Coefficient   : 0.048 
Contingency Coeff.: 0.048 
Cramer's V        : 0.048 
> oddsratio(Ctable)
log odds ratios for Sex and W.Hnd 

[1] -0.3750241
> oddsratio(Ctable, log = FALSE) # デフォルトでは対数表示
 odds ratios for Sex and W.Hnd 

[1] 0.6872727

2xJ表の解析

カテゴリカル変数が2つよりも大きい場合もchisq.test()は実行できる。

> Ctable <- xtabs(~ Treatment + Improved, data = Arthritis)
> Ctable
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21
> chisq.test(Ctable)

	Pearson's Chi-squared test

data:  Ctable
X-squared = 13.055, df = 2, p-value = 0.00146

ただし、このときは単にPlaceboとTreatedの分布が違う、程度のことしか言えない。どの値が帰無仮説からのずれが大きいかを調べる方法として、調整済み残差を用いる方法がある。

> Ctable <- xtabs(~ Treatment + Improved, data = Arthritis)
> result <- chisq.test(Ctable)
> v <- result$expected * (1 - rowSums(Ctable)/sum(Ctable)) %o%
+ (1 - colSums(Ctable)/sum(Ctable))
> (Ctable - result$expected) / sqrt(v)
         Improved
Treatment        None        Some      Marked
  Placebo  3.27419655 -0.09761768 -3.39563632
  Treated -3.27419655  0.09761768  3.39563632

この残差が1.96より大きければ帰無仮説からのずれが大きいと考える。

カテゴリカル変数に順序関係のある場合の解析

マンテル検定

カテゴリカル変数間に順序関係が有るとき、これに1, 2, 3...のようなスコアを与える検定手法としてマンテル検定がある。

manteltrend.test <- function(x, score = c(1:dim(x)[2])){
  ctotal <- colSums(x) # 列和
  rtotal <- rowSums(x) # 行和
  n <- sum(x) # データ総数
  expected <- (rtotal %o% ctotal) / n # 期待度数
  tscore <- (x - expected) %*% score
  v <- n * (score^2 %*% ctotal) - (score %*% ctotal)^2
  statistics <- n^2 * (n-1) * tscore[2]^2 / rtotal[1] / rtotal[2] / v
  p.value <- pchisq(statistics, 1, lower.tail = FALSE)
  list(statistics = statistics, p_value = p.value)
}

スコアは等間隔である必要は無いが、どのようなスコア間隔を用いるかは大勢に影響を与えない場合が多い。

> manteltrend.test(Ctable, score=0:2)
$statistics
         [,1]
[1,] 12.85902

$p_value
             [,1]
[1,] 0.0003358568

> manteltrend.test(Ctable, score=c(1, 2, 5))
$statistics
         [,1]
[1,] 12.60627

$p_value
             [,1]
[1,] 0.0003844559
Wilcoxonの順位和検定

Wilcoxonの順位和検定では、スコアを使わず、各データを順位に変換して検定を行う。Rではwilcox.test()で行えるが、あらかじめデータを1人1人のデータに戻しておく必要がある。

> x <- rep(col(Ctable)[1,], Ctable[1,])
> y <- rep(col(Ctable)[2,], Ctable[2,])
> wilcox.test(x, y)

	Wilcoxon rank sum test with continuity correction

data:  x and y
W = 517.5, p-value = 0.0003666
alternative hypothesis: true location shift is not equal to 0

 警告メッセージ: 
 wilcox.test.default(x, y): 
   タイがあるため、正確な p 値を計算することができません 
コクラン・アーミテージ検定

薬の投与量に応じて効果のありなしが分かれるようなデータがあったとして、効果ありの割合p_jが投与量c_jに対して、次のような関係にあったとする。
p_j = \alpha + \beta c_j
ここで\beta = 0であるとすると、投与量と効果ありの割合の間に関係がないということになる。このとき、ロジスティック回帰モデルのようにp_jを何らかの関数で変換したとしても、検定方法は変わらない。この検定はRではprop.trend.test()で実行できる。

> gd <- c(2, 3, 9, 15, 18)
> n <- rep(20, 5)
> prop.trend.test(gd, n, score = 1:5)

	Chi-squared Test for Trend in Proportions

data:  gd out of n ,
 using scores: 1 2 3 4 5
X-squared = 38.86, df = 1, p-value = 4.553e-10

I x J表の解析

独立性の検定

chisq.test()がこれまでと同様に使用できる。

> Ctable <- xtabs(Freq ~ Status + Alcohol, data = DanishWelfare)
> chisq.test(Ctable)

	Pearson's Chi-squared test

data:  Ctable
X-squared = 87.972, df = 4, p-value < 2.2e-16
線形関連の検定

行と列が共に順序カテゴリカルデータの場合、これらにスコアを与えて相関係数を求めることができる。相関係数をr、標本数をnとしたとき、
M^2 = (n - 1)r^2
は、標本数nが大きくなると近似的に自由度1のカイ二乗分布に近づく。このことを利用した検定を線形関連の検定と呼ぶ。
Rではcoinパッケージのlbl_test()関数で実行できる。

> library(coin)
> Ctable <- apply(jobsatisfaction, 1:2, sum) # 性別を合計
> Ctable
             Job.Satisfaction
Income        Very Dissatisfied A Little Satisfied Moderately Satisfied
  <5000                       2                  4                   13
  5000-15000                  2                  6                   22
  15000-25000                 0                  1                   15
  >25000                      0                  3                   13
             Job.Satisfaction
Income        Very Satisfied
  <5000                    3
  5000-15000               4
  15000-25000              8
  >25000                   8
> lbl_test(as.table(Ctable))

	Asymptotic Linear-by-Linear Association Test

data:  Job.Satisfaction (ordered) by
	 Income (<5000 < 5000-15000 < 15000-25000 < >25000)
Z = 2.7623, p-value = 0.005739
alternative hypothesis: two.sided

スコアは指定しなければ自動的に1, 2, 3, ...と設定されるが、指定したければlistで与える。

lbl_test(as.table(Ctable), 
         scores = list(Job.Satisfaction = 1:4,
                       Income = c(2.5, 10, 20, 30)))

なお、I x J表においてI=2の場合はマンテル検定と同じ結果が得られる。