この数日の半減期の変化を推定
この数日で「半減期が延びているのでは？」というツイートをいくつか見かけました。
そこで対数変換上での傾きが変わるモデルに変更し、推定し直しました。
- 11日からの経過日数が3日以下 ：-5
- 11日からの経過日数が13日以下 ：-5 + a1 + b1*(t - 3)
- 11日からの経過日数が17日以下 ：-5 + a1 + a2 + b1*(t - 3)
- 11日からの経過日数が17日より長い：-5 + a1 + a2 + b1*(17 - 3) + b2*(t - 17)
17日以降で半減期が変わったというモデルにしています。
このモデルで3/14～4/1の放射線量を推定すると次のようになります。
前半（3/28日まで）の半減期は約5日で後半（3/29以降）の半減期は約9日です。
ツイッターでも言われているように、I133の影響が薄れて来ているのかもしれません。
次にこの差が統計的に差があると言えるのか、Bootstrap法を使って検定します。
2,000回推定すると、（半減期2 - 半減期1）の分布は次のようになります。
| 2.5% | 50% | 97.5% | 片側P値 |
|---|---|---|---|
| -6.11日 | 3.55日 | 8.26日 | 0.052 |
半減期1より半減期2の方が平均的に3.55日長く、片側P値は0.052です。
「半減期1より半減期2の方が長い」という条件で検定するのであれば片側で十分ですので、有意といっても良いレベルだと思います。
統計学的な見地からも、現在のデータではここ数日は半減期が長くなっていると言ってもいいのかもしれません。
options(encoding="SJIS") HoushaSen <- read.csv("放射線.csv", as.is = T) plot(HoushaSen[, 3], HoushaSen[, 2], lwd = 2) plot(HoushaSen[, 3], log10(HoushaSen[, 2]), lwd = 2) #---ベースラインの補正 HoushaSenBase <- HoushaSen HoushaSenBase[, 2] <- HoushaSenBase[, 2] - 0.04 OptimHangenki <- function(Tvec, Hvec){ T <- Tvec H <- Hvec Hangenki <- function(Beta){ T1 <- 3 T2 <- 13 T3 <- 17 a1 <- Beta[1] a2 <- Beta[2] b1 <- Beta[3] b2 <- Beta[4] t <- T Pred <- ifelse(t <= T1, -5, ifelse(t <= T2, -5 + a1 + b1*(t - T1), ifelse(t <= T3, -5 + a1 + a2 + b1*(t - T1), ifelse(T3 < t , -5 + a1 + a2 + b1*(T3 - T1) + b2*(t - T3), NA)))) Resid <- log10(H) - Pred Resid%*%Resid } optim(c(0, 0, 0, 0), Hangenki) } (PredHangenki <- OptimHangenki(HoushaSenBase[, 3], HoushaSenBase[, 2])) T1 <- 3 T2 <- 13 T3 <- 17 Beta <- PredHangenki$par a1 <- Beta[1] a2 <- Beta[2] b1 <- Beta[3] b2 <- Beta[4] t <- HoushaSenBase[, 3] Pred <- ifelse(t <= T1, -5, ifelse(t <= T2, -5 + a1 + b1*(t - T1), ifelse(t <= T3, -5 + a1 + a2 + b1*(t - T1), ifelse(T3 < t , -5 + a1 + a2 + b1*(T3 - T1) + b2*(t - T3), NA)))) Resid <- log10(HoushaSenBase[, 2]) - Pred Resid%*%Resid (Half1 <- -(log10(20) - log10(10)) / b1) (Half2 <- -(log10(20) - log10(10)) / b2) plot(HoushaSenBase[, 3], log10(HoushaSenBase[, 2]), lwd = 2, xlab = "day (+11)", ylab = "log10()", main = "Sendai-shi Aoba-ku (baseline: -0.04μSv)") lines(t, Pred, lwd = 2) text(15, -0.6, paste("Estimated Half-period1: ", round(Half1, digits = 2), " day", sep = "")) text(15, -0.7, paste("Estimated Half-period2: ", round(Half2, digits = 2), " day", sep = "")) #bootstrap BootHangenki <- function(nboot, Housha, T1){ #---make data frame using T and D Data <- data.frame(Housha, T1) #---produce NULL data ResultBoot <- matrix(NA, nboot, 4) #---perform bootstrap SSI for(i in 1:nboot){ set.seed(i) BootID <- sample(1:nrow(Data), nrow(Data), replace = "T") BootSample <- Data[BootID, ] BootSample <- BootSample[order(as.numeric(rownames(BootSample))), ] T <- BootSample[, 2] H <- BootSample[, 1] PredHangenki <- OptimHangenki(T, H) ResultBoot[i, ] <- PredHangenki$par print(i) } a1 <- ResultBoot[, 1] a2 <- ResultBoot[, 2] b1 <- ResultBoot[, 3] b2 <- ResultBoot[, 4] hangen1 <- -(log10(20) - log10(10)) / b1 hangen2 <- -(log10(20) - log10(10)) / b2 print("median, 95% CI of Bootstrap hangenki") print(round(quantile(hangen1, probs = c(0.025, 0.5, 0.975), na.rm = T), digits =2)) print(round(quantile(hangen2, probs = c(0.025, 0.5, 0.975), na.rm = T), digits =2)) ResultBoot } Boot <- BootHangenki(2000, HoushaSenBase[, 2], HoushaSenBase[, 3]) Half1 <- -(log10(20) - log10(10)) / Boot[, 3] Half2 <- -(log10(20) - log10(10)) / Boot[, 4] HalfDiff <- Half2 - Half1 length(HalfDiff[HalfDiff < 0]) / length(HalfDiff) quantile(HalfDiff, probs = c(0.025, 0.5, 0.975))
