赤ちゃんに保湿剤を塗るとアトピー性皮膚炎が減るらしい。
赤ちゃんに毎日保湿剤 アトピー減
生後まもない赤ちゃんの皮膚に保湿剤を毎日塗ると、アトピー性皮膚炎になるリスクを30%減らすことができたとする研究成果を、国立成育医療研究センターのグループが発表しました。
特定の方法にアトピー性皮膚炎の予防効果があると証明できたのは、これが初めてだということです。この研究を行ったのは、国立成育医療研究センターの斎藤博久副研究所長らのグループです。
グループでは、アトピー性皮膚炎になった家族が1人以上いる、生後まもない赤ちゃん118人を無作為に半分に分けました。
そして、一方のグループの赤ちゃんには、一般の薬局などで売られている保湿剤を毎日、全身に塗り、もう一方のグループでは皮膚の乾燥している部分にのみワセリンを塗りました。
そして8か月後に調べたところ、保湿剤を塗った赤ちゃんのグループでは19人がアトピー性皮膚炎になったのに対し、ワセリンを塗ったグループは２８人で、保湿剤には発症のリスクを32%抑える効果があると証明できたとしています。
特定の方法にアトピー性皮膚炎の予防効果があると証明できたのはこれが初めてだということで、研究を行った国立成育医療研究センターの大矢幸弘医長は、「かゆみが出てしまってからでは遅いので、家族に患者がいるというハイリスクの赤ちゃんは、皮膚の状態が正常な段階で保湿剤を使うようお薦めできると思う」と話しています。
とあるのだが…
リスクが減るというのはたしかに32%で正しそう。
1-(19/59)/(28/59)
[1] 0.3214286
検定上は有意差がなかった。オッズ比としては半減するが有意ではない。
ato <- matrix(c(19, 28, 40, 31), nr=2) fisher.test(ato)
Fisher's Exact Test for Count Data data: ato p-value = 0.1321 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.2321089 1.1843236 sample estimates: odds ratio 0.5288016
…と言う話でおしまい、とせず、統計学とBioinformaticsをやるものとしては、この実験計画はどうすればよかったかを考えよう。
というわけで検出力検定をやる。リスク減とするかオッズ比とするかいろいろ混同しているが、カイ自乗検定用の検出力検定パッケージstatmodがあったので、リスク減だが考えているのはオッズ比ということにする(適当
(1) 32%のリスク減を検出力として設定したとして、2群でどれだけサンプルサイズが必要かを考える。
各々59人の時点では検出力が0.32程度しかない。検出力0.8を出そうとすると、各々160人くらいサンプルが必要そうである。
library(statmod) ato <- matrix(c(19, 28, 40, 31), nr=2) power.fisher.test(19/59, 28/59, 59, 59) n <- 59:200 d <- mapply(function(n) power.fisher.test(19/59, 28/59, n, n, nsim=1000), n) plot(n, d, xlab="Sample size (each group)", ylab="Power", ylim=c(0, 1))
(2) 各々59人のサンプルサイズとして、どれほどのEffect sizeを考えて実験を計画したのかを考える。
各々59人で検出力0.8になるEffect sizeを計算してみると、対照で28人のとき、11人アトピー性皮膚炎になれば検出力が0.8を超える。このとき、60%のリスク減となる。すげぇ。
1-(11/59)/(28/59)
[1] 0.6071429
res <- matrix(0, nr=60, nc=60) # 時間がかかる for(i in 0:59){ for(j in 0:59){ res[i+1, j+1] <- power.fisher.test(i/59, j/59, 59, 59, nsim=100) } print(i) # 進捗用 } cols <- rev(head(rainbow(5), -1)) image(0:59, 0:59, res, col=cols, xlab="", ylab="") abline(h=27, v=18, lty=3, lwd=2) legend("bottomright", legend=c("> 0.8", "> 0.5", "> 0.3", "< 0.3"), bg="white", pch=15, col=rev(cols), title="Power", cex=1.5)
というわけで、Effect sizeかサンプルサイズか、どちらにせよ設定を間違えているわけで、特定の方法にアトピー性皮膚炎の予防効果があると証明できたのはこれが初めてだ、というのはこの実験からは言えそうにない。
と思ったらカプランマイヤーだったらしい。
Journal of Allergy and Clinical Immunology, Volume 134, Issue 4, Pages 824–830.e6, October 2014.
原著を紹介しないNHKはクソだが読めない論文もクソとある筋から論文を入手した。
Methodの部分は
The primary outcome (cumulative rate of incidence of AD, eczema, or both by temporal observation) was analyzed by using the log-rank test. The significance level was set at .05. The Kaplan-Meier method was used to estimate the cumulative incidence of AD/eczema for each group, and the Cox regression model was applied to estimate the hazard ratio between groups.
とCox比例ハザードうんたらと書いてあったけど完全無視。
で、データを再現した。
ハザード率の推定にはbshazardを使った。
推定はハザード率をとして
になるから、
になる。
ハザード比(HR)はMSTが取れるならとしてもいいし、MSTが取れないなら平衡になったOSなどで
としても数学的には同じ(ログを取ると分子分母入れ替わることに注意)。
結果の部分で
During their first 32 weeks of life, 19 infants in the intervention group had AD/eczema compared with 28 infants in the control group. Calculation of cumulative incidence values for AD/eczema by using the Kaplan-Meier method showed that the intervention group maintained intact skin for a significantly longer period than the control group (P = .012, log-rank test; Fig 3). Cox regression analysis showed the risk of AD/eczema to be significantly lower in the intervention group (hazard ratio, 0.48; 95% CI, 0.27-0.86).
と、HR 0.48とほぼ同じような値。
保湿剤ではMSTが48週、対照では24週程度だから、保湿剤を使うとおよそ50%のリスク減となる。保湿剤が何でもよくて、副作用が許容できるなら安いし簡便だからとりあえずやってみてもいいと思う(手のひら返し
この論文はサンプルサイズ計算やpost-hoc analysisなどきちんと宣言してやっていたから好感が持てた。
library(survival) dat <- read.csv("clipboard") # 下のデータをコピペして使おう。 s0 <- Surv(time=dat$week, event=dat$cencor) ~ dat$group p <- pchisq(survdiff(s0)$chisq, 1, lower.tail=FALSE) plot(survfit(s0), col=c(2, 4), lwd=3, xlab="Weeks", ylab="Probability without AD/eczema") #abline(v=10*(0:3), lty=3) legend("bottomleft", legend=c("Intervention", "Control"), lwd=3, col=c(2, 4), lty=1, bty="n", cex=1.5, pch=3) legend("bottomright", legend=paste("p =", round(p, 3)), bty="n", cex=1.5)
library(bshazard) b <- bshazard(s0) hr <- print.bshazard(b)$rate # hazard rate log(2)/hr
# MST [1] 47.99497 24.67604
hr[1] / hr[2]
# Hazard ratio [1] 0.514138
week cencor group 0 0 1 5 0 1 6 0 1 6.5 1 1 7 0 1 7.2 0 1 7.5 1 1 7.8 1 1 8 0 1 8.8 1 1 9.5 1 1 9.5 0 1 11 1 1 11.6 1 1 12.3 1 1 12.3 1 1 12.7 1 1 12.9 1 1 13 1 1 13.3 1 1 13.8 1 1 15 0 1 15.2 1 1 16.3 0 1 20 0 1 21.5 0 1 22.5 0 1 23 1 1 23.4 0 1 25.2 0 1 25.4 0 1 26 0 1 27.5 1 1 30.5 0 1 31 0 1 31.1 0 1 31.2 0 1 31.3 0 1 32 0 1 32.5 0 1 33 1 1 33 1 1 33.1 0 1 33.2 0 1 33.3 0 1 33.2 0 1 33.4 0 1 33.5 0 1 33.5 0 1 33.6 0 1 33.7 0 1 34 0 1 33.9 0 1 34 0 1 34 0 1 34.1 0 1 34.2 0 1 34.3 0 1 34.3 0 1 0.5 0 2 4.5 0 2 4.7 0 2 5 1 2 5 0 2 5.1 1 2 5.2 1 2 6 0 2 6.1 1 2 6.2 1 2 7 1 2 7.5 0 2 7.7 0 2 8 1 2 8.2 0 2 9.9 1 2 10 1 2 10.1 1 2 10.1 0 2 10.2 1 2 10.3 0 2 10.6 1 2 11 1 2 11.2 1 2 11.4 1 2 12 1 2 12 0 2 12.7 1 2 12.7 0 2 12.8 0 2 13 1 2 13.1 1 2 13.2 0 2 14 1 2 16.5 1 2 17 0 2 18 1 2 20 1 2 21 0 2 23.8 0 2 24 1 2 24 1 2 24.5 1 2 26 1 2 28 1 2 31 0 2 31 0 2 31.2 0 2 31.3 0 2 31.4 0 2 31.5 0 2 31.7 0 2 31.8 0 2 32 0 2 32.2 0 2 32.4 0 2 32.5 0 2 33 0 2 34 0 2
