「統計学関連なんでもあり」の過去ログ--- 042

No.10686　p値について　　【初心者子】　2009/08/23(Sun) 18:32

いつもお世話になります。
初心者の質問にも回答していただき
非常にありがたいです。
今回の質問は
ログランクの検定でp値がp=0とでます。
どのくらいp値が低いとp=0と出るのでしょうか？

No.10688　追加：p値について　　【初心者子】　2009/08/23(Sun) 18:34

すいません。
解析ソフトはRです。

No.10689　Re: p値について　　【青木繁伸】　2009/08/23(Sun) 19:07

> 解析ソフトはRです。

どのライブラリの何という関数を使っているのですか。
実際に0と表示されたときのコンソールをコピーペーストしてくれるとはっきりするのですけど。

一般的に言えば，たとえば cat と sprintf の組み合わせだと，表示書式により四捨五入した結果によって，とても 0 とはいえないものも 0 として表示されます。
> x <- 0.035
> cat(sprintf("x = %.1f\n", x)) # 小数点以下1けた目まで表示せよ
x = 0.0
> cat(sprintf("x = %.2f\n", x)) # 小数点以下2けた目まで表示せよ
x = 0.04
> cat(sprintf("x = %.3f\n", x)) # 小数点以下3けた目まで表示せよ
x = 0.035
> cat(sprintf("x = %.4f\n", x)) # 小数点以下4けた目まで表示せよ
x = 0.0350
> cat(sprintf("x = %.5f\n", x)) # 小数点以下5けた目まで表示せよ
x = 0.03500
逆に，ほとんど 0 だろうというものでも，0 という表示にならないこともあります。
> y <- 1.23456789e-14
> cat(sprintf("y = %g\n", y))
y = 1.23457e-14
> old <- options(digits=20)
> cat(y)
1.23456789e-14
> old2 <- options(digits=3)
> cat(y)
1.23e-14
たとえば，カイ二乗分布の上側確率を求める pchisq 関数だと，
> pchisq(1482.5120154687307, 1, lower.tail=FALSE)
[1] 4.94e-324
> pchisq(1482.5120154687308, 1, lower.tail=FALSE)
[1] 0
のようになるので，自由度が1の場合には，カイ二乗値が 1482.5120154687307 と 1482.5120154687308 の間に，表示値が 0 になる境界があることがわかります。

No.10690　Re: p値について　　【初心者子】　2009/08/23(Sun) 19:36

下記にコピーいたします。
Rコマンダーで解析しております。
ご教授をお願いいたします。

> table(Ex169Mvial$event, Ex169Mvial$暴露)

C H
0 52 0
1 3 59

> Survival time description:
> tapply(Ex169Mvial$DAY, Ex169Mvial$暴露, summary, na.rm=TRUE)
$C
Min. 1st Qu. Median Mean 3rd Qu. Max.
2.690 9.690 9.690 9.442 9.690 9.690

$H
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.690 5.690 6.010 6.249 6.690 9.690

> Test results:
> survdiff(Surv(DAY , event)~暴露, data=Ex169Mvial)
Call:
survdiff(formula = Surv(DAY, event) ~ 暴露, data = Ex169Mvial)

N Observed Expected (O-E)^2/E (O-E)^2/V
暴露=C 55 3 39.9 34.2 118
暴露=H 59 59 22.1 61.9 118

Chisq= 118 on 1 degrees of freedom, p= 0

No.10691　Re: p値について　　【青木繁伸】　2009/08/23(Sun) 22:10

survdiff の出力結果は， Chisq で pchisq を使っているようですが，実際は pnorm を使っているようで（example(survdiff) で
> survdiff(Surv(jasa$futime, jasa$fustat) ~ offset(expect))
Call:
survdiff(formula = Surv(jasa$futime, jasa$fustat) ~ offset(expect))

Observed Expected        Z        p 
  75.000    0.587  -97.119    0.000 
というのが出る（ソースを読めばよいのだが面倒で。。。）で，pnorm がどのような引数をとると 0 を返すか，pchisq と同じくやって見ると，
> x <- 1400; pnorm(sqrt(x), lower.tail=FALSE)*2
[1] 2.101015e-306
> x <- 1410; pnorm(sqrt(x), lower.tail=FALSE)*2
[1] 0
これ以上正確に閾値を求めても意味がないので，これくらいにしておこう。
いずれにせよ，
Chisq= 118 on 1 degrees of freedom, p= 0 
なのだから，自由度1でカイ二乗値が118なら，正確には
> pchisq(118,1, lower.tail=FALSE)
[1] 1.73388e-27
で，0 に非常に非常に近いので，0と表示されようが1.73388e-27と表示されようが "< 0.001" と表示されようが，統計学的には同じこと。（注：統計学的には < 0.001 の表示が適切，1.73388e-27 でもよいけど，本来，このレベルの数値は近似値なので意味がないし。ただ，0 と表示するのはちょっと困るかもね。

No.10692　Re: p値について　　【青木繁伸】　2009/08/23(Sun) 22:27

survdiff の結果を表示するのは print.survdiff
そのソースを見ると，"Chisq" のラベルでの出力時に（別の形式で条件を書く方が適切だが，今回はそこまでする必要はないと思うので簡潔に），p は
format(signif(1 - pchisq(x$chisq, df), digits))
で出力される。示された例では，確かに，
> x$chisq <- 118
> df <- 1
> digits <- 10
> cat(format(signif(1 - pchisq(x$chisq, df), digits)))
0
となる。肝心な所は，
1 - pchisq(x$chisq, df)
となっているところ。ここは確かに，その結果は0になる。しかし，これと同じことを
pchisq(x$chisq, df, lower.tail=FALSE)
と書くと，
> x$chisq <- 118
> df <- 1
> digits <- 10
> cat(format(signif(pchisq(x$chisq, df, lower.tail=FALSE), digits)))
1.73388e-27
となる。
しかし，見かけは違っても，いずれもきわめて0に近いという意味ではどちらも正しい（繰り返すが，0 と表示するのはマズイとは思う）

なお，統計学を離れて数値計算的には 1 - pchisq(x$chisq, df) と chisq(x$chisq, df, lower.tail=FALSE) では，後者が適切。前者は，桁落ちが生じている。桁落ちとは，つまり，コンピュータのような近似計算！しかできない場合！に，似た数値の引き算をすると有効桁が極端に少なくなってしまうこと。コンピュータを使って数値計算をするときは，このようなことにいつも注意していないと不適切な答えを出してしまうという典型例ですね。（コンピュータは本質的に，近似計算しかできないんです。驚きましたか？）

結論。見かけではなく，意味を考えよう！

No.10693　Re: p値について　　【初心者子】　2009/08/24(Mon) 05:36

ありがとうございます。
ご教授をすべては理解できておりませんが
意味を考えたいと思います。

数値計算においては，コンピュータも統計も知らないことばかりです。
とても奥が深いです。