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

統計学備忘録（R言語のメモ）

since2016　ときどきTEXのメモ

PDFにページ数を挿入する場合

TeX

\LaTeX

PDFにページ数を挿入する場合
PDFのファイル名はNo1,No2,No3の三枚
同じフォルダに入れる

\documentclass[uplatex,11pt]{jsarticle}
\usepackage[top=15truemm,bottom=15truemm,left=15truemm,right=15truemm]{geometry}     %余白

%図
\usepackage[dvipdfmx]{graphicx}
%図の位置調整　改行後に[H]
\usepackage{here}  

% 章が進むごとに図番号をリセットする
\makeatletter% プリアンブルで定義開始
\renewcommand{\thefigure}{\arabic{figure}}
\@addtoreset{figure}{section}
\makeatother % プリアンブルで定義終了

%図とキャプションを狭くする
\setlength\abovecaptionskip{5pt}

%---------------------------------------------------%
\begin{document} 
%---------- 1 ----------%
\begin{figure}[H] %現在の位置に
  \centering
    \includegraphics[width=17cm]{No1.pdf} %図の挿入
\end{figure}
%---------- 2 ----------%
\newpage
\begin{figure}[H] %現在の位置に
  \centering
    \includegraphics[width=17cm]{No2.pdf} %図の挿入
\end{figure}
%---------- 3 ----------%
\newpage
\begin{figure}[H] %現在の位置に
  \centering
    \includegraphics[width=17cm]{No3.pdf} %図の挿入
\end{figure}

\end{document}

並べ替え検定と正規近似

ﾉﾝﾊﾟﾗﾒﾄﾘｯｸ

2017-07-11投稿, 2019.7.4更新

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

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

並べ替え検定( permutation\ test )

例）ｘ群とｙ群を比較します

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

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

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

並べ替え検定の考え方

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

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

取り出し方は、_6C_2

choose (6, 2)

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

ｘの平均は7、ｙの平均は12なので、ｙ－ｘ＝5となります．つまり15通り中、差が平均の差5より大きくなる確率を正確なｐ値として考えます．
f:id:yoshida931:20171019175722j:plain:w600

差が5以上になるのは2通りなので、片側検定のｐ値は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

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