用R分析2015台南登隔熱疫情?
許sir
[設定所需的函式庫(libraries)以及載入資料]
setwd("/media/hsusir/DATA/Rdata Practice/04Spatial/dengue fever")
dengue <- read.csv("./dengue-20151107-utf8.csv", fileEncoding="utf8")[Part 1].Data-distributed
1-1.敘述性統計
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 "七股區","下營區",..: 38 27 7 7 7 7 7 16 7 7 ...
## $ 里別 : Factor w/ 530 levels "一甲里","七甲里",..: 103 262 77 77 77 323 323 219 253 323 ...
## $ 道路名稱: Factor w/ 1489 levels "","30鄰五福街",..: 878 1310 296 296 296 183 183 1414 1357 183 ...
## $ 緯度座標: num 23 23 23 23 23 ...
## $ 經度座標: num 120 120 120 120 120 ...summary(dengue)## 確診日 區別 里別 道路名稱
## 2015/09/24: 712 北區 :3529 正覺里 : 332 金華路 : 690
## 2015/09/18: 642 中西區 :3455 永祥里 : 319 公園路 : 671
## 2015/09/21: 604 永康區 :2632 大豐里 : 299 西門路 : 627
## 2015/09/22: 599 北 區 :2179 勝利里 : 272 文賢路 : 430
## 2015/09/23: 584 安南區 :1831 五王里 : 269 長榮路 : 423
## 2015/09/10: 553 南區 :1824 成德里 : 266 府安路 : 343
## (Other) :18564 (Other):6808 (Other):20501 (Other):19074
## 緯度座標 經度座標
## Min. :22.63 Min. :118.8
## 1st Qu.:22.99 1st Qu.:120.2
## Median :23.00 Median :120.2
## Mean :23.01 Mean :120.2
## 3rd Qu.:23.01 3rd Qu.:120.2
## Max. :24.95 Max. :121.5
## 1-2.畫在地圖上來看
install.packages("ggmap")
install.packages("mapproj")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)
1-3.再畫一張比較大的圖
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)
1-4.訂出一個篩選資料的方形區域
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)
1-5.把實際的資料篩選出來
filter.idx1 <- dengue$緯度座標 > 22.8 & dengue$緯度座標 < 23.5
filter.idx2 <- dengue$經度座標 > 120 & dengue$經度座標 < 120.6
dengue.tn <- dengue[filter.idx1 & filter.idx2, ]1-6.把篩選好的資料畫在地圖上
- 這張圖就是臺南市全年的登革熱病例分佈地圖
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)
[Part 2].Data-ETL
- 從上面 dengue 的 summary 輸出中可以看到行政區的名稱有一些問題，我們將 dengue.tn$區別 的 levels 列出來看一下
levels(dengue.tn$區別)## [1] "七股區" "下營區" "中西區" "仁德區" "佳里區" "六甲區"
## [7] "北區" "北 區" "南化區" "南區" "南 區" "南 區"
## [13] "善化區" "大內區" "學甲區" "安南區" "安定區" "安平區"
## [19] "官田區" "將軍區" "山上區" "左鎮區" "後壁區" "新化區"
## [25] "新市區" "新營區" "東區" "東 區" "東山區" "柳營區"
## [31] "楠西區" "歸仁區" "永康區" "永康區 " "玉井區" "白河區"
## [37] "西港區" "關廟區" "鹽水區" "麻豆區" "龍崎區"- 部份的區別名稱有包含空白，而有些卻沒有，區別名稱不統一會造成程式將一個區域誤判為多個區域。這裡我們修正區別的名稱，統一將所有的空白去掉
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] "七股區" "下營區" "中西區" "仁德區" "佳里區" "六甲區" "北區"
## [8] "南化區" "南區" "善化區" "大內區" "學甲區" "安南區" "安定區"
## [15] "安平區" "官田區" "將軍區" "山上區" "左鎮區" "後壁區" "新化區"
## [22] "新市區" "新營區" "東區" "東山區" "柳營區" "楠西區" "歸仁區"
## [29] "永康區" "玉井區" "白河區" "西港區" "關廟區" "鹽水區" "麻豆區"
## [36] "龍崎區"[Part 3].整體趨勢分析
- 分析全臺南市本土登革熱病例數整體分佈與趨勢。找出疫情最嚴重的時段。
3-1.畫出每週登革熱的病例數統計圖
hist(as.Date(dengue.tn$確診日), breaks = "weeks",
freq = TRUE, main = "登革熱每週病例數", xlab = "日期",
ylab = "病例數", format = "%m/%d")
3-2.計算每個月的登革熱病例數
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 393barplot(table(dengue.tn$month), xlab = "月份", ylab = "病例數",
main = "登革熱每月病例數")
#使用 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("登革熱每週病例數")
- 從圖形上可以看出登革熱的疫情最嚴重的時期是在九月份前後
3-3.分析疫情最嚴重的區域。計算各個行政區的病例總數
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 5707barplot(dengue.region.summary, las = 2, horiz = TRUE,
main = "各行政區病例統計", xlab = "病例數")
pie(dengue.region.summary)
[Part 4].細部分析
- 將資料依據地理位置與時間區分，進行細部的分析
4-1.將最嚴重的五個行政區病例資料篩選出來?
dengue.top.reg <- dengue.tn[
dengue.tn$區別 == "北區" |
dengue.tn$區別 == "中西區" |
dengue.tn$區別 == "南區" |
dengue.tn$區別 == "東區" |
dengue.tn$區別 == "永康區", ]4-2.依據時間畫出這 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(區別 ~ .)
4-3.依照月份來畫圖?
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(區別 ~ .)
- 看起來這 5 個區域最嚴重的時間都在 9 月中旬附近
4-4.依據月份區分，畫出每個月的登革熱病例分佈地圖
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)
[Part 5].定點分析
- 假設某人居住在台南市，其住家的經緯度座標為 (22.997088, 120.201771)，而登革熱病媒蚊飛行活動範圍可遠至 400 ～ 800 公尺的地區，將此人住家 400 公尺以內的病例資料篩選出來，觀察每個月的疫情變化。
5-1.計算出兩點之間的距離，單位為公里
- 這是計算地球上兩點之間距離的函數，輸入兩點的經緯度，可以計算出兩點之間的距離，單位為公里。
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)
}5-2.篩選出 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, ]5-3.查看每個月的資料狀況
table(dengue.home$month)##
## 05 07 08 09 10 11
## 2 1 42 309 119 5barplot(table(dengue.home$month), xlab = "月份", ylab = "病例數",
main = "登革熱每月病例數（特定區域）")
5-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)
5-5.改用 jitter 的方式畫圖
- 由於經緯度資料的精確度不足，造成大量的資料點重疊，因此改用 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)
5-6.2 週之內新增的病例分佈地圖
- 假設今天是 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, ]5-7.依照時間決定顏色畫圖
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))