Libraries
> library(readr)
> library(dplyr)
> library(sp)
> library(magrittr)
> library(scanstatistics)
Counts
> JPN_popcov <- read_csv("JPN_popcov.csv", locale = locale(encoding = "SHIFT-JIS"))
> head(JPN_popcov, 5)
# A tibble: 5 x 4
week pref pop count
<dbl> <chr> <dbl> <dbl>
1 6 01.Hokkaido 52.9 0
2 6 02.Aomori-ken 12.6 0
3 6 03.Iwate-ken 12.4 0
4 6 04.Miyagi-ken 23.2 0
5 6 05.Akita-ken 9.81 0
> counts <- JPN_popcov %>%
+ filter(week >= 9 & week <= 11) %>%
+ df_to_matrix(time_col = "week", location_col = "pref", value_col = "count")
>
> counts[,1:5]
01.Hokkaido 02.Aomori-ken 03.Iwate-ken 04.Miyagi-ken 05.Akita-ken
9 24 0 0 1 0
10 7 0 0 0 1
11 24 0 0 0 0
Zones
> JPN_geo <- read_csv("JPN_geo.csv", locale = locale(encoding = "SHIFT-JIS"))
> head(JPN_geo, 5)
# A tibble: 5 x 3
pref lat lon
<chr> <dbl> <dbl>
1 01.Hokkaido 43.1 141.
2 02.Aomori-ken 40.8 141.
3 03.Iwate-ken 39.7 141.
4 04.Miyagi-ken 38.3 141.
5 05.Akita-ken 39.7 140.
> zones <- JPN_geo %>%
+ select(lon, lat) %>%
+ as.matrix %>%
+ spDists(x = ., y = ., longlat = TRUE) %>%
+ dist_to_knn(k = 10) %>%
+ knn_zones
> zones[[1]]
[1] 1
> zones[[2]]
[1] 1 2
> zones[[350]]
[1] 41 42 43 45 46
Baselines
> mod <- glm(count ~ offset(log(pop)) + 1 + I(week - 8),
+ family = poisson(link = "log"),
+ data = JPN_popcov %>% filter(week < 9))
> summary(mod)
Call:
glm(formula = count ~ offset(log(pop)) + 1 + I(week - 8), family = poisson(link = "log"),
data = JPN_popcov %>% filter(week < 9))
Deviance Residuals:
Min 1Q Median 3Q Max
-3.4720 -1.2999 -0.8976 -0.5937 6.0792
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.6825 0.1036 -25.897 < 2e-16 ***
I(week - 8) 0.7453 0.1151 6.472 9.64e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 393.77 on 140 degrees of freedom
Residual deviance: 346.00 on 139 degrees of freedom
AIC: 423.72
Number of Fisher Scoring iterations: 6
> ebp_baselines <- JPN_popcov %>%
+ filter(week >= 9 & week <= 11) %>%
+ mutate(mu = predict(mod, newdata = ., type = "response")) %>%
+ df_to_matrix(value_col = "mu")
> ebp_baselines[,1:5]
01.Hokkaido 02.Aomori-ken 03.Iwate-ken 04.Miyagi-ken 05.Akita-ken
9 7.617563 1.820088 1.788384 3.337548 1.413702
10 16.050575 3.835013 3.768211 7.032375 2.978739
11 33.819339 8.080557 7.939803 14.817554 6.276347
Calculation
> set.seed(1)
> poisson_result <- scan_eb_poisson(counts = counts,
+ zones = zones,
+ baselines = ebp_baselines,
+ n_mcsim = 999)
> print(poisson_result)
Data distribution: Poisson
Type of scan statistic: expectation-based
Setting: univariate
Number of locations considered: 47
Maximum duration considered: 3
Number of spatial zones: 360
Number of Monte Carlo replicates: 999
Monte Carlo P-value: 0.001
Gumbel P-value: 0
Most likely event duration: 1
ID of locations in MLC: 28
> prefectures <- as.character(JPN_geo$pref)
> prefectures[c(28)]
[1] "28.Hyogo-ken"