返り値： クラス
"power.htest"
のオブジェクトで，method
とnote
成分が加えら れた，与えられた引き数(
計算されたものを含む)
を成分に持つリストパラメータ
n, delta, power, sd
そしてsig.level
の内丁度一つが値NULL
と指定さ れる必要があり，そのパラメータは他のパラメータから計算される．最後の二つはNULL
でない既定値を持っており，これらを計算したければ明示的にNULL
と指定する必要があ ることを注意しよう．もしstrict=TRUE
ならば，両側検定の場合，検出力は真の効果の 反対方向への棄却確率を含む．これがないと真の差が0
ならば検出力は信頼係数の半分 になるであろう．注意：未知パラメータを得るために検出力方程式を解くために
uniroot()
が使われ，し たがって，特に不正な引き数が与えられたとき根をくくり出せないという，エラーを得る かも知れない．関連：
t.test(), uniroot()
．> power.t.test(n=20, delta=1) # 観測値数20の時の検出力
Two-sample t test power calculation # 二標本t検定
n = 20 # 二標本のそれぞれのデータ数
delta = 1 # 真の平均差
sd = 1 # 二標本の標準偏差
sig.level = 0.05 # 第一種の誤り確率
power = 0.8689528 # 計算された検出力
alternative = two.sided # 両側対立仮説
NOTE: n is number in *each* group
> power.t.test(power=.90, delta=1) # 検出力0.90を得るために必要な観測値数 Two-sample t test power calculation # 二標本t検定
n = 22.02110 # 必要な両群の観測値数
delta = 1 sd = 1 sig.level = 0.05
power = 0.9 alternative = two.sided NOTE: n is number in *each* group
> power.t.test(power=.90, delta=1, alt="one.sided")
Two-sample t test power calculation # 二標本t検定
n = 17.84713 # 必要な両群の観測値数
delta = 1 sd = 1 sig.level = 0.05
power = 0.9
alternative = one.sided # 片側対立仮説を指定
NOTE: n is number in *each* group
> power.t.test(n=50, delta=1, sd=1, sig.level=0.05, power=NULL, type="one.sample", alternative="one.sided") One-sample t test power calculation # 一標本t検定
n = 50 delta = 1
sd = 1 sig.level = 0.05
power = 1 # そのときの検出力
alternative = one.sided # 片側対立仮説を指定
準
1 − α
，帰無仮説p = q
，両側対立仮説p 6 = q
で検定する例である．モンテカルロ法 では，繰り返し数が異なれば異なった結果(
使った乱数系列にも依存)
を得るので，繰り 返し数を変えて何度も確かめてみる方が賢明である．また，計算時間は当然余計にかかる が，あれこれ文献を探すよりは早い．> foo <- function(N, n, p, m, q, alpha) {
X <- rbinom(N,n,p) # B(n,p)に従う標本N個
Y <- rbinom(N,m,q) # B(m,q)に従う標本N個
phat <- (X+Y)/(n+m) # 帰無仮説p=q下での共通比率の推定値
S <- sqrt((1/n+1/m)*phat*(1-phat)) # その時のX/n-Y/mの標準偏差推定値
z <- qnorm(1-alpha/2) # 両側検定棄却限界値
mean(abs(X/n - Y/m)/S > z) # N回の検定で棄却される標本比率 }
> bar <- function(n) { # m=0.2*nのケース．N=1e7で実行 x <- foo(1e7, n, p, ceiling(0.2*n), q, alpha) cat(n, ceiling(0.2*n), x, fill=TRUE) }
> p=0.3; q=0.5; alpha=0.05; beta <- 0.8
> set.seed(1) # 再現性のために乱数種を指定
> sapply(260:300, bar) # n=260:300で実行
260 52 0.7875916 261 53 0.7932245 262 53 0.7959385 263 53 0.7920496 264 53 0.7943363 265 53 0.7971905
266 54 0.8025416 # 標本検出力が初めて0.8を越える標本数
267 54 0.799305 268 54 0.8011001 (以下省略)
> bar <- function(n) { # 繰り返し数1e6で再実行
x <- foo(1e6, n, p, ceiling(0.2*n), q, alpha) cat(n, ceiling(0.2*n), x, fill=TRUE) }
> set.seed(1)
> sapply(260:300, bar) 260 52 0.787136 261 53 0.793089 262 53 0.795644 263 53 0.792638 264 53 0.794402 265 53 0.797296
266 54 0.802223 # <- 標本検出力が初めて0.8を越える標本数
267 54 0.798358 268 54 0.801 (以下省略)
参考：パッケージ
Hmisc
には成功確率とサイズが異なる二組の二項分布データの平均の 差の検定の検出力を計算する関数bpower()
がある．これによれば必要サイズは266
と54
である．> library(Hmisc) # パッケージHmiscを読み込む
> baz <- function(N) { # 標本検出力が初めて80%を越える標本数
power <- bpower(p1=P1, p2=P2,
n1=N, n2=ceiling(Nratio*N), alpha=Alpha)
if (power >= Power) # 要求パワーを越えた時点でストップ
{cat(N, ceiling(Nratio*N), power, fill=TRUE); stop()}
}
> P1 <- 0.3; P2 <- 0.5; Nratio <- 0.2; Alpha <- 0.05; Power <- 0.8
> sapply(100:10000, baz)) # Nを変えながら検出力が0.8を越えるまで繰り返し
266 54 0.8003697 # ループ中断によるエラーメッセージ．無視可能
第 4 章
時系列解析
R
の基本パッケージ(base
パッケージおよび付属パッケージstats)
には時系列デー タを処理し，解析する基本的な関数が用意されている．更にCRAN
にはより高度な時系 列解析用のアドオンパッケージがいくつも*1ある．時間の経過の中で観測されるデータは，データの値
x
i とともに，しばしば観測された時 間t
i の情報が重要なことがある．したがってデータは，観測時間{ t
i}
と対応する観測値{ x
i}
という対で表現される．こうしたデータの中で，特に{ t
i}
が秒・分・時間・日・週・月・四半期・年等の等間隔の間隔であるものを
(
離散)
時系列データ(discrete time series data)
と呼ぶ．時系列データを解析する統計理論を時系列解析(time series analysis)
と 総称し，その興味の中心は異なった時間での観測値間に想定される従属性構造のモデル化 にある．極端な場合として{ x
i}
が互いに独立な確率変数列の実現値とみなすことができ れば，時間情報は無意味になる．時系列解析は，従って，本質的に従属性を持つ確率変数 列の理論となり，「独立・同分布データの解析」という，統計学の基本的枠組から外れ，そ の理論・解析は一層困難となり，R
等の高機能の統計解析ソフトの重要性が一層ます．以下では，
base
およびstats
パッケージ中の時系列解析用の関数を解説する．R
の 時系列解析用の関数は，時系列オブジェクト(time-series object)
という，時間と対応す る観測値の対，および観測期間・周期等の情報を持つ特殊なオブジェクトを扱う．こうし たオブジェクトは"ts"
というクラスとされ，オブジェクト自身のクラス属性として記録 されている．時系列解析用の関数は，与えられた引数がクラス"ts"
であることをまず確 認し，それに基づいて処理を行うようになっている．時系列オブジェクトは
(
等間隔)
観測時間情報を持つため，時間に関する以下のような 幾つかの特殊な情報を持つ:
*1CRANにあるTask Viewと呼ばれる各種パッケージをテーマ別にまとめたもののうち，テーマFinance とEconometricsを参照．
自然な時間単位 年，月，週，一時間等 観測開始時間
start
観測終了時間
end
周波数
frequency
．単位時間内の観測時間の数．月別データなら，自然な時間単位「年」に対する周波数は
12
サンプリング比率
deltat
．自然な時間単位に対するサンプリング間隔を表す比率．月別データなら，自然な時間単位「年」に対する
deltat
は1/12
．周波数とサン プリング比率どちらか一方を与えれば良い周期
cycle
．各データの観測時間情報を表す自然な時間単位とその中での周波数を表す対．例えば年と月
(1993,5)
，年と第2
四半期(2002,Qtr2)
，月と日(12,13)
， 週と曜日(43,Fri)
等注意：時系列データに固有の特徴として，系統的にデータ値が欠損することがある．例え ば，株価等の経済時系列は，土日・祭日はデータが本来無い．
注意：