1 一標本 t 検定

このページは、市原清志先生の「バイオサイエンスの統計学」の問題を、統計ソフト”R”で解説するものです。詳細については、市原先生のテキスト（p28〜）を参照ください。

1.1 ［基本］

n組のデータにつき各ペアの差 \(d\) を求め、その平均値 \(\overline{d}\) を統計量とする。
＜帰無仮説＞
　　　\(d\sim N(0\:,\:\sigma_d^2)\)（各ペアに差がない）　とすると　\(\overline{d}\sim N\left(0\:,\:\dfrac{\sigma_d^2}{n}\right)\)　が成立。

差の母分散 \(\sigma_d^2\) が既知の場合は、 \(z=\dfrac{d}{\sigma_d/\sqrt{n}} \sim N(0,1)\)　によるZ検定となるが、
母分散が未知のため標本分散 \(s_d^2\) で代用し、\(\overline{d}\) を標準誤差 \(\dfrac{s_d}{\sqrt{n}}\) で割ると
　　　\(t=\dfrac{d}{s_d/\sqrt{n}} \sim t(n-1)\)：自由度 \((n-1)\) の t分布に従う。

1.2 ［例2-1］バセドウ病患者（p28）

6人のバセドウ病患者に薬を投与したとき、投与の前後で脈拍数に変化があるかどうか検定する。

pre<-c(98,88,100,96,107,114)
post<-c(86 ,73 ,95 ,92 ,99 ,116)

変化を見るにはRの個別推移図（スパゲティ・プロット）が分かりやすい。
preは前、postは後を表す。
xyplot を使うため、データを縦長形式 (long style) のデータフレームにする。

# 個別推移図（スパゲティ・プロット）
library("lattice");library("latticeExtra")#パッケージの読み込み
score <- c(pre, post)
prepost <- factor(c(rep("pre", 6), rep("post", 6)), levels = c("pre","post"))
df <- data.frame(score, prepost)
df$indiv <- factor(c(rep(1:6), rep(1:6)))#個体の識別番号
df

##    score prepost indiv
## 1     98     pre     1
## 2     88     pre     2
## 3    100     pre     3
## 4     96     pre     4
## 5    107     pre     5
## 6    114     pre     6
## 7     86    post     1
## 8     73    post     2
## 9     95    post     3
## 10    92    post     4
## 11    99    post     5
## 12   116    post     6

識別番号ごとに、変化を折れ線で示す。

xyplot(score ~ prepost, group = indiv, type = c("l"), data = df, xlab = "pre-post", ylab = "score")

薬の投与後、血圧は下がっているように見える。この変化が有意なものかどうかRで調べる。

帰無仮説「preとpostで変化がない」のもとで、
\(~~~~~d=pre-post\) は、\(d\sim N(0\:,\:\sigma_d^2)\)：0を平均とした正規分布に従い、
標本標準偏差（標準誤差）で標準化した \(~~~t=\dfrac{d}{s_d/\sqrt{6}} \sim t(5)\) 　は自由度5の t 分布に従う。
自由度 5の t 分布とは、このような調査を無限回行った場合の t値の分布であり、０を中心としている。

par(family = "HiraKakuProN-W3") 
curve(dnorm(x, 0, 1),xlim=c(-5, 5),ylim=c(0, 0.4), col=2,xaxp=c(-6,7,13))
curve(dt(x, 5),xlim=c(-5, 5),col=3, add=T) 
legend( "topleft",legend=c("標準正規分布","t分布"), col=c("red","green"),lty=c(1,1))

t検定の方が有意差が出にくい。
両側t検定のrコードは、t.testで”two.sided”をオプションとすればよい。

pre<-c(98,88,100,96,107,114);post<-c(86 ,73 ,95 ,92 ,99 ,116)
res<- t.test(pre,post,paired=TRUE,alternative ="two.sided")
res

## 
##  Paired t-test
## 
## data:  pre and post
## t = 2.8265, df = 5, p-value = 0.03683
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   0.6338082 13.3661918
## sample estimates:
## mean of the differences 
##                       7

t 値は、t = 2.8265 で 0 から離れているように思われる。
t 値が、この値かそれ以上に0から離れる確率は、p-value = 0.03683
この確率は基準とする0.05より小さく、t = 2.8265は異常値と考えられる。
よって、
帰無仮説「dの平均が 0 」がそもそも間違いであり、
対立仮説「 dの平均が 0 ではない」を真と判断する。

さらに、差の平均 d は95%の確率で区間（0.63 ,13.36）に入る。

今回の調査での差の平均 7 は、ほぼこの中央に位置している。

［注意］
　t 値や確率の計算はRが自動で行ってくれるが、その意味を理解し、結果の見方を知っておくことが重要である。

統計学備忘録様を参考にさせていただきました。

2 符号検定(sign test)

市原先生のテキスト（p42〜）を参照。

2.1 ［基本］

pairの差 \(X_1,X_2,\cdots ,X_n\) の符号の個数で少ないほうの個数を検定量 \(r\) とする。
　　　\(r=X_1,X_2,\cdots ,X_n\) のうち正の値の個数（負の個数が少なければ負の個数）

＜帰無仮説＞　「変化なし」
のもとでは、符号は均等に現れるので、\(r\sim B\left(n,\dfrac{1}{2}\right)\)　（二項分布に従う）

nが小さければ、直接２項確率を計算する。
大標本の場合は正規分布によって近似する　\(r\sim N\left(\dfrac{n}{2},\dfrac{n}{4}\right)\)
（連続補正をおこなう）

2.2 ［例題2-2］昇圧ホルモン値（p44）

2.2.1 データ

xyplot を使うため、データを縦長形式 (long style) のデータフレームにする。

#例2_2<-read.csv("例2_2.csv")
#lay<-例2_2[,1]#臥位
lay <- c(0.1,0.3 ,0.1, 1.7 ,0.4, 0.3, 0.2 ,0.8)
#stand<-例2_2[,2]#立位
stand <- c( 0.2, 0.5 ,0.8 ,2.3, 0.1, 0.7, 0.8 ,1.3)
score <- c(lay, stand)
posision <- factor(c(rep("lay", 8), rep("stand", 8)), levels = c("lay","stand"))
df <- data.frame(score, posision)
df$indiv <- factor(c(rep(1:8), rep(1:8)))#個体の識別番号
df

##    score posision indiv
## 1    0.1      lay     1
## 2    0.3      lay     2
## 3    0.1      lay     3
## 4    1.7      lay     4
## 5    0.4      lay     5
## 6    0.3      lay     6
## 7    0.2      lay     7
## 8    0.8      lay     8
## 9    0.2    stand     1
## 10   0.5    stand     2
## 11   0.8    stand     3
## 12   2.3    stand     4
## 13   0.1    stand     5
## 14   0.7    stand     6
## 15   0.8    stand     7
## 16   1.3    stand     8

各個体の変化をRの個別推移図で表示する。

each <- xyplot(score ~ posision, group = indiv, type = c("l"), data = df, xlab = "hormone", ylab = "score")
each

グラフで見ると、8人のうち1人だけホルモン値が下降している。
　　　 (検定量 r)＝(立位)ー(臥位) の符号の個数（少ないほう）とすると、\(r=1\)
帰無仮説「layとstandで差がない」のもとでは、\(r\sim B\left(8,\dfrac{1}{2}\right)\).

rの期待値は \(E(r)=4\) であるから、4前後が起こりやすいので、
　　　\(r=1\) は異常値と予想される。
\(P(r\leq 1)\) を計算し、有意水準より小さければ異常値と考えるのが片側検定の考え方である。

小標本なので直接二項確率を計算する。
　　　\({}_8 \rm{C}_1 (0.5)^7(0.5)^1+ {}_8 \rm{C}_0 (0.5)^8=0.035<0.5\)
有意水準より小さいので、片側検定で有意である。
また、 R では “binom.test” で処理することができる。
https://multivariate-statistics.com/2022/07/28/r-programming-binomial-test/ https://ukigorilla.com/post-1224/

Rを用いた検定ではbinom.test を用いる。

library(binom)
binom.test(1, 8, 0.5, alternative ="less")

## 
##  Exact binomial test
## 
## data:  1 and 8
## number of successes = 1, number of trials = 8, p-value = 0.03516
## alternative hypothesis: true probability of success is less than 0.5
## 95 percent confidence interval:
##  0.0000000 0.4706794
## sample estimates:
## probability of success 
##                  0.125

１回の試行で起こる確率が 0.5 の事象があり、試行を８回行うとき、
起きる回数が１回以下である確率が、0.03516 であることが読み取れる。
回数が多ければ正規分布に近似して信頼区間を求めるが、８回では正規近似はできない。binom.confint関数で様々な算出法での信頼区間が出力される。

binom.confint(1, 8)

## Warning in stats::prop.test(x[i], n[i]): Chi-squared approximation may be
## incorrect

##           method x n      mean         lower     upper
## 1  agresti-coull 1 8 0.1250000  0.0011332625 0.4921724
## 2     asymptotic 1 8 0.1250000 -0.1041723038 0.3541723
## 3          bayes 1 8 0.1666667  0.0005567934 0.3970442
## 4        cloglog 1 8 0.1250000  0.0065933250 0.4227058
## 5          exact 1 8 0.1250000  0.0031597235 0.5265097
## 6          logit 1 8 0.1250000  0.0172727791 0.5372756
## 7         probit 1 8 0.1250000  0.0117985066 0.4852142
## 8        profile 1 8 0.1250000  0.0121249351 0.4453946
## 9            lrt 1 8 0.1250000  0.0075891287 0.4453568
## 10     prop.test 1 8 0.1250000  0.0065599169 0.5332108
## 11        wilson 1 8 0.1250000  0.0224174915 0.4708882

数学的な導出については杉本典夫先生のサイトに詳しいが、難しそうなので立ち入らない。

2.3 ［例題2-3］（p45）

大標本なのでZ検定を行う。（連続補正あり）

3 ウイルコクスン符号順位検定

　　Wilcoxon signed-rank test　（市原p50）
正規分布を前提としたパラメトリックな方法で関連２群の差の検定を行うには、各ペアの差 \(d\) が正規分布に従う必要があります。
\(x_1,~x_2\) が独立に正規分布に従うならば、\(d=x_1-x_2\) も正規分布に従いますが、独立でないときは、\(d=x_1-x_2\)　が正規分布に従う保証はありません。
正規性の検定を別に行うこともできますが、このような場合は正規性を前提としないノンパラメトリックな検定がよく用いられています。
前述の符合検定もノンパラメトリックであり、ウイルコクスンもノンパラメトリックな検定方法です。ノンパラメトリックな方法でも、標本数が大きければ正規近似できる統計量を構成することができると言われています。
数理については当サイトWilcoxonで説明していますが、ここではシミュレーションで確かめてみましょう。

3.1 シミュレーション

例2-4を参考に8人の前後のデータを作り、差の符合別順位和stを計算します。

\(st\) の平均と分散は　 \(E(st)=\dfrac{n(n+1)}{4}\) 　 \(V(st)=\dfrac{n(n+1)(2n+1)}{24}\)　であることより、
　　　\(E(st)=\dfrac{8\cdot 9}{4}=18\) 　　 \(V(st)=\dfrac{8\cdot 9\cdot17}{24}=51\)

この操作を1000回行い st がどのような分布をするか確かめましょう。
また、理論曲線と比べてみましょう。

St = c()
X=c()
for (i in 1:1000) {
  x1<-rnorm(8,165,12);x2<-rnorm(8,165,15)
  x <- sign(x1-x2)*rank(abs(x1-x2))
  st <- sum(x[x>0])
  St <- c(St,st)
}
mean(St)

## [1] 18.082

hist(St,prob=T)
curve(dnorm(x,18,sqrt(51)),add=TRUE,col="red")

n=8 でも近似が良いことに驚かされます。

3.2 ［例題2-4］（p52）

https://data-science.gr.jp/implementation/ist_r_wilcoxon_signed_rank_test.html

前<-c(182,169,173,143,158,156,176,165)
後<-c(163,142,174,137,151,143,180,162)
x1<-前;x2<-後
x <- sign(x1-x2)*rank(abs(x1-x2))
st <- sum(x[x>0])
st

## [1] 32

ヒストグラムで st>32 である確率は小さいようですが、両側検定で判断すると

library(exactRankTests)

##  Package 'exactRankTests' is no longer under development.
##  Please consider using package 'coin' instead.

wilcox.exact(x=前,y=後,paired=T)

## 
##  Exact Wilcoxon signed rank test
## 
## data:  前 and 後
## V = 32, p-value = 0.05469
## alternative hypothesis: true mu is not equal to 0

ギリギリ有意にはならないようです。

p53参照