ベースライン補正＆Bootstrap法で半減期の信頼区間を推定・検定

先日の仙台市青葉区の放射線量の解析について、@tomoki_kitawakiさんからアドバイスを頂きました。

有難うございます。

ベースライン補正は対数変換前に行わないと加法として考慮されないとのことでした。

そのご指摘を受けベースライン値を0.04μSv/hとして、先にこれを引いて減衰モデルを適用しました。

すると、@tomoki_kitawakiさんの5.4日という結果とほぼ一致しました。

修正版の半減期の推定値は5.51日です。

14日より以前の放射線量はベースラインを引いているので0のはずです。

しかし0に対数変換をすると-Infになってしまうので、パラメータ推定の際には対数変換後の値を-5と設定しました。

次にこれをBootstrap法で2,000回計算しますと、半減期の分布は次のようなヒストグラムになりました。

各パラメータの信頼区間は次のようになっています。

パラメータ	95%下側限界	中央値	95%上側限界
半減期	4.73	5.54	7.04
15日増加量	0.259	0.303	0.358
25日増加量	0.0502	0.0690	0.0847

半減期の95%信頼区間は4.73日から7.04日。

つまりI133の半減期である8日と検定すると有意になります。

よって「I133の半減期を早める何らかの要因がある」ことが示唆されます。

この要因を探索するためには気象データ等をチェックする必要があると思われます。

取り急ぎ以上となります。

setwd("C:/Documents and Settings/Issei/My Documents/Dropbox/")
HoushaSen <- read.csv("iAnalysis/その他/放射線/放射線.csv", as.is = T)

#---ベースラインの補正
HoushaSenBase      <- HoushaSen
HoushaSenBase[, 2] <- HoushaSenBase[, 2] - 0.04


OptimHangenki <- function(Tvec, Hvec){
	T <- Tvec
	H <- Hvec
	Hangenki <- function(Beta){
		T1 <- 14
		T2 <- 24
		a1 <- Beta[1]
		a2 <- Beta[2]
		b1 <- Beta[3]
		t  <- T
		Pred  <- ifelse(t  <= T1, -5, 
			 ifelse(t  <= T2, -5 + a1 + b1*(t - T1),
			 ifelse(T2 <  t , -5 + a1 + a2 + b1*(t - T1), NA)))
		Resid <- log10(H) - Pred
		Resid%*%Resid
	}
	optim(c(0, 0, 0), Hangenki)
}


BootHangenki <- function(nboot, Housha, T1){
	#------bootstrap
	library(bootstrap)
	
	#---make data frame using T and D
	Data  <- data.frame(Housha, T1)
	
	#---produce NULL data
	ResultBoot      <- matrix(NA, nboot, 3)
	
	#---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]
	hangen <- -(log10(20) - log10(10)) / b1
	
	print("median, 95% CI of Bootstrap hangenki")
	print(round(quantile(hangen, probs = c(0.025, 0.5, 0.975), na.rm = T), digits =2))
	
	ResultBoot
}

Boot <- BootHangenki(2000, HoushaSenBase[, 2], HoushaSenBase[, 3])

library(ggplot2)
Temp <- data.frame(半減期 = -(log10(20) - log10(10)) / Boot[, 3])
ggplot(Temp, aes(半減期)) + geom_histogram()


#---25日、増加する前の対数放射線量
log10.25 <- -5 + a1 + Boot[, 3]*11

quantile(10^(-5 + Boot[, 1]), probs = c(0.025, 0.5, 0.975))
quantile(10^(log10.25 + Boot[, 2]) - 10^(log10.25), probs = c(0.025, 0.5, 0.975))
quantile(-(log10(20) - log10(10))/Boot[, 3], probs = c(0.025, 0.5, 0.975))