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

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

since2016　ときどきTEXのメモ

シンプソンのパラドクス

クロス表 定理・公式

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

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

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を操作しました．
上記のようなクロス表が作れたら、色々勉強できますよ！
医療関係者には最適な医療を提供してもらいたいですね！

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

𝜒^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

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

𝜒^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

九州理学療法士・作業療法士合同学会2017in宮崎で紹介したスクリプト

モンテカルロ法によるｔ検定の検定力分析
モンテカルロ法については以下を参照ください

yoshida931.hatenablog.com

tval <- numeric(length = 10000)　　#ｔ値(空のベクトルを設定)
cs <- 0
for(i in 1:10000){
 x <- rnorm(n = 30, mean = 10, sd = 3)  　# xの乱数30個生成
  y <- rnorm(n = 30, mean = 12, sd = 3) 　 # yの乱数30個生成
   res <- t.test(x, y, var.equal = T)　　　　　　　#xとｙのｔ検定の結果
    tval[i] <- res[[1]]　　　　　　　　　　　　　　　#ｔ値のみを抜き取る
     cs <- cs + ifelse(res[[3]] <  0.05, 1,0)　
}

#cs：ｐ値が0.05より小さい場合には1を足し、大きい場合には0を足す
#ｔ検定を1000回繰り返した結果、0.05より小さいP値の個数となる

cs/10000

データの要約、散布図、エラーバー付き平均値の一括処理

x<-iris$Sepal.Length[51:75]
y<-iris$Petal.Length[101:125]


sumxy <- function(x, y, TEST01){
 
  library(pwr)

  
  cat1 <- c("Min.", "1st Qu", "Median", "Mean", "3rd Qu.", "Max.")   
  xs <- summary( x ) 　　　　　　　　　　　　　　　　　　　　　　　
  ys <- summary(y)
  
  cat(cat1,"\n","x",xs,"\n","y",ys,"\n\n",file = "TEST01");
  tx <- t.test(x);　　　　　　　　　　　　　
  cat2<- c("tvalue","pvalue","95%CI");　　
  txt <- tx$statistic;　　　　　　　　　
  txp <- tx$p.value;　　　　　　　　　　
  txc1 <- tx$conf.int[1];　　
  txc2 <- tx$conf.int[2];　　　　　　
  x1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f",txt,txp,txc1,txc2); 
  
  ty <- t.test(y); 　　　
  tyt <- ty$statistic;
  typ <- ty$p.value;
  tyc1 <- ty$conf.int[1];
  tyc2 <- ty$conf.int[2];
  y1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f", tyt, typ, tyc1, tyc2);
  
  txy <- t.test(x,y,var.equal = T);
  txyt <- txy$statistic;　　　　　　
  txyp <- txy$p.value;　　　　　　　　　
  txyc1 <- txy$conf.int[1];　　　　　　　
  txyc2 <- txy$conf.int[2];　　　　　　　
  xy1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f",txyt,txyp,txyc1,txyc2);
  
  d <- (mean(x) - mean(y)) / sqrt((var(x) + var(y)) / 2);
  p <- pwr.t.test(n = length(x), d = d, sig.level = 0.05);
  
  xy <- c(x,y);
  zu1 <- pdf(file = "散布図.pdf");
  cate <- c(rep("1",length(x)),rep("2", length(y)));
  plot(xy ~ cate,xlim = c(0.3, 2.8),xaxp=c(1, 2, 1),xaxt="n",xlab = "", ylab = "") ;
  name<-c("X","Y")
  axis(side=1,at=c(1,2),labels=name);
  dev.off();
  zu1
  
  zu2 <- pdf(file = "エラーバー.pdf");
  se2 <- c(sqrt(var(x)/(length(x))), sqrt(var(y)/(length(y))));
  m2 <- c(mean(x), mean(y));
  plot(m2, xlim=c(0.5,2.5), type="b", ylim=c(5,7), xaxt="n",xlab = "", ylab = "");
  arrows(1:2,m2-se2,
         1:2,m2+se2,
         lwd=1.0,angle=90,length=0.1,code=3);
  name<-c("X","Y")
  axis(side=1,at=c(1,2),labels=name);
  dev.off();
  zu2
  
  zu3 <- pdf(file = "箱ひげ図.pdf");
  cate <- c(rep("1",length(x)),rep("2", length(y)));
  boxplot(xy ~ cate,xlim = c(0.3, 2.8),xaxp=c(1, 2, 1),xaxt="n",xlab = "", ylab = "") ;
  name<-c("X","Y")
  axis(side=1,at=c(1,2),labels=name);
  dev.off();
  zu3
  
  cat(cat2, "\n", "x", x1, "\n", "y", y1, "\n", "xy", xy1, "\n",
      "x:1標本t検定, y:1標本t検定, xy:独立2群t検定", "\n", "\n",
      "独立2群t検定の検定力", p[[4]], "\n", file = "TEST01", append =T);
}

sumxy(x, y, "TEST01")
#　出力されたファイルはDocumentsにあると思います