使用 R 分析 2015 年臺南市本土登革熱疫情狀況

這裡使用 R 來分析 2015 年臺南市本土登革熱疫情狀況。

今年台灣的登革熱疫情很嚴重，台南地區的病例更突破兩萬人，以下我們利用 R 來分析台南市的登革熱疫情。

登革熱資料來源

原始資料可以從臺南市政府資料開放平台的網站上取得，以下是我整理過的資料，後續分析所需要的資料都在這裡：

這兩個檔案是臺南市本土登革熱病例數，截至 2015 年 11 月 7 日的資料。兩個檔案內容相同，只差在檔案編碼，一個是 UTF-8，另外一個是 BIG-5，如果是 Windows 的使用者請下載 BIG-5 的檔案，而 Mac OS X 或 Linux 則可使用 UTF-8 的檔案。

資料整理

將資料讀進 R 中，並確認資料的正確性，排除有問題的資料再進行後續分析。

將臺南市本土登革熱病例數資料讀進 R 中：

dengue <- read.csv("dengue-20151107-utf8.csv")

檢查一下資料的狀況：

str(dengue)

'data.frame':	22258 obs. of  6 variables:
 $ 確診日  : Factor w/ 141 levels "2015/01/06","2015/01/19",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ 區別    : Factor w/ 41 levels "安定區","安南區",..: 12 9 6 6 6 6 6 2 6 6 ...
 $ 里別    : Factor w/ 530 levels " 安西里"," 仁愛里",..: 265 371 236 236 236 499 499 519 31 499 ...
 $ 道路名稱: Factor w/ 1489 levels ""," 建東街","?埕路",..: 764 1316 536 536 536 454 454 1341 1052 454 ...
 $ 緯度座標: num  23 23 23 23 23 ...
 $ 經度座標: num  120 120 120 120 120 ...

summary(dengue)

確診日           區別           里別
 2015/09/24:  712   北區   :3529   正覺里 :  332
 2015/09/18:  642   中西區 :3455   永祥里 :  319
 2015/09/21:  604   永康區 :2632   大豐里 :  299
 2015/09/22:  599   北　區 :2179   勝利里 :  272
 2015/09/23:  584   安南區 :1831   五王里 :  269
 2015/09/10:  553   南區   :1824   成德里 :  266
 (Other)   :18564   (Other):6808   (Other):20501
    道路名稱        緯度座標        經度座標
 金華路 :  690   Min.   :22.63   Min.   :118.8
 公園路 :  671   1st Qu.:22.99   1st Qu.:120.2
 西門路 :  627   Median :23.00   Median :120.2
 文賢路 :  430   Mean   :23.01   Mean   :120.2
 長榮路 :  423   3rd Qu.:23.01   3rd Qu.:120.2
 府安路 :  343   Max.   :24.95   Max.   :121.5
 (Other):19074

這裡的經緯度座標值有些問題，有些經緯度已經超出台南市的範圍，我們可以將資料以 ggmap 畫在地圖上來看：

library(ggmap)
library(mapproj)
map <- get_map(location = "Taiwan", zoom = 7,
  language = "zh-TW", maptype = "roadmap")
ggmap(map, darken = c(0.5, "white")) +
  geom_point(aes(x = 經度座標, y = 緯度座標),
  color = "red", data = dengue)

畫出來的圖形為：

dengue-map-outlier-1

再畫一張比較大的圖：

map <- get_map(location = "Tainan", zoom = 9,
  language = "zh-TW", maptype = "roadmap")
ggmap(map, darken = c(0.5, "white")) +
  geom_point(aes(x = 經度座標, y = 緯度座標),
  color = "red", data = dengue)

畫這張圖的時候，出現一行警告訊息，表示有三筆資料沒有畫在上面（這三筆可以從上面那張圖看出來）。

警告訊息：
Removed 3 rows containing missing values (geom_point).

這是畫出來的圖形：

dengue-map-outlier-2

我們透過圖形來訂出一個篩選資料的方形區域：

map <- get_map(location = "Tainan", zoom = 9,
  language = "zh-TW", maptype = "roadmap")
ggmap(map, darken = c(0.5, "white")) +
  geom_point(aes(x = 經度座標, y = 緯度座標),
  color = "red", data = dengue) +
  geom_rect(aes(xmin = 120, xmax = 120.6, ymin = 22.8, ymax = 23.5),
  alpha = 0.1)

畫出來的圖形為：

dengue-map-outlier-3

接著把實際的資料篩選出來：

filter.idx1 <- dengue$緯度座標 > 22.8 & dengue$緯度座標 < 23.5
filter.idx2 <- dengue$經度座標 > 120 & dengue$經度座標 < 120.6
dengue.tn <- dengue[filter.idx1 & filter.idx2, ]

把篩選好的資料畫在地圖上：

map <- get_map(location = c(lon = 120.246100, lat = 23.121198),
  zoom = 10, language = "zh-TW")
ggmap(map, darken = c(0.5, "white")) +
  geom_point(aes(x = 經度座標, y = 緯度座標),
  color = "red", data = dengue.tn)

畫出的圖形為

dengue-map-1

這張圖就是臺南市全年的登革熱病例分佈地圖。

從上面 dengue 的 summary 輸出中可以看到行政區的名稱有一些問題，我們將 dengue.tn$區別 的 levels 列出來看一下：

levels(dengue.tn$區別)

[1] "安定區"   "安南區"   "安平區"   "白河區"   "北　區"
 [6] "北區"     "大內區"   "東　區"   "東區"     "東山區"
[11] "官田區"   "關廟區"   "歸仁區"   "後壁區"   "佳里區"
[16] "將軍區"   "柳營區"   "六甲區"   "龍崎區"   "麻豆區"
[21] "南    區" "南　區"   "南化區"   "南區"     "楠西區"
[26] "七股區"   "仁德區"   "山上區"   "善化區"   "西港區"
[31] "下營區"   "新化區"   "新市區"   "新營區"   "學甲區"
[36] "鹽水區"   "永康區"   "永康區 "  "玉井區"   "中西區"
[41] "左鎮區"

部份的區別名稱有包含空白，而有些卻沒有，區別名稱不統一會造成程式將一個區域誤判為多個區域。這裡我們修正區別的名稱，統一將所有的空白去掉：

dengue.tn[dengue.tn$區別 == "北　區", ]$區別 <- "北區"
dengue.tn[dengue.tn$區別 == "東　區", ]$區別 <- "東區"
dengue.tn[dengue.tn$區別 == "南　區" | dengue.tn$區別 == "南    區", ]$區別 <- "南區"
dengue.tn[dengue.tn$區別 == "永康區 ", ]$區別 <- "永康區"

重新建立一次 factor，這樣可以將空的 levels 去掉：

dengue.tn$區別 <- factor(dengue.tn$區別)

然後再確認一次區別名稱：

levels(dengue.tn$區別)

[1] "安定區" "安南區" "安平區" "白河區" "北區"
 [6] "大內區" "東區"   "東山區" "官田區" "關廟區"
[11] "歸仁區" "後壁區" "佳里區" "將軍區" "柳營區"
[16] "六甲區" "龍崎區" "麻豆區" "南化區" "南區"
[21] "楠西區" "七股區" "仁德區" "山上區" "善化區"
[26] "西港區" "下營區" "新化區" "新市區" "新營區"
[31] "學甲區" "鹽水區" "永康區" "玉井區" "中西區"
[36] "左鎮區"

整體趨勢分析

分析全臺南市本土登革熱病例數整體分佈與趨勢。找出疫情最嚴重的時段。

畫出每週登革熱的病例數統計圖：

hist(as.Date(dengue.tn$確診日), breaks = "weeks",
  freq = TRUE, main = "登革熱每週病例數", xlab = "日期",
  ylab = "病例數", format = "%m/%d")

畫出的圖形為

dengue-hist-1

計算每個月的登革熱病例數：

dengue.tn$month <- format(as.Date(dengue.tn$確診日), "%m")
table(dengue.tn$month)

01    05    06    07    08    09    10    11
    2     2    11   202  3023 13054  5565   393

畫圖：

barplot(table(dengue.tn$month), xlab = "月份", ylab = "病例數",
  main = "登革熱每月病例數")

dengue-hist-3

也可以使用 ggplot2 來畫：

library(ggplot2)
library(scales)
ggplot(dengue.tn, aes(x=as.Date(確診日))) +
  stat_bin(binwidth=7, position="identity") +
  scale_x_date(breaks=date_breaks(width="1 month")) +
  theme(axis.text.x = element_text(angle=90)) +
  xlab("日期") + ylab("病例數") +
  ggtitle("登革熱每週病例數")

畫出的圖形為

dengue-hist-2

從圖形上可以看出登革熱的疫情最嚴重的時期是在九月份前後。

接下來分析疫情最嚴重的區域。計算各個行政區的病例總數：

dengue.region.summary <- sort(summary(dengue.tn$區別), decreasing = FALSE)
dengue.region.summary

東山區 龍崎區 後壁區 山上區 將軍區 楠西區
     3      3      5      6      7      8
白河區 大內區 南化區 左鎮區 六甲區 下營區
     9      9     10     10     14     15
鹽水區 柳營區 西港區 七股區 官田區 學甲區
    15     18     19     21     22     27
關廟區 安定區 麻豆區 玉井區 新市區 佳里區
    48     50     51     57     61     67
善化區 新營區 新化區 歸仁區 仁德區 安平區
    87    103    126    179    287    875
安南區 永康區   東區   南區 中西區   北區
  1831   2633   2964   3450   3455   5707

畫圖：

barplot(dengue.region.summary, las = 2, horiz = TRUE,
  main = "各行政區病例統計", xlab = "病例數")

畫出的圖形為

dengue-barplot-1

使用圓餅圖：

pie(dengue.region.summary)

畫出的圖形為

dengue-pie-chart-1

疫情最嚴重的五個區域依序為北區、中西區、南區、東區與永康區。

細部分析

將資料依據地理位置與時間區分，進行細部的分析。

將最嚴重的五個行政區病例資料篩選出來：

dengue.top.reg <- dengue.tn[
  dengue.tn$區別 == "北區" |
  dengue.tn$區別 == "中西區" |
  dengue.tn$區別 == "南區" |
  dengue.tn$區別 == "東區" |
  dengue.tn$區別 == "永康區", ]

依據時間畫出這 5 個行政區的疫情變化：

ggplot(dengue.top.reg, aes(x=as.Date(確診日))) +
  stat_bin(binwidth=7, position="identity") +
  scale_x_date(breaks=date_breaks(width="1 month")) +
  theme(axis.text.x = element_text(angle=90)) +
  xlab("日期") + ylab("病例數") +
  ggtitle("登革熱每週病例數") + facet_grid(區別 ~ .)

dengue-region-hist-1

依照月份來畫圖：

ggplot(dengue.top.reg, aes(x=as.Date(確診日))) +
  stat_bin(breaks=as.numeric(seq(as.Date('2015-1-1'),
    as.Date('2015-12-1'), '1 month')), position="identity") +
  scale_x_date(breaks=date_breaks(width="1 month")) +
  theme(axis.text.x = element_text(angle=90)) +
  xlab("日期") + ylab("病例數") +
  ggtitle("登革熱每月病例數") + facet_grid(區別 ~ .)

dengue-region-hist-2

看起來這 5 個區域最嚴重的時間都在 9 月中旬附近。

依據月份區分，畫出每個月的登革熱病例分佈地圖。

map <- get_map(location = c(lon = 120.246100, lat = 23.121198),
  zoom = 10, language = "zh-TW")
ggmap(map, darken = c(0.5, "white")) +
  geom_point(aes(x = 經度座標, y = 緯度座標),
  color = "red", data = dengue.tn) +
  facet_wrap(~ month)

dengue-map-wrt-time-1

定點分析

假設某人居住在台南市，其住家的經緯度座標為 (22.997088, 120.201771)，而登革熱病媒蚊飛行活動範圍可遠至 400 ～ 800 公尺的地區，將此人住家 400 公尺以內的病例資料篩選出來，觀察每個月的疫情變化。

這是計算地球上兩點之間距離的函數，輸入兩點的經緯度，可以計算出兩點之間的距離，單位為公里。

earthDist <- function (lon1, lat1, lon2, lat2){
  rad <- pi/180
  a1 <- lat1 * rad
  a2 <- lon1 * rad
  b1 <- lat2 * rad
  b2 <- lon2 * rad
  dlon <- b2 - a2
  dlat <- b1 - a1
  a <- (sin(dlat/2))^2 + cos(a1) * cos(b1) * (sin(dlon/2))^2
  c <- 2 * atan2(sqrt(a), sqrt(1 - a))
  R <- 6378.145
  d <- R * c
  return(d)
}

篩選出 400 公尺以內的病例資料：

home.pos <- c(22.997088, 120.201771) # (緯度, 經度)
home.dist <- earthDist(dengue.tn$經度座標, dengue.tn$緯度座標, home.pos[2],  home.pos[1])
home.idx <- home.dist <= 0.4;
dengue.home <- dengue.tn[home.idx, ]

查看每個月的資料狀況：

table(dengue.home$month)

05  07  08  09  10  11
  2   1  42 309 119   5

畫成圖形：

barplot(table(dengue.home$month), xlab = "月份", ylab = "病例數",
  main = "登革熱每月病例數（特定區域）")

dengue-hist-4

每個月的病例分佈：

map <- get_map(location = c(lon = home.pos[2], lat = home.pos[1]),
  zoom = 16, language = "zh-TW", color = "bw")
ggmap(map) +
  geom_point(aes(x = 經度座標, y = 緯度座標),
  color = "red", data = dengue.home, size = 5) +
  facet_wrap(~ month)

dengue-map-wrt-time-2

由於經緯度資料的精確度不足，造成大量的資料點重疊，因此改用 jitter 的方式畫圖：

map <- get_map(location = c(lon = home.pos[2], lat = home.pos[1]),
  zoom = 16, language = "zh-TW", color = "bw")
ggmap(map) +
  geom_jitter(aes(x = 經度座標, y = 緯度座標),
  size = 3, position = position_jitter(w = 0.0005, h = 0.0005),
  data = dengue.home, color = "red") +
  facet_wrap(~ month)

這樣可以比較清楚呈現資料的分布狀況：

dengue-map-wrt-time-3

假設今天是 2015 年 9 月 15 日，在地圖上畫出此人住家 400 公尺以內，2 週之內新增的病例分佈地圖。

篩選出此人住家 400 公尺以內，兩週之內新增的病例：

dengue.home$day.diff <- as.numeric(as.Date(dengue.home$確診日) - as.Date("2015/09/15"))
dengue.home.subset <- dengue.home[dengue.home$day.diff >= 0 & dengue.home$day.diff < 14, ]

依照時間決定顏色畫圖：

map <- get_map(location = c(lon = home.pos[2], lat = home.pos[1]),
  zoom = 16, language = "zh-TW", color = "bw")
ggmap(map) +
  geom_jitter(aes(x = 經度座標, y = 緯度座標, color = day.diff),
  size = 3, position = position_jitter(w = 0.0005, h = 0.0005),
  data = dengue.home.subset) +
  scale_colour_gradientn(colours=heat.colors(3))

dengue-map-wrt-time-4

登革熱資料來源#

資料整理#

整體趨勢分析#

細部分析#

定點分析#

參考資料#

登革熱資料來源

資料整理

整體趨勢分析

細部分析

定點分析

參考資料