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) | 勉強ノート | このブログの読者になる | 更新情報をチェックする
2009年04月14日

MeadowでYaTeXを使うための設定

;;YaTeX
(setq auto-mode-alist (cons (cons "\\.tex$" 'yatex-mode) auto-mode-alist))
(autoload 'yatex-mode "yatex" "Yet Another LaTeX mode" t)
;dvipdfmxをプリントコマンドに
(setq dviprint-command-format "dvipdfmx %s")

Meadowのインストール時にパッケージとしてYaTeXを選択していればこれを.emacsに書き足すだけでいいはず.

もちろん,LaTeX環境は別途準備してある前提で.

拡張子.texのファイルを読み込むとYaTeXが起動し,モードラインに「やてふ」と表示される.

ついでに[C-t l]をdvipdfmxに割り当てた..dvi作成後に実行で.pdf作成.

posted by Rion778 at 17:32 | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする
2009年04月12日

anthyやindian-1-columnがないってMeadowが怒る

なんか前回の設定で動くっちゃ動くんだけど起動時に

Searching for program:no such file or directory, anthy-agent

とか出てauto-autoloads.elなどの読み込みが止まる.止まってもなんか普通にESSとか動くんだけどIMEのON/OFFを切り替えるたびにミニバッファに上の警告が出てエラー音を発する(IMEのON/OFFはちゃんと切り替わる).

害は今のとこないんだけどなんとかしたかったのでなんとかした.

anthyを切る

IME使うようにしてるのでそもそもいらない.

Meadowのauto-autoloads.elをいくつか止めて快適に - わだいのたけひこのにっき

\Meadow\packages\pkginfo\anthyにあるauto-autoloads.elをリネームして読み込まないようにした.

これでanthy-agantないよ!って怒られなくなった.[-Aあ-\--]とかも出るようになった.起動時に読み込むもの増えた.

でも今度は

error: Invalid value for :char-spec :indian-1-column

とか起動時に言われる.ぐぬぬ…

indian-1-columnをコメントアウト

Meadow起動時のエラー - 教えて!goo

ググってるときに教えて!gooの「一般人」の適当な回答が引っ掛かったときの怒りを表現する言葉が欲しい.

とにかく,indian-1-columnというcharsetがもうないので,読み込まないようにしろと,そういうことらしい.

それでどこに書いてあるんだ,font-setup.elとか検索してもないぞ.と思ったらauto-autload.elの中だった.

Windows Meadow 2.10 多言語環境の設定 - NOX INSOMNIAE

\Meadow\packages\pkginfo\intlfontsの中のauto-autoloads.elを開くと,中に

(indian-1-column ("ind1c16-mule.bdf"));; MuleIndian-1

という記述があるのでコメントアウト.てか消していい気がする.僕はチキンなので消さない.

起動が重くなった

真面目にいろいろ読み込むせいか起動が遅いぞ.いろいろ読み込まないようにしたらいいのだろうか.これから調べる.

*パッケージ減らしたらいくらか早くなった

posted by Rion778 at 15:53 | Comment(0) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする
2009年04月08日

Meadow+ESSでR(Emacsとか分からないWindowsユーザー向け)[導入編]

2010/09/20追記 以下の記事の記述はやや古くややいい加減です。Emacs + ESSでR(for Windows users) - もうカツ丼でいいよなにちょっと内容を見なおしたものを上げてありますので、よろしければそちらを御覧ください。

Emacsのキーバインドを覚える手間を差し引いてもESSを導入する価値はあると思ったのでちょいとメモ.「Rはちょっとくらい使えるけど,Emacsとかよく分からない」というWindowsユーザー向け(というかそもそもMacやLinux使う人はEmacs使えるだろうしESSも自力でインストールできると思うので,RjpwikiのESSの項目なんかを参考にするといいと思います).

ESSとは?

ESS (Emacs Speaks Statistics) は Emacs や XEmacs から R などの統計解析アプリケーションを便利に使ってしまおう,という Lisp プログラムです.

ESSがどれだけれ能率向上をもたらしてくれるものか,などは「The R Book」p.55--56に叙述されています.(Rjpwiki ESS)

要はEmacsという非常に高機能なテキストエディタがあって,その上でRを使うとき便利なようにいろいろカスタムするプログラムがESSということ.

何を覚える必要があるのか?

  1. R
  2. Emacs
  3. ESS

3つだけ.簡単です.まずは準備を.(ただ最初に書きましたがRの使い方だけはどこかよそで予習しておきましょう.R自体はESSなくてもRguiで動くので.僕に聞くのはいろいろ間違いです><)

R,Meadow,ESSのインストールと設定

.exeクリックで一発!とはいかないのでちょっとがんばる.

Rのインストール

最新版はここにあるはずなので,「Download R 2.8.1 for Windows」などというリンクから.exeを落としてきて実行すればインストールが開始される.

途中の選択肢はESSを使う上ではどう選んでも関係ないはずだけれど,インストールしたフォルダは忘れないように.

Meadowのインストール

Windows用のEmacs実装としてMeadowというものがある.今回はこれにESSを入れて使う.Meadowのページから「ダウンロード」というリンクを辿ると「Netinstall packages」というものが見つかるので,「setup-ja.exe」を落とす.

実行する前に,Cドライブ直下に「Meadow」などの名前でフォルダを作っておき,そこへsetup-ja.exeを移動.

setup-ja.exeを実行してインストール.

  • ローカルパッケージディレクトリはC:\Meadow(作成したフォルダ.実行したフォルダへのパスが自動で入るはず)
  • インストール先ディレクトリもC:\Meadow(同じく変更の必要はない)
  • パッケージはデフォルトでいい気がする

最後にsetup.exeを実行する必要がある.実行すると「どこにある「.emacs」を読むようにする?」などと聞いてくる.画面には[C:\Meadow]などと出ているはず.何も入力せずにEnterを叩けばそこにある.emacsを読むようになる.書き換えてもいい.ただし場所は覚えておく

これでインストールは終了.デスクトップにショートカットが出現し,それをダブルクリックでMeadowが起動するようになっているはず.

.emacsの作成

Meadowは起動時に.emacsというファイルに書き込まれた設定を読み込む.なのでこれを作らないといけないのだが,「新規作成→テキスト→名前を.emacsに」しても「ファイル名を入力して下さい!」と怒られる.

なのでMeadowから作成する.Emacsの操作が分からない人はまじないの如く以下の作業を行う.

  1. Meadowを起動
  2. Ctlr+X,Ctlr+Fと連続してタイプ
  3. 一番下に「Find file: ~/bin/」と出る
  4. 続けて,「Find file: ~/bin/.emacs」と入力してEnter
  5. 次に,Ctlr+X,Ctlr+Wと連続してタイプ
  6. C:\Meadow\binの中に.emacsというファイルができているはず

完成した.emacsをインストール時に設定した「.emacsを読む場所」に合わせて移動すれば終わり.もちろんわかる人は最初からそこへ作成すればいい.

作ってしまえばこっちのもので,あとはメモ帳で開けば編集できるし,保存もできる.Windowsは難しい.ともかく,Meadowの扱いがわからないうちはメモ帳で編集できるので安心.

ESSのインストール

ESSはたぶんここに最新版があると思う.DownloadからLatest Released Versionって所のzipを落としてきて解凍する.

解凍してできたフォルダをフォルダごと\Meadow\site-lispの中へ放り込む.

.emacsの編集

ESSを使えるようにするため,早速.emacsを編集する必要がある.なんかよくわからなかったらとりあえず次のように書いて保存.

(require 'ess-site)
(setq ess-ask-for-ess-directory nil)
(setq ess-pre-run-hook
	'((lambda ()
	(setq default-process-coding-system '(sjis . sjis))
	)))

(set-language-environment "Japanese")
(setq default-input-method "MW32-IME")
(mw32-ime-initialize)

(defun ess:format-window-1 ()
	(split-window-horizontally)
	(other-window 1)
	(split-window)
	(other-window 1))
(add-hook 'ess-pre-run-hook 'ess:format-window-1)

(setq default-frame-alist
      (append (list '(foreground-color . "azure3")
	'(background-color . "black")
	'(border-color . "black")
	'(mouse-color . "white")
	'(cursor-color . "white")
	)
	default-frame-alist))

ここのをベースにちょこちょこいじってついでにIMEが動くようにしてあります.)

環境変数の設定

ここまで終わってもまだ動かない.環境変数というものを設定する必要がある.わかっている人はPathへRとMeadowの実行ファイルがあるフォルダ(bin)を,HomeへMeadow\binを追加.以下分からない人へ.

環境変数の設定は,

  • XPの場合:コントロールパネル(クラシック表示)→システム→詳細設定タブ など
  • Vistaの場合:コントロールパネル(クラシック表示)→システム→左側にある「システムの詳細設定」 など

で見つけられるはず.「ユーザー環境変数」と「システム環境変数」があるけど,ユーザー個人に設定を適用するか,コンピュータ全体に適用するか程度の違いなので以下好きな方へ設定.

  1. 変数からPATH(or Path)を選択.なければ「新規」で変数名にPATHと入れる.
  2. 変数値にC:\Meadow\bin (RunMW32.exeのあるフォルダ)を追加.新規に作成した場合はそのまま記入.すでに何か書いてある場合,末尾に「;」を記入したのちに続けて記入する.

同様にしてPATHへRの実行ファイルがあるディレクトリ(C:\Program Files\R\R-2.8.1\bin など)を追加し,HOMEへC:\Meadow\bin(.emacsのあるフォルダ)を追加.

確認

拡張子が.RなファイルをMeadowで開くとESSモードで起動.ツールバーに「R」とかのボタンが見えるはず.

M-x R (Alt+x rもしくはEsc x rと順番にキーを叩く)すると,画面が分割されて右下でRが起動.

.Rなファイルが表示されている場所(バッファ)にカーソルがある状態でC-c C-b (Ctrl+c Ctrl+bの順に叩く)すると中身がすべて実行される.結果はRのバッファなどに.

うまくいかなかったら…

どうしよう?よくわかりません><

手順確認→Google先生に聞く しても分からなかったら聞いてください.でもきっとわかりません.

使い方を覚える

RjpwikiにESSを取り扱ったページがあって, ショートカットなんかは一覧がある(ESSモード)ので参考に練習しましょう.

で終わらせたいところだけど,Windowsユーザーにとっての鬼門はMeadow(Emacs)の独特の操作や単語, ショートカットの表記法を覚えるところだと思う.一覧表とか見せられても量がすごいので「覚えられない!無理!」ということになる.

Emacsの操作方法を覚える

Meadow(Emacs)にはチュートリアルが入っている。

ツールバーの[Help]から[Emacs Tutorial]を選ぶとチュートリアルが始まる(ただのテキストだけど)ので, これをやる.

案外長くて丁寧にやると30分くらいかかるけど,覚えるべきことはカーソル操作,消去と挿入,ファイル操作程度なので, 実際それほど大変じゃない.

1回だけじゃ忘れるかもしれないけど,3回もやったら基本的な操作はほとんど覚えられるはず.

「こんな変な操作覚えて大丈夫か?」などという心配があるかもしれないが,大抵のものはMeadowで事足りる, というかMeadowの方が便利なことが多いし,他のエディタがEmacsのキーバインドをサポートしていることも多い(Visual C++とか).なので心配せず慣れてしまっていいと思う.

ESSの扱いを覚える

操作自体はさっきも書いたRjpwikiのESSのページを参考に覚えたらいい. ショートカットの部分を拡大印刷して手元において置くとかして.

それよりも操作に慣れることが重要.ちょっと慣れると急激にありがたさが分かってくる.

これは実際にいくつかプログラムを書いてみるのが一番いい.

RでProject Eulerなんかいかがですか?

posted by Rion778 at 16:46 | Comment(0) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする

広告


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

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

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


×

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