理学療法士がまとめたノート

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

テキストファイルに検定結果を出力

テキストファイルに検定結果を出力
投稿日2017.11.14

r-de-r様からコメントいただきましたので修正しております
素直に感動しております

テキスト形式のファイルに検定結果を書き込む練習をします
sprintf:書式指定変換した出力を文字列に格納
cat:文字列を表示する基本的な関数, file=でファイルに書き込み

備忘録
\n 改行
%.3f 小数点3桁まで出力
append = F:ファイルに上書き
append = T:ファイルに追加

Rのサンプルirisを使います

x <- iris$Sepal.Length[51:75]
y <- iris$Petal.Length[101:125]
#noteというファイルに、xとyの平均を書き込みます
cat("x平均", mean(x), "y平均", mean(y), file = "note.txt") #noteというファイル名

f:id:yoshida931:20170908160200j:plain

sprintfの使い方
xの平均を小数点以下1桁まで、項目名「Xの平均」と付いた文字列にして取り出します

sprintf("xの平均 %.1f", mean(x))
"xの平均 6.0"

t検定を実施して、strで中身を確認してみます

t <- t.test(x,y,var.equal = T)
str(t)

Two Sample t-test

data: x and y
t = 2.1954, df = 48, p-value = 0.033
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.03131358 0.71268642
sample estimates:
mean of x mean of y
6.012 5.640


List of 9
$ statistic : Named num 2.2
..- attr(*, "names")= chr "t"
$ parameter : Named num 48
..- attr(*, "names")= chr "df"
$ p.value : num 0.033
$ conf.int : atomic [1:2] 0.0313 0.7127
..- attr(*, "conf.level")= num 0.95

t値、p値、95%信頼区間のみ取り出してテキストファイル「note」に書き出してみます

tv <- sprintf("t値 %.4f\n", t[[1]]) 
pv <- sprintf("p値 %.4f\n", t[[3]]) 
ci <- sprintf("95%信頼区間 %.4f %.4f\n", t[[4]][1], t[[4]][2])     #ここがポイントです(%を小文字にする場合には%%と書きます)
cat(tv,pv, ci, sep="",file = "note.txt")

catのデフォルトでsep=" "となっているので改行した後にスペースが入るそうです.スペースなし(sep="")にして書き直しました.無駄なスペースが無くなってスッキリしました.

f:id:yoshida931:20171114133508j:plain:w300

参考と引用 Rプログラム (TAKENAKA's Web Page)

単回帰分析

投稿日2017.2.2
更新日2017.11.13

線形回帰分析

説明変数と目的変数との関係を直線でモデル化する回帰分析。

Xi <- c(35, 45, 55, 65, 75)
Yi <- c(114, 124, 143, 158, 166)
plot(Xi, Yi, pch = 16)
abline(lm(Yi~Xi), col = 2)

f:id:yoshida931:20171108084407j:plain:w300

どのようにして赤のラインを引いたのか?を考えていきます.

母回帰方程式 

Yi=\beta1 + \beta2 * Xi + ei\ \ (i=1,2,3,,,n)

これが大前提となる母集団のモデルです.一般的に知ることができないので、このモデルについて推定・検定することを回帰分析と呼んでいます.

母回帰係数 \beta1, \beta2
誤差項 ei (これが回帰分析のポイントになります)
 ・期待値=0
 ・分散は一定
 ・異なる誤差項は無相関

E(Yi)=\beta1+\beta2*Xi

  

最小二乗法

最小二乗法により、母回帰係数(\beta1, \beta2) を推定します. その推定値を(\hat{\beta1}, \hat{\beta2})とします

母回帰方程式より、誤差項ei=Yi-(\beta1 + \beta2 * Xi) となります

最小二乗法とは、S= \sum {ei}^2とおいた場合に、Sが最小になる(\hat{\beta1}, \hat{\beta2})を求める手法です.

S=\sum({Yi}^2-2β1Yi-2β2XiYi+{β1}^2+2β1β2Xi+{β2}^2{Xi}^2)

Sを最小にするために、\beta1, \beta2 でそれぞれ微分した二つの方程式を解くことで次式を求めます.

標本(偏)回帰係数
\hat{\beta2}=\frac{\sum{(Xi-\bar{Xi})(Yi-\bar{Yi})}}{\sum{(Xi-\bar{Xi})}^2}  \hat{\beta1}=\bar{Yi}-\hat{\beta2}\bar{Xi}

標本回帰方程式
Y=\hat{\beta1}+\hat{\beta2}\bar{X}

RのサンプルCO2(草の耐寒性実験)からデータを抜き出して回帰分析を行ってみます.

# データフレームから必要な部分だけ抜き出します.(CO2濃度95ml/Lと175ml/Lの二酸化炭素摂取量)
x <- subset(CO2, subset = conc == 95, select = c(uptake))     
y <- subset(CO2, subset = conc == 175, select = c(uptake)) 
x <- x$uptake
y <- y$uptake

summary(lm(formula=y~x) )

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.1062 -1.8958 -0.7628  2.0096 10.0376 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   1.4784     6.1126   0.242  0.81377   
x             1.6972     0.4876   3.480  0.00592 **
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.424 on 10 degrees of freedom
Multiple R-squared:  0.5478, Adjusted R-squared:  0.5026 
F-statistic: 12.11 on 1 and 10 DF,  p-value: 0.005917

plot(x,y)
abline(lm(y~x))

f:id:yoshida931:20171113125825j:plain:w300
標準化残差を確認してみます

ei<-y-(1.4784+1.6972*x)
eis<-ei/sqrt( (sum(ei^2) )/(length(x)-2) )  #標準化残差
plot(eis,xlab="CO2濃度95ml/L",ylab="CO2濃度175ml/L")
abline(h=0, lty=2, col=2)

f:id:yoshida931:20171113131309j:plain:w300

次にRの関数で求めた回帰モデル y = (1.4784+1.6972*x) を最小二乗法で求めてみます.
最小二乗法より標本回帰係数\hat{\beta2} を求めます

sum ( ( y - mean(y) ) * ( x - mean(x) ) ) / sum ( ( x - mean(x) )^2 )
= 1.697206

標本回帰係数\hat{\beta1} (y切片)を求めます

mean(y) - 1.697206*mean(x)
=1.478416

回帰モデル
y = 1.4784 + 1.6972*x

参考) 豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012

九州理学療法士・作業療法士合同学会2017in宮崎で紹介したスクリプト

モンテカルロ法によるt検定の検定力分析
モンテカルロ法については以下を参照ください

yoshida931.hatenablog.com

tval <- numeric(length = 10000)  #t値(空のベクトルを設定)
cs <- 0
for(i in 1:10000){
 x <- rnorm(n = 30, mean = 10, sd = 3)   # xの乱数30個生成
  y <- rnorm(n = 30, mean = 12, sd = 3)   # yの乱数30個生成
   res <- t.test(x, y, var.equal = T)       #xとyのt検定の結果
    tval[i] <- res[[1]]               #t値のみを抜き取る
     cs <- cs + ifelse(res[[3]] <  0.05, 1,0) 
}

#cs:p値が0.05より小さい場合には1を足し、大きい場合には0を足す
#t検定を1000回繰り返した結果、0.05より小さいP値の個数となる

cs/10000

データの要約、散布図、エラーバー付き平均値の一括処理

x<-iris$Sepal.Length[51:75]
y<-iris$Petal.Length[101:125]


sumxy <- function(x, y, TEST01){
 
  library(pwr)

  
  cat1 <- c("Min.", "1st Qu", "Median", "Mean", "3rd Qu.", "Max.")   
  xs <- summary( x )                        
  ys <- summary(y)
  
  cat(cat1,"\n","x",xs,"\n","y",ys,"\n\n",file = "TEST01");
  tx <- t.test(x);             
  cat2<- c("tvalue","pvalue","95%CI");  
  txt <- tx$statistic;         
  txp <- tx$p.value;          
  txc1 <- tx$conf.int[1];  
  txc2 <- tx$conf.int[2];      
  x1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f",txt,txp,txc1,txc2); 
  
  ty <- t.test(y);    
  tyt <- ty$statistic;
  typ <- ty$p.value;
  tyc1 <- ty$conf.int[1];
  tyc2 <- ty$conf.int[2];
  y1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f", tyt, typ, tyc1, tyc2);
  
  txy <- t.test(x,y,var.equal = T);
  txyt <- txy$statistic;      
  txyp <- txy$p.value;         
  txyc1 <- txy$conf.int[1];       
  txyc2 <- txy$conf.int[2];       
  xy1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f",txyt,txyp,txyc1,txyc2);
  
  d <- (mean(x) - mean(y)) / sqrt((var(x) + var(y)) / 2);
  p <- pwr.t.test(n = length(x), d = d, sig.level = 0.05);
  
  xy <- c(x,y);
  zu1 <- pdf(file = "散布図.pdf");
  cate <- c(rep("1",length(x)),rep("2", length(y)));
  plot(xy ~ cate,xlim = c(0.3, 2.8),xaxp=c(1, 2, 1),xaxt="n",xlab = "", ylab = "") ;
  name<-c("X","Y")
  axis(side=1,at=c(1,2),labels=name);
  dev.off();
  zu1
  
  zu2 <- pdf(file = "エラーバー.pdf");
  se2 <- c(sqrt(var(x)/(length(x))), sqrt(var(y)/(length(y))));
  m2 <- c(mean(x), mean(y));
  plot(m2, xlim=c(0.5,2.5), type="b", ylim=c(5,7), xaxt="n",xlab = "", ylab = "");
  arrows(1:2,m2-se2,
         1:2,m2+se2,
         lwd=1.0,angle=90,length=0.1,code=3);
  name<-c("X","Y")
  axis(side=1,at=c(1,2),labels=name);
  dev.off();
  zu2
  
  zu3 <- pdf(file = "箱ひげ図.pdf");
  cate <- c(rep("1",length(x)),rep("2", length(y)));
  boxplot(xy ~ cate,xlim = c(0.3, 2.8),xaxp=c(1, 2, 1),xaxt="n",xlab = "", ylab = "") ;
  name<-c("X","Y")
  axis(side=1,at=c(1,2),labels=name);
  dev.off();
  zu3
  
  cat(cat2, "\n", "x", x1, "\n", "y", y1, "\n", "xy", xy1, "\n",
      "x:1標本t検定, y:1標本t検定, xy:独立2群t検定", "\n", "\n",
      "独立2群t検定の検定力", p[[4]], "\n", file = "TEST01", append =T);
}

sumxy(x, y, "TEST01")
# 出力されたファイルはDocumentsにあると思います

分散分析(繰り返し vs 反復)

月元敬;対応あり, 繰り返し, 反復: 心理データ解析基礎において区別すべき 3 つの 「R」. 岐阜大学教育学部研究報告. 人文科学, 2015, 64.1: 67-73.

「繰り返し」と「反復」は区別なく混同している参考書も多いようですが、今回は上記の文献をもとに違いを整理しておきます.

対応 relatedness
反復 repetition
繰り返し replication    繰り返し数 number of replication

まず、これらの言葉の整理をするためには、被験者も因子として考えることが重要です

繰り返し

実験内容としては「水準の設定を最初から施してサンプルを一つとる」ことになります. 「繰り返し」とは、この実験自体を繰り返すことを意味しています.

# 検査A,B,Cを同じ条件で、それぞれ別の被験者に実施する実験.
# ランダムに選択した被験者に対して、それぞれの検査を4回繰り返します(繰り返し数4回)   
# 完全無作為化法

検査A 被験者8,  被験者4,  被験者5,  被験者10
検査B 被験者3,  被験者1,  被験者12, 被験者7
検査C 被験者11, 被験者2,  被験者6,  被験者9

一元配置分散分析(1要因3水準、繰り変え数4回) つまり3×4=12回の実験をランダムな順序で実施することになります. これを「対応のない」という表現をする場合もあります

反復

 繰り返しとは全く概念が異なることに注意しておきましょう!
反復」はフィッシャーの3原則のうちの1つです. 反復測定とは、「ある条件下において、複数の異なる個体から、それぞれ1つだけの観測値を繰り返して測定する方法」というのが正しいでしょう。

乱塊法 

# 例)繰り返しのある一元配置分散分析
#     無作為抽出した繰り返し数3回の例・・・実験順序の系統誤差が大きい

検査α 被験者6→被験者8→被験者9
検査β 被験者4→被験者1→被験者2
検査γ 被験者5→被験者3→被験者7

# 次のように書き換えると偏りが見えてきます

検査順序
β → β → γ → β → γ → α → γ → α → α

# 次に1日に3名実施して場合、1日を1ブロック(同じ条件で実験できる1領域)としたら・・・
# 偏りのある採取方法になっています

1日目 2日目 3日目
  β   β   γ
  β   γ   α
  γ   α   α

# 採取日による偏りを解消し、.

1日目 α → β → γ
2日目 β → γ → α
3日目 β → α → γ

実験順序
     1日目 2日目 3日目
      α   β   β
検査要因  β   γ   α
      γ   α   γ
# 繰り返しのなし二元配置分散分析

無作為に抽出されたサンプルに対して、制御因子(再現性あり)ブロック因子(再現性なし)の2因子を用いる特殊な二元配置法が乱塊法と呼ばれています.検査項目が制御因子となります.1日目、2日目、3日目には再現性がないのでブロック因子となります.そのブロックが順次増えていくことを「反復」といいます.

参考)フィッシャーの三原則 randomization・replication・ブロック化(local control)

# 乱塊法
# 被験者もブロック因子として考えられます

 検査A
        A1   A2   A3   A4
      A  4    2    1    3
被験者要因 B   1    2    4    3
      C   2    3    1    4   (数字は検査順番)

# 繰り返しのない二元配置分散分析とよばれている配置になります
# 一要因被験者内計画で対応のある一元配置分散分析です

しかし!!! 被験者A,B,Cに検査A1を実施して、次に検査A2、A3、A4の実験を行うことを「反復」と呼ぶ場合があるようです. どちらでも概念的にも実用的にも問題は内容ですが誤用と考えられます.

対応のある、対応のない

繰り返しがある場合には、別々の被験者からのデータ収集になるため「対応がない」と表現されています.また繰り返しがない場合には、上述の例のように被験者AのA1、A2、A3、A4は同一被験者のデータになるので、「対応がある」という表現が使用されているようです.

まとめ

対応のない1要因で分類される多群の検定
繰り返しのある一元配置分散分析


対応のある1要因で分類される多群の検定
繰り返しのない二元配置分散分析
反復測定による一元配置分散分析


対応のない2要因で分類される多群の検定
繰り返しのある二元配置分散分析


対応のある2要因で分類される多群の検定
反復測定による二元配置分散分析


分散分析の基本

まだ理解できていない.なので書き直し・・・ 一元配置分散分析

言葉の整理

要因(factor), 因子(factor):実験結果に影響を与える要素.それぞれの分野で使い分ける場合もあるので注意.このブログでは要因と因子の区別をせず「要因」で統一します.
〇元配置分散分析:〇は要因の数
水準:要因を分ける条件( 要因に含まれる項目 )
:=水準
平方和:偏差の二乗の合計
変動:=平方和
水準内変動:= 群内変動 = 誤差変動 = 残差変動
平均平方和:平方和を自由度で割った値(分散)
効果:目的変数への影響
主効果:水準によって効果が異なる場合、その要因に主効果があるという
交互作用:2つの要因が目的変数に対して互いに影響を及ぼしている場合
   相殺効果:一方の要因の効果が高いと他方の要因の効果がそれをうち消すこと
   相乗効果:2つの要因により効果が大きくなること
要因効果:主効果と交互作用の総称

分散分析とは

  • 分散分析は、3群以上のそれぞれの群の母平均が等しいという帰無仮説のもとに検定を行う分析方法です.t検定を繰り返して使用することはできません.例えば3群の比較を5%有意水準でそれぞれt検定を行った場合、帰無仮説は「3通り全ての検定において差がみられない」となり、対立仮説は「少なくとも一組には差がみられる」となります.つまり、有意水準は0.05ではなく1- 0.953 = 0.142 で検定したことになります.つまり、第一種の過誤が生じる危険性が14.2%となります.

  • 分散分析では、各群の母平均の推定を行うのではなく、各群が及ぼす効果(要因の効果)の分散を推定します.そのために群間の平均の変動、各群の郡内変動から平均平方和 (分散) を求め、要因の効果の有無をF分布から推測します.


一元配置分散分析(対応なし)

対応のない一要因で分類される多群の検定(要因分散分析).
実験A,B,Cのそれぞれにn人が含まれており、全部で(3×n)人の被験者からなる実験です.
群間に差があるかどうかを知りたい

        1   2   3  ・・・  n
実験A  x11 x12 x12 ・・・  x1n
実験B  x21 x22 x23 ・・・  x2n
実験C  x31 x32 x33 ・・・  x3n
実験D  x41 x42 x43 ・・・  x4n  
  
サンプル(samp)  繰り返し数 5
実験A   10   6  10   9  10
実験B   10   5   5  12   4
実験C    5   4  11   4   6  
実験D    9   5   2   3   1

上記サンプル(samp)を使用して考えていきます.
実験結果のバラつきは分散として得られ、分散の原因として次の2点が考えられます.

  • 原因1)実験A、B、Cの効果によるもの(要因の効果
  • 原因2)被験者の問題、測定の誤差、偶然などによるもの(誤差項

重要 データの構造を表すモデル
誤差項(測定値 - 各群の平均 )は互いに独立で、正規分布に従うことを仮定しています.
測定値 = 全平均 + (総平均 - 各群の平均) + (測定値 - 各群の平均 )より

全変動 = 水準間変動 + 水準内変動  
∑∑(測定値-全平均)^2 = n*∑(総平均-各群の平均)^2+∑(測定値-各群の平均 )^2 

  ポイント:∑∑2(総平均-各群の平均)(測定値-各群の平均) = 0

分散分析の帰無仮説
実験A,B,C,Dの母平均は等しい.実験A,B,C,Dの母分散は0である.
したがいまして群間の差に着目して検定することになります

F分布
全ての群の母平均が等しい場合、統計量FはF分布に従います.

F=\frac{群間平方和}{群間自由度}\div\frac{群内平方和}{群内自由度}

サンプル(samp)の場合、群間平方和の自由度は3、群内平方和の自由度は5*4-4=16となります.したがって自由度(3,4)のF分布に従うことになります

curve(df(x, 3, 16), from = 0, to = 5, type = "l", ylim = c(0, 0.8))
f:id:yoshida931:20171031142612j:plain:w300

検定の結果として帰無仮説が棄却された場合には、要因の効果が認められる(群間に差がある)と推測します.

Rを使って一元配置分散分析(対応なし)

data <- c(10,10,5,9,6,5,4,5,10,5,11,2,9,12,4,3,10,4,6,1)
exp <- factor(c(rep(c("A","B","C","D"),5)))
summary(aov(data~exp))

            Df Sum Sq Mean Sq F value Pr(>F)  
exp          3  66.15   22.05   2.579 0.0898 .
Residuals   16 136.80    8.55                 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

結果、有意水準0.05で群間の差は認められませんでした.

F値の算出方法
一元配置分散分析 (対応なし) F値の算出方法 - 統計学備忘録 since2016


一元配置分散分析(対応あり)

対応のある1要因で分類される多群の検定(反復測定分散分析)
正確には繰り返しのない二元配置分散分析になります
サンプル:被験者4名(S1,S2,S3,S4)が5種類(testA,B,C,D,E)の検査を受けた結果
同じ被験者が繰り返して試験を受けることになる(繰り返し数は5回)
群間に点数差があるかどうかを知りたい
  ・被験者間に差があるか
  ・検査間に差があるか

サンプル(samp)  
data <- c(10,10,8,9,6,5,4,5,10,12,11,9,9,8,7,3,8,4,6,2)
data <- array(data, dim = c(4,5))
name <- list(Sub=c("s1","s2","s3","s4"), TEST=c("A","B","C","D","E"))
dimnames(data)<-name
data

    TEST
Sub   A B  C D E
  s1 10 6 10 9 8
  s2 10 5 12 8 4
  s3  8 4 11 7 6
  s4  9 5  9 3 2

Rを使って一元配置分散分析(対応あり)

data <- c(10,10,8,9,6,5,4,5,10,12,11,9,9,8,7,3,8,4,6,2)
sub <- factor(c(rep(c("s1","s2","s3","s4"),5)))
test <- factor(c(rep("A",4),rep("B",4),rep("C",4),rep("D",4),rep("E",4)))
summary(aov(data~sub+test))
  
            Df Sum Sq Mean Sq F value   Pr(>F)    
sub          3   24.2   8.067   3.681 0.043466 *  
test         4   99.7  24.925  11.373 0.000475 ***
Residuals   12   26.3   2.192                     
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

対応ありを無視したら…

summary(aov(data~sub))

#結果
            Df Sum Sq Mean Sq F value Pr(>F)
sub          3   24.2   8.067   1.024  0.408
Residuals   16  126.0   7.875 


summary(aov(data~test))

            Df Sum Sq Mean Sq F value  Pr(>F)   
test         4   99.7  24.925   7.403 0.00168 **
Residuals   15   50.5   3.367                   
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

自由度、平方和は同じ値で、もちろん平均平方も同じ値になっています.
しかし対応ありを無視したら、残差が大きなり、F値が小さくなっていることが確認できました.
対応ありを考慮すると有意差が出やすくなることが理解できました.


二元配置分散分析(対応なし)

対応のない2要因で分類される多群の検定 (要因分散分析)
2つの要因の組み合わせによる目的変数への影響を分析します

#サンプル:テストA、Bで各ポイント(p1,p2,p3,p4,p5)を獲得した人数
data <- c(10,10,8,9,6,5,4,5,10,12,11,9,9,8,7,3,8,4,6,2,4,10,8,8,4,6,2,10,12,13)
data <- array(data, dim = c(5,3,2))
name <- list(point=c("p1","p2","p3","p4","p5"), ","=c("","",""),TEST=c("A","B"))
dimnames(data)<-name
data

#二元配置分散分析(対応なし)のデータ配列のイメージ
#繰り返し数 3
TEST = A
point         
   p1 10  5 11
   p2 10  4  9
   p3  8  5  9
   p4  9 10  8
   p5  6 12  7

TEST = B
point        
   p1 3  4  6
   p2 8 10  2
   p3 4  8 10
   p4 6  8 12
   p5 2  4 13

分析はこのページで確認
二元配置分散分析 (1) - 統計学備忘録 since2016


二元配置分散分析(対応あり)

対応のある2要因で分類される多群の検定 (反復測定分散分析)
2つの要因の組み合わせによる目的変数への影響を分析します

#例)5人の被験者がAMとPMに試験(t1,t2,t3)を受けた時の点数
# 時刻(AM,PM)と試験(t1,t2,t3)
data <- c(10,10,13,9,12,9,11,9,10,12,7,9,9,8,7,3,8,4,6,2,4,10,8,8,6,9,8,10,12,10)
data2 <- array(data, dim = c(5,3,2))
name <- list(sub=c("s1","s2","s3","s4","s5"), "test"=c("t1","t2","t3"),"time"=c("AM","PM"))
dimnames(data2)<-name
data2  

#二元配置分散分析(対応あり)のデータ配列のイメージ
, , time = AM

    test
sub  t1 t2 t3
  s1 10  9  7
  s2 10 11  9
  s3 13  9  9
  s4  9 10  8
  s5 12 12  7

, , time = PM

    test
sub   t1 t2 t3
  s1  3  4  9
  s2  8 10  8
  s3  4  8 10
  s4  6  8 12
  s5  2  6 10

分析方法はちょっと特別な方法になります

#準備
sub2 <- factor(c(rep(c("s1","s2","s3","s4","s5"),6)))
tes2 <- factor(rep(c(rep("t1",5),rep("t2",5),rep("t3",5)),2))
tim2 <- factor(c(rep("AM",15),rep("PM",15)))

#被験者による効果、交互作用による効果を調べます
sub2
sub2:poi2
sub2:tes2
sub2:poi2:tes2

summary(aov(data~tes2*tim2+Error(sub2+sub2:tes2+sub2:tim2+sub2:tes2:tim2)))

#結果
Error: sub2
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  4  19.53   4.883               

Error: sub2:tes2
          Df Sum Sq Mean Sq F value Pr(>F)
tes2       2  8.267   4.133   2.988  0.107
Residuals  8 11.067   1.383               

Error: sub2:tim2
          Df Sum Sq Mean Sq F value Pr(>F)  
tim2       1  45.63   45.63   11.75 0.0266 *
Residuals  4  15.53    3.88                 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Error: sub2:tes2:tim2
          Df Sum Sq Mean Sq F value  Pr(>F)   
tes2:tim2  2  81.07   40.53   11.47 0.00447 **
Residuals  8  28.27    3.53                   
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

交互作用の効果が有意になっています.したがって、単純効果の検定を行います.さらに単純効果が有意な場合には多重比較を実施します.次回、投稿する予定です.



参考)
山田 剛史, 杉澤 武俊, 村井 潤一郎; Rによるやさしい統計学, オーム社 , 2008
小野 滋;読めば必ずわかる分散分析の基礎  http://elsur.jpn.org/resource/anova.pdf

感度と特異度

感度と特異度
数か所間違いがありましたので再投稿します  
投稿日2017.7.31
更新日2017.10.26

      検査陽性   検査陰性  
罹患    ×(+)     ×(-)
健常    ○(+)(-)  

感 度 (陽性反応的中度) = \frac{×(+)}{×(+)\ +\ ×(-)} 罹患、検査陽性

特異度 (陰性反応的中度) = \frac{○(-)}{○(+)\ +\ ○(-)} 罹患なし、検査陰性

偽陰性 = \frac{×(-)}{×(+)\ +\ ×(-)}  罹患、検査陰性

偽陽性 = \frac{○(+)}{○(+)\ +\ ○(-)}  罹患なし、検査陽性

感 度 = 1 - 偽陰性
特異度 = 1 - 偽陽性

Rを使用してROC曲線を描いてみます

# 準備
install.packages("ROCR")  
library(ROCR) 

# サンプルを用意します(参考文献)  
--------------------------------------------
         01234点
    罹患    1     0     6    11    12
罹患なし     9     2    11     8     0    
--------------------------------------------

# 検査の各点数の該当者数(0は含めない)
x <-c (rep(0, 10),rep(1, 2),rep(2, 17),rep(3, 19),rep(4, 12))
# xに罹患=1、罹患なし=0を当てはめます
y <-c (1, rep(0, 9), rep(0, 2), rep(1, 6), rep(0, 11), rep(1, 11), rep(0, 8), rep(1, 12))
# 次の操作はベクトルxがソートされていない場合に使用
pd <- prediction(x,y)
# 次の操作で偽陽性率vs感度を求めます
# tpr : True positive rate(感度)、fpr : False positive rate(偽陽性率)
roc <- performance (pd, "tpr", "fpr")
# 一気に計算してくれます!
"x.values":偽陽性率
  [[1]]
[1] 0.0000000 0.0000000 0.2666667 0.6333333 0.7000000 1.0000000
"y.values":感度
  [[1]]
[1] 0.0000000 0.4000000 0.7666667 0.9666667 0.9666667 1.0000000
# rocの散布図がROC曲線となります
plot( roc )

f:id:yoshida931:20171026124248j:plain

aucの算出 (area under the curve) 

auc.tmp <- performance (pd, "auc")
auc <- as.numeric (auc.tmp @ y.values)

 auc = 0.8327778

サンプルデータ
出典)柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010,p12 

参考ブログ  k-lab http://www.nemotos.net/?p=836
ほとんどそのまま使用しております.
ありがとうございました.

 

Fisherの直接法

Fisherの正確確率検定やFisherの直接確率検定、他にFisherの正確検定などと呼ばれています(統一してくれれば良いのにといつも思います).もともとカイ二乗検定は近似法でP値を求めています.一つのセルに度数が4以下が存在する場合には近似法の精度が落ちるため、直接P値を計算します.

袋の中にN個の玉が入っています.赤玉M個、白玉N-M個.

f:id:yoshida931:20171020153149j:plain:w300

個数 残り 総数
赤玉 x M-x M
白玉 n-x (N-M)-(n-x) N-M
n N-n N

帰無仮説H_0: 赤玉がx個出る確率 = 白玉がn-x個出る確率
対立仮説H_1: 赤玉がx個出る確率 < 白玉がn-x個出る確率

上記のような2×2表が生起する確率を直接求めることになります.周辺和が固定されているので、その確率は超幾何分布に従います.したがってN個入りの袋からn個取り出したとき、赤玉の数がx個・白玉がn-x個となる確率は次のようになります.

P(X=x)=\frac{_MC_x\times_{N-M}C_{n-x}}{_NC_n}

Fisherの直接法ではこの確率を使用します. それでは参考・引用文献の例題を参考にして勉強していきます.

例) 小規模な二群並行ランダム化試験

有効 無効
対照群 7 8 15
試験群 12 2 14
合計 19 10 29

帰無仮説H_0: 対照群が有効となる確率 = 試験群が有効となる確率
対立仮説H_1: 対照群が有効となる確率 < 試験群が有効となる確率

上記のような実験結果から試験群の有効率が対照群の有効率より大きいことを証明します.まずは何も考えずRを使用してカイ二乗検定を行ってみます.

da <- c(7,12,8,2)
md <- matrix(da,2,2)
ch <- chisq.test ( md, correct=F )
  
    Pearson's Chi-squared test

data:  md
X-squared = 4.8871, df = 1, p-value = 0.02706

Warning message:
In chisq.test(md, correct = F) :  カイ自乗近似は不正確かもしれません 

次に直接法で計算してみます.
対立仮説H_1: 対照群が有効となる確率 < 試験群が有効となる確率の方向を示すデータの確率を求めます.

     [有効] [無効]
[対照群]    a   b
[試験群]    c   d

c>12 となる確率が対立仮説の方向を強く示すデータとなります

a,b,c,d = 7,8,12,2
a,b,c,d = 6,9,13,1
a,b,c,d = 5,10,14,0

それぞれが生起する確率を求めます

c1 <- choose(15,7)
c2 <- choose(14,12)
c3 <- choose(29,19)
c1*c2/c3
 = 0.02923538

c21 <- choose(15,6)
c22 <- choose(14,13)
c23 <- choose(29,19)
c21*c22/c23
= 0.003498251

c31 <- choose(15,5)
c32 <- choose(14,14)
c33 <- choose(29,19)
 = 0.000149925

p-value = c31*c32/c33+c21*c22/c23+c1*c2/c3
        = 0.03288356

#Rの関数を使用
mx=matrix(c(7, 8, 12, 2), nrow=2, byrow=T)
fisher.test(mx,alternative = "less")  #ここでは割合の小さい対照群が上段にあるので"less"

    Fisher's Exact Test for Count Data

data:  mx
p-value = 0.03288
alternative hypothesis: true odds ratio is less than 1
95 percent confidence interval:
 0.00000 0.86222
sample estimates:
odds ratio 
 0.1565941 

p-value=0.03288356
c=12以上になる確率は <0.05 ということになりました.したがって、「a,b,c,d = 7,8,12,2」というサンプルは、帰無仮説のもとで起こりえない(有意水準0.05)組合せで、対照群が有効となる確率 < 試験群が有効となる確率という結果になります.

参考・引用
柳川 堯 , 荒木 由布子;バイオ統計の基礎―医薬統計入門,近代科学社 ,2010, 181-183