Latest update on 2024年3月1日 (金) at 19:06:18.
6:25に目が覚めたが、テレビ体操は布団の中で腕だけ動かした。疲れが抜けていない。
今日は2限の全学のオンライン会議と午後の来年度の卒業研究の学生との面談が時間が決まった業務だが、合間に締め切り過ぎの仕事と締切間近の仕事とメール処理を進めなくては。
いろいろ済んだところで17:00を回った。メールボックスを見に行ったら、朝倉書店から佐藤俊哉『宇宙怪人しまりす 統計よりも重要なことを学ぶ』、ISBN 978-4-254-12297-8(Amazon | honto | e-hon)(以下しまりす本3と呼ぶ)が届いていた。アマゾンで予約購入していたものは既に自宅に届いていて、読み始めているので、この謹呈分はラボに置いて院生や学部生に読ませよう。『シンプル衛生公衆衛生学2024』南江堂、ISBN 978-4-524-21022-0(Amazon | honto | e-hon)も執筆者分が1冊届いていた。
しまりす本3のpp.12「回復割合の差は20%」のp値を計算しようとしたが合わない。何かぼくのミスだと思うが原因がわからない。
x <- c(58, 62); n <- c(80, 100); X <- cbind(x, n-x) colnames(X) <- c("回復", "未回復") rownames(X) <- c("ヨクナール", "経過観察") print(X) s <- c(colSums(X)[1]/sum(X), colSums(X)[2]/sum(X)) # 帰無仮説:回復割合が処置と独立の場合の期待度数E print(E <- rbind(ヨクナール=rowSums(X)[1]*s, 経過観察=rowSums(X)[2]*s)) E/rowSums(E) # 確認:回復割合にヨクナールと経過観察で差はない 1-pchisq(sum((X-E)^2/E), 1) # 帰無仮説の下でのp値 Z <- c(0.1, -0.1) # ヨクナール群の回復割合が経過観察群より20%良い場合の期待度数EE print(EE <- rbind(ヨクナール=rowSums(X)[1]*(s+Z), 経過観察=rowSums(X)[2]*(s-Z))) EE/rowSums(EE) # 確認:回復割合がヨクナールの方が20%良い 1-pchisq(sum((X-EE)^2/EE), 1) # ヨクナールの回復割合が20%高い仮説の下でのp値
というRコードを実行すると、 差がないという帰無仮説(しまりす本3では「ゼロ仮説」)の下でのp値は本と同じく13.8%になるのだけれども、最後に表示される「回復割合の差が20%」という仮説の下でのp値は16.4%になってしまって、しまりす本3に書かれている17.2%とは合わない。しまりす本3ではpp.13-14に17.2%の算出方法が書かれていて、上のコードで使ったカイ二乗適合度検定とは違う方法(下記)なのだけれども、どちらでも一致するはずなんだがなあ。自分の頭の回転が鈍くて嫌になる。
rdp <- function(a, b, N1, N0, np=0) { RD <- a/N1-b/N0 seRD <- sqrt(a*(N1-a)/N1^3+b*(N0-b)/N0^3) Z <- (RD-np)/seRD return(list(est=RD, se=seRD, p=2*pnorm(-abs(Z)))) }
とRで関数定義すると、確かにrdp(58, 62, 80, 100, np=0.2)は17.2%になるのだが、今度はrdp(58, 62, 80, 100)(rdp(58, 62, 80, 100, np=0)と同じ。Rでは、関数定義の際に引数に初期値を与えておくと、実行時には引数ごと省略可能)が13.2%になってしまって本と合わない。
たぶん自分の初歩的な勘違いかコーディングミスだと思うが、わからないので、読者諸賢からのコメントを期待して統計相談メモにも書いておく。
青春アカペラ甲子園 全国ハモネプリーグ~2023春~の配信販売が始まったのでmoraで購入した。神戸大学ふなでが「輝く月のように」という曲をカバーしているが、相変わらず素晴らしいリードボーカルの声の伸びだと思う。今度の土曜夜に生放送で行われるハモネプにも応募はしていたが50チームから25チームが選ばれる段階で残れなかったのは残念。その分、全国大会に出場するハヰカラ使節団には実力を発揮して爪痕を残してきて欲しい。
▼前【1566】(代休ということになっているが仕事(2024年2月26日) ) ▲次【1568】(会議2つか3つ(2024年2月28日) ) ●Top