この月の記事一覧です。

2009年04月27日

回帰係数(傾き)の多重比較

*お決まりの「作法」として"回帰係数の多重比較"というモノがあるわけではない,と思います.そのまま鵜呑みにしてえらい人に怒られても知りませんよ.

共分散分析

共分散分析を使うと,共変量でもって調整したうえで分散分析をすることができる.たとえば次のようなデータ.

放牧がその後の種子生産に与える影響を調べたデータで,それぞれ初期の根の直径と種子生産量によって散布図が描いてある(統計学:Rを用いた入門書より.データはココ).

このとき,単純に種子生産量の平均値のみを2つの処理間で比べると,「放牧した方が種子生産量が多くなる」という結果が得られる.

しかし散布図から判断すると,それは放牧区に初期の個体サイズが大きなサンプルが集中していた結果のように見える.たとえば回帰直線を2本引いたら,明らかに放牧区の線は「下」にくるだろう.

それで,2本の回帰直線は違うのかとかそういうことをやるのが共分散分析.

これをやるにはそれぞれの回帰直線が平行という仮定が必要になる.要するにその回帰直線を平均値とし,回帰直線からの偏差を分散として分散分析をするんだから.

なので共分散分析は「平行性の検定(回帰の同質性の検定)」「要因効果の検定(共分散分析)」という2つのステップを踏む.

Rで

Rでは次のように行う.まず平行性の検定.

summary.aov(lm(Fruit~Root*Grazing))

ここでFruitは種子生産量,Rootは根の直径,Grazingはグループを表すベクトル.これにより次のような分散分析表が得られる.

             Df  Sum Sq Mean Sq  F value    Pr(>F)    
Root          1 16795.0 16795.0 359.9681 < 2.2e-16 ***
Grazing       1  5264.4  5264.4 112.8316 1.209e-12 ***
Root:Grazing  1     4.8     4.8   0.1031      0.75    
Residuals    36  1679.6    46.7                 

ここでRoot:Grazingとなっている部分は交互作用を表しており,ここが有意でなければ各回帰直線は平行であるということがいえる.この例ではp値が0.75なので問題ない.

共分散分析は次のように行う.

summary.lm(lm(Fruit~Root + Grazing))

モデル式の書き方がちょっと違う.これにより次のような結果が得られる.

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -127.829      9.664  -13.23 1.35e-15 ***
Root              23.560      1.149   20.51  < 2e-16 ***
GrazingUngrazed   36.103      3.357   10.75 6.11e-13 ***

ここでGrazingUngrazedとなっている部分のEstimateが2つの回帰式がどれだけ離れているかを表しており,ここのp値が有意なら2つの回帰式は有意に離れている(つまり切片が有意に異なる)ということが言える.

(もうちょっと説明をしたいところだけど正直イメージがうまく描けてないのでパスです.青木先生のところなんかをよーく読んだら分かるかもしれません.僕は流し読んだのでわかりません><)

回帰式の多重比較

共分散分析では平行性の検定を行う.で,ここで棄却されなければ次へ進める.

逆に言うと,平行性の検定が棄却されるということは切片はともかくとして2つの回帰式の傾きが有意に異なるということ.

というわけで,群が2つしかない場合にそれらの回帰式の傾きを比べる方法については解決だ.上に書いた方法で平行性を調べたらいい.また,「どこかに差がある」ということが分かればいいんだったらやっぱ同じ方法でいい.

でもやっぱ「どことどこに差があるの?」ってことが気になる.

というわけで,群の中から2つを選び総当たりで平行性の検定を行うことにしよう.そのとき問題になるのは検定の多重性というやつ.

Holmの方法による多重比較

たとえば群が3つあったら対比較は3回やらないといけない.危険率が5%として,すべての群が等しいのにどこかに差がある判定をしてしまう確率は1-0.95^3で0.14くらいになる(実際はちょっと違うらしいけど).

それで,全体として有意水準5%で検定をしたいときに,対比較の回数に応じて1回1回の対比較における有意水準を低く調整してやることで何とかしようという多重比較法がある.ボンフェローニ(Bonferroni)の方法をはじめとする有意水準調整法というやつだ.

今回はホルム(Holm)の方法というやつにする.これはそれぞれの対比較で得られたp値を小さい順に並べて,i番目のp値に対する有意水準をα/(k-i+1)にするという方法.kは対比較の回数,αは全体としてそれ以下におさめたい有意水準.ボンフェローニの方法だとすべての有意水準をα/kに調整してしまうので過剰に保守的になってしまうけど,ホルムの方法ではもうちょっとゆるい感じになっている.ゆるいといっても全体としての有意水準はα以下に保たれる.

で,やることは決まった.対比較を繰り返して,Holmの方法で有意水準を調整するんだ.

Rで

関数にしてみました.

pairwise.covar <- function(
                           data.f,       #データフレーム
                           g.i=1,        #グループ変数の列
                           i.i=2,        #独立変数の列
                           d.i=3         #従属変数の列
                           )
{
  group <- data.f[[g.i]]                  #グループ変数
  g.name <- levels(group)               #水準名
  independent <- data.f[[i.i]]            #従属変数
  dependent <- data.f[[d.i]]              #独立変数
  g.n <- length(levels(group))          #グループ数
  k <- g.n*(g.n-1)/2                    #組み合わせパターン数
  pvals <- numeric(k)                   #p値用ベクトル
  cmb <- combn(g.name,2)                #全組み合わせ
  cmb.name <- character(k)              #組み合わせの名前(出力用)
  sign <- character(k)                  #アスタリスク用
  for(i in 1:k){
    pvals[i] <- round(anova(
                            lm(dependent[group==cmb[1,i] | group==cmb[2,i]]~
                               independent[group==cmb[1,i]|group==cmb[2,i]]*
                               group[group==cmb[1,i]|group==cmb[2,i]]
                               )
                            )$Pr[3],7)          #p値の計算と格納
    cmb.name[i] <- paste(cmb[1,i],":",cmb[2,i]) #組み合わせ名格納
    if(pvals[i] < 0.001){
      sign[i] <- "***"                  #0.1%以下は***
    }else if(pvals[i] < 0.01){
      sign[i] <- "**"                   #1%以下は**
    }else if(pvals[i] < 0.05){
      sign[i] <- "*"                    #5%以下は*
    }else if(pvals[i] < 0.1){
      sign[i] <- "."                    #10%以下は.
    }else{
      sign[i] <- " "                    #それ以外は" "
    }
  }
  o <- order(pvals)
  cbind(cmb.name[o],
        p.value=(pvals*(k-order(pvals)+1))[o],
        sign[o])                        #結果のベクトル
}

アスタリスクとか飾りだけどさ,ほら,見やすさのために…

使い方は

pairwise.covar(データフレーム, 
				グループ変数の列,
				独立変数の列,
				従属変数の列)

列の指定は省略可.デフォルトでそれぞれ1, 2, 3.

使ってみる.まずデータの作成.

## 仮データの作成
set.seed(5)
x <- 1:20
y1 <- x+rnorm(20)
y2 <- x*2+rnorm(20)
y3 <- x+10+rnorm(20)
group <- as.factor(rep(1:3, rep(20,3)))
independent <- rep(x, 3)
dependent <- c(y1, y2, y3)
mydata <- data.frame(group, independent, dependent)

xが独立変数,y1~3が従属変数.プロットするとこうなる.

凡例忘れた...○がy1,△がy2,+がy3.

見て分かるようにy1とy3は切片が異なるけど傾きは同じ.y2だけ傾きが違う.

これに対してさっきのを使うと

> pairwise.covar(mydata)
             p.value         
[1,] "1 : 2" "0"        "***"
[2,] "2 : 3" "0"        "***"
[3,] "1 : 3" "0.639056" " "  

となる.ちゃんと計算できてる?

ちなみに有意水準を減らす代わりにp値を増やしているので場合によってはp値が1を超えるとかいうことになるかもしれない.あとp値が考えている有意水準より高いからって傾きが等しい=平行というのにはちょっと無理があるかも.ボンフェローニよりマシとはいえホルムの方法も保守的なのには違いがないので.

posted by Rion778 at 21:44 | TrackBack(0) | 勉強ノート | このブログの読者になる | 更新情報をチェックする

広告


この広告は60日以上更新がないブログに表示がされております。

以下のいずれかの方法で非表示にすることが可能です。

・記事の投稿、編集をおこなう
・マイブログの【設定】 > 【広告設定】 より、「60日間更新が無い場合」 の 「広告を表示しない」にチェックを入れて保存する。


×

この広告は1年以上新しい記事の投稿がないブログに表示されております。