<< 目次へ
最終更新日: 2009.9.8
経時測定データ
事前・事後データの解析法では2時点において繰り返しデータがとられたが、薬物試験などでは3日後、1週間後、3週間後といったように2時点以上にわたって繰り返しデータがとられることがある。そのような場合、ある時点(例えば3日後から3週間後の変化)に関心があるならば、同じように差をとって1標本データとして扱うことができる。しかしながら、そのようなことは他の時点において観測されたデータを切り捨てる(無視する)ことになるので好ましい方法とはいえない。そこで、ここではいわゆる対応ありの分散分析を紹介する。
3日後 | 1週間後 | 3週間後 |
2 2 3 2 1 4 |
2 4 2 3 3 5 |
6 7 8 6 6 10 |
分散分析についての詳しい解説は分散分析の基本原理を参照のこと。ここでは分散分析表を作成するための方法のみを紹介する。
まず時点要因による平方和SST、個体要因による平方和SSI、誤差による平方和SSE、全平方和SSYは次式によって表される。
ここに、
である。上から順に全平均、個体iの平均(行平均)、時点jの平均(列平均)となっている。
時点要因についての検定統計量は:
によって計算され、このF値は自由度p - 1, (p - 1)(n - 1)のF分布に従う。この検定統計量を用いて、以下の帰無仮説:
について検定を行うことになる。
続いて個体要因についての検定統計量は:
によって計算され、このF値は自由度p - 1, (p - 1)(n - 1)のF分布に従う。この検定統計量を用いて、以下の帰無仮説:
について検定を行う。
以下に分散分析表を載せる。
SS | DF | V | F | |
T | SST | p - 1 | VT = SST / (p - 1) | VT / VE |
I | SSI | n - 1 | VI = SSI / (n - 1) | VI / VE |
E | SSE | (p - 1)(n - 1) | VE = SSE / {(n - 1)(p - 1)} | |
Y | SSY | np - 1 |
■ ソフトウェアによる実行
関数paired.anova()を定義する。
paired.anova <- function(dat){ # 列に時点要因、行に個体要因となる行列を指定する n <- nrow(dat) # 行数の取得 p <- ncol(dat) # 列数の取得 ydd <- mean(dat) # 全平均 yid <- apply(dat, 1, mean) # 行平均 ydj <- apply(dat, 2, mean) # 列平均 SSY <- sum((dat - ydd)^2) # 全平方和 SST <- n * sum((ydj - ydd)^2) # 時点の平方和 SSI <- p * sum((yid - ydd)^2) # 個体の平方和 SSE <- SSY - (SST + SSI) # 誤差平方和 DFT <- p - 1 # 時点の自由度 DFI <- n - 1 # 個体の自由度 DFE <- (p - 1)*(n - 1) # 誤差の自由度 DFY <- n*p - 1 # 全自由度 VT <- SST / (p - 1) # 時点の分散 VI <- SSI / (n - 1) # 個体の分散 VE <- SSE / ((n-1)*(p-1)) # 誤差の分散 FT <- VT / VE # 時点に対するF値 FI <- VI / VE # 個体に対するF値 pT <- 1 - pf(FT, DFT, DFE) # 時点のp値 pI <- 1 - pf(FI, DFI, DFE) # 個体のp値 # 分散分析表 result <- matrix(c( SST, DFT, VT, FT, pT, SSI, DFI, VI, FI, pI, SSE, DFE, VE, NA, NA, SSY, DFY, NA, NA, NA), ncol=5, byrow=T) # 分散分析表の行ラベルと列ラベルを変更 rownames(result) <- c("T", "I", "E", "Y") colnames(result) <- c("SS", "DF", "V", "F", "p") result }
この関数を用いて解析してみる。引数にはデータ行列(行:個体要因、列:時点要因)を指定する。なお、1.505478e-06というのは0.000001505478ということで、6.133402e-03は0.006133402である。このe(またはEと記されることもある)は浮動小数点数といって、コンピュータ上での非常に小さな、あるいは非常に大きな値を表現するときに用いる表記法の1つである。
> thr.day <- c(2, 2, 3, 2, 1, 4) > one.week <- c(2, 4, 2, 3, 3, 5) > thr.week <- c(6, 7, 8, 6, 6, 10) > dat <- cbind(thr.day, one.week, thr.week) > dat thr.day one.week thr.week [1,] 2 2 6 [2,] 2 4 7 [3,] 3 2 8 [4,] 2 3 6 [5,] 1 3 6 [6,] 4 5 10 > paired.anova(dat) SS DF V F p T 80.111111 2 40.0555556 68.018868 1.505478e-06 I 19.111111 5 3.8222222 6.490566 6.133402e-03 E 5.888889 10 0.5888889 NA NA Y 105.111111 17 NA NA NA
参考文献