カテゴリカルデータ解析読書メモ(第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()に渡すことで、結果を視覚的に判断することもできる。

ポアソン回帰分析

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

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

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)であれば、sprayBであればというようにやるのだが、まとめて求めるには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")

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

オフセットによる調整法

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

> 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で日本語を表示するために指定

圧倒的に福岡の係数が大きいが、これをもって福岡は交通事故が多いと判断して良いだろうか？
ここで、人口当たりの交通事故発生件数で解析することを試みる。
交通事故発生件数の平均を、人口をとすると、次のようなモデルを解析することになる。

真数の商は対数の差であるので、右辺にを含めて解析することでを解析することができる。の様な項をオフセットと呼び、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")

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

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

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

> 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

そのように作成したので当然ではあるが、過分散が確認できる。
このようなデータの解析には負の二項分布モデルを使用することができ、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))

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

練習問題解答

カテゴリカルデータ解析第7章練習問題解答 · GitHub

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

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

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

ロジット変換

比率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

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

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

理解が怪しいので保留。

練習問題解答

カテゴリカルデータ解析第6章練習問題解答 · GitHub

カテゴリカルデータ解析読書メモ(第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)

しかし、層別にプロットすると、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)

条件付き独立

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 

また、テキストに記述のある「マンテル・ヘンセル統計量を用いた別の均一性の検定法(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

層別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

練習問題解答

カテゴリカルデータ解析第5章練習問題解答 · GitHub

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

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

オッズとオッズ比

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

独立性のカイ二乗検定

2つの二項割合、について、帰無仮説をとする。このとき、割合の推定量の差を分散のルートで補正した量の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 値を計算することができません 

コクラン・アーミテージ検定

薬の投与量に応じて効果のありなしが分かれるようなデータがあったとして、効果ありの割合が投与量に対して、次のような関係にあったとする。

ここでであるとすると、投与量と効果ありの割合の間に関係がないということになる。このとき、ロジスティック回帰モデルのようにを何らかの関数で変換したとしても、検定方法は変わらない。この検定は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

線形関連の検定

行と列が共に順序カテゴリカルデータの場合、これらにスコアを与えて相関係数を求めることができる。相関係数を、標本数をとしたとき、

は、標本数が大きくなると近似的に自由度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の場合はマンテル検定と同じ結果が得られる。

練習問題解答

カテゴリカルデータ解析第4章練習問題解答 · GitHub

*1:ミルクティー - Wikipedia

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

第3章の内容は大体知ってる事なのでささっと。

二項検定

> binom.test(405, 765)

	Exact binomial test

data:  405 and 765
number of successes = 405, number of trials = 765, p-value =
0.1116
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.4933327 0.5652633
sample estimates:
probability of success 
             0.5294118 

pのデフォルトは0.5なのでテキストの例では省略できる。
ついでに95%信頼区間が分かる。

多項分布

多項分布用に*multinom()関数シリーズが用意されている。

> p <- c(0.382, 0.219, 0.305, 0.094)
> x <- c(39, 16, 27, 15)
> dmultinom(x, prob = p)
[1] 0.0001047239

カイ二乗検定

> x <- c(62, 52, 53, 38, 46, 49)
> chisq.test(x)

	Chi-squared test for given probabilities

data:  x
X-squared = 6.36, df = 5, p-value = 0.2727

chisq.test()のpの初期値は等確率なのでこちらもテキストの例だと省略出来る。

練習問題解答

https://gist.github.com/nozma/d5217d65876ddc527ad54c206671d92a

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

気付いたら買って本棚に飾ってあったので読んでいる。

カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)

準備

install.packages("vcd")
library(vcd)
data("Arthritis")

Arthritisは腱鞘炎に関する臨床試験のデータで、処置と経過、性別と年齢のデータが個票として入っている。

個票形式のデータの集計とプロット

xtabs()は表形式の集計を行う。集計の仕方を式の形で与えることができる。

> xtabs(~ Improved, data = Arthritis)
Improved
  None   Some Marked 
    42     14     28 

このままbarplot()に与えれば棒グラフが作成できる。

IMP2 <- IMP[order(IMP, decreasing = TRUE)]
barplot(IMP2)

クロス集計をしたい場合は、条件を+で繋ぐ。

> xtabs(~ Treatment + Improved, data = Arthritis)
         Improved
Treatment None Some Marked
  Placebo   29    7      7
  Treated   13    7     21

3変数以上は配列になる。

> xtabs(~ Treatment + Improved + Sex, data = Arthritis)
, , Sex = Female

         Improved
Treatment None Some Marked
  Placebo   19    7      6
  Treated    6    5     16

, , Sex = Male

         Improved
Treatment None Some Marked
  Placebo   10    0      1
  Treated    7    2      5

I()を用いると、特定の変数の組み合わせをまとめて行や列に指定できる。表形式に抑えたい場合に便利だ。

> xtabs(~ I(Sex:Treatment) + Improved, data = Arthritis)
                Improved
I(Sex:Treatment) None Some Marked
  Female:Placebo   19    7      6
  Female:Treated    6    5     16
  Male:Placebo     10    0      1
  Male:Treated      7    2      5

集計結果を比率で表示したい場合はprop.table()を用いる。margin=引数によって行和または列和を1にするように制御できる。

> prop.table(Ctable, margin = 1) # margin 1:行和を1 2:列和を1
         Improved
Treatment      None      Some    Marked
  Placebo 0.6744186 0.1627907 0.1627907
  Treated 0.3170732 0.1707317 0.5121951

prop.table()で比率に変換した表をbarplot()に渡すことで、帯グラフが作成できる。

> Ctable.p <- prop.table(Ctable, margin = 1)
> barplot(t(Ctable.p), horiz = T, axes = F,
+ xlab = "Improved(None Some Marked)")
> axis(1, seq(0, 1, by = 0.2), paste(seq(0, 100, by = 20), "%"))

多変数間のクロス集計結果を図示するにはモザイクプロットが便利である。

> Ctable3 <- xtabs(~ Treatment +  Sex + Improved, data = Arthritis)
> mosaicplot(Ctable3, col = TRUE)

なお、個票形式のデータであればxtabsにより集計する必要はなく、直接mosaicplot()に渡せばよい。

mosaicplot(~ Treatment + Sex + Improved, data = Arthritis,
           col = TRUE)

集計データの扱い

集計データとは、変数の組み合わせとその度数をまとめた形式のデータである。

> data("DanishWelfare")
> head(DanishWelfare)
  Freq Alcohol Income  Status         Urban
1    1      <1   0-50   Widow    Copenhagen
2    4      <1   0-50   Widow SubCopenhagen
3    1      <1   0-50   Widow     LargeCity
4    8      <1   0-50   Widow          City
5    6      <1   0-50   Widow       Country
6   14      <1   0-50 Married    Copenhagen

集計データをxtabs()により集計する際には、式の左辺に度数を指定する必要がある。

> xtabs(Freq ~ Alcohol, data = DanishWelfare)
Alcohol
  <1  1-2   >2 
2339 2083  722 

モザイクプロットを作成したい場合、xtabs()による集計を省略することはできないので注意する。

Ctable4 <- xtabs(Freq ~ Status + Alcohol, data = DanishWelfare)
mosaicplot(Ctable4, col = TRUE)

表形式→集計形式→個票

既に表形式になってしまっているデータを集計形式に戻すのは簡単で、単にdata.frame()に渡せばよい。

> dtable <- matrix(c(2, 18, 3, 17, 9, 11, 15, 5, 18, 2), nr = 2)
> dimnames(dtable) <- list(effect = list("あり", "なし"),
+                          dose = list("1", "2", "3", "4", "5"))
> class(dtable) <- "table"
> dose_effect <- data.frame(dtable)
> head(dose_effect)
  effect dose Freq
1   あり    1    2
2   なし    1   18
3   あり    2    3
4   なし    2   17
5   あり    3    9
6   なし    3   11

集計形式から個票形式に戻すにはapply()などを駆使する。

> DW <- data.frame(lapply(DanishWelfare,
+ function(i) rep(i, DanishWelfare[,"Freq"])))
> DW2 <- DW[-1]
> head(DW2)
  Alcohol Income Status         Urban
1      <1   0-50  Widow    Copenhagen
2      <1   0-50  Widow SubCopenhagen
3      <1   0-50  Widow SubCopenhagen
4      <1   0-50  Widow SubCopenhagen
5      <1   0-50  Widow SubCopenhagen
6      <1   0-50  Widow     LargeCity

カテゴリーの編集

カテゴリー変数を編集したい場合は、直接書き換えてしまうことができる。

> levels(DanishWelfare$Income)
[1] "0-50"    "50-100"  "100-150" ">150"   
> levels(DanishWelfare$Income) <- c("0-100", "0-100",
+ "100<", "100<")
> xtabs(Freq ~ Income, data = DanishWelfare)
Income
0-100  100< 
 2042  3102 

量的変数を大きさで区切ってカテゴリー変数にしたい場合は、cut()が便利だ。

> Arthritis$CatAge <- cut(Arthritis$Age, 
+                         breaks = 2:8 * 10,
+                         labels = 2:7 * 10,
+                         right = TRUE)
> head(Arthritis)
  ID Treatment  Sex Age Improved CatAge
1 57   Treated Male  27     Some     20
2 46   Treated Male  29     None     20
3 77   Treated Male  30     None     20
4 17   Treated Male  32   Marked     30
5 36   Treated Male  46   Marked     40
6 23   Treated Male  58   Marked     50

練習問題解答

https://gist.github.com/nozma/845a4e793b7dadf6302296a2441a4033

Rで複数条件抽出&集計

このようなデータフレームがあって、条件c1とc2に基づいて何らかの集計をしたいとする。
また、データフレームはdfというオブジェクトに代入されているとする。

tapply

> with(df, tapply(v, list(c1, c2), mean))
         1        2
a 1.689910 2.563432
b 1.780409 2.568772

by

> with(df, by(v, list(c1, c2), mean))
: a
: 1
[1] 1.68991
----------------------------------------------------------------------------- 
: b
: 1
[1] 1.780409
----------------------------------------------------------------------------- 
: a
: 2
[1] 2.563432
----------------------------------------------------------------------------- 
: b
: 2
[1] 2.568772

aggregate

> with(df, aggregate(v, list(c1, c2), mean))
  Group.1 Group.2        x
1       a       1 1.689910
2       b       1 1.780409
3       a       2 2.563432
4       b       2 2.568772

dplyr(group_by, summarize)

> library(dplyr)
> df %>% group_by(c1, c2) %>% summarize(mean(v))
Source: local data frame [4 x 3]
Groups: c1 [?]

      c1    c2  mean(v)
  <fctr> <int>    <dbl>
1      a     1 1.689910
2      a     2 2.563432
3      b     1 1.780409
4      b     2 2.568772

まとめ

Rでやるべき。

Excelで複数条件抽出&集計

このようなデータがあって、条件1と条件2に基づいて何らかの集計をしたいとする。

AVERAGEIFS

平均値を計算したいのであれば、AVERAGEIFS関数があるので、例えばこのように入力する。

=AVERAGEIFS(平均対象範囲,条件範囲1,条件1,条件範囲2,条件2)

最初の例で言えば、条件がこのように入力されているとして、

入力は次のようになる。

=AVERAGEIFS(C:C,A:A,F3,B:B,G3)

AVERAGEIFSは高速に動作するので、列全体を範囲指定してもあまり問題が無い。
同様の集計関数にはCOUNTIFS、SUMIFSなどがあり、いずれも高速である。なお、Excel2007以降で使える。

配列数式

平均、カウント、合計以外の集計をしたい場合、悪名高き配列数式を使うという手段がある。
配列数式は配列を渡して配列が返ってくる仕組みで、一人で使っている分には便利だが、他人にファイルを引き継ごうとすると配列数式の説明を漏れなく行う必要があり、大変面倒くさいという代物である。
AVERAGEやVAR.Sなどのセル範囲を引数として受け取る関数は、セル範囲の代わりに配列を渡しても良いので、IFを組み合わせて配列を作成し、それを渡せば良い。

=AVERAGE(IF(IF(A1:A10000=F3,TRUE,FALSE)*
            IF(B1:B10000=G3,TRUE,FALSE),
            C1:C10000,""))

IF文中のTRUE、FALSEは1、0でも良い。
この式をCtrl+Shiftを押下しながら確定すると、配列数式に変換されて式の両端に{}が付く。見た目には{}以外の変化が無いので、大変分かりにくい。
そして配列数式の最大の欠点として、「遅い」というものがある。上記例ではわざわざセル範囲を1万行に制限してあるが、これをしないと1秒くらいの計算時間がかかってしまう。配列数式を使う場合、範囲指定は狭ければ狭いほど良い。

名前付き範囲

この配列数式の欠点を補う方法として、名前付き範囲を使用するという方法がある。
Excelにおける名前とは変数名のようなもので、指定のセル範囲を指定の名前で表記できるというものである。例えば、前述の式を次のように記述できる。

=AVERAGE(IF(IF(条件1=F3,1,0)*
            IF(条件2=G3,1,0),
            値,""))

これで大分可読性が向上する。
名前はリボンの[数式]から[名前の定義]で定義する。
ここで、次のように入力する。

このようにOFFSET関数で参照範囲を定義することで、参照範囲が入力状況に対応して自動的に伸縮し、必要最低限の参照範囲となり、計算時間を節約できる。
なお、COUNTAで参照する列は、定義から想像出来ると思うが、全ての行に値が入っているような列でなければならない。MATCH関数等と組み合わせることで空白を許容することも可能であるが、長くなるのでここでは説明しない*1。
なお、名前には拡大倍率を40%以下にするとどこに名前を設定しているか一目で確認できるという便利機能があるが、上記の方法で設定した名前にはこの手段が使えない。

まとめ

Rでやるべき。

*1:詳しい方法についてはExcel Hacks 第2版― プロが教える究極のテクニック140選に記述があるが、多分適当にググっても出てくるだろう。

グループ変数に応じて複数の折れ線グラフを書く

緑色の本のp.119に載っている図と似たものを書こうとして、

d <- data.frame(q = rep(c("q0.1", "q0.3", "q0.8"), c(9, 9, 9)),
                p = c(dbinom(0:8, 8, 0.1), dbinom(0:8, 8, 0.3), dbinom(0:8, 8, 0.8)),
                y = rep(0:8, 3))

というようなデータを準備したものの、baseパッケージではどうプロットしたら良いものか若干悩んだ。

baseの例

# base
matplot(reshape(d, timevar = "q", idvar = "y", direction = "wide")[, -1], 
        ylab = "p", xlab = "y", xaxt = "n", type = "b", pch = 16, lty = 1)
axis(1, at = 1:9, labels = 0:8)
legend("topright", inset = 0.05,
       legend = c("q = 0.1", "q = 0.3", "q = 0.8"),
       col = 1:3,
       pch = 1,
       box.lty = 0)

linesで重ねる方法もあるけどいずれにせよ面倒だし使い回ししにくい。もっと良い方法はないものか。

続いてlatticeで描いてみた。

# lattice
library(lattice)
xyplot(p ~ y, data = d, group = q, type = "b", 
       auto.key = list(text =c("q = 0.1", "q = 0.3", "q = 0.8")))

凡例の編集をしようと思わなければ書式が素直で一番理解しやすい様に思う。

ggplot2

# ggplot2
# install.packages("ggplot2")
library(ggplot2)
ggplot(data = d, aes(x = y, y = p, col = q)) + geom_point() + geom_line() +
  scale_color_hue(label = c(q0.1 = "q = 0.1", q0.3 = "q = 0.3", q0.8 = "q = 0.8"))

latticeとさほど手間は変わらないけど、ggplot2はどうも覚えにくい気がする。初心者にaes()とかどう教えたら良いのかと毎回悩む。

RStudioのキーバインド変更

Windowsでのお話。

Tools -> Modify Keyboard Shortcuts...

で出来るようになっているという話だけど、Ctrl+hをBackSpaceにしたいとかCtrl+zをPgUpにしたいとか思ってもどうもうまく出来ないので結局AutoHotkeyを使うことにした。ほんの少しゴニョゴニョしたいだけならAutoHotkeyが一番手軽な気がする。

1. AutoHotkeyを入手&インストール(AutoHotkey)
2. 拡張子.ahkなファイルに下記のスクリプトを記述して保存後ダブルクリック

#IfWinActive ahk_class Qt5QWindowIcon
^h::Send {BS}
^z::Send {PgUp}
#IfWinActive

ahk_classはAutoHotkeyに同梱のActive Window Info (Window Spy) で調べる。