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

統計学備忘録 since2016

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

勉強会資料

サンプルをコピーします

グループ ROMT
A   18
A   8
B   23
B   27
A   14
B   15
A   15
B   22
B   15
A   16
A   13
A   23
B   20
A   20
B   20
B   20
B   18
B   25
A   14
A   20

#コピーしているデータをRに読み込ませます
data01 <- read.table("clipboard",header = T)

これでdata01の中にコピーしたデータが入りました
入ってるかどうか確認します

data01

次のサンプルをdata02の中に読み込ませます

一次判定 歩行
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できる
非該当   できない
非該当   できない
非該当   できない
非該当   できない
非該当   できない
非該当   できない
非該当   できない
非該当   できない
非該当   できない
非該当   できない
要支援   できる
要支援   できる
要支援   できる
要支援   できる
要支援   できる
要支援   できる
要支援   できる
要支援   できる
要支援   できる
要支援   できる
要支援   できない
要支援   できない
要支援   できない
要支援   できない
要支援   できない
要支援   できない
要支援   できない
要支援   できない
要支援   できない
要支援   できない
要支援   できない
要支援   できない
要支援   できない
要支援   できない
要支援   できない

#コピーしているデータをRに読み込ませます
data02 <- read.table("clipboard",header = T)

これでdata02の中にコピーしたデータが入りました
入ってるかどうか確認します

data02

エラーバー付きのグラフ ggplot2

使用するパッケージ

install.packages("ggplot2")
library(ggplot2)

データセットを作成します

#乱数を使用のため、実施毎にグラフは異なります
dat <- c(c(rnorm(15)),c(rnorm(15,3,1.2)),c(rnorm(15,10,2)),c(rnorm(15,17,1.5)))
pre_post <- c(rep("前", 30),rep("後", 30))
treat <- c(rep("A",15),rep("B",15),rep("A",15),rep("B",15))
da01 <- data.frame(dat, pre_post, treat);head(dataf)

            dat pre_post treat
1  0.0004020798       前     A
2  0.3236960969       前     A
3  1.0026311892       前     A
4  0.3526106439       前     A
5 -0.1307876208       前     A
6 -1.5476749466       前     A

列pre_postの配列を確認します

levels(da01$pre_post)
[1] "後" "前"

このままではx軸が後→前の配列になるので入れ替えます

da01 <- transform(dataf, pre_post= factor(pre_post, levels = c("前", "後")))
levels(da01$pre_post)
[1] "前" "後"

二群比較

ggplot(da01, aes(x = pre_post, y = dat)) + geom_point()+theme_bw()

f:id:yoshida931:20190825213846p:plain:w400

標準誤差をエラーバーとしたグラフを描きます

#position_jitterdodgeで図の重なり具合を調整しておきます
pos_1 <- position_jitterdodge(
  jitter.width  = 0,#横にずらす
  jitter.height = 0,#縦にずらす
  dodge.width   = 0.5#カテゴリで分ける
)

(g <- ggplot(data = da02,
       aes(
         x      = pre_post,
         y      = dat,
         colour = treat,
       )) +
  geom_jitter(alpha = 0.4, position = pos_1) +  #alpha重なりの濃さ
  stat_summary(aes(x = as.numeric(pre_post) + 0.07), fun.y = "mean", geom = "point", size = 3, position = pos_1)+#平均値のプロット
  stat_summary(aes(x = as.numeric(pre_post) + 0.07),fun.data = "mean_se",#mean_seで標準誤差、mean_cl_normalで95%信頼区間
                  # 標準偏差   fun.ymin = function(x) mean(x) - sd(x),fun.ymax = function(x) mean(x) + sd(x),
               geom = "errorbar",
               width = 0.1,#ひげの横幅
               lwd = 1,#ひげの棒の幅
               #fun.args = list(conf.int = 0.95),#信頼区間のときに使用する
               position = pos_1
  ) +
  theme_bw()) #背景の設定

f:id:yoshida931:20190825195621p:plain:w400

単回帰分析

fitA <- lm(dat ~ pre_post , data =da01[da01$treat=="A",])
summary(fitA)

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.2681     0.5873  -0.456    0.652    
pre_post後    4.0898     0.8306   4.924 3.42e-05 ***

fitB <- lm(dat ~ pre_post , data =da01[da01$treat=="B",])
summary(fitB) 

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.0459     0.6852   4.445 0.000126 ***
pre_post後    1.0633     0.9690   1.097 0.281873  

単回帰の結果を書き込んでみます

(g2 <- g + 
geom_signif(stat = "identity",
            data = data.frame(x = c(0.87, 1.12),
                              xend = c(1.87, 2.12),
                              y = c(9, 12),
                              annotation = c("**", "N.S.")),
            aes(x = x, xend = xend, y = y, yend = y, annotation = annotation),
            col = c("red","#20B2AA")) ) +
  ylim (-5,13)

f:id:yoshida931:20190825223325p:plain:w400

重回帰分析

fit <- lm(dat ~ pre_post + treat + pre_post * treat, data =da01)
summary(fit)

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -0.2681     0.6381  -0.420 0.676005    
pre_post後          4.0898     0.9025   4.532 3.12e-05 ***
treatB              3.3140     0.9025   3.672 0.000539 ***
pre_post後:treatB  -3.0265     1.2763  -2.371 0.021189 *  

#分散分析表
anova(fit)

Response: dat
               Df Sum Sq Mean Sq F value    Pr(>F)    
pre_post        1  99.58  99.579 16.3016 0.0001658 ***
treat           1  48.64  48.640  7.9627 0.0065965 ** 
pre_post:treat  1  34.35  34.349  5.6231 0.0211892 *  
Residuals      56 342.08   6.109  

この記事で使用したデータフレーム
コピーしてRにペースト
da01 <- read.table("clipboard", header=T)
これで再現できます

> da01
    X         dat pre_post treat
1   1 -0.29733004       前     A
2   2 -2.63812280       前     A
3   3 -0.90097072       前     A
4   4  1.06843016       前     A
5   5 -3.03846846       前     A
6   6  2.14097694       前     A
7   7  2.47494865       前     A
8   8 -0.02154341       前     A
9   9  0.56411223       前     A
10 10  2.81826063       前     A
11 11 -1.22557353       前     A
12 12 -1.96770187       前     A
13 13 -2.41952982       前     A
14 14 -2.04788520       前     A
15 15 -1.49396856       前     A
16 16  4.42978119       前     B
17 17  1.14121396       前     B
18 18  2.81408173       前     B
19 19  3.88617124       前     B
20 20  0.27836700       前     B
21 21  7.04775580       前     B
22 22  6.76956940       前     B
23 23 -1.01243943       前     B
24 24  2.28982618       前     B
25 25  2.56697055       前     B
26 26  3.92271476       前     B
27 27  0.56754099       前     B
28 28  2.55950816       前     B
29 29  5.45781910       前     B
30 30  5.97792826       前     B
31 31  3.57321119       後     A
32 32  3.84212737       後     A
33 33  1.38192000       後     A
34 34  1.89198329       後     A
35 35  5.90437940       後     A
36 36  3.16262753       後     A
37 37  1.29383105       後     A
38 38  7.08585415       後     A
39 39  1.53787689       後     A
40 40  5.82484000       後     A
41 41  4.45300261       後     A
42 42  3.94823267       後     A
43 43  2.22314781       後     A
44 44  5.21597500       後     A
45 45 -1.10466187       後     A
46 46  1.03430878       後     B
47 47  3.26798971       後     B
48 48  4.46817078       後     B
49 49  4.30615587       後     B
50 50  3.13762549       後     B
51 51  3.12041534       後     B
52 52  5.40843117       後     B
53 53  2.01223446       後     B
54 54 -2.57111800       後     B
55 55  4.01156402       後     B
56 56  6.58882849       後     B
57 57  6.94067008       後     B
58 58  3.12482569       後     B
59 59  3.81695957       後     B
60 60  8.53560520       後     B

独立している2群を比較するグラフ ggplot2

データの準備

dat <- c(c(rnorm(15)),c(rnorm(15,3,1.2)),c(rnorm(15,10,2)),c(rnorm(15,17,1.5)))
pre_post <- c(rep("前", 30),rep("後", 30))
treat <- c(rep("A",15),rep("B",15),rep("A",15),rep("B",15))
dataf <- data.frame(dat, pre_post, treat)
head(dataf)
          dat pre_post treat
1 -0.27386197       前     A
2 -0.03189246       前     A
3 -0.33039790       前     A
4 -0.82833272       前     A
5  0.33606814       前     A
6  0.51831329       前     A

パッケージのインストールを忘れずに

install.packages("ggplot2")
library(ggplot2)

stripchartで簡単に作図してみます

stripchart(dat~pre_post, data = dataf, 
           vertical=TRUE,
           pch=16,
           xlim=c(0.5,2.5))

f:id:yoshida931:20190815173750p:plain:w400

ggplot2で2群を比較するグラフの編集に挑戦してみます

(g1 <- ggplot(dataf, aes(x = pre_post, y = dat)) + geom_point())

f:id:yoshida931:20190815174001p:plain:w400

x軸の「前」「後」の順を変えてみます
前後を数字(0,1)にした列pre_post2を追加します(前=0、後=1で設定します)

# dplyrパッケージのmutateを使用します
library(dplyr)
dataf2 <- mutate(dataf, pre_post2 = if_else(pre_post == "前" , true = 0, false = 1))
head(dataf2, n=3)

          dat pre_post treat pre_post2
1 -0.27386197       前     A         0
2 -0.03189246       前     A         0
3 -0.33039790       前     A         0

#pre_post2という新しい列が追加されました
#データセット名もdataf2に変更しました

reorder で前後の順を変更します

(g2 <- ggplot(dataf2, aes(x = reorder(pre_post, pre_post2), y = dat)) + geom_point())

f:id:yoshida931:20190815174633p:plain:w400

このグラフに回帰直線を描き加えてみます

fit <- lm(dat~pre_post2 , data = dataf2)
(g2 + geom_abline(slope=fit$coef[2], intercept=fit$coef[1]))   

切片がずれる!!!
f:id:yoshida931:20190815174905p:plain:w400

仕方がないので、曝露を数字「0、1」に戻してみます

(g3 <- ggplot(dataf2, aes(x = pre_post2, y = dat)) +
    geom_point() +
    geom_abline(slope=fit$coef[2], intercept=fit$coef[1])
  )    

f:id:yoshida931:20190815175702p:plain:w400

バックの色や軸を編集

(g3 +
    xlim(-0.5,1.5) +
    theme(panel.background = element_blank(),
          axis.line = element_line(colour="black"))
)

f:id:yoshida931:20190815180055p:plain:w400

軸を文字にしたい場合にはパワーポイントで編集するしかないか・・・?

最後にt検定で確認してみます

t.test(dat ~ pre_post2, data = dataf2)

    Welch Two Sample t-test

data:  dat by pre_post2
t = -13.911, df = 45.549, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -12.751126  -9.526776
sample estimates:
mean in group 0 mean in group 1 
       1.899348       13.038299 

fit$coef[1]  #非曝露群の平均
(Intercept) 
   1.899348 
fit$coef[2]+fit$coef[1]  #曝露群の平均
pre_post2 
  13.0383 

Rで作る分割表

投稿日:2018/06/28, 最終更新日2019/08/02
次のデータフレームを分割表にして解析してみます

x <- data.frame(c(rep("あり",15),rep("なし",15)),
           c(rep("あり",9),rep("なし",6),rep("あり",5),rep("なし",10)),
           c(rep("有効",12),rep("無効",3),rep("有効",4),"無効","有効","有効",rep("無効",8)))
colnames(x) <- c("筋トレ","歩行練習","効果")

2変数の分割表の作成

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

      歩行練習
筋トレ あり なし
  あり    9    6
  なし    5   10

効果=有効のみの分割表の作成

x2 <- xtabs(~筋トレ+歩行練習, data=x[x$効果=="有効",])

      歩行練習
筋トレ あり なし
  あり    9    3
  なし    4    2

周辺合計

addmargins(x2)

      効果
      歩行練習
筋トレ あり なし Sum
  あり    9    6  15
  なし    5   10  15
  Sum    14   16  30

列の入れ替え

x3 <- xtabs(~筋トレ+効果, data=x)

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

# 列を入れ替え
x3[c("あり","なし"),c("有効","無効")]   

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

#または行列を直接作ることも可能
xx <- matrix(c(12,6,3,9),2,2)
colnames(xx) <- paste(c("有効","無効"))
rownames(xx) <- paste(c("あり","なし"))

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

転置

x3_2

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

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

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

歩行練習で層別

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 
------------------------------------------------------

ロジスティック回帰で基準を変更する場合

glm(formula = 効果 ~ 筋トレ + 歩行練習, family = binomial, data = x)
Coefficients:
 (Intercept)    筋トレなし  歩行練習なし  
       3.496        -1.738        -3.319  

x$筋トレ2 <- relevel(x$筋トレ2, ref="なし")
glm(formula = 効果 ~ 筋トレ2 + 歩行練習, family = binomial, data = x)
Coefficients:
 (Intercept)    筋トレあり  歩行練習なし  
       1.758         1.738        -3.319  
#xに1列追加されているので注意

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

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