December 8, 2013
By

(This article was first published on Wiekvoet, and kindly contributed to R-bloggers)

Fukushima is still a topic which gets headlines, and somewhere in a comment was a link to actual and historical radiation data: Fukushima prefecture radioactivity measurement map.  It takes a bit of time to load but then you have a map of the Fukushima region with current radiation levels at measurement stations. Somewhere on that website there is data from the past year, this post examines the decrease at four of the stations.

### Data

If you look at the map, there are loads of blue dots, indicating low radiation, and a line of yellow, orange and red dots, striking north east from the reactor. It is an interactive map. You can click on any dot and see some details for that region. And at the bottom of that is a link to a detail page. On that page there are, among other things, graphs such as this one.
One of these plots contains 'transit of air dose rates in the most recent 365 days. It is a dynamic plot, the link to it contains the data, just like the plot above. It is pretty easy to extract the data from the link. Since getting the links was manual work, I decided to restrict my data to four dots. To discriminate between locations I copied some header and a number. For example 【相双地方】室原公民館の測定値 (Google translate: Measurements [phase] of bi-district community center Murohara) with number 704. The complete data with location information are at the bottom of the page.
There may be a more simple way to get the data, but since I don't speak Japanese, I cannot say either way.

### Analysis

Since the data looks pretty choppy, I wanted a robust model. To quote the robust task view:  "robust depends on robustbase, the former providing convenient routines for the casual user where the latter will contain the underlying functionality, and provide the more advanced statistician with a large range of options for robust modeling." Regarding robust I am certainly the casual user. Proportional decrease by using a linear model on log transformed data seemed an appropriate modelling strategy.
The plot below contains the data and the predictions. For me, the predictions seem just like I would make them when ignoring some of the crazy jumps.
The data themselves, show a clear correlation. The same drop mid January, a jump in September, though I do not know an explanation.
Half lives can be estimated from the models. Two of them are around 2. Cesium-134 has a half live of 2.04 years and has apparently been released, it is easy to speculate Cesium-134 is part of the source of this radiation.
Lower CI     mean Upper CI
community center Murohara         3.171937 4.463582 7.529786
meeting place to the power before 1.922496 2.074042 2.251525
Sugiuchi meeting place            1.764994 1.817172 1.872530
township village center           2.602469 2.854075 3.159538

### code

#【相双地方】室原公民館の測定値
#704
# Measurements [phase] of bi-district community center Murohara
#【相双地方】杉内集会所の測定値
#628
# Measured value of the bi-phase region] Sugiuchi meeting place
# 【相双地方】前乗集会所の測定値
#741
# Measurement of the meeting place to the power [bi-phase region] before
# 【相双地方】町区集落センターの測定値
#1376
# Measured value of the bi-phase region] township village center

ys <- gsub('.*t:','',ys)
ys <- gsub('&chxl.*','',ys)
dose <- scan(textConnection(ys),sep=',')
data.frame(dose=dose,station=i,relday = -(length(dose):1))
}

locations <-data.frame(station=c(704,628,741,1376),
loc=factor(c('community center Murohara','Sugiuchi meeting place','meeting place to the power before', 'township village center')))

r1$doseus <- r1$dose/100
r1$day <- as.Date("05-12-2013",format='%d-%m-%Y')+r1$relday
Sys.setlocale("LC_TIME",'C')

r1 <- merge(r1,locations)
r1$status <- 'Observed' library(robust) models <- lapply(levels(r1$loc),function(st) {
r2 <- r1[r1$loc==st,] lmrob(log10(doseus) ~ relday,data=r2) }) topred <- r1[r1$station==r1$station[1],c('relday','day')] preds <- lapply(1:nlevels(r1$loc),function(i) {
preds <- topred
preds$doseus <- 10^predict(models[[i]],preds) preds$status <- 'Predicted'
preds$loc <- levels(r1$loc)[i]
preds
}
)
preds <- do.call(rbind,preds)
preds <- merge(preds,locations)
all <- rbind(subset(r1,select=-dose),preds)

library(lattice)
png('fuk1.png')
xyplot(doseus ~ day | loc,groups= status,data=all,type='l',
scales=list(y=list(log=10) ),
ylab=expression(mu * "Sv per hour" ))
dev.off()

halflives <-     t(sa <- sapply(la,function(x) {
point <- coef(x)[2]
sd <- sqrt(vcov(x)[2,2])
yy <- c(log10(.5)/c(0.025=point-1.96*sd,mean=point,0.975=point+1.96*sd))/365
names(yy) <- c('Lower CI','mean','Upper CI')
yy
}))
rownames(halflives) <- levels(r1\$loc)
halflives