myscratchpad

# 4種類の棒グラフ
barplot(1:4, angle = 45 * 0:3, density = 15, col = "black")

# 8種類の棒グラフ
barplot(1:8, angle = 45 * 0:7, , density = c(rep(15, 4), rep(25, 4)), col = "black")

reshape関数とstack関数を使うやり方がある
reshape関数はもともと入っている変数名を使える
stack関数は並べ替えるだけ

reshape関数
これはもともとaov用
spfp.q1 <- data.frame(
  a = factor(c(rep(1,20), rep(2,20))),
  b = factor(rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5)),2)),
  s = factor(c(rep(1:5, 4), rep(6:10, 4))),
  result = c(3,3,1,3,5, 4,3,4,5,7, 6,6,6,4,8, 5,7,8,7,9,
             3,5,2,4,6, 2,6,3,6,4, 3,2,3,6,5, 2,3,3,4,6))
comment(spfp.q1) <- "2way, mixed (1within, 1between), balanced"
## データ解析テクニカルブックより、松井先生のコード。感謝

# 横長にする
spfp.q1.2 <- reshape(spfp.q1, idvar = c("a","s"), timevar = "b", direction = "wide")
colnames(spfp.q1.2) <- c("a", "s", "b1", "b2", "b3", "b4")

# reshape関数で元の縦長にする
spfp.q1.3 <- reshape(spfp.q1.2, varying = list(3:6), idvar = "s", timevar = "b", direction = "long")
## varyingで縦に並べる変数 (反復測定変数) を指定する。idvarで被験者要因。timevar は反復測定変数の水準名の変数
spfp.q1.3 # "b"と"s"が入れ替わっているが気にしない

stack関数
stack関数を使うと反復測定要因を縦に並べたデータフレームを簡単に作れる

personality-projectの説明がわかりやすい
recall.data <- read.table("http://personality-project.org/r/datasets/recall1.data",header=TRUE)
sums <- recall.data[,9:12]   #縦に並べる列をとりだす。
sums # データチェック
stackeds <- stack(sums) #stack関数で並べかえる
stackeds #データチェック
後は変数名を入れてるだけ
numcases <- 27 #被験者数
numvariables <- 4 #反復要因の水準数
## 後は変数、水準名を整理しているだけ
recall.df=data.frame(recall=stackeds,
 subj=factor(rep(paste("subj", 1:numcases, sep=""), numvariables)),
 time=factor(rep(rep(c("short", "long"), c(numcases, numcases)), 2)),
 study=factor(rep(c("d45", "d90"), c(numcases*2, numcases*2))))

Quick-Rの解説が詳しい。以下の数式の画像ファイルは同サイトから

基本は、どの関数もサンプルサイズ (n)、効果量等 (d, f, r...) 、有意水準 (第1種の過誤. sig.level) 、検定力 (1 - 第2種の過誤. power) の引数をもつ。いずれかを指定しないことで、望む値が得られる。

対応のないt検定で、被験者20、効果量が1、有意水準を.05にだったときの検定力を知りたい (南風原, 2002, p. 166-177)
効果量は以下で計算



pwr.t.test(n = 20, d = 1, sig.level = 0.05, type = "two.sample") # power = を 指定しない
     Two-sample t test power calculation

              n = 20
              d = 1
      sig.level = 0.05
          power = 0.868953
    alternative = two.sided

 NOTE: n is number in *each* group

対応のないt検定で、効果量が1のとき、有意水準を.05、検定力.80にしたいときの被験者数を知りたい
pwr.t.test(d = 1, sig.level = 0.05, power = 0.80, type = "two.sample") # n = を指定しない
     Two-sample t test power calculation

              n = 16.71473
              d = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

 NOTE: n is number in *each* group

  type = "one.sample"とすると1標本の平均値を調べられる。pairedと同じく対応のあるデータで用いる。one.sampleは差得点が0かどうかの検定を行うものだが、出力される結果はpairedと同じ。ただし、0じゃない値との検定もできるはできる。
つまり、0と異なるかどうかのt検定=pairedのt検定ということ
alternativeで両側、片側の指定

以下、それぞれの関数の覚書。というかQuick-Rのパクリ
1要因の分散分析
pwr.anova.test(k = , n = , f = , sig.level = , power = )
k: 群の数、n: 各群のサンプルサイズ、f: 効果量。



多元配置分散分析は川端先生のサイトに解説あり。感謝

相関係数
pwr.r.test(n = , r = , sig.level = , power = )
n: サンプルサイズ、r: 相関係数
## 南風原 (2002), p. 144-145と同じ結果になるやつ -> pwr.r.test(n = 20, r =.4 , sig.level =.05 )

線形モデル
pwr.f2.test(u =, v = , f2 = , sig.level = , power = )
u, v: 自由度、f2: 効果量



比率の検定
pwr.2p.test(h = , n = , sig.level =, power = )
h: 効果量



カイ二乗検定
pwr.chisq.test(w =, N = , df = , sig.level =, power = )
w: 効果量、N: 全体のサンプルサイズ、df: 自由度



cohen.ES関数で、dとかf2とかの"大きい"とかいわれる数値を知ることができる。
別に自分の実験の効果量を調べられるわけではない
cohen.ES(test="r",size="medium")

帰無分布で有意になる数値が対立仮説の分布で得られる確率を検定力という
帰無分布とはH = 0とする普通の検定 (フィッシャー式) で使われる分布
対立分布とは (そんな呼び方はないかもしれないが) 対立仮説としてH = 5くらいを考えたとしたら母集団の数値が5のサンプリング分布。通常、実際に実験や調査を行って算出した標本の平均値やら相関やらを用いる

検定力が低いとは、結局サンプルサイズが小さいか分散が大きくてサンプリング分布の標準誤差が小さい、つまり母数の推定精度が低いということ。ゆえに、標本の数値が帰無仮説を棄却しても、対立仮説の分布上ではそれより小さい (有意にならない) 数値が出てくる、すなわち第2種の過誤を犯す可能性が高いということ。

検定力を上げるとはサンプリング分布の標準誤差を小さくして、幅を狭くする。そうするともう一回サンプリングして棄却域より小さい数値が得られる確率は小さくなる。
サンプリング分布の標準誤差は大体分母にサンプルサイズを入れるので、サンプルサイズを大きくすれば検定力は高くなる。

なお、原理的には例えば平均値なら標準偏差を小さくすれば標準誤差 (分子だから) も小さいので検定力は高いし、また群Aと群Bの差が大きければ分布の幅が広くても対立分布で棄却域より小さい値が得られる確率は小さくなる。


検定力分析の使用法

• 実験・調査で得られた数値を母数とするサンプリング分布で、棄却域より大きい数値が得られる確率 (第2種の過誤を犯さない確率) を調べる。高いと検定力が高いということ
• 実験・調査をやる前に、母集団の平均値とかの差はこのくらいにだろうと予想する。サンプリングの数値にはばらつきがあるので、きちんと想定した母集団の数値が得られるサンプルサイズを決める。ぶっちゃけていえば、サンプリング分布の幅を小さくして有意になる (帰無仮説を棄却し、かつ対立仮説の面積を占める) サンプルサイズを決める (いつも思うんだけど、最初から平均値差とか相関を予測しろって、先行研究とか予備調査から調べておけってことか？仮にそうしたとしても、そんなに簡単に調査の結果なんか予測できないだろう。みんなできるの？)

pwrパッケージで検定力分析をしよう、 と思ったらなぜかCRANから削除されていた。
http://cran.r-project.org/web/packages/pwr/index.html
アーカイブのgzというのはLinux用のものだし、困ったと思ってたら
昔のパッケージやらR本体やらは全部補完されていた。便利だのう
http://cran.us.r-project.org/bin/windows/contrib/2.8/pwr_1.1.zip
http://cran.us.r-project.org/bin/windows/contrib/

参考は
http://www.statmethods.net/stats/power.html
川端先生のサイト