Hatena::ブログ(Diary)

アイアナ:データ分析や人工知能(AI)などの技術雑記

2009 | 11 |
2010 | 02 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2011 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2012 | 01 | 03 | 04 | 05 | 10 | 12 |
2013 | 01 | 02 | 04 |
2014 | 03 | 12 |
2016 | 03 |
2017 | 10 |

2011-06-30

小技メモ(ggplot2に直線を引く)

@triadsouさんが紹介していたページの記事です。

http://stackoverflow.com/questions/6520758/adding-line-with-points-to-a-plot-in-ggplot2

library(ggplot2)

addlinetoplot <- function(dataset, varx, vary) { 
  list(
    geom_line(data=dataset, aes_string(x=varx, y=vary)), 
    geom_point(data=dataset, aes_string(x=varx, y=vary))
  )
}

df1 <- data.frame(c1 = c(1:10), c2 = c(1:10))
c1 <- c(1:10)
csq <- c1^2
df2 <- data.frame(c1 = c(1:10), csq)

pltbase <- ggplot() + geom_line(data = df1, aes(x=c1, y=c2))
pltbase + addlinetoplot(df2, varx = "c1", vary = "csq")

f:id:isseing333:20110630150908j:image

2011-06-25

おしゃれStatistics I

統計学勉強会、「おしゃれStatistics」についての記事です。

Googleグループを作っておりますので、どなた様もお気軽にご登録下さい。

略称は「おしゃスタ」ですww


第1回は一ヶ月ほど前に開催しておりましたので、その時のレポートを致します。

開催日:2011/05/05 14:00〜16:00

場所:新宿歌舞伎町ルノアール

人数:7人

f:id:isseing333:20110625202600p:image

歌舞伎町ルノアールの掲示板。Statisticsの綴りが間違っています笑



内容:教科書『Statistics』の1〜2章


1. コントロールされた実験(ランダム化試験)

  • 新薬の効果は治療群と対照群の「比較」によって測定する
  • そもそも対照群が無い場合
  • ワクチンの効果と流行の効果を分離できない
  • 1916年から流行したポリオに対するソークワクチンの例

f:id:isseing333:20110625203547p:image

  1. 学年によって群分けしている
    • ポリオは接触感染するので学年2の発症率が学年1、3より高いといった可能性がある
    • その場合はワクチンに不利なバイアスとなる
  2. 同意のとれた子供だけ治療している

  • 門脈大静脈シャント、門脈体循環シャント
    • 門脈は肝臓に入って肝静脈を通って心臓へ血液を流す
    • 肝硬変などで肝臓のろ過が悪くなると肝臓→心臓の流れが悪くなる
    • 食道の血管圧が高くなり静脈瘤の危険性が出てくる
    • シャントによって門脈と肝静脈の間にバイパスを作って回避する
    • 手術時間が長く、危険も多い。
    • リスクベネフィットはどうなのか?
    • ランダム化試験と対照の無い研究の比較
      • ランダム化試験では「効果なし」
      • 対照の無い研究、対照はあってもランダム化してない試験では「効果あり」
    • ランダム化していない試験では選択バイアスが存在する
      • 症状の重い患者を対照としており、軽い患者に手術をしている
      • 手術に有利なバイアスとなっている
    • 「適格基準」を定めて同等の患者同士で比較をするべき

  • 既存対照、ヒストリカルコントロール
    • 対照群の性質が治療群と同等であるか?
    • 既存対照試験では効果があるが、ランダム化試験では効果なしという治療法も多い
    • 既存対照試験では交絡が見られる

  • サマリ
    1. 効果を測定するときは治療群と対照群の比較を行う
    2. 2群の治療以外の性質が似ている場合は、結果の差異は治療の効果である
    3. 治療以外の性質が異なっている場合は交絡の可能性がある
    4. 治療群と対照群を似た集団にするためにランダム化試験が行われる
    5. 治療されることへの意識の強さに影響されないようにするため対照群はプラセボ投与
    6. 対象者と評価者による先入観のバイアスを防ぐため、二重盲見を行う



2. 観察研究

  • サマリ
    1. 観察研究では研究者が対象者を割り付けることはできない
      • 喫煙研究の例では、ある対象者に喫煙を勧めることは不可能である
      • もともと喫煙している人を治療群、していない人を対照群とする
    2. 観察研究は「関連(assiciation)」を見つけることはできるが、それは必ずしも因果ではない
    3. 治療効果は治療群を決定する因子によって交絡していることが多く、しばしば結果は誤解を招く。
    4. 研究→対照群があるか→同時期か→対照群が制御されているか→ランダムか
    5. 観察研究の場合、2群が均一であり比較可能かを確認する
      • 交絡因子があるかどうか
    6. 交絡因子を調整する
      • 似たような集団に分けて比較する(マッチング、層別)
      • 傾向スコアという調整法もある
    7. 観察研究の弱点は未知の交絡要因であり、ランダム化試験はこれを最小限にできる

f:id:isseing333:20110625204302p:image

メモ:治療・対照群は介入A・B群、介入・非介入群などのように適宜読み替える


  • 広告の効果はどうやって測定する?
    • アマゾンはランダム化試験をやったが社会的に叩かれた
    • 観察研究のように、既にあるデータから考えられる交絡要因を調整した解析をするのが現実的なのではないか?
      • 一般的な交絡要因:性、年齢
      • 時間帯、曜日、趣味、嗜好、社会情勢など
    • そもそも調整して得られる効果は「全体的な効果」
      • 経営的には意味の無い値かも
      • 「どんな層にどれくらいウケているか」といった記述的な分析の方が必要とされているのかも?


以上が第1回目の内容です。第1回目は本の内容を1時間、質問で1時間といった感じでした。

本の内容をもうちょっと少なくして、意見交換の時間を増やしてもいいかなと思っています。

統計について質問できる場も少ないと思いますので、内容に関わらすいろいろ質問して下さい。

統計と社内システムをどうやって結びつけるか?」などのように漠然とした疑問でも結構です!


あと日本語の本だと、これも良いかなと思っています。

おしゃスタで使うかどうか検討中です。。


第2回目はGoogleグループでも告知しますが7/7、19:00〜21:00、新大久保で行う予定です。

(あまり広い部屋を取ってないのですぐ埋まってしまうかも、、、)

これから月一で開催できればなと思っています!

2011-06-16

時系列データの解析(厚労省公開の医療費データ)

この本に沿って時系列データの解析方法をまとめました。


Rによる時系列分析入門

Rによる時系列分析入門



サンプルデータを使っても面白くないので、厚労省が公開している医療費のデータを使いました。

厚労省の医療費データベース


例によってこのデータはエクセルで公開されていて、そのまま解析できる状態じゃありません

今回は入院の総医療費だけを扱ったので、その部分だけ加工してcsvにしました。

一応、加工したデータはダウンロードのページに置いてます。



それでは、解析していきます。


まずはデータ読み込みと加工

Iryouhi <- read.csv("医療費.csv", as.is = T)
Nyuin <- ts(Iryouhi[, 2], frequency = 12, start = c(1984, 4))

read.csv関数csvファイルを読み込み、ts関数で時系列データの形に変換します。

時系列プロットにしてイメージを掴みます。

plot(Nyuin, yaxt = "n", ylab = "入院医療費(1ヵ月毎)", xlab = "年月", main = "入院医療費の経時変化")
axis(2, c(5*10^11, 7.5*10^11, 10^12), c("5000億", "7500億", "1兆"))

f:id:isseing333:20110616122819j:image


2010年では一ヵ月に約1兆円かかっています。

単純計算で年間約12兆、入院は全体の3分の1だと考えると全体は36兆。

大体計算は合いそうですね。


ギザギザなっているのは周期的な季節変動。

全体的に医療費が上がっているのは分かるのですが、季節変動を取り除くとどうなるかやってみます。


NyuinTS <- decompose(Nyuin)
plot(NyuinTS)

f:id:isseing333:20110616122820j:image


decompose()関数で勝手に季節変動を取り除いてくれます。

上から元のデータ、季節変動を除いたトレンド、季節変動、ランダムな成分となっています。


トレンドだけを抽出するとこのようになります。

plot(NyuinTS$trend)

f:id:isseing333:20110616122821j:image


移動平均を行っても、同じような結果になります。

NyuinMA12 <- filter(Nyuin, rep(1, 12))
plot(NyuinMA12)


このトレンドは一見するとただ単に毎年上がっているように見えるのですが、2000年までの医療費の上がり方とそれ以降の上がり方が違うように見えます。その事をさらに確認するべく、自己相関プロットを描いてみます。

acf(Nyuin)

f:id:isseing333:20110616130031j:image


右肩上がりで上昇しているトレンドを持つ時系列データでは、自己相関プロットこのようにLag(間隔)が大きくなるにつれて下がっていくような図になります。

入院医療費が全体的に上昇するトレンドを持っているので、このような図が得られたのです。


これを、入院医療費の期間を区切って描いてみます。

par(mfrow=c(3, 2))
acf(window(Nyuin, start = c(1985, 4), end = c(1990, 3)), main = "1985-1990")
acf(window(Nyuin, start = c(1990, 4), end = c(1995, 3)), main = "1990-1995")
acf(window(Nyuin, start = c(1995, 4), end = c(2000, 3)), main = "1995-2000")
acf(window(Nyuin, start = c(2000, 4), end = c(2005, 3)), main = "2000-2005")
acf(window(Nyuin, start = c(2005, 4), end = c(2010, 3)), main = "2005-2010")
par(mfrow=c(1, 1))

f:id:isseing333:20110616143620j:image


時系列を5年ごとに区切って自己相関プロットを描いています。

結果を見ると明らかに1985年から2000年までの3つのプロットと2000年から2010年までの2つのプロットの傾向が違うことが分かります。


2000年から2005年の間は自己相関が低く、この期間の時系列に上昇傾向のトレンドが弱いことが分かります。この時期を思い返してみると、小泉元首相構造改革の時期とほぼ重なっているように思えます。


構造改革小泉さんは「医療費を抑制する」ことを打ち出していましたが、こんなところでそれがデータに現れているのではないでしょうか?


また、2005年から2010年の自己相関もあまり高くありません。

1985年から2000年までの医療費上昇率より、ここ5年間の医療費上昇率の方が低いことが伺えます。


というより、過去の上昇率が異常に高かったのでは?

さらに議論をするためには、多方面からの分析が必要でしょう。



さて、時系列分析に戻ります。

時系列データにモデルを当てはめ、将来の予測を行うこともできます。

NyuinSARIMA <- arima(Nyuin, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 1), period = 12))
nseason         <- 120
NyuinPredSARIMA <- predict(NyuinSARIMA, n.ahead = nseason)
ts.plot(Nyuin, NyuinPredSARIMA$pred, col = c(1, "blue"), main = "SARIMA予測", 
	ylab = "入院医療費(1ヵ月毎)", xlab = "年月")
text(1990, c(5*10^11, 7.5*10^11, 10^12, 1.25*10^12), c("5000億", "7500億", "1兆", "1.25兆"))

f:id:isseing333:20110616150009j:image


これは季節ARIMA(SARIMA)モデルで将来の予測をしたものです。

ARIMAモデルはAR(自己相関), 差分の回数, MA(移動平均)の3つをモデル化したもので、それぞれの次数によって予測結果は変わります。ですので、今回の予測結果が全てではありませんが、どちらにせよこのまま医療費が上がっていくことは予想されるでしょう。

最後に、その他の時系列解析のコードを載せておきます。

#------ホルト・ウィンタース法(汎用的な時系列モデル)
NyuinHW1 <- HoltWinters(Nyuin)
NyuinHW2 <- HoltWinters(Nyuin, seasonal = "mult")
NyuinHW3 <- HoltWinters(Nyuin, gamma = 0, beta = 0)
plot(NyuinHW1)
plot(NyuinHW2)
plot(NyuinHW3)

nseason      <- 1200
NyuinPred1    <- predict(NyuinHW1, n.ahead = nseason)
NyuinPred2    <- predict(NyuinHW2, n.ahead = nseason)
NyuinPred3    <- predict(NyuinHW3, n.ahead = nseason)
ts.plot(Nyuin, NyuinPred1, col = c(1, "blue"), main = "HW1予測")
ts.plot(Nyuin, NyuinPred2, col = c(1, "blue"), main = "HW2予測")
ts.plot(Nyuin, NyuinPred3, col = c(1, "blue"), main = "HW3予測")


#------ARIMA、SARIMA(→信頼区間は?)

#---order:AR(自己相関), 差分の回数, MA(移動平均)
#---予測のプロットを見ると、このARIMAはちょっとおかしい

NyuinARIMA  <- arima(Nyuin, order = c(1, 1, 1))
NyuinSARIMA <- arima(Nyuin, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 1), period = 12))

plot(Nyuin - NyuinARIMA$residuals, main = "ARIMA")
plot(Nyuin - NyuinSARIMA$residuals, main = "SARIMA")

#---残差の検定(Ljung-Box検定)→有意でなければホワイト・ノイズ
plot(NyuinARIMA$residuals, main = "残差(ARIMA)")
Box.test(NyuinARIMA$residuals)
Box.test(NyuinSARIMA$residuals)

#---lagを変化させて、「s次まで自己相関が無い」ことをチェックする
tsdiag(NyuinARIMA)
tsdiag(NyuinSARIMA)


nseason         <- 120
NyuinPredARIMA  <- predict(NyuinARIMA,  n.ahead = nseason)
NyuinPredSARIMA <- predict(NyuinSARIMA, n.ahead = nseason)
ts.plot(Nyuin, NyuinPredARIMA$pred, col = c(1, "blue"), main = "ARIMA予測")
ts.plot(Nyuin, NyuinPredSARIMA$pred, col = c(1, "green"), main = "SARIMA予測", 
	ylab = "入院医療費(1ヵ月毎)", xlab = "年月")
text(1990, c(5*10^11, 7.5*10^11, 10^12, 1.25*10^12), c("5000億", "7500億", "1兆", "1.25兆"))



#------単位根検定(修正ディッキー・フラー検定、ADF検定)
#---Call: lm()とあるように、線形回帰モデルを利用している
library(urca)
summary(ur.df(Nyuin))
summary(ur.df(Nyuin, type = "trend"))
summary(ur.df(Nyuin, type = "trend", lag = 12))
summary(ur.df(Nyuin, type = "drift", lag = 12))




#------スペクトル
#---ペリオドグラム(高速フーリエ変換、spans=で修正Daniell平滑化を行う)
spectrum(Nyuin, method = "pgram")


#---ARモデル
#---order=でARの次数を指定できる、指定しなければAICが最小の次数を選ぶ
spectrum(Nyuin, method = "ar")


#---フィルター後のスペクトル
NyuinFilter <- filter(Nyuin, rep(1, 12))
NyuinFilter <- NyuinFilter[!is.na(NyuinFilter)]
spectrum(NyuinFilter, method = "ar")


#---※スルツキー効果:差分をとることで見せかけの周期性が見られる現象

2011-06-10

小技メモ(hwriterでcssjavascriptを使う)

hwriterはRを使って簡単にHTMLレポートを作るのに適したパッケージです。

cssjavascriptにも対応しています。


htmlファイルを作るディレクトリに.cssファイル、.javascriptを保存しておく。


  • 利用するjavascriptcssを指定(日本語を扱う場合は「charset = "shift-jis"」を指定しておく)。
p <- openPage(".html", charset = "shift-jis", link.javascript = , link.css = )

  • cssクラスの指定(row.bgcolorは先頭行の色を指定)。

テーブルオブジェクトはそのままHTMLにできないのでデータフレームにしておきます。

hwrite(テーブル, p, table.class = "cssクラス名", row.bgcolor = '#ffdc98', br = T)
hwrite(テーブルじゃないオブジェクト, p, class = "cssクラス名", br = T)

javascriptはopenPageで読み込むだけで有効になるはずです。

最後に次のコードでクローズ。


closePage(p)

お試しあれ。

CRAN Task View: Clinical Trial Design, Monitoring, and Analysis (のまとめと補足)

タスクビューの臨床試験部分を半日本語訳しましたので、記事にしておきます。

ソースはこちら→http://cran.r-project.org/web/views/ClinicalTrials.html


タスクビュー:臨床試験デザイン、モニタリング、解析

メンテナンス者: Ed Zhang

連絡先: Ed.Zhang.jr at gmail.com

バージョン: 2010-07-07


デザインとモニタリング

blockrand

  • ブロックランダム化試験の割付を行う
  • PDFのランダム化カードを作る
male <- blockrand(n=100, id.prefix='M', block.prefix='M',stratum='Male')
female <- blockrand(n=100, id.prefix='F', block.prefix='F',stratum='Female')
my.study <- rbind(male,female)
plot.blockrand(my.study,'mystudy.pdf',
	top=list(text=c('My Study','Patient: %ID%','Treatment: %TREAT%'),
	col=c('black','black','red'),font=c(1,1,4)),
	middle=list(text=c("My Study","Sex: %STRAT%","Patient: %ID%"),
	col=c('black','blue','green'),font=c(1,2,3)),
	bottom="Call 123-4567 to report patient entry",
	cut.marks=TRUE)

experiment

  • 前回参照
  • ブロックランダム化試験の割付を行う
    • randomize(, block= )

GroupSeqgsDesign

  • 群逐次デザイン(wald先生が教科書)
    • メモ:ワルド・スコア・尤度比検定(3種の検定があるけど3つともあんまり結果は変わらない)
  • α消費関数に基づく
    • αエラーを「消費する」という概念を利用(αエラーの概念は多重検定とかだけに応用されているわけじゃない)
    • メモ:適応デザイン、adaptive design

Hmisc ldBand()関数ldbounds

  • 群逐次デザインのLan-DeMets中止境界を描く(ld98プログラムを利用、Department of Biostatistics, University of Wisconsin written)

PwrGSD

  • 群逐次デザインの検出力を計算

seqmon

  • 不明
  • Armitage-McPherson and Rowe Algorithm(Schoenfeld (2001))

デザインと解析

clinfun

  • フェーズII試験、最適・ミニマックス2段階フェーズII試験(サイモンデザイン)、毒性のモニタリングのための繰り返し有意性検定、群逐次フェーズIII試験

DoseFinding

  • 用量発見試験(薬物動態試験は?)
  • 多重対比試験、非線形用量反応モデル、MCP(multiple comparison procedures)デザインへの準備
  • 正規性で等分散なエンドポイントに適用できる

MCPMod

  • 多重比較手順(multiple comparison procedures)を組み合わせた用量発見試験
    • メモ:ゲートキーピングのパッケージはあるかな?

特殊なデザインの解析

bifactrial

  • 2要因、3要因臨床試験
  • ブートストラップ法を使う解析もある

ClinicalRobustPriors

MChtest

  • モンテカルロ仮説検定
    • メモ:試験の設定が複雑な場合は、数式でサンプルサイズを計算するのではなくてシミュレーションで出すこともある(製薬業界)
  • 2つの異なる逐次境界を比較(?)、Besag and Clifford (1991)
  • P値の信頼区間も計算

speff2trial

  • 2群ランダム化試験
  • IPW(inverse probability weighting)
  • 質的、2値エンドポイント

一般的な解析

Base

chisq.test, prop.test, binom.test, t.test, wilcox.test, kruskal.test, mcnemar.test, cor.test, power.t.test, power.prop.test, power.anova.test, lm, glm, nls, anova

※基本的なサンプルサイズ設計はBaseのpower.t.test, power.prop.test, power.anova.testなどで十分


asypow

  • 漸近尤度比を使った検出力の計算

binomSamSize

  • 2項分布でのサンプルサイズ

coin

  • 条件付推量
    • ランダム化をやっていてもランダムサンプリングの仮定がフィットしない時など
    • 並べ替え検定かな?

epibasix

  • マッチング(matched pair)の解析
    • 関数:diffdetect, n4means, n4props

epicalc

  • サンプルサイズ設計
    • 関数:n.for.2means, n.for.2p, n.for.equi.2p, n.for.noninferior.2p, n.for.cluster.2p

HH

  • ae.dotplot()関数で有害事象の確認ができる

Hmisc

  • 200ほどの有用関数(osakaRのスライドシェア、林さん)
  • データ解析、高水準グラフ、サンプルサイズ設計、SASデータセットの変換、欠測の補完、高度な表の作成、クラスタリング文字列の複製、Texコードの作成、再コーディング、ブートストラップ繰り返し解析

multcomp

  • 線形モデルの同時検定、同時信頼区間(GLM、GLIM、GLMM、生存時間解析)

survival

ssanv


メタアナリシス

copas

meta

metafor

  • メタアナリシスの固定効果、変量効果
  • 2×2表のデータにはMH、Petoの方法が利用できる

meta

2011-06-01

ディリクレ過程混合モデル(The Dirichlet Process Mixture (DPM) Model)は難しいです。。。

今日は言語処理の分野で良く使われていると言われるディリクレ過程混合モデル(DPM)についてです。

記事にはしているものの、自分でもぶっちゃけ良く分かっていませんorz

あまり慣れない分布なのでイメージが付きにくいですね。。。


ディリクレ分布については多項分布の共役事前分布になっているという程度の知識しかありません。

ベータ分布とか出てこられても、、、ねぇ。



参考にさせて頂いたのは以下の2つの資料。

http://d.hatena.ne.jp/ryamada22/20050528/1117267472

http://biocomp.bioen.uiuc.edu/journal_club_web/dirichlet.pdf


下のURLにあるPDF資料に沿ってホワイトボードで解説したのが次の画像です。(我ながら汚いですが、、)


f:id:isseing333:20110601193602j:image


f:id:isseing333:20110601193603j:image



PDFの方の資料に沿って説明しているのですが、ポリアの壺のところで断念。

とりあえずRに関数はあるので分析&クラスタリングはできるのですが、内部的になにをやっているのがイメージを付けにくいですね。


あとはRでの乱数発生。

#---ガンマ分布
#---ある制限時間内に2匹、3匹、5匹釣ることが期待されている状況
n <- 100
mean(rgamma(n,2))
mean(rgamma(n,3))
mean(rgamma(n,5))

#---ディリクレ分布
library(MCMCpack)
rdirichlet(n,c(2,3,5))

そのうち分かるようになるかなと思いつつ、とりあえず記録しておきます。

CRAN Task View: Design of Experiments & Analysis of Experimental Data (の日本語訳)

twitterで#RtransというRのドキュメント翻訳プロジェクトに関わっているのですが、Task Viewの方に浮気してしまいました。。。

いやでもこっちのドキュメントも重要だよね?と自分を納得させつつ、とりあえず訳します。



実験計画法の部分の訳です(ソースは→http://cran.r-project.org/web/views/ExperimentalDesign.html

どこかwikiでRドキュメントの訳をまとまているところがあったら転載して頂いて構いません。

(@yokkunsが運営していたような。。。

→無事転載されました!http://rwiki.tkul.jp/index.php?CRANTaskView



以下、和訳です。

英語苦手なんで誤訳があったらすみません。

また、分野によって若干訳し方が違うかもしれませんがご了承下さい。





タスクビュー:実験計画(Design of Experiments, DoE)と実験データの解析


メンテナンス者:Ulrike Groemping

連絡先:groemping at bht-berlin.de

バージョン:2010-09-06


このタスクビューは実験計画と実験から得られたデータの解析に関するRパッケージの情報を集めています。改良に関する示唆はお気軽にご連絡下さい。またこの分野に属すと思われる新しいパッケージや大きなパッケージのアップデートがあればご連絡下さい。このwebに詳細が書いてあります。http://prof.beuth-hochschule.de/groemping/DoE


実験計画は様々な分野で適用されており、いろいろなフィールドの需要のために手法が調整されている。このタスクビューは最も一般的なパッケージのセクションから始まり、農学・工業のような専門的なセクション、専門的な実験計画の目的のために作られたパッケージに関するセクションへと続いていく。もちろん、分野への分割はいつも明確なわけではないが、より専門的なセクションは一般的な文脈で適用することができる。



一般的な目的のための実験計画

一般的な目的の実験計画のために作られたパッケージはほとんどない。まず最初に、baseパッケージにある一般化線形モデルは、実験が計画されたデータを解析するためにとても重要である(特にlm()関数、aov()関数、線形モデルオブジェクトから作られるメソッド関数)。これらはKuhnert and Venables (2005, p. 109 ff)に詳しく書いてあるが、Vikneswaran (2005)は実験計画のための特殊な利用方法を指摘した(contrasts()関数、多重比較関数、model.tables()、replications()、plot.design()のなどような便利な関数を使って)。granovaは単純な構造の実験(一元配置、二元配置、対応のあるデータ)を少し変わったグラフィックで表示する。AlgDesignは量的変数がある場合と無い場合の完全要因実験、混合デザイン(つまりデザインは要因を足すと1=100%となる)、正確または近似的なD-、A-、I-最適デザインを作る。注意:AlgDesignの作者・メンテナンス者であるBob Wheelerはこの仕事から退きたがっており、引き継いでくれる人を探している。興味のある方はBobに連絡して下さい。conf.designはブロックに交互作用のあるデザインを作ることが出来(conf.design関数)、既存のデザインを組み合わせることが出来る(例えば、工業実験における田口の配列内、外デザインに有用である)。daeは実験計画やR因子に関するいくつかの有用な関数がある。例えば、いくつかの入れ子構造を扱うことのできる実験方法(Bailey 1981)、いくつかの因子を組み合わせて1つにしたり1つの因子をいくつかに分けたりする関数がある。さらにaov()関数によって返される処理された後のオブジェクトの機能を提供することができる。例えば、2段階実験のイェーツの効果を抽出できる。blockToolsはブロックのサイズが小さい場合に、均一なブロックで終わるように割り付けを行う。



農学・植物育種実験のための実験計画

agricolae農学や植物育種のために特化した実験計画のための豊富な機能を提供する。この実験計画は他の目的にも有用である。これは格子デザイン、要因デザイン、乱塊デザイン(ブロックランダム化)、完全ランダム化デザイン、(グレコ・)ラテン方格デザイン、不完全ブロックデザイン、αデザインの計画を行う。例えば治療比較の目的やノンパラメトリック検定などの実験データの解析のための機能もあるが、特殊なタイプの実験のための限られた可能性もある。



工業(産業)実験のための実験計画

いくつかのパッケージは工業実験を扱う。工業実験は意図的な交絡を非常に分別し、誤差のための余分な自由度を持つ。

DoE.baseFrF2DoE.wrapperなどの一連のパッケージは主要な文法やアウトプットの構造が同じであり、中心的な特徴は実験計画のオブジェクトのためのdesignクラスである。4つめのRcmdrPlugin.DoE(βバージョン)パッケージはこれら3つのパッケージの機能を、コマンドラインプログラミングを使えない・使いたくないユーザーのためにRコマンダー(Rcmdrパッケージ、Fox 2005)に統合している。

DoE.baseは完全要因実験や主効果の実験のための直行配列を作る(Kuhfeld 2009で144もの実行がリストされ、???)。このパッケージは他の2つのためのインフラを作る。それは実験計画データフレームのためのdesignクラスや、このメソッドのためのいくつかのクラスである。直行配列の質を評価するためのいくつかの実験機能もある。FrF2はプラケット・バーマン(Plackett-Burman)型スクリーニングのような2段階の部分要因実験を行う。部分要因実験は最高分解能最小収差(maximum resolution minimum aberration)デザインがデフォルトであり、デザインの併合分類(incorporated catalogue)によって様々な方法でカスタマイズできる(Chan, Sun and Wu 1993によって分類されたデザイン、Block and Mee 2005、とXu 2009に分類されたさらに大きなデザインを含む。加えてFrF2.catlg128はとても大きい完全分類(分解能IVの、特殊な目的のための23因子128実行デザイン)を作ることができる)。解析の視点からは、FrF2は簡単なグラフィカル解析ツールである(正規、半正規効果プロットBsMDで改良された)、主効果プロットMinitabソフトのような交互作用プロット行列、3因子の組み合わせの立方プロット)。また2段階因子の部分要因のための別の構造を示すことができる。DoE.wrapperrsmlhsDiceDesignAlgDesignの機能を他の2つのパッケージの文脈に継承させ、便利な機能や他の可能性を付け加えることができる。特にこのパッケージは自動的にCCDs(central composite designs、中央複合デザイン)のキューブ部分を選ぶ、または既にあるキューブをCCDに付け足すことができる。

BHH2はBox, Hunter and Hunter第2版に書いてあり、いくつかのデータセットがある。おおくの因子と関係の定義リストから2段階の完全・部分要因実験を作る(ffDesMatrix()関数FrF2パッケージより快適ではない)。2段階要因実験データを解析するためのいくつかの関数がある。anovaPlot()関数は残差と比較した効果の大きさを評価し、lambdaPlot()はBox-Cox変換した効果の統計的有意差を評価する。BsMDはBox and Meyer(1986)が提案したベイジアンチャートや、2段階部分要因実験の中で効果が発現しているかを評価する効果プロット(正規、半正規、レンス(Lenth))を作成する。要因実験の計画や解析を行うツールから離れて、Rはrsmパッケージによって曲面最適化surface optimization)を行う。このパッケージは1次・2次反応曲面モデルの逐次最適化を行い、急上昇するように近づく最適化と線形モデルの応答関数の視覚化を提供している。また反応曲面の探索コーディングが促進されている。



コンピュータ実験のための実験計画

質的因子を使ったコンピュータ実験は特殊なタイプの実験デザインを必要とする。多くの異なる段階の因子が含まれることがよくあり、繰り返しが有益ではない。さらに調査している現象を適切に表現するために、線形や2次モデルを仮定するには実験の範囲がとても広い。その結果、実験空間を可能な限り埋めること(空間充填デザイン、space-filling designs)が理想的である。例えいくつかの要因が無関係であると判明しても、それぞれの実行(run)は追加の情報を生み出す。lhsパッケージはこの目的のためにラテン超立方デザイン(latin hypercube designs)を作る。さらに、実行しているフォローアップ実験を強調したコンピュータ実験の解析方法も提供する。同じような姿勢のもう1つのDiceDesignパッケージは、さらに空間充填デザインを実行する方法やコンピュータ実験の質の指標を追加する。DiceKrigingコンピュータ実験からメタモデルを生成するクリック手法を提供する。



その他の目的の実験計画

さらにいくつかのパッケージは特殊な状況での実験計画を扱う。desirailityは多基準解析を単純にするために、いくつかのターゲットの基準を望ましい関数に組み合わせる。experimentは臨床実験のツールを含む。例えばランダム化ツールや臨床試験のためのいくつかの特別な解析である。DoseFindingは用量発見実験のデザインと解析を行う関数を提供する(例えば薬剤のフェーズII臨床試験)。これはMCPModパッケージ(メンテナンスはされていない、Bornkamp, Pinheiro and Brets 2009で紹介されている)に、用量発見デザインのための特殊な最適化デザイン(MED最適デザイン、D-最適デザイン、両方の混合、Detteら 2008)の便利さを組み合わせている。idDesign連鎖平衡デザインを、crossdesは計量センサー(sensometrics)に利用できるいくつかのタイプのクロスオーバーデザインを(ラテン方格、多直交ラテン方格、ヨーデン方格も含む)利用できる。SensoMineRは例えば三角テスト(triangle test)のような計量センサー研究のための特殊なデザインを含んでいる。





  • Atkinson, A.C. and Donev, A.N. (1992). Optimum Experimental Designs . Oxford: Clarendon Press.
  • Bailey, R.A. (1981). A unified approach to design of experiments. Journal of the Royal Statistical Society, Series A 144 , 214-223.
  • Ball, R.D. (2005). Experimental Designs for Reliable Detection of Linkage Disequilibrium in Unstructured Random -Population Association Studies. Genetics 170 , 859-873.
  • Block, R. and Mee, R. (2005). Resolution IV Designs with 128 Runs. Journal of Quality Technology 37 , 282-293.
  • Bornkamp B., Pinheiro J. C., and Bretz, F. (2009). MCPMod: An R Package for the Design and Analysis of Dose-Finding Studies . Journal of Statistical Software 29 (7), 1-23.
  • Box G. E. P, Hunter, W. C. and Hunter, J. S. (2005). Statistics for Experimenters (2nd edition). New York: Wiley.
  • Box, G. E. P and R. D. Meyer (1986). An Analysis for Unreplicated Fractional Factorials. Technometrics 28 , 11-18.
  • Box, G. E. P and R. D. Meyer (1993). Finding the Active Factors in Fractionated Screening Experiments. Journal of Quality Technology 25 , 94-105.
  • Chen, J., Sun, D.X. and Wu, C.F.J. (1993). A catalogue of 2-level and 3-level orthogonal arrays. International Statistical Review 61 , 131-145.
  • Collings, B. J. (1989). Quick Confounding. Technometrics 31 , 107-110.
  • Daniel, C. (1959). Use of Half Normal Plots in Interpreting Two Level Experiments. Technometrics 1 , 311-340.
  • Derringer, G. and Suich, R. (1980). Simultaneous Optimization of Several Response Variables. Journal of Quality Technology 12 , 214-219.
  • Dette, H., Bretz, F., Pepelyshev, A. and Pinheiro, J. C. (2008). Optimal Designs for Dose Finding Studies. Journal of the American Statisical Association 103 , 1225-1237.
  • Federov, V.V. (1972). Theory of Optimal Experiments . Academic Press, New York.
  • Fox, J. (2005). The R Commander: A Basic-Statistics Graphical User Interface to R . Journal of Statistical Software 14 (9), 1-42.
  • Groemping, U. (2009). Design of Experiments in R . Presentation at UseR! 2009 in Rennes, France.
  • Hoaglin D., Mosteller F. and Tukey J. (eds., 1991). Fundamentals of Exploratory Analysis of Variance . Wiley, New York.
  • Jones, B. and Kenward, M.G. (1989). Design and Analysis of Cross-Over Trials . Chapman and Hall, London.
  • Johnson, M.E., Moore L.M. and Ylvisaker D. (1990). Minimax and maximin distance designs. Journal of Statistical Planning and Inference , 26 , 131-148.
  • Kuhfeld, W. (2009). Orthogonal arrays. Website courtesy of SAS Institute Inc., accessed August 4th 2010. URL http://support.sas.com/techsup/technote/ts723.html .
  • Kuhnert, P. and Venables, B. (2005) An Introduction to R: Software for Statistical Modelling & Computing . URL http://CRAN.R-project.org/doc/contrib/Kuhnert+Venables-R_Course_Notes.zip . (PDF document (about 360 pages) of lecture notes in combination with the data sets and R scripts)
  • Kunert, J. (1998). Sensory Experiments as Crossover Studies. Food Quality and Preference 9 , 243-253.
  • Lenth, R.V. (1989). Quick and Easy Analysis of Unreplicated Factorials. Technometrics 31 , 469-473.
  • Lenth, R.V. (2009). Response-Surface Methods in R, Using rsm . Journal of Statistical Software 32 (7), 1-17.
  • Mee, R. (2009). A Comprehensive Guide to Factorial Two-Level Experimentation . New York: Springer.
  • Myers, R. H. and Montgomery, C. M. (1995). Response Surface Methodology: Process and Product Optimization Using Designed Experiments . New York: Wiley.
  • Plackett, R.L. and Burman, J.P. (1946). The design of optimum multifactorial experiments. Biometrika 33 , 305-325.
  • Rosenbaum, P. (1989). Exploratory Plots for Paired Data. The American Statistician 43 , 108-109.
  • Sacks, J., Welch, W.J., Mitchell, T.J. and Wynn, H.P. (1989). Design and analysis of computer experiments. Statistical Science 4 , 409-435.
  • Santner T.J., Williams B.J. and Notz W.I. (2003). The Design and Analysis of Computer Experiments . Springer, New York.
  • Sen S, Satagopan JM and Churchill GA (2005). Quantitative Trait Locus Study Design from an Information Perspective. Genetics 170 , 447-464.
  • Stein, M. (1987). Large Sample Properties of Simulations Using Latin Hypercube Sampling. Technometrics 29 , 143-151.
  • Stocki, R. (2005). A Method to Improve Design Reliability Using Optimal Latin Hypercube Sampling. Computer Assisted Mechanics and Engineering Sciences 12 , 87-105.
  • Vikneswaran (2005). An R companion to "Experimental Design". URL http://CRAN.R-project.org/doc/contrib/Vikneswaran-ED_companion.pdf . (The file accompanies the book "Experimental Design with --Applications in Management, Engineering and the Sciences" by Berger and Maurer, 2002.)
  • Xu, H. (2009). Algorithmic Construction of Efficient Fractional Factorial Designs With Large Run Sizes. Technometrics 51 , 262-277.