車輪を発明する

ファーマコメトリクス・モデリング研究所

rpois() 関数の非再現性について

R

set.seed() 指定してあっても、ある条件下で rpois() によるポアソン乱数発生結果が再現しない。
set.seed() で種を指定した上で、平均パラメータベクトルを与えてポアソン乱数を発生させる。
次に、また同じ種を指定した上で、別の平均パラメータベクトル（ベクトルのサイズは同じ。平均値が異なる）を与えると、結果の乱数が途中から異なってくる。初めのうちは同じなのだ。途中から以降はすべて異なってくる。場合がある。いや、この現象は再現できた。
考えるに、rpois() の中で、別の乱数発生ルーチンが呼ばれており、それが呼ばれる回数が平均パラメータの値次第で異なっているのだろうか。もしそうだとすると、この現象は説明がつく。
今のところ回避する方法が見つからない。こういうことが起こる（起こりうる）と意識した上で臨むしかない。

iPhone 死んだ

日記

20013 年から使ってきていた iPhone がとうとう死にました。2017 年ころからは sim なしになったので wifi のみで、ほぼ music 端末として使っていた。
3 月後半から調子が悪くなったので新たに iPod touch を購入し、それでもこの iPhone は music 端末としてだましだまし使ってきたのだけれど、今朝、とうとう最期を迎えた。リンゴループからどうしても抜け出さない。バッテリが消耗しきったのではないかと思う。ケーブルを挿しても充電を始めないのだ。充電開始のときのブルッという震えが起こらない。バッテリは一度交換してあるのだ。やっぱり、こういう交換バッテリはオリジナルのものよりも持ちが悪いのですね。
iPad を music 端末として使うことにした。でも、ライブラリの同期がまったく進まない。散々ため込んできたアルバムがリストアップされてこない。検索して聴くか。

KMMC plot

モデル評価 R 論文


https://www.page-meeting.org/pdf_assets/4280-2012-06%20PAGE%20KMMC.pdf

Survival model の共変量の効果を視覚的に評価する方法。うまい方法です。
シミュレーションが必要になる。R コードを書いておきましょう。

library(survival)
library(ggplot2)

set.seed(2901)

summary(fit.1 <- survreg(Surv(day, e1c0) ~ age, ds, dist="exponential"))
#fit.1$coefficients

# observed KMMC ----
summary(ds)

ds.kmmc <- data.frame(
  day=seq(7, 365, by=7)
)
ds.kmmc$age.mean <- NA

for (d in seq(nrow(ds.kmmc))) {
  age.ge <- ds$age[ds$day >= ds.kmmc$day[d]]
  ds.kmmc$age.mean[d] <- mean(age.ge)
}

g1 <- ggplot(ds.kmmc, aes(day, age.mean)) +
  geom_line()
#g1

# simulation ----
nsim <- 100

.sim.kmmc <- ds.kmmc
.sim.kmmc$isim <- NA
sim.kmmc <- data.frame()
for (isim in seq(nsim)) {
  ds.sim <- ds
  ds.sim$day <- round(rexp(nrow(ds), rate=exp(-(fit.1$coefficients[1]+fit.1$coefficients[2]*ds$age))), 0)

  for (d in seq(nrow(.sim.kmmc))) {
    age.ge <- ds.sim$age[ds.sim$day >= .sim.kmmc$day[d]]
    .sim.kmmc$age.mean[d] <- mean(age.ge)
  }
  .sim.kmmc$isim <- isim
  sim.kmmc <- rbind(sim.kmmc, .sim.kmmc)
}

res <- as.data.frame(matrix(unlist(by(sim.kmmc, sim.kmmc$day, FUN=function(d1){
  quantile(d1$age.mean, probs=c(0.05, 0.5, 0.95))
})), ncol=3, byrow=T))
names(res) <- c("P05","P50","P95")
res$day <- ds.kmmc$day

g2 <- ggplot(ds.kmmc, aes(day)) +
  geom_ribbon(aes(ymin=P05, ymax=P95), res, fill="green3", alpha=0.5) +
  geom_line(aes(y=P05), res, colour="black") +
  geom_line(aes(y=P95), res, colour="black") +
  geom_line(aes(y=age.mean), colour="blue", size=2) +
  ylab('Mean')
g2
ggsave('KMMC.png', dpi=600)

疫学、統計学への嫌悪感

統計

できの悪い統計学者や疫学者ほど嫌いなものはない。
彼らは自分たちを正義の人だと考えている。裁判官である。自分たちこそは科学を知っている。科学の論理をよく理解している。
世の者どもはどうしてこの論理を理解していないのかね。どうして勉強しようとしないのか。まあ、難しいから仕方なかろう。私たちが正しい道筋を教えて進ぜるから、君たちはそれに従いたまえ。
そういったニュアンスを彼らの文章の節々から感じる。統計学者のいう検定の原理や p 値の解釈、信頼区間の解釈も同じだ。最近では「因果関係とは何か」「因果と相関は違うのだ」という説明にこのニュアンスを感じることが多い。
統計学は科学の文法である、とは誰の言葉だったか。
科学の基礎を仕切る文法裁判官にでもなっているつもりなのだろう。人の間違いを指摘するのは心地よいものだ。
一方で、できの悪い科学者たちも自分たちの不出来・不勉強を自覚して裁判官の託宣にすがろうとする。裁判官のお言葉がこうだったから、とそれを金科玉条にして他の考えに耳を貸さない。
統計を勉強した医者の書いた統計の本で勉強してはならない、という言葉もあった。統計を勉強し、疫学を知った医者ほどたちの悪いものはないということを、今、ある本を読んでいてつらつら思った。
もちろん、統計学者、疫学者全員に当てはまることではない。立派な人もいると思う。いや、「いる」と断言はしない。なぜなら知らないから。だから「いると思う」と文末を濁す。
しかし、多くの統計学者はやはり、自分たちこそは科学の論理を支えているという自負があって、それが例の偉そうな口ぶりにつながってしまうのだと思う。