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

統計学備忘録 since2016

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

シンプソンのパラドクス

最近、統計学の学習が進んでおりません
ちょっとお仕事が忙しくて…

今日は少しだけ勉強しておきます.
観察研究には重要な概念です.

x <- c(110,70,90,120)
x <- matrix(x,2,2)
rownames(x)<-c("治療A","治療B")
colnames(x)<-c("改善","変化なし")
addmargins(x) 

        改善  変化なし  Sum
治療A  110       90     200
治療B   70      120     190
Sum    180      210     390

これだけ見て、「治療Aが良い!」と結論付けるのは早いのです.
患者さんの症状別にみた次の表を見てみましょう!

y <- c(3,10,40,93,17,31,42,20,90,29,8,7)
dat <- array(y, dim = c(2,2,3))
name <- list("治療"=c("A","B"),"変化"=c("改善","変化なし"),"患者"=c("重度","中等度","軽度"))
dimnames(dat)<-name

dat

, , 患者 = 重度

    変化
治療 改善 変化なし
   A    3       40
   B   10       93

, , 患者 = 中等度

    変化
治療 改善 変化なし
   A   17       42
   B   31       20

, , 患者 = 軽度

    変化
治療 改善 変化なし
   A   90        8
   B   29        7

#患者の重症度によって分類してみたら・・・必ずしも治療Aが良いとは言えません.
#症状が重度、または中等度の患者さんには治療Aはあまり効果がないようです.
#症状が中等度の患者さんには治療Bが効果的のようです.
#観察研究では交絡因子の影響をブロックすることが重要になります.

今日は、久しぶりにRを操作しました.
上記のようなクロス表が作れたら、色々勉強できますよ!
医療関係者には最適な医療を提供してもらいたいですね!

カイ二乗検定後の残差分析

2017.12.5更新

残差=実測値-期待値

標準化残差= \frac{残差}{\sqrt{期待値 }}

標準化残差の分散=\sqrt{(1-\frac{列計}{n})(1-\frac{行計}{n})}

調整済残差=\frac{(実測値-期待値)}{\sqrt{期待値(1-\frac{列計}{n})(1-\frac{行計}{n})}}

#各質問に関するyesの回答者数の割合について有意水準5%で検定します.
xxx <- c(20,10,2,5,410,350,200,120)
xxx <- matrix(xxx,4,2)

   yes  no
Q1  20 410
Q2  10 350
Q3   2 200
Q4   5 120

rname <- c("Q1","Q2","Q3","Q4")  #行名
cname <- c("yes","no")           #列名
rownames(xxx)<-paste(rname)   #行名のペースト
colnames(xxx)<-paste(cname)   #列名ののペースト

#周辺和
addmargins(xxx)             

    yes   no  Sum
Q1   20  410  430
Q2   10  350  360
Q3    2  200  202
Q4    5  120  125
Sum  37 1080 1117

#期待値
xe <- chisq.test(xxx)$expected

         yes       no
Q1 14.243509 415.7565
Q2 11.924799 348.0752
Q3  6.691137 195.3089
Q4  4.140555 120.8594

#残差=実測値-期待値
re <- xxx-xe

          yes         no
Q1  5.7564906 -5.7564906
Q2 -1.9247986  1.9247986
Q3 -4.6911370  4.6911370
Q4  0.8594449 -0.8594449

#調整済残差
asr <- re/sqrt(xe*outer(1-rowSums(xxx)/sum(xxx), 1-colSums(xxx)/sum(xxx)))

1-rowSums(xxx)/sum(xxx)
       Q1        Q2        Q3        Q4 
0.6150403 0.6777081 0.8191585 0.8880931 

1-colSums(xxx)/sum(xxx)
       yes         no 
0.96687556 0.03312444 

outer(1-rowSums(xxx)/sum(xxx), 1-colSums(xxx)/sum(xxx))
         yes         no
Q1 0.5946674 0.02037287
Q2 0.6552594 0.02244870
Q3 0.7920243 0.02713417
Q4 0.8586755 0.02941759

asr 
          yes         no
Q1  1.9779359 -1.9779359
Q2 -0.6885780  0.6885780
Q3 -2.0377875  2.0377875
Q4  0.4557999 -0.4557999

#累積分布
 pnorm(asr,lower.tail=TRUE )
          yes         no
Q1 0.97603203 0.02396797
Q2 0.24554445 0.75445555
Q3 0.02078559 0.97921441
Q4 0.67573307 0.32426693

Q1にyesと回答した人の割合は、他の質問項目に比べて有意に大きく、Q3にyesと回答した人の割合は、他の質問項目に比べて有意に小さいことがわかりました.

# r-de-r様からアドバイスいただきました
# 調整済み残差を直接算出するRの関数です
xxx <- c(20,10,2,5,410,350,200,120)
xxx <- matrix(xxx,4,2)
cht<-chisq.test(xxx)
cht$stdres

           [,1]       [,2]
[1,]  1.9779359 -1.9779359
[2,] -0.6885780  0.6885780
[3,] -2.0377875  2.0377875
[4,]  0.4557999 -0.4557999

カイ自乗近似は不正確かもしれません

このコメントが出る場合は以下のようなときです

期待度数が 0になるセルがある場合、もしくは期待度数が 5未満になるセルが全体の 20% を超える場合

R: Fisher's Exact Test for Count Data

that is if no cell has expected counts less than 1 and more than 80% of the cells have expected counts at least 5, otherwise the exact calculation is used.  (カイ二乗分布近似は、期待度数が1未満になるセルがない場合、また80%以上のセルに少なくとの5個の期待度数がある場合に使用される.そうでない場合には、正確検定が使用される.)

カイ自乗近似は不正確かもしれません
このコメントが出た場合には、フィッシャーの正確確率検定を実施します   

yoshida931.hatenablog.com

統計検定2級

統計検定2級をこれから目指す人のために
私がお勧めする本です

確率が苦手のまま検定に望んではいけません!
この本を読んでから統計学に入りましょう!

高校生が感動した確率・統計の授業 (PHP新書)

高校生が感動した確率・統計の授業 (PHP新書)

次にお勧めするのは、簡単な統計の本です.
嫌にならない程度に慣れることが重要だと思います.
Rを使いながら、この本で勉強してみてください.
Amazonでも4つ星がついてました.

Rによるやさしい統計学

Rによるやさしい統計学

でも、やはり基礎をしっかり学習する必要があります.
ボロボロになるまで見たのが、これです.

改訂版 日本統計学会公式認定 統計検定2級対応 統計学基礎

改訂版 日本統計学会公式認定 統計検定2級対応 統計学基礎

上記の本をある程度理解できる人は、この本をみるべきでしょう!
未だに完全理解しておりませんが・・・

統計学入門 (基礎統計学?)

統計学入門 (基礎統計学?)

Rを使いながら、視覚的に統計を学習した方が良いと思います.
非常に良い本が出版されました.これです.

Rとグラフで実感する生命科学のための統計入門

Rとグラフで実感する生命科学のための統計入門

最後に、検定前に毎日のように取り組んだ問題集が次の2冊です.

統計学演習

統計学演習


日本統計学会公式認定 統計検定 2級 公式問題集[2014〜2016年]

日本統計学会公式認定 統計検定 2級 公式問題集[2014〜2016年]

私は理学療法士ですが、将来的に統計学が重要視されることを想定して学習しています.
次に紹介する本が、それを示唆しています.

バイオ統計の基礎―医薬統計入門 バイオ統計シリーズ1

バイオ統計の基礎―医薬統計入門 バイオ統計シリーズ1

統計ソフトに入力するだけではなく、統計解析の意味を理解しましょう!!!

決定係数(寄与率)

決定係数 {R}^2
回帰係数がどの程度よく当てはまっているか、つまりXがYをうまく説明できているか、ということを明らかにすることは重要になります.決定係数は、その当てはまりを示す基準としてして使用されています.

観測値の平方和(全変動) Sy
回帰による平方和  Sr
残差による平方和 Se

全変動 Sy= 回帰による平方和Sr + 残差平方和Se

決定係数 {R}^2=\frac{Sr}{Sy}

観測値の平方和(全変動) は回帰方程式で説明できる変動と説明できない変動に分けられます.決定係数は回帰による平方和(回帰方程式で説明できる変動)の割合です.ピアソンの相関係数の二乗と同じ値になります.残差二乗和が小さく、決定係数が1に近いほどよいモデルと言えます.(重回帰の場合には、必ずしもそうではありません)

サンプルは回帰分析のページと同じものを使います

x <- subset(CO2, subset = conc == 95, select = c(uptake))     
y <- subset(CO2, subset = conc == 175, select = c(uptake)) 
x <- x$uptake
y <- y$uptake

#回帰分析の結果
Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   1.4784     6.1126   0.242  0.81377   
x             1.6972     0.4876   3.480  0.00592 **

#残差を可視化
e<-y-(1.4784+1.6972*x)
plot(e,xlab="x",ylab="残差", xaxt="n")
abline(h=0, lty=2, col=2)

回帰残差のプロット
f:id:yoshida931:20171128182430j:plain:w400
あまり良いモデルではないようです・・・

平方和の分解から決定係数を算出

全変動を求めてみます

sum((y-mean(y))^2)
=432.7167

残差による変動を求めてみます

yの予測値をyeとします
ye <- 1.6972*x + 1.4784
sum((y-ye)^2)
=195.682

決定係数を求めてみます  

(432.7167-195.682)/432.7167
=0.5477826
#やはり、このモデルではxがyを十分に説明できていません

決定係数は相関係数の二乗になります
相関係数rを求めてみます

cov <- sum((x-mean(x))*(y-mean(y)))
sxy <- sqrt(sum((x-mean(x))^2))*sqrt(sum((y-mean(y))^2))
r = cov/sxy
r^2=(cov/sxy)^2=0.5477825
#決定係数と同じ値になりました

誤差項の分散は回帰分析のあてはまりの良さを表します
推定値の標準誤差 (Residual standard error) を求めてみます
残差 e= yi -1.6972*x - 1.4784

e <- y - 1.6972*x - 1.4784
s = sqrt(sum(e^2)/(length(x)-2))
  =4.423596
#全ての残差が±2sには入っていません
#上の図の中に±2sを書き加えてみたら、1個だけはみ出ておりました・・・

f:id:yoshida931:20171128183358j:plain:w400

続きは下記のページへ
回帰係数の区間推定 - 統計学備忘録 since2016

予測値と残差の相関は0になります

cor(e,ye)
= 3.885081e-06

f:id:yoshida931:20171128184453j:plain:w400

次回は回帰係数の検定を勉強します

テキストファイルに検定結果を出力

テキストファイルに検定結果を出力
投稿日2017.11.14

r-de-r様からコメントいただきましたので修正しております
素直に感動しております

テキスト形式のファイルに検定結果を書き込む練習をします
sprintf:書式指定変換した出力を文字列に格納
cat:文字列を表示する基本的な関数, file=でファイルに書き込み

備忘録
\n 改行
%.3f 小数点3桁まで出力
append = F:ファイルに上書き
append = T:ファイルに追加

Rのサンプルirisを使います

x <- iris$Sepal.Length[51:75]
y <- iris$Petal.Length[101:125]
#noteというファイルに、xとyの平均を書き込みます
cat("x平均", mean(x), "y平均", mean(y), file = "note.txt") #noteというファイル名

f:id:yoshida931:20170908160200j:plain

sprintfの使い方
xの平均を小数点以下1桁まで、項目名「Xの平均」と付いた文字列にして取り出します

sprintf("xの平均 %.1f", mean(x))
"xの平均 6.0"

t検定を実施して、strで中身を確認してみます

t <- t.test(x,y,var.equal = T)
str(t)

Two Sample t-test

data: x and y
t = 2.1954, df = 48, p-value = 0.033
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.03131358 0.71268642
sample estimates:
mean of x mean of y
6.012 5.640


List of 9
$ statistic : Named num 2.2
..- attr(*, "names")= chr "t"
$ parameter : Named num 48
..- attr(*, "names")= chr "df"
$ p.value : num 0.033
$ conf.int : atomic [1:2] 0.0313 0.7127
..- attr(*, "conf.level")= num 0.95

t値、p値、95%信頼区間のみ取り出してテキストファイル「note」に書き出してみます

tv <- sprintf("t値 %.4f\n", t[[1]]) 
pv <- sprintf("p値 %.4f\n", t[[3]]) 
ci <- sprintf("95%信頼区間 %.4f %.4f\n", t[[4]][1], t[[4]][2])     #ここがポイントです(%を小文字にする場合には%%と書きます)
cat(tv,pv, ci, sep="",file = "note.txt")

catのデフォルトでsep=" "となっているので改行した後にスペースが入るそうです.スペースなし(sep="")にして書き直しました.無駄なスペースが無くなってスッキリしました.

f:id:yoshida931:20171114133508j:plain:w300

参考と引用 Rプログラム (TAKENAKA's Web Page)

単回帰分析

投稿日2017.2.2
更新日2017.11.13

線形回帰分析

説明変数と目的変数との関係を直線でモデル化する回帰分析。

Xi <- c(35, 45, 55, 65, 75)
Yi <- c(114, 124, 143, 158, 166)
plot(Xi, Yi, pch = 16)
abline(lm(Yi~Xi), col = 2)

f:id:yoshida931:20171108084407j:plain:w300

どのようにして赤のラインを引いたのか?を考えていきます.

母回帰方程式 

Yi=\beta1 + \beta2 * Xi + ei\ \ (i=1,2,3,,,n)

これが大前提となる母集団のモデルです.一般的に知ることができないので、このモデルについて推定・検定することを回帰分析と呼んでいます.

母回帰係数 \beta1, \beta2
誤差項 ei (これが回帰分析のポイントになります)
 ・期待値=0
 ・分散は一定
 ・異なる誤差項は無相関

E(Yi)=\beta1+\beta2*Xi

  

最小二乗法

最小二乗法により、母回帰係数(\beta1, \beta2) を推定します. その推定値を(\hat{\beta1}, \hat{\beta2})とします

母回帰方程式より、誤差項ei=Yi-(\beta1 + \beta2 * Xi) となります

最小二乗法とは、S= \sum {ei}^2とおいた場合に、Sが最小になる(\hat{\beta1}, \hat{\beta2})を求める手法です.

S=\sum({Yi}^2-2β1Yi-2β2XiYi+{β1}^2+2β1β2Xi+{β2}^2{Xi}^2)

Sを最小にするために、\beta1, \beta2 でそれぞれ微分した二つの方程式を解くことで次式を求めます.

標本(偏)回帰係数
\hat{\beta2}=\frac{\sum{(Xi-\bar{Xi})(Yi-\bar{Yi})}}{\sum{(Xi-\bar{Xi})}^2}  \hat{\beta1}=\bar{Yi}-\hat{\beta2}\bar{Xi}

標本回帰方程式
Y=\hat{\beta1}+\hat{\beta2}\bar{X}

RのサンプルCO2(草の耐寒性実験)からデータを抜き出して回帰分析を行ってみます.

# データフレームから必要な部分だけ抜き出します.(CO2濃度95ml/Lと175ml/Lの二酸化炭素摂取量)
x <- subset(CO2, subset = conc == 95, select = c(uptake))     
y <- subset(CO2, subset = conc == 175, select = c(uptake)) 
x <- x$uptake
y <- y$uptake

summary(lm(formula=y~x) )

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.1062 -1.8958 -0.7628  2.0096 10.0376 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   1.4784     6.1126   0.242  0.81377   
x             1.6972     0.4876   3.480  0.00592 **
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.424 on 10 degrees of freedom
Multiple R-squared:  0.5478, Adjusted R-squared:  0.5026 
F-statistic: 12.11 on 1 and 10 DF,  p-value: 0.005917

plot(x,y)
abline(lm(y~x))

f:id:yoshida931:20171113125825j:plain:w300
標準化残差を確認してみます

ei<-y-(1.4784+1.6972*x)
eis<-ei/sqrt( (sum(ei^2) )/(length(x)-2) )  #標準化残差
plot(eis,xlab="CO2濃度95ml/L",ylab="CO2濃度175ml/L")
abline(h=0, lty=2, col=2)

f:id:yoshida931:20171113131309j:plain:w300

次にRの関数で求めた回帰モデル y = (1.4784+1.6972*x) を最小二乗法で求めてみます.
最小二乗法より標本回帰係数\hat{\beta2} を求めます

sum ( ( y - mean(y) ) * ( x - mean(x) ) ) / sum ( ( x - mean(x) )^2 )
= 1.697206

標本回帰係数\hat{\beta1} (y切片)を求めます

mean(y) - 1.697206*mean(x)
=1.478416

回帰モデル
y = 1.4784 + 1.6972*x

参考) 豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012