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

No.04591　３群以上の共分散分析　　【おすぎ】　2007/10/28(Sun) 14:58

3群のデータを使って共分散分析をしようとしています。
Rなどのソフトを使って解析しようとしても，2群の群間比較はできるのですが，3群以上のデータでは解析できなさそうです。
もし，3群以上でできる方法があれば教えてください。
それとも，この解析は統計的に間違っているのでしょうか？

No.04592　Re: 3群以上の共分散分析　　【波音】　2007/10/29(Mon) 10:59

できますよ。

例えば，glm()を使用して

> result <- glm(Y ~ X + G) #Xは連続型データでGはカテゴリカル型データ
> summary(aov(result)) #分散分析表の出力

としてあげれば解析してくれます。2群間の比較なら変数Gは

> G <- as.factor(c(1,1,1,2,2,2)) #カテゴリカル型データとして指定する

としますし，3群なら

> G <- as.factor(c(1,1,1,2,2,2,3,3,3))

としてあげるだけです。

例題としてhttp://www.oup.com/uk/orc/bin/9780199252312/01resources/datasets/excel/separate/にあるDatasheets Chapter 06というExcelワークシートにあるLeprosyという名前のSheetにあるデータを同じように解析すると次のようになります。
> dat <- excel.w(3) #青木先生の所から拝借した関数を使用して取り込み
Read 90 items
# データフレームの作成
> Leprosy <- data.frame(TREATMT=as.factor(dat[,1]), BACBEF=dat[,2], BACAFTER=dat[,3])
> attach(Leprosy) # 作成したデータフレームにパスを通す
> result <- glm(BACAFTER ~ BACBEF + TREATMT)
> summary(aov(result))
            Df Sum Sq Mean Sq F value    Pr(>F)    
BACBEF       1 587.48  587.48 40.1800 1.033e-06 ***
TREATMT      2  83.35   41.67  2.8502     0.076 .  
Residuals   26 380.15   14.62                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

No.04593　Re: 3群以上の共分散分析　　【波音】　2007/10/29(Mon) 11:04

一応，グラフを描くのなら次のようにしてできます（もっと少ない行で描けるのかもしれませんが，，，）。


# plot()を用いて処置1のデータをプロットする
plot(BACBEF[TREATMT==1], BACAFTER[TREATMT==1],
xlim=c(0,25), ylim=c(0,20),
xlab="BACBEF", ylab="BACAFTER",
pch=18)

# points()を用いて処置2のデータを図に追加
points(BACBEF[TREATMT==2], BACAFTER[TREATMT==2], pch=17)

# 同じく処置3のデータをプロット
points(BACBEF[TREATMT==3], BACAFTER[TREATMT==3], pch=8)

x <- seq(0, 25, 0.1) # 直線の描写につかうx軸の値
model1 <- -1.6026 + 0.8831*x + 0 # モデル1のy軸の値
points(x, model1, type="l") # モデル1の直線を描く
model2 <- -1.6026 + 0.8831*x + 0.8641
points(x, model2, type="l", lty=2)
model3 <- -1.6026 + 0.8831*x + 3.9060
points(x, model3, type="l", lty=3)

# 直線の凡例
legend(17, 5, # 凡例を書く座標
        c("model1", "model2", "model3"), # 凡例のラベル
        lty=c(1, 2, 3), # 凡例のラインのスタイル
)

# 点の凡例
legend(0, 20,
c("treatmt1", "treatmt2", "treatmt3"),
pch=c(18, 17, 8) # 点のスタイル
)

No.04594　Re: 3群以上の共分散分析　　【青木繁伸】　2007/10/29(Mon) 15:29

一応，凡例を描く前までを
plot(BACBEF, BACAFTER,
xlim=c(0,25), ylim=c(0,20),
xlab="BACBEF", ylab="BACAFTER",
pch=c(18, 17, 8)[TREATMT]) # 群別に別記号で
p <- coefficients(result) # 係数を取り出す
abline(p[1], p[2], lty=1) # abline で
abline(p[1]+p[3], p[2], lty=2)
abline(p[1]+p[4], p[2], lty=3)

No.04596　Re: 3群以上の共分散分析　　【波音】　2007/10/29(Mon) 18:13

私のように群のデータごとにpoints()を使用して描かなくても，最初のplot()のpch=でこのように指定すれば簡単だったのですね。直線を描く場合も係数を取り出してabline()を利用した方が簡単ですね。勉強になりました。