9月2日（火）午前

mychisq.test()の検証 （やり直し）

先生のご指摘をいただいて、同じデータセットを用いてそれぞれのp値を出した方がいいことに気づく。

どうしてこんなことに気が付かなかったのかと思う。

というわけで、やり直し。

s <- function(){
sample(1:3,20,replace=TRUE)
}
s.t <- function(){
table(s())}#テーブル作成

tmp.list<-list()
cnt<-1
for(i in 1:100){
tmp.list[[cnt]]<-s.t()
cnt<-cnt+1
}

x<-c()
for(i in 1:100)x<-c(x,mychisq.test(tmp.list[[i]]))
y<-c()
for(i in 1:100)y<-c(y,chisq.test(tmp.list[[i]])$p.value)
#データをベクトルにする

var.test(x, y)

t.test(x,y,paired=TRUE)

> var.test(x, y)

        F test to compare two variances

data:  x and y
F = 0.9076, num df = 99, denom df = 99, p-value = 0.6305
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.6106748 1.3489140
sample estimates:
ratio of variances 
         0.9076055 

> t.test(x,y,paired=TRUE)

        Paired t-test

data:  x and y
t = -17.5323, df = 99, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.03771064 -0.03004266
sample estimates:
mean of the differences 
            -0.03387665 

結果的としては、（xの平均値）=（yの平均値）という帰無仮説は棄却されることになった。

しかし、

sample estimates:
mean of the differences 
            -0.03387665

ここを見ると、2つのペアの差の平均はとても小さいことが分かる。

すなわち、xとyの値は常に同じような値をとってはいるものの、ほとんど常にxのほうが少しだけyよりも小さいゆえに、帰無仮説が棄却されたという可能性も考えられる。

この場合は別に関数の作成に失敗したとは言えない。

ここへきて、母分散と母平均の検定をして帰無仮説が採択されればＯＫ…という考えが浅はかだったと痛感。

そこで気付いたのだが、

> x-y
  [1] -2.278622e-02 -4.185318e-02 -4.185318e-02 -5.557578e-02 -2.225453e-02
  [6] -6.568243e-02  2.151072e-02 -3.374049e-02 -1.873449e-02 -6.568243e-02
 [11] -3.559429e-02 -3.333014e-02 -3.374049e-02 -4.185318e-02 -3.374049e-02
 [16] -4.185318e-02 -9.096699e-02 -2.193640e-02  2.202013e-07 -2.278622e-02
 [21] -3.333014e-02 -2.378328e-02 -2.278622e-02 -6.568243e-02 -3.333014e-02
 [26]  2.202013e-07 -1.873449e-02 -3.333014e-02 -3.333014e-02 -3.374049e-02
 [31] -9.096699e-02 -4.185318e-02 -6.568243e-02 -2.378328e-02 -2.225453e-02
 [36]  4.601890e-03  2.202013e-07 -6.568243e-02 -4.185318e-02 -6.568243e-02
 [41] -2.278622e-02 -3.374049e-02 -3.333014e-02 -4.185318e-02 -3.374049e-02
 [46] -2.630728e-02 -9.096699e-02 -3.333014e-02 -2.378328e-02 -1.873449e-02
 [51] -4.185318e-02 -4.185318e-02 -1.873449e-02 -4.185318e-02 -3.559429e-02
 [56] -2.278622e-02 -3.333014e-02  4.601890e-03 -3.333014e-02 -4.185318e-02
 [61] -3.333014e-02 -3.374049e-02 -2.225453e-02 -3.559429e-02 -4.185318e-02
 [66] -4.185318e-02 -3.374049e-02 -3.333014e-02 -1.873449e-02 -8.485686e-03
 [71] -4.185318e-02 -6.568243e-02 -6.568243e-02 -3.374049e-02 -3.374049e-02
 [76] -4.185318e-02 -4.185318e-02 -3.374049e-02 -3.333014e-02 -2.378328e-02
 [81] -2.278622e-02 -5.557578e-02 -2.278622e-02 -3.374049e-02 -6.568243e-02
 [86] -4.185318e-02 -1.643381e-02 -3.374049e-02 -2.278622e-02  4.601890e-03
 [91] -2.278622e-02 -2.378328e-02 -3.374049e-02 -1.439835e-02 -3.374049e-02
 [96] -2.378328e-02 -3.333014e-02 -3.374049e-02 -4.185318e-02 -1.873449e-02

これを眺めるだけで、関数の作成には成功したと言ってよい気がする。

• 9.096699e-02はさすがに気になるが、全体としては間違いなく対応関係が見られる。

ただし、おそらく近似方法の問題か何かで、ほとんどの場合でmychisq.test()で求めた値のほうが小さくなるようである。

生存分析　2

いろいろ考えていたが、やはり生存曲線を描くための関数をつくってみるのが最も数学的理論の理解にもRの勉強に役立つと思ったので、そうしようと思う。

しばらくの間、様々なテーマについてガリガリ関数の自作をすることになりそうである。

RのパッケージMASSに、gehanという生存時間データがあるので、まずはそれをカプラン・マイヤー法で生存曲線にすることから始める。

成果は午後にまとめようと思う。

以上