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

統計学備忘録 since2016

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

対応のある一元配置分散分析、多重比較

ID   1回目    2回目    3回目
1  93     109    124 
2  119    136    132 
3  115    121    118 
4  113    161    122 
5  123    125    111 
6  116    102    115 
7  104    115    127 
8  113    126    113 
9  111    124    132 
10 115    93     138 
11 89     105    123 
12 100    146    137 
13 126    110    127 
14 138    100    132 
15 92     134    124 
16 104    122    142 
17 102    104    134 
18 119    97     132 
19 119    118    139 
20 117    109    120 

準備として、上述のデータセットを以下のように変更します

ID   検査順   結果
1  1回目    93 
1  2回目    109 
1  3回目    124 
2  1回目    119 
2  2回目    136 
2  3回目    132
3  1回目    115 
3  2回目    121 
3  3回目    118 
・・・

一要因被験者内計画で対応のある一元配置分散分析です
また、要因1=ID, 要因2=検査順となっているので、 繰り返しのない二元配置分散分析ともよばれます
(1回目、2回目、3回目にそれぞれ繰り返して測定したら、繰り返しのある二元配置分散分析になります)

それでは上記のようなデータセットの変換を実施します

dataA <- read.table("clipboard",header = T)    

IDの配列はそれぞれ3倍必要になります

ID <- c(rep(dataA[,1],each=3))

検査順の配列は1回目、2回目、3回目を20回繰り返します

re <- c("1回目","2回目","3回目")
times <- c(rep(re,20))

最後の「結果」の配列を整えます
やや面倒ですが、元データの操作はしなことが基本ですので 辛抱してR上で展開します

#2列目~4列目のデータのみを取り出して行列に変換
zx <- data.matrix(dataA[,2:4])
#行を基準にデータを整列させるために、行列を転置
tzx <- t(zx)
#ベクトルに変換
outcome <- as.vector(tzx)

ID, 検査順, 結果を列にしたデータフレームを作成します

dataB <- data.frame(ID, times, outcome)
head(dataB)

  ID times outcome
1  1 1回目      93
2  1 2回目     109
3  1 3回目     124
4  2 1回目     119
5  2 2回目     136
6  2 3回目     132

これで準備が整いました.
分散分析を実施します

summary(aov(dataB$outcome~dataB$ID+dataB$times))

            Df Sum Sq Mean Sq F value  Pr(>F)   
dataB$ID     1      2     2.2   0.012 0.91306   
dataB$times  2   2491  1245.5   6.967 0.00199 **
Residuals   56  10011   178.8                   
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
> 

有意差があったので多重比較を実行します
テューキーHSD(Tukey's‘Honest Significant Difference’method)
全ての群の組み合わせについて母平均の差の検定を行う
ボンフェローニ(Bonferroni)
全ての対比を行う.対比較の数に応じて有意水準を調整. 対比較の数が多くなると検出力が低くなる.
テューキー(Tukey's multiple comparison test)
全ての群の組み合わせについて母平均の差の検定
シェッフェ(Scheffe's multiple comparison test)
全ての対比を行い、その中で有意なものを見つけ出す検定。
ィッシャーの最小有意差 Fisherʼs LSD
多重性が考慮されていないため3群のみに限定される.
Dunnettの方法
対照群のみとの比較

#ボンフェローニ(Bonferroni)
pairwise.t.test (dataB$outcome,dataB$times, p.adjust.method="bonferroni")

Pairwise comparisons using t tests with pooled SD 

data:  dataB$outcome and dataB$times 

      1回目  2回目 
2回目 0.3881 -     
3回目 0.0013 0.0941

P value adjustment method: bonferroni 
#テューキー(Tukey's multiple comparison test)
TukeyHSD (aov(dataB$outcome ~ dataB$times), ordered = TRUE)
Tukey multiple comparisons of means
    95% family-wise confidence level
    factor levels have been ordered

Fit: aov(formula = dataB$outcome ~ dataB$times)

$`dataB$times`
             diff        lwr    upr     p adj
2回目-1回目  6.45 -3.6360015 16.536 0.2806492
3回目-1回目 15.70  5.6139985 25.786 0.0012081
3回目-2回目  9.25 -0.8360015 19.336 0.0785692

ボンフェローニ, テューキーより1回目と3回目には有意差がみられた

ダミー変数の作成 makedummies

下のようなカテゴリカルデータをダミー変数に変更します

treat    outcome
NO  NO
NO  YES
YES NO
NO  YES
YES YES
YES YES
NO  YES
YES NO
NO  NO
NO  YES
NO  NO
NO  NO
YES YES
YES YES
YES NO
YES NO
NO  NO
NO  NO
YES YES
NO  NO
YES NO
YES NO
YES YES
NO  YES
NO  YES
YES YES
NO  NO
NO  NO
YES NO
NO  NO
NO  YES
NO  NO
NO  YES
YES NO
NO  YES

サンプルデータをコピーしてRにペーストします

s_data <- read.table("clipboard", header=T )

makedummiesをインストールします

install.packages("makedummies")
library(makedummies)

一瞬でダミー変数に代わります

d_data <- makedummies(s_data, basal_level = TRUE)

  treat NA. outcome NA..1
1      1   0       0     1
2      1   0       1     0
3      0   1       1     0
4      1   0       0     1
5      0   1       0     1
6      1   0       1     0
7      1   0       1     0
8      0   1       1     0
9      1   0       1     0

treat YESを1にしたい場合=2列目を使用
treat YESを0にしたい場合=1列目を使用
outcome YESを1にしたい場合=4列目を使用
outcome YESを0にしたい場合=3列目を使用

#treat YESを1にしたい場合=2列目   vs   outcome YESを0にしたい場合=3列目
d_data2 <- xtabs(~d_data[,2]+d_data[,3])

           d_data[, 3]
d_data[, 2]  0  1
          0  9 11
          1  7  8

name <- list(treat=c("NO","YES"),outcome=c("YES","NO"))
dimnames(d_data2)<-name
d_data2

     outcome
treat YES NO
  NO    9 11
  YES   7  8

パッケージEpiで分割表の解析を一瞬で!

install.packages("Epi")
library(Epi) 
twoby2(d_data2)

2 by 2 table analysis: 
------------------------------------------------------ 
Outcome   : YES 
Comparing : NO vs. YES 

    YES NO    P(YES) 95% conf. interval
NO    9 11    0.4500    0.2532   0.6638
YES   7  8    0.4667    0.2409   0.7070

                                    95% conf. interval
             Relative Risk:  0.9643    0.4664   1.9935
         Sample Odds Ratio:  0.9351    0.2440   3.5836
Conditional MLE Odds Ratio:  0.9369    0.2004   4.3980
    Probability difference: -0.0167   -0.3178   0.2850

             Exact P-value: 1 
        Asymptotic P-value: 0.922 
------------------------------------------------------

makedummiesの詳しい解説は、
github.com

正規分布の色塗り

plot(dnorm, -4, 4, xaxt="n")
xvals <- seq(-4, -1.96, length=10)        # -4以上-1.33以下 領域をx軸方向に10個の多角形(台形)に等分割
dvals <- dnorm(xvals)    # 対応するグラフの高さ
polygon(c(xvals,rev(xvals)),
        c(rep(0,10),rev(dvals)),col=5)       # 塗りつぶす

xvals2 <- seq(1.96, 4, length=10)        
dvals2 <- dnorm(xvals2)   
polygon(c(xvals2,rev(xvals2)),
        c(rep(0,10),rev(dvals2)),col=5) 

name<-c("-1.96","0","1.96")
axis(side=1,at=c(-1.96,0,1.96),labels=name)  

[f:id:yoshida931:20180719112607p:plain:500]

ヒストグラムのビン調整

head(iris)

  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

変数Speciesがsetosaのデータを抽出
列Speciesがsetosatである行を全て抜き出す

iris_setosa <- iris[iris$Species=="setosa",]

パッケージMASSをインストール

install.packages("MASS")
library("MASS")

ヒストグラムのビン調整(0.1, 0.2, 0.3, 0.4)

par(mfrow = c(2, 2))
truehist(iris_setosa$Sepal.Length, h=0.1, prob=F) 
truehist(iris_setosa$Sepal.Length, h=0.2, prob=F) 
truehist(iris_setosa$Sepal.Length, h=0.3, prob=F) 
truehist(iris_setosa$Sepal.Length, h=0.4, prob=F) 
par(mfrow = c(1, 1))   

f:id:yoshida931:20180716223044p:plain:w600

ノンパラメトリック 相関係数

     [,1] [,2] [,3]
[1,]   27   14    5
[2,]   10   17   26
[3,]    5   12   50

上記の分割表から行の順序スコアと列の順序スコアを算出してデータセットを作成します

#RANK 
xr1 <- c(rep(23.5,46),rep(73,53),rep(133,67))
xc1 <- c(rep(21.5,27),rep(64,14),rep(125,5),rep(21.5,10),rep(64,17),
         rep(125,26),rep(21.5,5),rep(64,12),rep(125,50))

データのイメージ

xr1_xc1 <- data.frame(xr1, xc1)

相関係数を求めます

cor_r_p <- cor.test(xr1,xc1,method = "pearson") #ピアソン
cor_r_s <- cor.test(xr1,xc1,method = "spearman") #スピアマン

上記の分割表から表スコアを0,1,2としてデータセットを作成します

#TABLE 
xr2 <- c(rep(0,46),rep(1,53),rep(2,67))
xc2 <- c(rep(0,27),rep(1,14),rep(2,5),rep(0,10),rep(1,17),
         rep(2,26),rep(0,5),rep(1,12),rep(2,50))

データセットにイメージ

xr2_xc2 <- data.frame(xr2, xc2)

表スコアから相関係数を求めます

cor_t_p <- cor.test(xr2,xc2,method = "pearson") #ピアソン
cor_t_s <- cor.test(xr2,xc2,method = "spearman") #スピアマン

Rで作る分割表

次のデータを分割表にして解析してみます
まずはコピーして

筋トレ    歩行練習    効果
あり  あり  有効
あり  あり  有効
あり  あり  有効
あり  あり  有効
あり  あり  有効
あり  あり  有効
あり  あり  有効
あり  あり  有効
あり  あり  有効
あり  なし  有効
あり  なし  有効
あり  なし  有効
あり  なし  無効
あり  なし  無効
あり  なし  無効
なし  あり  有効
なし  あり  有効
なし  あり  有効
なし  あり  有効
なし  あり  無効
なし  なし  有効
なし  なし  有効
なし  なし  無効
なし  なし  無効
なし  なし  無効
なし  なし  無効
なし  なし  無効
なし  なし  無効
なし  なし  無効
なし  なし  無効

Rにペーストします

x <- read.table("clipboard",header=T)

2変数の分割表の作成

x2 <- xtabs(~筋トレ+歩行練習, data=x)

       効果
筋トレ 無効 有効
  あり    3   12
  なし    9    6

周辺合計

addmargins(x2)

      効果
筋トレ 無効 有効 Sum
  あり    3   12  15
  なし    9    6  15
  Sum    12   18  30

行や列の入れ替え

x3 <- cbind(x2[,2],x2[,1])
colnames(x3) <- paste(c("有効","無効"))
x3

     有効 無効
あり   12    3
なし    6    9

#行の入れ替えはrbindで!

#または行列を直接作って解析
xx <- matrix(c(12,6,3,9),2,2)
colnames(xx) <- paste(c("有効","無効"))
rownames(xx) <- paste(c("あり","なし"))

     有効 無効
あり   12    3
なし    6    9

転置

t(x2)   # aperm(x2)でも同じ

      筋トレ
効果   あり なし
  無効    3    9
  有効   12    6

歩行練習で層別

x4 <- xtabs(~筋トレ+, data=x)

, , 歩行練習 = あり

      効果
筋トレ 無効 有効
  あり    0    9
  なし    1    4

, , 歩行練習 = なし

      効果
筋トレ 無効 有効
  あり    3    3
  なし    8    2

コクラン-マンテル-ヘンツェル検定で共通オッズ比の検定

mantelhaen.test(x4) 

    Mantel-Haenszel chi-squared test with continuity correction

data:  x4
Mantel-Haenszel X-squared = 1.4761, df = 1, p-value = 0.2244
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 0.02148238 1.42558684
sample estimates:
common odds ratio 
            0.175 

ついでに2変数の分割表解析

xx <- matrix(c(12,6,3,9),2,2)
colnames(xx) <- paste(c("有効","無効"))
rownames(xx) <- paste(c("あり","なし"))

install.packages("Epi")  
library(Epi)            
twoby2(xx)

2 by 2 table analysis: 
------------------------------------------------------ 
Outcome   : 有効 
Comparing : あり vs. なし 

     有効 無効    P(有効) 95% conf. interval
あり   12    3        0.8    0.5302   0.9341
なし    6    9        0.4    0.1918   0.6519

                                   95% conf. interval
             Relative Risk: 2.0000    1.0240   3.9063
         Sample Odds Ratio: 6.0000    1.1717  30.7246
Conditional MLE Odds Ratio: 5.6137    0.9423  44.9348
    Probability difference: 0.4000    0.0504   0.6398

             Exact P-value: 0.0604 
        Asymptotic P-value: 0.0315 
------------------------------------------------------

Rstudioの小ネタ (パッケージやファイルの保存方法)

Rstudioを閉じても、PC再起動してもファイルの読み込みやインストールしたパッケージは残せます (ただしパッケージは休んでいますので、起動するときにはlibraryで起こしましょう)

データ処理する前に必ず行う作業は以下の通りです
まずRstudioを起動させた状態で・・・
FileメニュからNewProject → NewProject → NewProjectと進んでいきます
Projectの名前を適当に書いてください(例:test001)、次に場所をしてしてください(as you like).

以下のように一つのフォルダの中に「---.Rhistory」と「---.Rproj」が入ってたら OK
データファイルも同じフォルダに収納します
f:id:yoshida931:20180604094238p:plain

次回から「---.Rproj」をダブルクリックしてRstudioを起動させてください
前回までに取り込んだデータは右上に表示され、そのまま使用可能です.
また右下に同じフォルダにあるデータファイルが表示されます.ここから直接データを取り込むこともできます.
またインストールしたパッケージも記憶されていますので、library( )で起動させるだけで使用できます.
f:id:yoshida931:20180604094641p:plain