「統計学関連なんでもあり」の過去ログ--- 048

No.22179　ロジスティック回帰、Tukey法　　【ぐれいぷ】　2016/10/19(Wed) 09:35

初めて質問します。よろしくお願いします。
生物系の実験をしていて，以下のようなデータの解析をしたいと考えています。

> a <- read.table("clipboard",header=T)
> a
   med tree  a  N
1    A    P  2 16
2    A    Q  1 15
3    A    R  1 16
4    A    S  0 17
5    B    P 10 16
6    B    Q  9 15
7    B    R  5 16
8    B    S  8 17
9    C    P  5 16
10   C    Q  4 15
11   C    R  4 16
12   C    S  0 17

（4本の樹(P,Q,R,S)の葉に3種の薬剤(A,B,C)を処理して，処理した葉N枚のうち，a枚で障害が出た，というデータです。）

　薬剤の効果を比較したいと思い，ロジスティック回帰後にRのmultcompのパッケージを用い，Tukey法で検定しました。

> med <- a$med
> r <- glm(cbind(a$a,a$N-a$a)~med,family=binomial)
> mul <- glht(r,linfct=mcp(med="Tukey"))
> summary(mul)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: glm(formula = cbind(a$a, a$N - a$a) ~ med, family = binomial)

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)    
B - A == 0   2.7081     0.5737   4.720  < 0.001 ***
C - A == 0   1.3412     0.6026   2.225  0.06470 .  
C - B == 0  -1.3669     0.3988  -3.428  0.00169 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

　薬剤の効果は比較できるのですが，「1つの薬剤に関しては単純な4反復になっていて，上の解析では同じ樹を使っている効果が除外されているのでは？」と質問を受けました。

　薬剤の効果と樹の効果については，パッケージcarを使って，

> r2 <- glm(cbind(a$a,a$N-a$a)~a$med*a$tree,family=binomial)
> Anova(r2)
Analysis of Deviance Table (Type II tests)

Response: cbind(a$a, a$N - a$a)
             LR Chisq Df Pr(>Chisq)    
a$med          36.264  2  1.335e-08 ***
a$tree          7.878  3    0.04861 *  
a$med:a$tree    8.477  6    0.20523    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>

と解析してみたのですが，上記のTukeyで，樹の効果が分かるような解析方法はありますでしょうか？

No.22180　Re: ロジスティック回帰、Tukey法　　【太郎】　2016/10/19(Wed) 14:38

　単純に　a/N（障害発生葉の比率）をデータとして，薬剤3水準と樹4水準の乱塊法（繰り返しのない2元配置法）を適用するというのはいかがですか？
　
　比率型データなので角変換してから分散分析して，薬剤，樹の多重比較（Tukeyなど）をしてもよいと思います。

No.22181　Re: ロジスティック回帰、Tukey法　　【ぐれいぷ】　2016/10/19(Wed) 23:43

太郎様
　回答していただきありがとうございます。私も少し前まで，比率の検定→変換という方法を用いておりましたが，変換して分散分析するよりも一般化線形モデルを用いたほうがいいと聞きましたので，上記の方法を用いました。
　上記で樹の効果が分かる方法がなければ，角変換してみます。

No.22183　Re: ロジスティック回帰、Tukey法　　【鈴木康弘】　2016/10/20(Thu) 16:43

　すみません，multcomp というパッケージのことは初めて聞いたし，ダウンロードもしてないのですが，
> r <- glm(cbind(a$a,a$N-a$a)~a$med+a$tree,family=binomial)
　とやって glht コマンドを使っても，ご希望のようにいかないのでしょうか？

No.22184　Re: ロジスティック回帰、Tukey法　　【ぐれいぷ】　2016/10/21(Fri) 10:28

鈴木康弘様
　回答ありがとうございます。はい，多重比較の結果は出るのですが，その中にtreeの効果は出ません。

No.22185　Re: ロジスティック回帰、Tukey法　　【荒】　2016/10/21(Fri) 20:15

同じ樹から繰り返しデータを取っているため一般化線形混合モデルの適用ではないでしょうか。
ただしこの種の解析をしたことがないため，望みの結果が得られているかどうかは自信がありません。

> library(lme4)
> r <- glmer(cbind(a, N - a) ~ med + (1|tree), data = a, family = binomial)
> summary(r)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(a, N - a) ~ med + (1 | tree)
   Data: a

     AIC      BIC   logLik deviance df.resid 
    51.4     53.3    -21.7     43.4        8 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7409 -0.2627  0.3620  0.4520  0.7766 

Random effects:
 Groups Name        Variance Std.Dev.
 tree   (Intercept) 0.1319   0.3632  
Number of obs: 12, groups:  tree, 4

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.7576     0.5521  -4.994 5.90e-07 ***
medB          2.7654     0.5808   4.761 1.92e-06 ***
medC          1.3595     0.6053   2.246   0.0247 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr) medB  
medB -0.849       
medC -0.805  0.766

> library(multcomp)
> mul <- glht(r, linfct = mcp(med = "Tukey"))
> summary(mul)

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: glmer(formula = cbind(a, N - a) ~ med + (1 | tree), data = a, 
    family = binomial)

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)    
B - A == 0   2.7654     0.5808   4.761  < 1e-04 ***
C - A == 0   1.3595     0.6053   2.246  0.06171 .  
C - B == 0  -1.4059     0.4065  -3.458  0.00151 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

No.22186　Re: ロジスティック回帰、Tukey法　　【ぐれいぷ】　2016/10/21(Fri) 22:22

荒様
　回答ありがとうございます。Rを使って解析していただいて助かります。一般化線形混合モデルは今まで扱ったことがないので知りませんでした。勉強してみます。
　今後ともよろしくお願いいたします。

No.22187　Re: ロジスティック回帰、Tukey法　　【荒】　2016/10/22(Sat) 20:02

後で気が付いたのですが，このような対応のあるデータにTukey法を適用できないのではと思いました。

ですので，この場合であれば
A-B, A-C, B-C の3通りの組み合せで解析を行い，その後P値をHolm法で調整するのがいいのではないでしょうか。

combn関数で各組み合わせが配列として得られるため，
以下のプログラムで検定結果が得られるはずです。
このままでは結果の表示が貧弱なため，もう少し改変が必要ですが。
comb <- combn(unique(a$tree), 2)
n <- length(comb)
p <- numeric(n) # P値を入れる変数

for (i in seq(n)) {
  tmp <- subset(a, tree %in% comb[, i])
  r <- glmer(cbind(a, N - a) ~ med + (1|tree), data = tmp, family = binomial)
  p[i] <- summary(r)$coefficients[2, 4] # P値
}

p.adjust(p) # デフォルトではHolm法で補正

No.22188　Re: ロジスティック回帰、Tukey法　　【ぐれいぷ】　2016/10/24(Mon) 14:40

荒様
　たびたびお答え下さりありがとうございます。compなど，今まで使ったことが
ないため，勉強して使ってみます。