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

統計学備忘録 since2016

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

対応のあるデータのグラフ

id <- 1:20
group <- c(rep("A",10),rep("B",10))
x <- c(rnorm(10,20,2), rnorm(10,15,3))
y <- c(rnorm(10,25,1), rnorm(10,20,0.5))
data <- data.frame(id, group, x, y)

x <- c(1,2)
t_data <- t(data[,3:4])
mg <- ifelse(data$group == "A", 1, 2)
matplot(x, t_data, type = "l", lty = 1, pch = 1, col = mg, xlab = "", ylab = "",xaxt="n", xlim = c(0.9,2.1))
name<-c("前","後")
axis(side=1,at=c(1, 2),labels=name)  #指定した場所に挿入
legend(1.8, 15, c("治療A","治療B"), col = c(1, 2), lty = 1)

f:id:yoshida931:20190711140353p:plain:w500

makedummiesを使ったダミー変数の作成 

2018-07-19投稿、最終更新日2019-07-07

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

治療 効果
NO  NO
YES 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 YES
YES NO
NO  NO
NO  NO
YES YES
NO  NO
YES YES
YES YES
YES YES
NO  YES
NO  NO
YES YES
NO  NO
NO  NO
YES NO
NO  NO
NO  YES
NO  NO
NO  YES
YES NO
YES YES

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

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

summary(s_data)
  治療     効果   
 NO :18   NO :17  
 YES:17   YES:18  

makedummiesをインストールします

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

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

d_data <- makedummies(s_data, basal_level = TRUE)

   治療_NO 治療_YES 効果_NO 効果_YES
1        1        0       1        0
2        0        1       0        1
3        0        1       1        0
4        1        0       0        1
5        0        1       0        1

colnames(d_data) <- c("治療なし","治療あり","効果なし","効果あり")

   効果なし 効果あり 治療なし 治療あり
1         1        0        1        0
2         0        1        0        1
3         0        1        1        0
4         1        0        0        1
5         0        1        0        1
#1列目と3列目を使用した場合の分割表を作成
#行=治療(0=あり、1=なし)
#列=効果(0=あり、1=なし)

d_data2 <- xtabs(~d_data[,2]+d_data[,3])

           d_data[, 3]
d_data[, 1]  0  1
          0 12  5
          1  6 12

name <- list(治療=c("あり","なし"), 効果=c("あり","なし"))
dimnames(d_data2)<-name
d_data2

      効果
治療   あり なし
  あり   12    5
  なし    6   12

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

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

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

     あり なし    P(あり) 95% conf. interval
あり   12    5     0.7059    0.4581   0.8720
なし    6   12     0.3333    0.1580   0.5712

                                   95% conf. interval
             Relative Risk: 2.1176    1.0289   4.3584
         Sample Odds Ratio: 4.8000    1.1471  20.0850
Conditional MLE Odds Ratio: 4.5683    0.9465  25.7201
    Probability difference: 0.3725    0.0427   0.6073

             Exact P-value: 0.0437 
        Asymptotic P-value: 0.0317 
------------------------------------------------------

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

オッズ比の信頼区間が1をまたいでいるのでパラメーターは有意になるはず!
ロジスティック回帰で確認してみます

fit <- glm(効果あり~治療あり, data = d_data, family = binomial)
summary(fit)

Call:
glm(formula = 効果あり ~ 治療あり, family = binomial, data = d_data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4823  -0.8346  -0.8346   0.9005   1.5645  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -0.8755     0.5323  -1.645   0.1000  
治療あり      1.5686     0.7303   2.148   0.0317 *
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 48.492  on 34  degrees of freedom
Residual deviance: 43.512  on 33  degrees of freedom
AIC: 47.512

Number of Fisher Scoring iterations: 4

並べ替え検定と正規近似

2017-07-11投稿, 2019.7.4更新

正確なp値 ( exact\ p\ value )

特定の確率分布をもとに推定を行うのではなく、母集団の未知のパラメータやサンプリング誤差が入らないため計算上も正しいp値が得られる.

並べ替え検定( permutation\ test )

例)x群とy群を比較します

x <- c(5, 9)           #平均 = 7  
y <- c(6, 12, 14, 16)  #平均 = 12  

xとyは同じ母集団からのサンプリングと考えます.それぞれのグループへの割付の際にたまたま差が生じました.なおxとyは正規分布には従いません.

帰無仮説H_0:xとyに差はない
対立仮説H_1:xよりyが大きい

並べ替え検定の考え方

もしH_0が正しいと考えるとき、このサンプリングの平均差がどの程度大きいのかを考えます.(  5 , 6 , 9 , 12 , 14 , 16 )  がマークされている、6個の同質の玉が袋に入っていると考えます.
f:id:yoshida931:20171019175236j:plain:w350

帰無仮説は、どのように取り出してもx(2個)とy(4個)が示す増加量が等しいということになります.

取り出し方は、_6C_2

choose (6, 2)  

の15通りあります.帰無仮説が正しいとすると、15通り全てが「xとyは等しい」ということになります.

xの平均は7、yの平均は12なので、y-x=5となります.つまり15通り中、差が平均の差5より大きくなる確率を正確なp値として考えます.
f:id:yoshida931:20171019175722j:plain:w600

差が5以上になるのは2通りなので、片側検定のp値は2/15 =0.1333333 、両側検定のp値は 4/15 = 0.2666667 となります.

Rでは以下のような計算式になります
Package ‘coin’, March 8, 2019, Version 1.3-0, Date 2019-03-04より
https://cran.r-project.org/web/packages/coin/coin.pdf

#サンプルデータセット  
diffusion <- data.frame(
  pd = c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46,
         1.15, 0.88, 0.90, 0.74, 1.21),
  age = factor(rep(c("At term", "12-26 Weeks"), c(10, 5)))
)
# ソート
(sortlist <- order(diffusion$pd))
( diffusion <- diffusion[sortlist,] )
diffusion$no <- c(1:15);diffusion
#セット
> diffusion
     pd         age no
9  0.73     At term  1
14 0.74 12-26 Weeks  2
1  0.80     At term  3
2  0.83     At term  4
12 0.88 12-26 Weeks  5
13 0.90 12-26 Weeks  6
4  1.04     At term  7
11 1.15 12-26 Weeks  8
15 1.21 12-26 Weeks  9
6  1.38     At term 10
5  1.45     At term 11
10 1.46     At term 12
8  1.64     At term 13
3  1.89     At term 14
7  1.91     At term 15

Rのパッケージ"coin"を使って検定します

#install.packages("coin")
library(coin)
#正確なP値を求めます:並び替え検定、Exact Wilcoxon-Mann-Whitney test  
(wt <- wilcox_test(pd ~ age, data = diffusion, distribution = "exact", conf.int = TRUE))

    Exact Wilcoxon-Mann-Whitney Test

data:  pd by age (12-26 Weeks, At term)
Z = -1.2247, p-value = 0.2544
alternative hypothesis: true mu is not equal to 0
95 percent confidence interval:
 -0.76  0.15
sample estimates:
difference in location 
                -0.305 

期待値や分散など

statistic(wt, type = "linear")# sum of the ranks for age = "12-26 Weeks"
12-26 Weeks 30

expectation(wt)# 平均順位の差の期待値
12-26 Weeks        40 

covariance(wt)# 平均順位の差の分散
            12-26 Weeks
12-26 Weeks    66.66667

pvalue(wt)
[1] 0.2544123

confint(wt)
95 percent confidence interval:
 -0.76  0.15 
sample estimates:
difference in location 
                -0.305 

Mann-Whitney の U 検定, wilcoxonの順位和検定(Wilcoxon rank-sum test)
正規近似を行って検定する方法.なお,この正規近似は m,n が 7 より大きければかなり正確であることも示されている.
ノンパラメトリック検定

(wt2 <- wilcox_test(pd ~ age, data = diffusion,  conf.int = TRUE))

    Asymptotic Wilcoxon-Mann-Whitney Test

data:  pd by age (12-26 Weeks, At term)
Z = -1.2247, p-value = 0.2207
alternative hypothesis: true mu is not equal to 0
95 percent confidence interval:
 -0.7599742  0.1499702
sample estimates:
difference in location 
            -0.3038399 

コンピュータの発達により、並べ替え検定や正確確率検定を行うのは無理ではなくなったので、わざわざp値の近似値を求める従来のパラメトリック検定やノンパラメトリック検定よりも、直観的かつ、わかりやすい結果が得られるといえるだろう. 水本 篤(2010)より引用

参考
柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010,p162-164
水本 篤:統計数理研究所共同研究リポート 238『言語コーパス分析における数理データの統計的処理手法の検討』(2010) pp. 1–14

Rで簡単(二元配置分散分析:データには対応なし、繰り返し数は5)

投稿2016.11.17 最終更新日2019.6.5

検定は帰無仮説が大切です!
商品のの主効果:帰無仮説カップ麺でもインスタント麺でも評価の母平均は等しい
スープの主効果:帰無仮説=スープが違っても評価の母平均は等しい
商品とスープの主効果:帰無仮説=商品とスープの組合せと評価は関係ない

商品 スープ 評価
カップ とんこつ 10
カップ とんこつ 11
カップ とんこつ 11
カップ とんこつ 9
カップ とんこつ 9
カップ 醤油 11
カップ 醤油 8
カップ 醤油 10
カップ 醤油 8
カップ 醤油 8
カップ みそ 7
カップ みそ 3
カップ みそ 5
カップ みそ 2
カップ みそ 3
インスタント とんこつ 11
インスタント とんこつ 11
インスタント とんこつ 10
インスタント とんこつ 9
インスタント とんこつ 10
インスタント 醤油 9
インスタント 醤油 7
インスタント 醤油 9
インスタント 醤油 8
インスタント 醤油 8
インスタント みそ 4
インスタント みそ 4
インスタント みそ 2
インスタント みそ 3
インスタント みそ 3

コピーしてRに読み込みます

x<-read.table("clipboard",header=T)
summary(aov(x$評価~x$商品*x$スープ))

                Df Sum Sq Mean Sq F value   Pr(>F)    
x$商品           1   1.63    1.63   1.077    0.310    
x$スープ         2 231.67  115.83  76.374 3.93e-11 ***
x$商品:x$スープ  2   1.67    0.83   0.549    0.584    
Residuals       24  36.40    1.52                     
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

スープによる主効果あり!
次の式でも同じ結果になります

summary(aov(x$評価~x$商品+x$スープ+x$商品:x$スープ))

ベクトルを使用した方法

評価<-c(x$評価)
商品<-factor(c(rep("カップ",15),rep("インスタント",15)))
スープ<-factor(rep(c(rep("とんこつ",5),rep("醤油",5),rep("みそ",5)),2))
summary(aov(x$評価~x$商品*x$スープ))

検証するために図を描きましょう!

par(mfrow=c(1,2)) 
interaction.plot(商品,スープ,評価)
interaction.plot(スープ,商品,評価)
par(mfrow=c(1,1)) 

f:id:yoshida931:20190606143049p:plain

 

 

相関行列と共分散行列

相関行列

x1 <- c(93, 89, 115, 90, 75)
x2 <- c(121, 136, 121, 161, 125)
x3 <- c(101, 115, 118, 122, 106)

#データフレームにすること
dat <- data.frame(x1, x2, x3)       
 
cor(dat)

           x1         x2        x3
x1  1.0000000 -0.2220449 0.3981155
x2 -0.2220449  1.0000000 0.6767614
x3  0.3981155  0.6767614 1.0000000

散布図

pairs(dat)

f:id:yoshida931:20190604112501p:plain:w400

共分散行列

これもデータフレームで    
var(dat)

       x1     x2    x3
x1 207.80 -54.15 49.80
x2 -54.15 286.20 99.35
x3  49.80  99.35 75.30

参考

install.packages("psych")
library(psych)
x <- rnorm(100,5,3)
y <- rnorm(100,15,8)
z <- rnorm(100,2,0.5)
dat2 <- data.frame(x,y,z)
pairs.panels(dat2)

f:id:yoshida931:20190604113119p:plain:w400

対応のある一元配置分散分析(繰り返しのない二元配置分散分析)、多重比較

投稿日2018.7.21, 最終更新日2019.6.4

以下のデータセットを使用します

x1 <- c(93, 89, 115, 90, 75)
x2 <- c(121, 136, 121, 161, 125)
x3 <- c(101, 115, 118, 122, 106)

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

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

times <- c(rep("1回目",5),rep("2回目",5),rep("3回目",5))
ID <- c(rep(LETTERS[1:5],3))
data <- c(x1, x2, x3)

イメージ

dataf <- data.frame(ID,times,data)

   ID times data
1   A 1回目   93
2   B 1回目   89
3   C 1回目  115
4   D 1回目   90
5   E 1回目   75
6   A 2回目  121
7   B 2回目  136
8   C 2回目  121
9   D 2回目  161
10  E 2回目  125
11  A 3回目  101
12  B 3回目  115
13  C 3回目  118
14  D 3回目  122
15  E 3回目  106

結果

summary(aov(data~times+ID))

            Df Sum Sq Mean Sq F value Pr(>F)
times        2  792.1   396.1   2.473  0.146       #反復による変動
ID           4 1040.4   260.1   1.624  0.259        #患者さんによる変動
Residuals    8 1281.2   160.2                         #水準内変動

有意差があったので多重比較を実行します
テューキー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 comparisons using t tests with pooled SD 
data:  dataf$data and dataf$times 

      1回目  2回目 
2回目 0.0017 -     
3回目 0.1216 0.1119

P value adjustment method: bonferroni 

テューキー(Tukey's multiple comparison test)

TukeyHSD (aov(dataf$data ~ dataf$times), ordered = TRUE)
  Tukey multiple comparisons of means
    95% family-wise confidence level
    factor levels have been ordered

Fit: aov(formula = dataf$data ~ dataf$times)

$`dataf$times`
            diff       lwr      upr     p adj
3回目-1回目 20.0 -3.243611 43.24361 0.0950704
2回目-1回目 40.4 17.156389 63.64361 0.0015343
2回目-3回目 20.4 -2.843611 43.64361 0.0879623

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

#注意事項
相関係数より
           x1         x2        x3
x1  1.0000000 -0.2220449 0.3981155
x2 -0.2220449  1.0000000 0.6767614
x3  0.3981155  0.6767614 1.0000000

x2, x3に相関がみられます。反復測定なのでよくあることです。分散分析の前提はiidですので、
x2とx3が「独立しているか?」という点で疑問があります.
このような場合には、混合効果モデルなどを使用した解析が必要になります.

散布図に直線を描く

f:id:yoshida931:20190524000831p:plain:w400

x <- rnorm(20)  #乱数を使用してるので毎回変わります
y <- c(1:20)
plot(x,y)

xs <- seq(min(x), max(x), length=1000)

# xの影響を含むモデル
fit <- glm(y~x,family = poisson)
lines(xs, exp(fit$coef[1] + xs*fit$coef[2]))

# 切片のみのモデル
fit.null <- glm(y~1,family = poisson)
lines(xs, exp(fit.null$coef[1] + xs*0))