薬指使えないだけで、こうも作業効率が下がるとは、、、

R 統計

前回のおさらいから
odds ratioの比較に、ただただ区間推定して比較してみたが、vcdパッケージには、woolf検定（層ごとにオッズ比は同じという帰無仮説を検定）をしてくれるパッケージもある。

library(vcd)
woolf_test(data2)

         Woolf-test on Homogeneity of Odds Ratios (no 3-Way assoc.)

data:  data2 
X-squared = 11.1711, df = 1, p-value = 0.0008308

と層ごとに違うんじゃないかという結果が返される。
woolf_testってアンダーバーが関数に含まれていることに気付かず、何回か

エラー：  関数 "woolf.test" を見つけることができませんでした

という、表示をしてしまった。


さて、いわゆる３次元表を一般化線形モデルを使って、Rでやるためには２つの方法がある。

  • ポワソン分布を誤差分布として、解析
  • （おそらく最も）注目したいカテゴリを目的変数とみなし、その変数を２値応答数として(cbind関数を使って）解析

まずはポワソン分布を誤差分布とした解析から。

データをmatrix構造からdata frame構造にする（関数があるかもしれませんが、とりあえず自力でやった）

count <- c(45, 13, 27, 21, 19, 17, 39, 11)
pop <- c(rep("tokyo", 4), rep("osaka", 4))
spi <- c(rep(c("saku", "saku", "ume", "ume"), 2))
state <- c(rep(c("ok", "out"), 4))
data3 <- data.frame(count, pop, spi, state)             
data3
  count   pop  spi state
1    45 tokyo saku    ok
2    13 tokyo saku   out
3    27 tokyo  ume    ok
4    21 tokyo  ume   out
5    19 osaka saku    ok
6    17 osaka saku   out
7    39 osaka  ume    ok
8    11 osaka  ume   out
mosaicplot(xtabs(count ~ pop + spi + state, data3),
		color = TRUE, main = NULL, las=1)


これを単純にポワソン分布を誤差構造にした、３次の交互作用を含めたモデル（いわゆる飽和モデル）を当てはめると、

fitglm <- glm(count ~ pop * spi * state, data = data3,
               family = poisson)
fitglm

Call:  glm(formula = count ~ pop * spi * state, family = poisson, data = data3)

Coefficients:
             (Intercept)                  poptokyo                    spiume  
                  2.9444                    0.8622                    0.7191  
                stateout           poptokyo:spiume         poptokyo:stateout  
                 -0.1112                   -1.2299                   -1.1305  
         spiume:stateout  poptokyo:spiume:stateout  
                 -1.1544                    2.1448  

Degrees of Freedom: 7 Total (i.e. Null);  0 Residual
Null Deviance:      41.49 
Residual Deviance: -1.599e-14   AIC: 55.32

となる。すべてのパラメータを使って推定しているため、自由度は0である。これを高次の交互作用を取り除いたモデルと逸脱度分析すると、

fitglm2 <- update(fitglm, ~. - pop:spi:state)
anova(fitglm, fitglm2, test = "Chi")
Analysis of Deviance Table

Model 1: count ~ pop * spi * state
Model 2: count ~ pop + spi + state + pop:spi + pop:state + spi:state
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1         0       0.00                          
2         1      11.54 -1   -11.54 0.0006812 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

交互作用が重要な説明を与えている、つまり東京と大阪で、種の罹病率が異なる（共通オッズ比が異なるという前回の結果と一致）のではないかと考察することができる。

（続く）