2017-12-27

正規分布

下記のサイトに移転しました

2017-12-19

クラスター分析

グラフ

Rのサンプルattitudeを使用して、クラスタ分析（階層的方法）を勉強します．
attitudeは、管理者態度のデータです．無作為に選ばれた35名の雇用者よりアンケート．好意的な割合が数値化されています．
rating全般的評価、 complaints雇用者からの苦情処理、 privileges特別な特権の付与、 learning学習機会の付与、 raises業績による昇進機械、 critical批判的態度、 advance進歩

# xxはdata.frame

(result<-hclust(dist(xx),method="single"))
plot(result, main = "最近隣法")

(result<-hclust(dist(xx),method="complete"))
plot(result, main = "最遠隣法")

(result<-hclust(dist(xx),method="average"))
plot(result, main = "群平均法")

(result<-hclust(dist(xx),method="centroid"))
plot(result, main = "重心法")

(result<-hclust(dist(xx),method="median"))
plot(result, main = "メディアン法")

(result<-hclust(dist(xx),method="ward.D"))
plot(result, main = "ウォード法 ward.D")

(result<-hclust(dist(xx),method="ward.D2"))
plot(result, main = "ウォード法 ward.D2")

実際にウォード法を実行してみます　まずは変数クラスタを求めてみます

xx <- t(attitude)
(result<-hclust(dist(xx),method="ward.D"))
plot(result, main = "ウォード法 ward.D")

f:id:yoshida931:20171219114804p:plain:w400
進歩、批判的態度を除くと2グループに管理者の評価は分かれるようです進歩、批判的態度、向上意欲をくみ取る能力、雇用者の全般的な評価といった分類になりそうです（これかは私の勝手な見解です）

つぎにサンプルクラスタをもとめてみます

xx <- attitude 
(result<-hclust(dist(xx),method="ward.D"))
plot(result, main = "ウォード法 ward.D")

f:id:yoshida931:20171219114910p:plain:w700
雇用者が数グループに分けられそうです．
こんなデータが病院経営戦略にも役立つのかも…？

引用・参考間瀬茂；R 基本統計関数マニュアル(https://cran.r-project.org/doc/contrib/manuals-jp/Mase-Rstatman.pdf)

2017-12-12

シンプソンのパラドクス

クロス表定理・公式

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

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

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

ｶｲ二乗検定後の残差分析

𝜒^2分布クロス表

2017.12.5更新

残差＝実測値-期待値

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

標準化残差の分散= $(1-\frac{列周辺和}{n})(1-\frac{行周辺和}{n})$

標準化残差は、近似的に正規分布N(0 , 標準化残差の分散) に従う

$調整済（標準化）残差=\frac{残差}{\sqrt{期待値×標準化残差の分散}}$

$調整済（標準化）残差=\frac{(実測値-期待値)}{\sqrt{期待値(1-\frac{列計}{n})(1-\frac{行計}{n})}}$ $\approx N(0, 1 )$

#各質問に関するyesの回答者数の割合について有意水準５％で検定します．
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

2017-12-01

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

𝜒^2分布クロス表

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

期待度数が 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