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

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

since2016　ときどきTEXのメモ

LaTeXでプレゼン

プログラム TeX

無料ソフトのみで統計からプレゼンまで！
spss, word, power pointを使用しないで以下のような２枚のスライドを作ってみました．

使ったもの

  • Hatena Blog
  • R
  • LeTeX

f:id:yoshida931:20180206114953p:plain:w500

f:id:yoshida931:20180205173355p:plain:w500
画像は オッズ比の信頼区間 - 統計学備忘録 since2016 からのコピペです

%LaTexの記載は以下の通りです．
%タイプセットはpLaTex(ptex2pdf)を選択します

\documentclass[slide,papersize]{jsarticle}
\usepackage[dvipdfmx]{color}
\usepackage{pxfonts}
\def \sheet #1{
\section*{\centering \large \bfseries #1}
}

\title{\LaTeX!でプレゼンに挑戦}

\usepackage[dvipdfmx]{graphicx}  %図
%http://www.wankuma.com/seminar/20090912nagoya09/5.pdf

\begin{document}

\maketitle

\begin{center}
   \includegraphics[width=60mm,scale=0.8]{WS000004.jpg} %絵も小さくしておきます
\end{center}
  
$p'A = \frac{a}{a+b}, p'B = \frac{c}{c+d}, \frac{\frac{p’A}{1-p’A}}{\frac{p’B}{1-p’B}}=\frac{p’A(1 - p'B)}{p’B(1 - p'A)}=\frac{\frac{a}{b}}{\frac{c}{d}}$

\end{document}

全て、以下のページに書かれてます
http://www.wankuma.com/seminar/20090912nagoya09/5.pdf

オッズ比の信頼区間

クロス表 比率 区間推定

オッズ比、見込み比（odds ratio）または交差積比（cross-product）

前提
比は分母が小さくなると、数値が大きくなりすぎて正規近似の精度が悪くなります．比の対数であれば高い精度で正規近似することが可能になります． したがって、比の対数を考えていくことになります． 

xm <- matrix(c("a","b","c","d"), nrow=2, byrow=T)
name <- list("暴露"=c("あり(A)","なし(B)"),"疾病"=c("あり","なし"))
dimnames(xm)<-name
xm

         疾病
暴露      あり なし
  あり(A) "a"  "b" 
  なし(B) "c"  "d"


p'A = \frac{a}{a+b} 　 p'B = \frac{c}{c+d}

オッズ比＝\frac{\frac{p’A}{1-p’A}}{\frac{p’B}{1-p’B}}=\frac{p’A(1 - p'B)}{p’B(1 - p'A)}=\frac{\frac{a}{b}}{\frac{c}{d}}

暴露なし(A群)が、暴露あり(B群)より疾病に罹患する可能性が大きいか、小さいかを表す尺度．
p'A, p'B が小さい場合には、リスク比と同じ値になります．


オッズ比の信頼区間

比の対数を考えていきます．
標本のオッズ比＝標本オッズ比( \frac{\frac{a}{b}}{\frac{c}{d}} =\frac{ad}{cb}=or)

母集団のオッズ比を母オッズ比(OR)とします．

標本オッズ比の対数 log\ or、　母リスク比の対数 log OR

正規近似で信頼区間を求めていくためには、以下の式が必要になります．

Z = \frac{ log\ or - logOR}{ log\ or の標準誤差}

したがって、以下の式が95％信頼区間を求める式になります

log\ or - 1.96(log\ orの標準誤差) \leqq logOR \leqq log\ or + 1.96(log\ orの標準誤差)

以下、log\ orの標準誤差 = SE とします

log\ or - 1.96SE \leqq logOR \leqq log\ or + 1.96SE

or\times exp(-1.96\times SE) \leqq OR \leqq or\times exp(1.96\times SE)


log\ orの分散

log\ or=log\ a + log\ d - log\ b - log\ c より

V(log\ or)=\frac{1}{a} + \frac{1}{d} + \frac{1}{b} + \frac{1}{c} 　デルタ法(delta method)によって近似的に求めたもの(参考webより)

#Rで95％信頼区間を求めてみます
rr <- a*(c+d)/(a+b)/c      #  = (a/(a+b))/(c/(c+d)) 
exp(log(rr)+c(1, -1)*qnorm ( c(0.025,0.975) )*sqrt(b/a/(a+b)+d/c/(c+d)))

例）
クロス表イメージ
(x <- matrix(c(3, 5, 6, 8 ),2,byrow=T))
     [,1] [,2]
[1,]    3    5
[2,]    6    8

a<-3; b<-5; c<-6; d<-8　     #参考webより
( or <-  a*d/(b*c) )        #オッズ比
[1] 0.8
or*exp(qnorm ( c(0.025, 0.975) )*sqrt( 1/a + 1/b + 1/c + 1/d) )
[1] 0.1348801  4.7449559

#Ｒの関数を使えば瞬間で色々出力してくれます（やっぱり凄い）
install.packages("Epi")  
library(Epi)            
twoby2(x)

2 by 2 table analysis: 
------------------------------------------------------ 
Outcome   : Col 1 
Comparing : Row 1 vs. Row 2 

      Col 1 Col 2    P(Col 1) 95% conf. interval
Row 1     3     5      0.3750    0.1254   0.7152
Row 2     6     8      0.4286    0.2065   0.6837

                                    95% conf. interval
             Relative Risk:  0.8750    0.2972   2.5763 　　　　#リスク比と95％信頼区間
         Sample Odds Ratio:  0.8000    0.1349   4.7450        #オッズ比と95％信頼区間
Conditional MLE Odds Ratio:  0.8081    0.0884   6.3780　#Fisherの正確検定によるオッズ比と95％信頼区間
    Probability difference: -0.0536   -0.3956   0.3312

             Exact P-value: 1 　　　　　　#Fisherの正確検定によるp値
        Asymptotic P-value: 0.8059 　　#正規近似によるp値
------------------------------------------------------


参考web1　統計学入門−第3章
参考web2　R による統計処理
参考文献）柳川 堯 ; 観察データの多変量解析―疫学データの因果分析,近代科学社 ,2016

リスク比の信頼区間

クロス表 区間推定 比率 二項分布

ポイント：比の対数をとり、正規近似する

リスク比は疫学における指標の1つです．一般的には相対危険度（相対リスク，relative risk，RR）として利用されています．

xm <- matrix(c("a","b","c","d"), nrow=2, byrow=T)
name <- list("暴露"=c("あり(A群)","なし(B群)"),"疾病"=c("あり","なし"))
dimnames(xm)<-name
xm

            疾病
暴露         あり なし
  あり(A群)  "a"  "b" 
  なし(B群)  "c"  "d"

前提
比は分母が小さくなると、数値が大きくなりすぎて正規近似の精度が悪くなります．比の対数であれば高い精度で正規近似することが可能になります． したがって、比の対数を考えていくことになります．

定義と言葉の整理

標本比率　\frac{a}{a+b}=p'A の母比率をPA

標本比率　\frac{c}{c+d}=p'B の母比率をPB

log\frac{a}{a+b}=log(p'A),　log\frac{c}{c+d}=log(p'B)


リスク差（過剰絶対リスク） = p'A-p'B

リスク比 RR = \frac{p'A}{p'B}　Ａ群のＢ群における疾病の頻度の比

過剰リスク（過剰相対リスク） = \frac{p'A-p'B}{p'B}　

オッズ比=\frac{\frac{a}{b}}{\frac{c}{d}}


リスク比の信頼区間

比の対数を考えていきます．
標本のリスク比＝標本リスク比( \frac{p'A}{p'B} =rr)、母集団のリスク比を母リスク比(\frac{PA}{PB} =RR)とします．

標本リスク比rrの対数 log\ rr、　母リスク比RRの対数 log RR

正規近似で信頼区間を求めていくためには、以下の式が必要になります．

Z = \frac{ log\ rr - logRR}{ log\ rr の標準誤差}

したがって、以下の式が95％信頼区間を求める式になります

log\ rr - 1.96(log\ rrの標準誤差) \leqq logRR \leqq log\ rr + 1.96(log\ rrの標準誤差)

以下、log\ rrの標準誤差 = SE とします

log\ rr - 1.96SE \leqq logRR \leqq log\ rr + 1.96SE

rr\times exp(-1.96\times SE) \leqq RR \leqq rr\times exp(1.96\times SE)


log\ rrの分散

log\ rr=log\ a - log(a+b) - log\ c + log(c+d) より

V(log\ rr)=\frac{1}{a} - \frac{1}{a+b} + \frac{1}{c} - \frac{1}{c+d} 　デルタ法(delta method)によって近似的に求めたもの

「矢野様より指摘いただき、修正しております」

#Rで95％信頼区間を求めてみます
rr <- a*(c+d)/(a+b)/c      #  = (a/(a+b))/(c/(c+d)) 
exp(log(rr)+qnorm ( c(0.025,0.975) )*sqrt(b/a/(a+b)+d/c/(c+d)))

または
exp(log(rr)+c(-1,1)*qnorm (0.975)*sqrt(b/a/(a+b)+d/c/(c+d)))

例）
クロス表イメージ
matrix(c(76, 399, 129, 332),2,byrow=T)
     [,1] [,2]
[1,]   76  399
[2,]  129  332

a<-76; b<-399; c<-129; d<-332　#参考web2より
( rr <-  a*(c+d)/(a+b)/c )    #リスク比
[1] 0.5717829
rr*exp(qnorm ( c(0.025,0.975) )*sqrt(b/a/(a+b)+d/c/(c+d)))
[1] 0.4440632 0.7362370

または
rr*exp(c(-1,1)*qnorm (0.975)*sqrt(b/a/(a+b)+d/c/(c+d)))
[1] 0.4440632 0.7362370


参考web1　R -- 相対危険度（対応のない場合）
参考web2　R による統計処理
参考文献）柳川 堯 ; 観察データの多変量解析―疫学データの因果分析,近代科学社 ,2016

有意差とは・・・？

統計検定 t分布

乱数を発生させて、set.seed( )で記憶させてシミュレーションしてみます
乱数なので再現できませんが、set.seed( )を使用することで再確認できます

set.seed(1)     #もう一度確かめたいときはset.seed( )で乱数を記憶させておきます．( )の中は何でもＯＫ．
x20 <- rnorm(20,3.25,2.25)　　#平均3.25,標準偏差2.25のデータ(乱数)20個
set.seed(2)
y20 <- rnorm(20,2.12,1.95)　　#平均2.12,標準偏差1.95のデータ(乱数)20個
t.test(x20,y20)　　#独立した2群の差の検定を実施してみます

        Welch Two Sample t-test

data:  x20 and y20
t = 1.808, df = 37.999, p-value = 0.07852
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1409049  2.4959643
sample estimates:
mean of x mean of y 
 3.678679  2.501149 

♯ p-value =  0.07852　危険率0.05でも有意差は認められませんでした・・・

次に同じ平均値、同じ標準偏差で100個ずつ用意して、検定してみます

set.seed(3)
x100 <- rnorm(100,3.25,2.25)
set.seed(4)
y100 <- rnorm(100,2.12,1.95)
t.test(x100,y100)

        Welch Two Sample t-test

data:  x100 and y100
t = 3.6836, df = 196.81, p-value = 0.0002971
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.4491119 1.4841008
sample estimates:
mean of x mean of y 
 3.274830  2.308224 

P = 0.0002971  危険率0.01でも有意差あり！

統計量(この場合はｔ値)を考えれば当然の結果なのですが、
まだこのことに気づかずに有意差のみで議論している理学療法士も多いようです.
独立した2群であれば、それぞれのサンプルを増やすことで有意差を出すことが可能になります！

では、箱ひげ図で確認してみます

boxplot(x20,y20,x100,y100,xaxt="n")   
name<-c("x20","y20","x100","y100")
axis(side=1,at=c(1,2,3,4),labels=name)　

f:id:yoshida931:20181217225224p:plain

x20とy20に有意差がなくて、x100とy100に有意差がある・・・ようには見えませんね！
平均値と標準偏差が等しい正規分布からの乱数ですので、同じような箱ひげ図になるのは当然です．
x20とy20に有意差がなくて、x100とy100に有意差があるという結果の理解に苦しみます．
ここに数字のマジックが隠れています．

どのように有意差検定を実施しているのか、理解することが必要になります．
以下のページで、n増加→統計量大→ｐ値小の構図が理解できると思います．

yoshida931.hatenablog.com

さて次はどうでしょうか？

# x20,y20の平均値は同じで、標準偏差を1/10にしてみましょう
set.seed(5)     
x202 <- rnorm(20,3.25,0.225)　　
set.seed(6)
y202 <- rnorm(20,2.12,0.195)　
t.test(x202,y202)

    Welch Two Sample t-test
data:  x202 and y202
t = 15.657, df = 37.994, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.8960426 1.1621642
sample estimates:
mean of x mean of y 
 3.186874  2.157771 

# p-value = 2.2e-16 　危険率0.01でも有意差あり！

箱ひげ図で確認してみましょう．
f:id:yoshida931:20180124195732p:plain:w400

平均値が同じでも、ばらつき（分散、標準偏差）に違いがあれば統計結果は異なります．

有意差に一喜一憂するのではなく、サンプル数や分散、またそのデータがもつ性質をよく考えて、差を考察するべきです．Ｐ値よりも信頼区間や効果量などが比較の参考になるでしょう．
yoshida931.hatenablog.com

平均

統計検定 基本

算術平均 arithmetin mean (=相加平均)

1回　10人
2回　15人
3回　8人
平均回数は
(10*1+15*2+3*8)/(10+15+8)

度数分布からの平均
真の平均の近似値なので多少のズレが生じます
最初と最後の階級が少ない場合には無視して求めます（年収平均など）
無視できない数であれば、最初の下限は0にして、最後の上限は直線の範囲と同じ範囲にします


幾何平均 geometoric mean (＝相乗平均)

そもそも「幾何」とは、「いくばく」と読むので、わずか・すこしという意味があります 例1）利回りの例で考えてみます

10%, 15%, -5%, 13% と年々変化する利回り
平均は
 ( 1.1 + 1.15 + 0.95 + 1.13 ) / 4 
=1.0825 ではありません！
正解は・・・
( 1.1 * 1.15 * 0.95 * 1.13 )^(1/4)
= sqrt( sqrt( 1.1 * 1.15 * 0.95 * 1.13 ) )
=1.079501 
年平均利率は7.95％が正解です

奇数の場合は対数で計算します

例2）極端な例ですが、消費税で考えます

100円の品物の買う場合に
昨年5%になり、翌年15％に引きあがりました
この2年間の平均上昇率は何％でしょうか？
( 5 + 15 ) / 2 = 10% ではありません
(1.05 * 1.15 ) ^ (1/2) =1.098863
平均9.88％が正解です


調和平均 harmonic mean

距離÷時間＝速度を例にとって考えてみます

片道100kmを往復します
往路は10km/h, 復路は15km/hかかりました
平均速度は(10+15)/2=12.5km/hではありません
往路の時間は100/10時間、復路の時間は100/15時間で合計50/3時間になります
したがって平均時速は
(100*2)/( (100/10)+(100/15) ) = 12 km/h

移動平均

時系列 
x <- runif(500,-1,1)    #一様分布の乱数500個
plot(y2,type = "l")

f:id:yoshida931:20180122151901p:plain:w600

二乗平均平方根（RMS）

xr <-sqrt( x^2)

f:id:yoshida931:20180122152328p:plain:w600

移動平均

install.packages("TTR")
library(TTR)
x5 <- SMA(xr, 5)　# 移動平均間隔5
plot(x5,type = "l")
x10 <- SMA(xr, 10)　# 移動平均間隔10
plot(x10,type = "l")

f:id:yoshida931:20180122152720p:plain:w600

f:id:yoshida931:20180122152747p:plain:w600