BWTE2
Manuscript: Temporal and spatial variation in commercial poultry exposure to an avian influenza natural host.
1 Load Data
1.1 Poultry Data
USDA Census Data (1997-2012): https://quickstats.nass.usda.gov/
USDA.df = read.csv("./Poultry_Census_all.csv", stringsAsFactors=FALSE)
dim(USDA.df)## [1] 49855 21head(USDA.df)## Program Year Period Week.Ending Geo.Level State State.ANSI
## 1 CENSUS 2012 END OF DEC NA COUNTY ALABAMA 1
## 2 CENSUS 2012 END OF DEC NA COUNTY ALABAMA 1
## 3 CENSUS 2012 END OF DEC NA COUNTY ALABAMA 1
## 4 CENSUS 2012 END OF DEC NA COUNTY ALABAMA 1
## 5 CENSUS 2012 END OF DEC NA COUNTY ALABAMA 1
## 6 CENSUS 2012 END OF DEC NA COUNTY ALABAMA 1
## Ag.District Ag.District.Code County County.ANSI Zip.Code Region
## 1 BLACK BELT 40 AUTAUGA 1 NA NA
## 2 BLACK BELT 40 DALLAS 47 NA NA
## 3 BLACK BELT 40 ELMORE 51 NA NA
## 4 BLACK BELT 40 ELMORE 51 NA NA
## 5 BLACK BELT 40 GREENE 63 NA NA
## 6 BLACK BELT 40 HALE 65 NA NA
## watershed_code Watershed Commodity
## 1 0 NA TURKEYS
## 2 0 NA CHICKENS
## 3 0 NA CHICKENS
## 4 0 NA CHICKENS
## 5 0 NA CHICKENS
## 6 0 NA CHICKENS
## Data.Item Domain Domain.Category Value
## 1 TURKEYS - INVENTORY TOTAL NOT SPECIFIED (D)
## 2 CHICKENS, BROILERS - INVENTORY TOTAL NOT SPECIFIED (D)
## 3 CHICKENS, BROILERS - INVENTORY TOTAL NOT SPECIFIED (D)
## 4 CHICKENS, PULLETS, REPLACEMENT - INVENTORY TOTAL NOT SPECIFIED (D)
## 5 CHICKENS, ROOSTERS - INVENTORY TOTAL NOT SPECIFIED (D)
## 6 CHICKENS, PULLETS, REPLACEMENT - INVENTORY TOTAL NOT SPECIFIED (D)
## CV....
## 1 (D)
## 2 (D)
## 3 (D)
## 4 (D)
## 5 (D)
## 6 (D)unique(USDA.df$State)## [1] "ALABAMA" "ALASKA" "ARIZONA" "ARKANSAS"
## [5] "CALIFORNIA" "COLORADO" "CONNECTICUT" "DELAWARE"
## [9] "FLORIDA" "GEORGIA" "HAWAII" "IDAHO"
## [13] "ILLINOIS" "INDIANA" "IOWA" "KANSAS"
## [17] "KENTUCKY" "LOUISIANA" "MAINE" "MARYLAND"
## [21] "MASSACHUSETTS" "MICHIGAN" "MINNESOTA" "MISSISSIPPI"
## [25] "MISSOURI" "MONTANA" "NEBRASKA" "NEVADA"
## [29] "NEW HAMPSHIRE" "NEW JERSEY" "NEW MEXICO" "NEW YORK"
## [33] "NORTH CAROLINA" "NORTH DAKOTA" "OHIO" "OKLAHOMA"
## [37] "OREGON" "PENNSYLVANIA" "RHODE ISLAND" "SOUTH CAROLINA"
## [41] "SOUTH DAKOTA" "TENNESSEE" "TEXAS" "UTAH"
## [45] "VERMONT" "VIRGINIA" "WASHINGTON" "WEST VIRGINIA"
## [49] "WISCONSIN" "WYOMING"unique(USDA.df$Data.Item)## [1] "TURKEYS - INVENTORY"
## [2] "CHICKENS, BROILERS - INVENTORY"
## [3] "CHICKENS, PULLETS, REPLACEMENT - INVENTORY"
## [4] "CHICKENS, ROOSTERS - INVENTORY"
## [5] "DUCKS - INVENTORY"
## [6] "CHICKENS, LAYERS - INVENTORY"Organize Pultry Data First, finding total for each county, poultry class, and year. Second, averaging counties across years. Last, convert each poultry class to a seperate column.
USDA.df$Value2 = ifelse(USDA.df$Value == " (D)", 0, USDA.df$Value)
USDA.df$Value2 = as.numeric(gsub(",","",USDA.df$Value2))
USDA.piv = as.data.frame(
USDA.df %>%
group_by(State.ANSI, County.ANSI, Data.Item, Year) %>%
summarise(InvSum = sum(Value2)) %>%
group_by(State.ANSI, County.ANSI, Data.Item) %>%
summarise(InvSum = mean(InvSum)))
USDA.piv$Data.Item = ifelse(USDA.piv$Data.Item == "TURKEYS - INVENTORY", "Turkey",
ifelse(USDA.piv$Data.Item == "CHICKENS, BROILERS - INVENTORY", "Broilers",
ifelse(USDA.piv$Data.Item == "CHICKENS, PULLETS, REPLACEMENT - INVENTORY", "Pullets",
ifelse(USDA.piv$Data.Item == "CHICKENS, ROOSTERS - INVENTORY", "Roosters",
ifelse(USDA.piv$Data.Item == "DUCKS - INVENTORY", "Ducks",
ifelse(USDA.piv$Data.Item == "CHICKENS, LAYERS - INVENTORY", "Layers", NA))))))
USDA.piv = dcast(USDA.piv, State.ANSI + County.ANSI ~ Data.Item, value.var = "InvSum", fill=FALSE)Join to County Uisng USDA Jurisdicational Boundaries: https://www.nass.usda.gov/Publications/AgCensus/2012/Online_Resources/Ag_Atlas_Maps/
Counties = readOGR(dsn = "./Bounds",
layer = "CoUS_GCS12",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\LabWork\Patuxent\BWTE_2\BWTE2\Bounds", layer: "CoUS_GCS12"
## with 3070 features
## It has 8 fieldsLL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
Counties = spTransform(Counties, LL84)
Counties$State.ANSI = as.integer(substr(Counties$atlas_stco, 1, 2))
Counties$County.ANSI = Counties$cntyn
Counties@data = merge(Counties@data,
USDA.piv, all.x = TRUE,
by = c("State.ANSI", "County.ANSI"))
Counties@data[is.na(Counties@data)] = 0
#State boundaries for mapping
States = readOGR(dsn = "./Bounds",
layer = "StUS_GCS12",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\LabWork\Patuxent\BWTE_2\BWTE2\Bounds", layer: "StUS_GCS12"
## with 48 features
## It has 7 fieldsStates$States = as.integer(States$atlas_st)
States = spTransform(States, LL84)Creating Raster and Point Version for Mapping
tmp.ras = raster(res=0.2)
extent(tmp.ras) = extent(States)
Domain.r = rasterize(States,
tmp.ras,
field = 0,
background = NA)
Grd.pnt = rasterToPoints(Domain.r, sp = TRUE)
Grd.pnt@data = Grd.pnt@data %>%
mutate(Long = Grd.pnt@coords[,1],
Lat = Grd.pnt@coords[,2]) %>%
select(-layer)Poultry Census Maps
USs_df = fortify(States)## Regions defined for each PolygonsCounties@data$id = rownames(Counties@data)
USc_df = fortify(Counties, region = "id")
USc_df = left_join(USc_df, Counties@data, by = "id")
USc_df$brksB = cut(USc_df$Broilers, right = F,
breaks = c(0, 100, 1000, 10000, 100000, 1000000, 10000000, 40000000),
labels = c("[0 - 100)", "[100 - 1 k)", "[1 K - 10 k)", "[10 k - 100 k)", "[100 k - 1 m)",
"[1 m - 10 m)", "[10 m - 40 m)"))
Br.plt = ggplot(USc_df,
aes(long,lat, group=group, fill = brksB)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.5) +
scale_fill_viridis(name="Broiler", discrete=T, direction = -1) +
xlab(" ") +
ylab("Latitude") +
scale_x_longitude(xmin=-140, xmax=-50, step=10) +
scale_y_latitude(ymin=-10, ymax=90, step=5) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x=element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold", colour="transparent"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_blank())
USc_df$brksB = cut(USc_df$Layers, right = F,
breaks= c(0, 100, 1000, 10000, 100000, 1000000, 10000000),
labels = c("[0 - 100)", "[100 - 1 k)", "[1 K - 10 k)", "[10 k - 100 k)", "[100 k - 1 m)",
"[1 m - 10 m)"))
La.plt = ggplot(USc_df,
aes(long,lat, group=group, fill = brksB)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.5) +
scale_fill_viridis(name="Layer", discrete=T, direction = -1) +
xlab(" ") +
ylab("Latitude") +
scale_x_longitude(xmin=-140, xmax=-50, step=10) +
scale_y_latitude(ymin=-10, ymax=90, step=5) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.text.x = element_blank(),
axis.ticks.x=element_blank(),
axis.text.y = element_blank(),
axis.ticks.y=element_blank(),
axis.title.x = element_text(size=16, face="bold", colour="transparent"),
axis.title.y = element_text(size=16, face="bold", colour="transparent"),
plot.title = element_blank())
USc_df$brksB = cut(USc_df$Pullets, right = F,
breaks= c(0, 100, 1000, 10000, 100000, 1000000, 2000000, 3000000),
labels = c("[0 - 100)", "[100 - 1 k)", "[1 K - 10 k)", "[10 k - 100 k)", "[100 k - 1 m)",
"[1 m - 2 m)", "[2 m - 3 m)"))
Pu.plt = ggplot(USc_df,
aes(long,lat, group=group, fill = brksB)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.5) +
scale_fill_viridis(name="Pullet", discrete=T, direction = -1) +
xlab(" ") +
ylab("Latitude") +
scale_x_longitude(xmin=-140, xmax=-50, step=10) +
scale_y_latitude(ymin=-10, ymax=90, step=5) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x=element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold", colour="transparent"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_blank())
USc_df$brksB = cut(USc_df$Roosters, right = F,
breaks= c(0, 100, 1000, 10000, 100000, 300000),
labels = c("[0 - 100)", "[100 - 1 k)", "[1 K - 10 k)", "[10 k - 100 k)", "[100 k - 300 k)"))
Ro.plt = ggplot(USc_df,
aes(long,lat, group=group, fill = brksB)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.5) +
scale_fill_viridis(name="Rooster", discrete=T, direction = -1) +
xlab("Longitude") +
ylab("Latitude") +
scale_x_longitude(xmin=-140, xmax=-50, step=10) +
scale_y_latitude(ymin=-10, ymax=90, step=5) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.text.x = element_blank(),
axis.ticks.x=element_blank(),
axis.text.y = element_blank(),
axis.ticks.y=element_blank(),
axis.title.x = element_text(size=16, face="bold", colour="transparent"),
axis.title.y = element_text(size=16, face="bold", colour="transparent"),
plot.title = element_blank())
USc_df$brksB = cut(USc_df$Turkey, right = F,
breaks= c(0, 100, 1000, 10000, 100000, 1000000, 2000000, 3000000, 4000000),
labels = c("[0 - 100)", "[100 - 1 k)", "[1 K - 10 k)", "[10 k - 100 k)", "[100 k - 1 m)",
"[1 m - 2 m)", "[2 m - 3 m)", "[3 m - 4 m)"))
Tu.plt = ggplot(USc_df,
aes(long,lat, group=group, fill = brksB)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.5) +
scale_fill_viridis(name="Turkey", discrete=T, direction = -1) +
xlab(" ") +
ylab("Latitude") +
scale_x_longitude(xmin=-140, xmax=-50, step=10) +
scale_y_latitude(ymin=-10, ymax=90, step=5) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
#axis.text.x = element_blank(),
axis.ticks.x=element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_blank())
USc_df$brksB = cut(USc_df$Ducks, right = F,
breaks= c(0, 100, 1000, 10000, 100000, 1000000),
labels = c("[0 - 100)", "[100 - 1 k)", "[1 K - 10 k)", "[10 k - 100 k)", "[100 k - 1 m)"))
du.plt = ggplot(USc_df,
aes(long,lat, group=group, fill = brksB)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.5) +
scale_fill_viridis(name="Duck", discrete=T, direction = -1) +
xlab("Longitude") +
ylab(" ") +
scale_x_longitude(xmin=-140, xmax=-50, step=10) +
scale_y_latitude(ymin=-10, ymax=90, step=5) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size=12, face="bold"),
#axis.text.x = element_blank(),
#axis.ticks.x=element_blank(),
axis.text.y = element_blank(),
axis.ticks.y=element_blank(),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold", colour="transparent"),
plot.title = element_blank())
grid.arrange(Br.plt, La.plt, Pu.plt, Ro.plt, Tu.plt, du.plt, ncol = 2)
1.2 ARGOS Telemetry Data
Removing “dead” locations.
Selecting “plausible” records from the DAF_Filter.
Formating date, adding local time (America/Chicago) and getting Hour, Day, Week, season.
Argos_sp = read.csv("C:/Users/humph173/Documents/LabWork/Patuxent/Douglas_080218/BWTeal_2013-2017_Ramey_argos_tabular_diag_gis_filtered.csv",
stringsAsFactors=FALSE,
header = TRUE, sep=",")[,1:35] %>%
filter(Fate == "alive",
DAF_Filter == 0) %>%
mutate(ID = Animal_ID,
ObsDate = as.POSIXct(Location_Datetime_GMT,
tz = "GMT",
format = "%Y/%m/%d %H:%M:%S"),
ObsDateD = as.Date(Location_Datetime_GMT),
LocT = ymd_hms(format(ObsDate, tz="America/Chicago", usetz=TRUE)),
LocHR = hour(LocT),
LocDY = day(LocT),
LocWK = week(LocT),
LocMN = month(LocT),
LocYR = year(LocT))
Argos_sp$Season = ifelse(Argos_sp$LocMN == 12 | Argos_sp$LocMN == 1, "Winter",
ifelse(Argos_sp$LocMN == 3, "Spring",
ifelse(Argos_sp$LocMN >= 4 & Argos_sp$LocMN <= 8, "Breeding",
ifelse(Argos_sp$LocMN >= 9 & Argos_sp$LocMN <= 11, "Fall",
ifelse(Argos_sp$LocMN == 2 & Argos_sp$Day <= 15, "Winter",
ifelse(Argos_sp$LocMN == 2 & Argos_sp$Day >= 16, "Spring",
NA))))))
#Convert to points
Locs = cbind(Argos_sp$Longitude, Argos_sp$Latitude)
Sp.pnts = SpatialPointsDataFrame(Locs, Argos_sp)
proj4string(Sp.pnts) = proj4string(States)
Sp.pnts$US = is.na(over(Sp.pnts, States))[,1]
Sp.pnts = subset(Sp.pnts, US == "FALSE")
nProj = nProj = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0
+x_0=0 +y_0=0 +k=1 +units=km +nadgrids=@null +no_defs"
Sp.pntsp = spTransform(Sp.pnts, nProj)
Countiesp = spTransform(Counties, nProj)
Statesp = spTransform(States, nProj)Int.set = Sp.pnts@data %>%
mutate(Week = LocWK,
ID = Animal_ID,
Source = "Argos") %>%
select(Long, Lat, Week, ID, Source)1.3 eBird Data
Cleaning eBird data using custom “CleanBird” function. Selecting blue-winged teal observations for years 2013-2017.
MnYr = min(Sp.pnts$LocYR)-1
MxYr = max(Sp.pnts$LocYR)+1
CleanBird("C:/Users/humph173/Documents/LabWork/Patuxent/eBird/ebd_buwtea_relFeb-2018/ebd_buwtea_relFeb-2018.txt", "BWTE", MnYr, MxYr)## Parsed with column specification:
## cols(
## .default = col_character(),
## `LAST EDITED DATE` = col_datetime(format = ""),
## `TAXONOMIC ORDER` = col_double(),
## `SUBSPECIES COMMON NAME` = col_logical(),
## `SUBSPECIES SCIENTIFIC NAME` = col_logical(),
## `BREEDING BIRD ATLAS CODE` = col_logical(),
## `BREEDING BIRD ATLAS CATEGORY` = col_logical(),
## `BCR CODE` = col_logical(),
## `USFWS CODE` = col_logical(),
## `ATLAS BLOCK` = col_logical(),
## LATITUDE = col_double(),
## LONGITUDE = col_double(),
## `OBSERVATION DATE` = col_date(format = ""),
## `TIME OBSERVATIONS STARTED` = col_time(format = ""),
## `DURATION MINUTES` = col_double(),
## `EFFORT DISTANCE KM` = col_double(),
## `EFFORT AREA HA` = col_double(),
## `NUMBER OBSERVERS` = col_double(),
## `ALL SPECIES REPORTED` = col_double(),
## `HAS MEDIA` = col_double(),
## APPROVED = col_double()
## # ... with 3 more columns
## )## See spec(...) for full column specifications.eB.set = BWTE %>%
mutate(Week = week(Date),
ID = 1:dim(BWTE)[1],
Source = "eBird") %>%
select(Long, Lat, Week, ID, Source)Join and Convert Observations to Points
BWTE.pnt = rbind(Int.set, eB.set)
BWTE.pnt = SpatialPointsDataFrame(BWTE.pnt[,1:2], BWTE.pnt)
proj4string(BWTE.pnt) = proj4string(States)
BWTE.pnt$Domain = is.na(over(BWTE.pnt, States))[,1]
BWTE.pnt = subset(BWTE.pnt, Domain == FALSE)
BWTE.pntp = spTransform(BWTE.pnt, nProj)1.4 Major Flyways
USFWS Flyways.
Flyways = readOGR(dsn = "./Flyway",
layer = "Flyway",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\LabWork\Patuxent\BWTE_2\BWTE2\Flyway", layer: "Flyway"
## with 4 features
## It has 1 fieldsFlyways = spTransform(Flyways, proj4string(Counties))
Counties$Flyways = over(Counties, Flyways)[,1]
Counties$FlywayI = as.integer(factor(Counties$Flyways))View Flyways
NAmer.reg = readOGR(dsn = "C:/Users/humph173/Documents/LabWork/Patuxent/Patux_level1/Arc/NAmer/States",
layer = "NA_regions",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\LabWork\Patuxent\Patux_level1\Arc\NAmer\States", layer: "NA_regions"
## with 127 features
## It has 36 fields
## Integer64 fields read as strings: FID_NA_reg FID_Cuba FID_ FID_NA_r_1 FID_0 MEXST0_ MEXST0_ID TOTPOP MALE FEMALE FID_NA_sta FID_Bahama OBJECTID TOTPOP_CYNAmer = gUnaryUnion(NAmer.reg)
Paths = readOGR(dsn = "C:/Users/humph173/Documents/LabWork/Patuxent/Patux_level1/LInes",
layer = "Paths",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\LabWork\Patuxent\Patux_level1\LInes", layer: "Paths"
## with 42 features
## It has 2 fields
## Integer64 fields read as strings: Id#PathsCut = crop(Paths, USs_df)
NA2 = fortify(NAmer)
Pth2 = fortify(Paths)
FlywaysC = crop(Flyways, States)
FlywaysC@data$id = rownames(FlywaysC@data)
Flyw_df = fortify(FlywaysC, region = "id")
Flyw_df = left_join(Flyw_df, FlywaysC@data, by = "id")
myCol = c("blue", "yellow", "darkred", "green")
names(myCol) = levels(factor(Flyw_df$NAME))
colScale = scale_fill_manual(name = "Flyway", values = myCol,
labels = c("Atlantic", "Central", "Mississippi", "Pacific"))
Fly.plt = ggplot() +
geom_polygon(data = NA2,
aes(long,lat, group=group),
fill = "gray98", col="darkgray") +
geom_path(data=Pth2,
aes(long, lat,
group=group, fill=NULL),
# col = tmp.df$AI,
col= "red",
size=0.25) +
geom_polygon(data = USs_df,
aes(long,lat, group=group),
fill = "transparent", col="darkgray") +
#ShpScale +
#colScale +
xlab("Longitude") +
ylab("Latitude") +
#ggtitle("Historic A.I. Outbreaks") +
scale_x_longitude(xmin=-150, xmax=80, step=10) +
scale_y_latitude(ymin=20, ymax=60, step=10) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
#legend.direction = "horizontal",
legend.position = c(0.18,0.25),
strip.text = element_text(size=16, face="bold"),
strip.background = element_blank(),
#legend.position=c(0.65, 0.05),
#legend.key.size = unit(1.5,"line"),
legend.text = element_text(size=10, face="bold"),
legend.title = element_text(size=18, face="bold"),
axis.text.x = element_text(size=12, face="bold"),
axis.text.y =element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y =element_text(size=16, face="bold"),
plot.title = element_blank())
Fly.plt + geom_polygon(data = Flyw_df,
aes(long,lat, group=group,
fill = factor(Flyw_df$NAME)),
col ="black", alpha=0.1) + colScale
2 Match to Counties
2.1 Weekly Occurrences
By County. If a point is within the county, setting “OBS” (observations) to 1. This is the full data set (Argos + eBird)
Week.lvls = levels(factor(BWTE.pntp$Week))
for(i in 1:length(Week.lvls)){
Tmp.poly = Countiesp
Tmp.pnt = subset(BWTE.pntp, Week == i)
Tmp.poly$OBS = over(Tmp.poly, Tmp.pnt)[,1]
Tmp.poly$OBS = ifelse(is.na(Tmp.poly$OBS) == FALSE, 1, 0)
Tmp.poly$Week = i
if(i == 1){RT.poly.df = Tmp.poly@data}
else{RT.poly.df = rbind(RT.poly.df, Tmp.poly@data)}
}Weekly Occurrence by County (Argos, Mississippi and central Flyways)
ARGOS data only.
BWTE.pntpM = subset(BWTE.pntp, Source == "Argos")
MissFlyway = subset(Counties, FlywayI == 2 | FlywayI == 3)
MissFlywayp = spTransform(MissFlyway, proj4string(Countiesp))
Week.lvls = levels(factor(BWTE.pntpM$Week))
for(i in 1:length(Week.lvls)){
Tmp.poly = MissFlywayp
Tmp.pnt = subset(BWTE.pntpM, Week == i)
Tmp.poly$OBS = over(Tmp.poly, Tmp.pnt)[,1]
Tmp.poly$OBS = ifelse(is.na(Tmp.poly$OBS) == FALSE, 1, 0)
Tmp.poly$Week = i
if(i == 1){RT.poly.df2 = Tmp.poly@data}
else{RT.poly.df2 = rbind(RT.poly.df2, Tmp.poly@data)}
}Weekly Occurrence by County (eBird, Mississippi and central Flyways)
eBird data only.
BWTE.pntpM = subset(BWTE.pntp, Source == "eBird")
MissFlyway = subset(Counties, FlywayI == 2 | FlywayI == 3)
MissFlywayp = spTransform(MissFlyway, proj4string(Countiesp))
Week.lvls = levels(factor(BWTE.pntpM$Week))
for(i in 1:length(Week.lvls)){
Tmp.poly = MissFlywayp
Tmp.pnt = subset(BWTE.pntpM, Week == i)
Tmp.poly$OBS = over(Tmp.poly, Tmp.pnt)[,1]
Tmp.poly$OBS = ifelse(is.na(Tmp.poly$OBS) == FALSE, 1, 0)
Tmp.poly$Week = i
if(i == 1){RT.poly.df3 = Tmp.poly@data}
else{RT.poly.df3 = rbind(RT.poly.df3, Tmp.poly@data)}
}Weekly Counts by County (Argos + eBird, entire U.S.)
eB.set.cnt = eB.set
eB.set.cnt$Count = 1
BWTE.counts = BWTE %>%
mutate(Source = "eBird",
Week = week(Date),
ID = 1:dim(BWTE)[1]) %>%
select(Long, Lat, Week, ID, Source, Count)
BWTE.counts = rbind(BWTE.counts, eB.set.cnt)
BWTE.counts = SpatialPointsDataFrame(BWTE.counts[,1:2], BWTE.counts)
proj4string(BWTE.counts) = proj4string(States)
BWTE.counts$Domain = is.na(over(BWTE.counts, States))[,1]
BWTE.counts = subset(BWTE.counts, Domain == FALSE)
BWTE.counts$Count[is.na(BWTE.counts$Count)] = 1
BWTE.countsp = spTransform(BWTE.counts, nProj)
Countiesp$Region = Counties$atlas_stco
Week.lvls = levels(factor(BWTE.countsp$Week))
for(i in 1:length(Week.lvls)){
Tmp.poly = Countiesp
Tmp.pnt = subset(BWTE.countsp, Week == i)
Tmp.pnt$ERegion = over(Tmp.pnt, Tmp.poly)[,"Region"]
Count.frm = as.data.frame(
Tmp.pnt@data %>%
group_by(ERegion) %>%
summarise(BSum = sum(Count)))
Tmp.poly$BSum = with(Count.frm,
BSum[match(
Tmp.poly$Region,
ERegion)])
Tmp.poly$BSum[is.na(Tmp.poly$BSum)] = 0
Tmp.poly$Week = i
if(i == 1){count.poly.df = Tmp.poly@data}
else{count.poly.df = rbind(count.poly.df, Tmp.poly@data)}
}Match to Flyways.
RT.poly.df$FlywayI = with(Counties@data,
FlywayI[match(
RT.poly.df$atlas_stco,
atlas_stco)])
count.poly.df$FlywayI = with(Counties@data,
FlywayI[match(
count.poly.df$atlas_stco,
atlas_stco)])
RT.poly.df2$FlywayI = with(MissFlyway@data,
FlywayI[match(
RT.poly.df2$atlas_stco,
atlas_stco)])
RT.poly.df3$FlywayI = with(MissFlyway@data,
FlywayI[match(
RT.poly.df3$atlas_stco,
atlas_stco)])3 Phase I Analysis
Comparing ARGOS and eBird Temporal Signatures across the Central and Mississippi Flyways.
3.1 Divide ARGOS and eBird
Subset Mississippi and Central Flyways and assign spatial and temporal indicies.
RT.poly.df2 = RT.poly.df2 %>% #ARGOS
filter(FlywayI == 2 | FlywayI == 3)
RT.poly.df3 = RT.poly.df3 %>% #eBird
filter(FlywayI == 2 | FlywayI == 3)Add Integer-type County Identifier
MissFlyway$Region = as.integer(factor(MissFlyway$atlas_stco))
RT.poly.df2$Region = as.integer(factor(RT.poly.df2$atlas_stco))
RT.poly.df3$Region = as.integer(factor(RT.poly.df3$atlas_stco))
MFMod.argos = RT.poly.df2
MFMod.ebird = RT.poly.df3 Assign Temporal and Spatial Indices
MFMod.argos$Week1 = MFMod.argos$Week2 = as.integer(as.factor(MFMod.argos$Week))
MFMod.argos$Region1 = MFMod.argos$Region2 = as.integer(as.factor(MFMod.argos$Region))
MFMod.ebird$Week1 = MFMod.ebird$Week2 = as.integer(as.factor(MFMod.ebird$Week))
MFMod.ebird$Region1 = MFMod.ebird$Region2 = as.integer(as.factor(MFMod.ebird$Region))Add Space-Time Interaction Index
MFMod.argos$Region.Wk = paste("ID", MFMod.argos$Region, "W", MFMod.argos$Week, sep="")
MFMod.argos$ID.Region.Wk = as.integer(as.factor(MFMod.argos$Region.Wk))
MFMod.ebird$Region.Wk = paste("ID", MFMod.ebird$Region, "W", MFMod.ebird$Week, sep="")
MFMod.ebird$ID.Region.Wk = as.integer(as.factor(MFMod.ebird$Region.Wk))3.2 Spatial Neighborhood Graph
Identifying each county’s neigbors.
MFnb = poly2nb(MissFlyway,
queen = T)
nb2INLA("JMF", MFnb)
JMF = inla.read.graph("JMF")
summary(MFnb)## Neighbour list object:
## Number of regions: 1930
## Number of nonzero links: 11394
## Percentage nonzero weights: 0.3058874
## Average number of links: 5.903627
## 1 region with no links:
## 2480
## Link number distribution:
##
## 0 1 2 3 4 5 6 7 8 9 10 11
## 1 5 16 53 138 401 733 446 116 18 2 1
## 5 least connected regions:
## 1236 1293 1773 2989 3055 with 1 link
## 1 most connected region:
## 1632 with 11 linksView Spatial Neighborhood Graph
plot(MissFlyway, col='gray', border='blue')
xy = coordinates(MissFlyway)
plot(MFnb, xy, col='red', lwd=0.5, cex=0.25, add=TRUE)
3.3 ARGOS Model
Model Priors and Matricies (ARGOS)
Sr = length(unique(MFMod.argos$Region1))
Tm = length(unique(MFMod.argos$Week1))
g = inla.read.graph("JMF")
Q.xi = matrix(0, g$n, g$n)
for (i in 1:g$n){
Q.xi[i,i]=g$nnbs[[i]]
Q.xi[i,g$nbs[[i]]]=-1}
Q.Leroux = diag(Sr)-Q.xi
lunif = "expression:
a=1; b=1;
beta = exp(theta)/(1+exp(theta));
logdens = lgamma(a+b)-lgamma(a)-lgamma(b)+(a-1)*beta+(b-1)*(1-beta);
log_jacobian = log(beta*(1-beta));
return(logdens+log_jacobian)"
sdunif = "expression: logdens = -log_precision/2;
return(logdens)"
D1 = diff(diag(Tm),differences=1)
Q.gammaRW1 = t(D1)%*%D1
D2 = diff(diag(Tm),differences=2)
Q.gammaRW2 = t(D2)%*%D2Spatial temporal model
formula.T1 = OBS ~ f(Region1,
model="generic1",
Cmatrix=Q.Leroux,
constr=TRUE,
hyper=list(prec=list(prior=sdunif),
beta=list(prior=lunif))) +
f(Week1,
model="rw1",
constr=TRUE,
hyper=list(prec=list(prior=sdunif)))
#theta1 = Model.st$internal.summary.hyperpar$mean
theta1 = c(-2.622310, 1.528891, 1.834132)
Model.stA = inla(formula.T1,
family = "binomial",
verbose = TRUE,
control.predictor = list(
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta1),
control.inla = list(strategy = "gaussian",
int.strategy = "eb"),
control.compute = list(cpo = TRUE, waic=TRUE, dic=TRUE),
data = MFMod.argos)
#save(list=c("MFMod.argos", "Q.Leroux", "g", "sdunif", "lunif"), file="./HPC010219/STFMA2.RData")st.dA = as.data.frame(Model.stA$summary.random$Week1)[,c(1:6)]
names(st.dA) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
st.dA$Data = "ARGOS"3.4 eBird Model
Model Priors and Matricies (eBird)
Sr = length(unique(MFMod.ebird$Region1))
Tm = length(unique(MFMod.ebird$Week1))
g = inla.read.graph("JMF")
Q.xi = matrix(0, g$n, g$n)
for (i in 1:g$n){
Q.xi[i,i]=g$nnbs[[i]]
Q.xi[i,g$nbs[[i]]]=-1}
Q.Leroux = diag(Sr)-Q.xi
lunif = "expression:
a=1; b=1;
beta = exp(theta)/(1+exp(theta));
logdens = lgamma(a+b)-lgamma(a)-lgamma(b)+(a-1)*beta+(b-1)*(1-beta);
log_jacobian = log(beta*(1-beta));
return(logdens+log_jacobian)"
sdunif = "expression: logdens = -log_precision/2;
return(logdens)"
D1 = diff(diag(Tm),differences=1)
Q.gammaRW1 = t(D1)%*%D1
D2 = diff(diag(Tm),differences=2)
Q.gammaRW2 = t(D2)%*%D2Spatial temporal model
formula.T1 = OBS ~ f(Region1,
model="generic1",
Cmatrix=Q.Leroux,
constr=TRUE,
hyper=list(prec=list(prior=sdunif),
beta=list(prior=lunif))) +
f(Week1,
model="rw1",
constr=TRUE,
hyper=list(prec=list(prior=sdunif)))
#theta1 = Model.st$internal.summary.hyperpar$mean
theta1 = c(-2.622310, 1.528891, 1.834132)
Model.stE = inla(formula.T1,
family = "binomial",
verbose = TRUE,
control.predictor = list(
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta1),
control.inla = list(strategy = "gaussian",
int.strategy = "eb"),
control.compute = list(cpo = TRUE, waic=TRUE, dic=TRUE),
data = MFMod.ebird)
#save(list=c("MFMod.ebird", "Q.Leroux", "g", "sdunif", "lunif"), file="./HPC010219/STFME2.RData")st.dE = as.data.frame(Model.stE$summary.random$Week1)[,c(1:6)]
names(st.dE) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
st.dE$Data = "eBird"3.5 Compare
Correlation
cor.test(st.dA$Mean, st.dE$Mean)##
## Pearson's product-moment correlation
##
## data: st.dA$Mean and st.dE$Mean
## t = 7.1114, df = 51, p-value = 3.629e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5379147 0.8195936
## sample estimates:
## cor
## 0.7056155Time Effect
Here, seasons are estimated for each data set seperately by finding where the horizontal line is crossed.
MModels.df = rbind(st.dA, st.dE)
MModels.df$Data2 = ifelse(MModels.df$Data == "ARGOS", "Telemetry",
ifelse(MModels.df$Data == "eBird", "eBird", NA))
XXX = MModels.df %>% filter(Data == "ARGOS")
Myzeros = FindZero(XXX$Mean, XXX$ID)
dateRanges <- data.frame(
from=c(1, Myzeros[1], Myzeros[2], Myzeros[3], Myzeros[4]),
to=c( Myzeros[1], Myzeros[2], Myzeros[3], Myzeros[4], 53),
Season = as.factor(c("Overwinter", "Migration", "Breeding", "Migration", "Overwinter")),
Data2 = "Telemetry")
XXX = MModels.df %>% filter(Data == "eBird")
Myzeros = FindZero(XXX$Mean, XXX$ID)
dateRanges1 <- data.frame(
from=c(1, Myzeros[1], Myzeros[2], Myzeros[3], Myzeros[4]),
to=c( Myzeros[1], Myzeros[2], Myzeros[3], Myzeros[4], 53),
Season = as.factor(c("Overwinter", "Migration", "Breeding", "Migration", "Overwinter")),
Data2 = "eBird")
dateRanges = rbind(dateRanges, dateRanges1)
yCol = c("darkgreen", "tan", "darkred")
names(myCol) = levels(factor(dateRanges$Season))
colScale = scale_fill_manual(name = "Season", values = myCol,
labels = c("Breeding", "Migration", "Overwinter"))
Lplot = ggplot() + geom_rect(data = dateRanges,
aes(xmin = from, xmax = to,
ymin = -Inf, ymax = Inf,
fill = Season), alpha=0.1) +
colScale
Lplot + geom_smooth(data = MModels.df,
aes(ID, Mean),
col = "black",
linetype= "solid",
method = "loess",
span = 0.15,
se = FALSE,
lwd = 1.1) +
geom_smooth(data = MModels.df, aes(ID, Q025),
col = "darkgrey",
method = "loess",
span = 0.15,
se = FALSE,
linetype= "dashed",
lwd = 0.7) +
geom_smooth(data = MModels.df, aes(ID, Q975),
col = "darkgrey",
method = "loess",
span = 0.15,
se = FALSE,
linetype= "dashed",
lwd = 0.7) +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
facet_wrap(~Data2, ncol = 1) +
xlab(" ") +
ylab("Temporal Effect (logit)") +
theme_classic() +
scale_x_continuous(breaks = seq(1, 53, 13),
limits =c(1,53)) +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "bottom",
legend.title=element_text(size=18),
legend.text=element_text(size=16),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
#panel.border = element_rect(colour = "black"),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))## Warning: Removed 3 rows containing missing values (geom_rect).
BWTE Occurrence Probability ARGOS Telemetry Series. Showing every third week.
MFMod.argos$Fit = Model.stA$summary.fitted.values$mean
MFMod.argos$Data = "ARGOS"
Mod.dat10 = MFMod.argos
Mod.dat10 = Mod.dat10 %>% filter(Week <= 53)
Week.levls = seq(1, 53, 3)
for(i in 1:length(Week.levls)){
tmp.wk = Mod.dat10 %>% filter(Week == Week.levls[i])
MissFlyway$Fit = with(tmp.wk,
Fit[match(
MissFlyway$Region,
Region)])
MissFlyway$Week = paste("Week", Week.levls[i], sep = " ")
MissFlyway$WeekOrd = i
MissFlyway@data$id = rownames(MissFlyway@data)
USc_df = fortify(MissFlyway, region = "id")
USc_df = left_join(USc_df, MissFlyway@data, by = "id")
if(i == 1){USc_df2 = USc_df}
else{USc_df2 = rbind(USc_df2, USc_df)}
}
pOrd = arrange(USc_df2, WeekOrd)
USc_df2$Week = factor(USc_df2$Week, levels=unique(pOrd$Week))
MyPick = rev(viridis(100, option = "C"))
StatesM = crop(States, MissFlyway)
USsM_df = fortify(StatesM)## Regions defined for each Polygonsggplot(USc_df2,
aes(long,lat, group=group, fill = Fit)) +
geom_polygon(col="transparent") +
geom_polygon(data=USsM_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.1) +
scale_fill_gradientn(name = "Probability",
colours=MyPick,
breaks=seq(0, 1, 0.5),
limits=c(0,1)) +
xlab("") +
ylab("") +
facet_wrap(~Week, ncol=4) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.position=c(0.8,0.1),
legend.direction = "vertical",
legend.title=element_text(size=18, face="bold"),
legend.text=element_text(size=14),
strip.background = element_blank(),
strip.text = element_text(size=12, face="bold"),
plot.title = element_blank(),
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank())
4 Phase II Analysis
Modeling BWTE Occurrence probabilty across the entire United States.
4.1 Assign spatial Temporal Indicies
Assign Integer-type county identifier and copy data.
Counties$Region = as.integer(factor(Counties$atlas_stco))
RT.poly.df$Region = as.integer(factor(RT.poly.df$atlas_stco))
Mod.dat = RT.poly.df Temporal and Spatial Indices
Mod.dat$Week1 = Mod.dat$Week2 = as.integer(as.factor(Mod.dat$Week))
Mod.dat$Region1 = Mod.dat$Region2 = as.integer(as.factor(Mod.dat$Region))Space-Time Interaction Index
Mod.dat$Region.Wk = paste("ID", Mod.dat$Region, "W", Mod.dat$Week, sep="")
Mod.dat$ID.Region.Wk = as.integer(as.factor(Mod.dat$Region.Wk))4.2 Spatial Neighborhood Graph
Identifying each county’s neigbors.
nb = poly2nb(Counties,
queen = T)
nb2INLA("J", nb)
J = inla.read.graph("J")
summary(nb)## Neighbour list object:
## Number of regions: 3070
## Number of nonzero links: 18016
## Percentage nonzero weights: 0.1911532
## Average number of links: 5.868404
## 3 regions with no links:
## 1190 1833 2908
## Link number distribution:
##
## 0 1 2 3 4 5 6 7 8 9 10 11 13 14
## 3 11 25 80 292 665 1073 657 211 41 9 1 1 1
## 11 least connected regions:
## 177 194 1184 1236 1293 1842 2276 2845 2885 2895 2989 with 1 link
## 1 most connected region:
## 2758 with 14 linksView Spatial Neighborhood Graph
plot(Counties, col='gray', border='blue')
xy = coordinates(Counties)
plot(nb, xy, col='red', lwd=0.5, cex=0.25, add=TRUE)
4.3 Models
Model Priors and Matricies (eBird)
Sr = length(unique(Mod.dat$Region1))
Tm = length(unique(Mod.dat$Week1))
g = inla.read.graph("J")
Q.xi = matrix(0, g$n, g$n)
for (i in 1:g$n){
Q.xi[i,i]=g$nnbs[[i]]
Q.xi[i,g$nbs[[i]]]=-1}
Q.Leroux = diag(Sr)-Q.xi
lunif = "expression:
a=1; b=1;
beta = exp(theta)/(1+exp(theta));
logdens = lgamma(a+b)-lgamma(a)-lgamma(b)+(a-1)*beta+(b-1)*(1-beta);
log_jacobian = log(beta*(1-beta));
return(logdens+log_jacobian)"
sdunif = "expression: logdens = -log_precision/2;
return(logdens)"
D1 = diff(diag(Tm),differences=1)
Q.gammaRW1 = t(D1)%*%D1
D2 = diff(diag(Tm),differences=2)
Q.gammaRW2 = t(D2)%*%D2Spatial Model Only
formula.T1 = OBS ~ f(Region,
model="generic1",
Cmatrix=Q.Leroux,
constr=TRUE,
hyper=list(prec=list(prior=sdunif),
beta=list(prior=lunif)))
#theta1 = Model.space$internal.summary.hyperpar$mean
theta1 = c(-2.147839, 2.057896)
Model.space = inla(formula.T1,
family = "binomial",
verbose = TRUE,
control.predictor = list(
compute = TRUE,
link = 1),
#control.mode = list(restart = TRUE, theta = theta1),
control.inla = list(strategy = "gaussian",
int.strategy = "eb"),
control.compute = list(cpo = TRUE, waic=TRUE, dic=TRUE),
data = Mod.dat)summary(Model.space)##
## Call:
## c("inla(formula = formula.T1, family = \"binomial\", data = Mod.dat, ", " verbose = TRUE, control.compute = list(cpo = TRUE, waic = TRUE, ", " dic = TRUE), control.predictor = list(compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.8752 117.4967 2.6310 122.0029
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -1.8218 0.0144 -1.85 -1.8218 -1.7936 -1.8218 0
##
## Random effects:
## Name Model
## Region Generic1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for Region 0.1168 0.0042 0.1096 0.1164 0.1260 0.1151
## Beta for Region 0.8843 0.0263 0.8229 0.8880 0.9256 0.8958
##
## Expected number of effective parameters(std dev): 2360.85(0.00)
## Number of equivalent replicates : 33.81
##
## Deviance Information Criterion (DIC) ...............: 66181.27
## Deviance Information Criterion (DIC, saturated) ....: 66181.25
## Effective number of parameters .....................: 2452.64
##
## Watanabe-Akaike information criterion (WAIC) ...: 65801.81
## Effective number of parameters .................: 2006.18
##
## Marginal log-Likelihood: -34483.52
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computeddic.m1 = c(Model.space$dic$dic, Model.space$waic$waic)
dic.m1## [1] 66181.27 65801.81Spatial and Temporal Model
formula.T1 = OBS ~ f(Region1,
model="generic1",
Cmatrix=Q.Leroux,
constr=TRUE,
hyper=list(prec=list(prior=sdunif),
beta=list(prior=lunif))) +
f(Week1,
model="rw1",
constr=TRUE,
replicate = FlywayI,
hyper=list(prec=list(prior=sdunif)))
#theta1 = Model.st$internal.summary.hyperpar$mean
theta1 = c(-2.622310, 1.528891, 1.834132)
Model.st = inla(formula.T1,
family = "binomial",
verbose = TRUE,
control.predictor = list(
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta1),
control.inla = list(strategy = "gaussian",
int.strategy = "eb"),
control.compute = list(cpo = TRUE, waic=TRUE, dic=TRUE),
data = Mod.dat)
#save(list=c("Mod.dat", "Q.Leroux", "g", "sdunif", "lunif"), file="./HPC010219/STF.RData") #save to rin on HPCCsummary(Model.st)##
## Call:
## c("inla(formula = formula.T1, family = \"binomial\", data = Mod.dat, ", " verbose = TRUE, control.compute = list(cpo = TRUE, waic = TRUE, ", " dic = TRUE), control.predictor = list(compute = TRUE, ", " link = 1), control.inla = list(strategy = \"gaussian\", ", " int.strategy = \"eb\"), control.mode = list(restart = TRUE, ", " theta = theta1), num.threads = 8)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 3.2894 169.3773 9.8672 182.5339
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -2.7671 0.0159 -2.7983 -2.7671 -2.7359 -2.7671 0
##
## Random effects:
## Name Model
## Region1 Generic1 model
## Week1 RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for Region1 0.0738 0.0067 0.0648 0.0725 0.0899 0.0679
## Beta for Region1 0.8021 0.0600 0.6530 0.8150 0.8828 0.8451
## Precision for Week1 5.3061 0.7999 4.0132 5.2042 7.1411 4.9683
##
## Expected number of effective parameters(std dev): 2746.16(0.00)
## Number of equivalent replicates : 59.25
##
## Deviance Information Criterion (DIC) ...............: 93574.36
## Deviance Information Criterion (DIC, saturated) ....: NULL
## Effective number of parameters .....................: 2884.47
##
## Watanabe-Akaike information criterion (WAIC) ...: 93127.00
## Effective number of parameters .................: 2366.17
##
## Marginal log-Likelihood: -49145.66
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computeddic.m2 = c(Model.st$dic$dic, Model.st$waic$waic)
dic.m2## [1] 93574.36 93127.00st.d = as.data.frame(Model.st$summary.random$Week1)[,c(1:6)]
names(st.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")4.4 Results
Plotting 1 week to check code.
Mod.dat$Fit = Model.st$summary.fitted.values$mean
Mod.dat10 = Mod.dat %>% filter(Week == 1)
Counties$Fit = with(Mod.dat10,
Fit[match(
Counties$Region,
Region)])
Counties@data$id = rownames(Counties@data)
USc_df = fortify(Counties, region = "id")
USc_df = left_join(USc_df, Counties@data, by = "id")
ggplot(USc_df,
aes(long,lat, group=group, fill = Fit)) +
geom_polygon(col="#2b2b2b") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.5) +
scale_fill_viridis(name="Residence", discrete=F,
direction = -1, option = "C") +
xlab("Longitude") +
ylab("Latitude") +
scale_x_longitude(xmin=-140, xmax=-50, step=10) +
scale_y_latitude(ymin=-10, ymax=90, step=10) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_blank())
4.4.1 Fit Estimates
Matching model fitted values by week and county for mapping. Also, comparing estimated values to actual observations.
Mod.dat$Fit = Model.st$summary.fitted.values$mean
Qk.frm = as.data.frame(
Mod.dat %>%
group_by(Week) %>%
summarise(Mean = mean(Fit)))
Threshold.df = Mod.dat %>%
select(OBS, Fit, Week, FlywayI) %>%
mutate(ID = 1:length(OBS))
Threshold.df$FlywayN = ifelse(Mod.dat$FlywayI == 1, "Atlantic",
ifelse(Mod.dat$FlywayI == 2, "Central",
ifelse(Mod.dat$FlywayI == 3, "Mississippi",
ifelse(Mod.dat$FlywayI == 4, "Pacific", NA))))
Mod.dat$FlywayN = ifelse(Mod.dat$FlywayI == 1, "Atlantic",
ifelse(Mod.dat$FlywayI == 2, "Central",
ifelse(Mod.dat$FlywayI == 3, "Mississippi",
ifelse(Mod.dat$FlywayI == 4, "Pacific", NA))))
Flyway.lvls = levels(factor(Threshold.df$FlywayN))
for(i in 1:length(Flyway.lvls)){
tmp.th = Threshold.df %>% filter(FlywayN == Flyway.lvls[i])
th.wks = levels(factor(tmp.th$Week))
for(j in 1:length(th.wks)){
tmp.thw = tmp.th %>% filter(Week == th.wks[j])
NSMod.df = as.data.frame(cbind(tmp.thw$ID, tmp.thw$OBS, tmp.thw$Fit))
OptT.atl = optimal.thresholds(NSMod.df, opt.methods = "MaxSens+Spec")
NSMod.atl = presence.absence.accuracy(NSMod.df, threshold = OptT.atl[1,2])
NSMod.atl$TSS = NSMod.atl$sensitivity + NSMod.atl$specificity - 1
NSMod.atl$MSE = mean((tmp.thw$Fit - tmp.thw$OBS)^2)
NSMod.atl$FlywayN = Flyway.lvls[i]
NSMod.atl$Week = as.integer(th.wks[j])
if(j == th.wks[1]){WeekV.set = NSMod.atl}
else{WeekV.set = rbind(WeekV.set, NSMod.atl)}
}
if(i == 1){Fly.set = WeekV.set}
else{Fly.set = rbind(Fly.set, WeekV.set)}
}
#Match by Flyway and Week
Fly.set.mrg = Fly.set %>% select(FlywayN, Week, threshold)
Mod.dat2 = merge(Mod.dat, Fly.set.mrg, by = c("FlywayN", "Week"))
Mod.dat2$Fit.t = ifelse(Mod.dat2$Fit < Mod.dat2$threshold, 0, Mod.dat2$Fit)Sensitivity and Specificity
myShapes = c(0:3)
names(myShapes) = levels(factor(Fly.set$FlywayN))
ShpScale = scale_shape_manual(name = "Flyway", values = myShapes,
labels = c("Atlantic", "Central", "Mississippi", "Pacific"))
SenSpec = ggplot(Fly.set, aes(sensitivity, specificity, col = Week, shape=factor(Fly.set$FlywayN))) +
geom_point(size = 4) +
ShpScale +
scale_color_gradient(low="blue", high="red",
breaks = seq(0, 50, 10),
labels = paste(c(1, seq(10,50,10))),
limit = c(0,53)) +
xlim(0, 1) +
ylim(0, 1) +
xlab("Sensitivity") +
ylab("Specificity") +
theme_bw() +
theme(legend.position = c(0.125, 0.5),
legend.title=element_text(size=20),
legend.text=element_text(size=16),
panel.grid.major = element_line(colour="darkgray", size=0.25),
panel.grid.minor = element_line(colour="lightgray", size=0.1),
strip.background = element_blank(),
#panel.border = element_rect(colour = "black"),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
SenSpec2 = ggplot(Fly.set, aes(AUC, MSE, col = Week, shape=factor(Fly.set$FlywayN))) +
geom_point(size = 3) +
ShpScale +
scale_color_gradient(low="blue", high="red",
breaks = seq(0, 50, 10),
labels = paste(c(1, seq(10,50,10))),
limit = c(0,53)) +
xlim(0, 1) +
ylim(0, 1) +
xlab("AUC Receiver Operating Characteristics") +
ylab("Mean Squared Error") +
theme_bw() +
theme(legend.position = "none",
legend.title=element_blank(),
legend.text=element_blank(),
panel.grid.major = element_line(colour="darkgray", size=0.25),
panel.grid.minor = element_line(colour="lightgray", size=0.1),
strip.background = element_blank(),
#panel.border = element_rect(colour = "black"),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
grid.arrange(SenSpec, SenSpec2, ncol=1)
Map Series
Showing estimated BWTE occurrence probability for every third week.
Mod.dat10 = Mod.dat2 %>% filter(Week <= 53)
Week.levls = seq(1, 53, 3)
for(i in 1:length(Week.levls)){
tmp.wk = Mod.dat10 %>% filter(Week == Week.levls[i])
Counties$Fitm = with(tmp.wk,
Fit[match(
Counties$Region,
Region)])
Counties$Week = paste("Week", Week.levls[i], sep = " ")
Counties$WeekOrd = i
Counties@data$id = rownames(Counties@data)
USc_df = fortify(Counties, region = "id")
USc_df = left_join(USc_df, Counties@data, by = "id")
if(i == 1){USc_df2 = USc_df}
else{USc_df2 = rbind(USc_df2, USc_df)}
}
pOrd = arrange(USc_df2, WeekOrd)
USc_df2$Week = factor(USc_df2$Week, levels=unique(pOrd$Week))
USc_df2 %>%
group_by(Week) %>%
summarise(Mean = mean(Fitm))## # A tibble: 18 x 2
## Week Mean
## <fct> <dbl>
## 1 Week 1 0.0775
## 2 Week 4 0.0677
## 3 Week 7 0.0844
## 4 Week 10 0.178
## 5 Week 13 0.430
## 6 Week 16 0.482
## 7 Week 19 0.371
## 8 Week 22 0.223
## 9 Week 25 0.140
## 10 Week 28 0.0948
## 11 Week 31 0.109
## 12 Week 34 0.223
## 13 Week 37 0.346
## 14 Week 40 0.289
## 15 Week 43 0.205
## 16 Week 46 0.0882
## 17 Week 49 0.0721
## 18 Week 52 0.0654ggplot(USc_df2,
aes(long,lat, group=group, fill = Fitm)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.1) +
scale_fill_gradientn(name = "Occurrence Probability",
colours=MyPick,
breaks=seq(0, 1, 0.5),
limits=c(0,1)) +
xlab("") +
ylab("") +
facet_wrap(~Week, ncol=4) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.position=c(0.7,0.1),
legend.direction = "horizontal",
legend.title=element_text(size=20, face="bold"),
legend.text=element_text(size=16),
strip.background = element_blank(),
strip.text = element_text(size=12, face="bold"),
plot.title = element_blank(),
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank())
4.5 Domestic Poultry Exposure
Exposure = BWTE Occurrence Probability x COmmercial Poultry Abundance.
Broiler Exposure. Showing Every three weeks.
Mod.dat10 = Mod.dat2 %>% filter(Week <= 53)
Mod.dat10$Broiler = Mod.dat10$Fit.t*Mod.dat10$Broilers
Mod.dat10$Broiler = Mod.dat10$Broiler/1000000
Week.levls = seq(1, 53, 3)
for(i in 1:length(Week.levls)){
tmp.wk = Mod.dat10 %>% filter(Week == Week.levls[i])
Counties$Broiler = with(tmp.wk,
Broiler[match(
Counties$Region,
Region)])
Counties$Week = paste("Week", Week.levls[i], sep = " ")
Counties$WeekOrd = i
Counties@data$id = rownames(Counties@data)
USc_df = fortify(Counties, region = "id")
USc_df = left_join(USc_df, Counties@data, by = "id")
if(i == 1){USc_df2 = USc_df}
else{USc_df2 = rbind(USc_df2, USc_df)}
}
pOrd = arrange(USc_df2, WeekOrd)
USc_df2$Week = factor(USc_df2$Week, levels=unique(pOrd$Week))
USc_df2$Broiler.cut = USc_df2$Broiler
USc_df2$Broiler.cut[USc_df2$Broiler.cut >= 1] = 1
USc_df2$Broiler.cut = ifelse(USc_df2$Broiler.cut == 0, 0, ceiling(USc_df2$Broiler.cut*100000)/100000)
UniqueVals = length(unique(USc_df2$Broiler.cut))
MyPick = rev(viridis(UniqueVals, option = "C"))
MyPick = MyPick[25:length(MyPick)]
USc_df2$brksB = cut(USc_df2$Broiler.cut, right = F,
breaks = c(0, 0.0001, 0.001, 0.01, 0.1, 1, 10, 20, 40),
labels = c("[0 - 100)", "[100 - 1 k)", "[1 K - 10 k)", "[10 k - 100 k)", "[100 k - 1 m)",
"[1 m - 10 m)", "[10 m - 20 m)", "[20 m - 40 m)"))
ggplot(USc_df2,
aes(long,lat, group=group, fill = brksB)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.1) +
scale_fill_viridis(name="Broiler (millions)", discrete=T,
direction = -1, option = "C") +
xlab("") +
ylab("") +
facet_wrap(~Week, ncol=4) +
#scale_x_longitude(xmin=-140, xmax=-50, step=10) +
#scale_y_latitude(ymin=-10, ymax=90, step=10) +
#coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
#legend.position = "bottom",
legend.position=c(0.7,0.1),
legend.direction = "horizontal",
legend.title=element_text(size=20, face="bold"),
legend.text=element_text(size=16),
strip.background = element_blank(),
strip.text = element_text(size=12, face="bold"),
plot.title = element_blank(),
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank())
MS Figure Broiler
#Get
Mod.dat10 = Mod.dat2 %>% filter(Week <= 53)
Mod.dat10$Broiler = Mod.dat10$Fit.t*Mod.dat10$Broilers
Mod.dat10$Broiler = Mod.dat10$Broiler/1000000
Week.levls = c(17, 29, 38, 50)
for(i in 1:length(Week.levls)){
tmp.wk = Mod.dat10 %>% filter(Week == Week.levls[i])
Counties$Broiler = with(tmp.wk,
Broiler[match(
Counties$Region,
Region)])
Counties$Fitt = with(tmp.wk,
Fit[match(
Counties$Region,
Region)])
Counties$Week = paste("Week", Week.levls[i], sep = " ")
Counties$WeekOrd = i
Counties@data$id = rownames(Counties@data)
USc_df = fortify(Counties, region = "id")
USc_df = left_join(USc_df, Counties@data, by = "id")
if(i == 1){USc_df2 = USc_df}
else{USc_df2 = rbind(USc_df2, USc_df)}
}
pOrd = arrange(USc_df2, WeekOrd)
USc_df2$Week = factor(USc_df2$Week, levels=unique(pOrd$Week))
USc_df2$Broiler.cut = USc_df2$Broiler
USc_df2$Broiler.cut[USc_df2$Broiler.cut >= 10] = 10
USc_df2$Broiler.cut = ifelse(USc_df2$Broiler.cut == 0, 0, ceiling(USc_df2$Broiler.cut*100000)/100000)
UniqueVals = length(unique(USc_df2$Broiler.cut))
MyPick = rev(viridis(UniqueVals, option = "C"))
MyPick = c(MyPick[1:2], MyPick[8:length(MyPick)])
USc_df2WK1 = USc_df2 %>% filter(Week == "Week 17")
USc_df2WK2 = USc_df2 %>% filter(Week == "Week 29")
USc_df2WK3 = USc_df2 %>% filter(Week == "Week 38")
USc_df2WK4 = USc_df2 %>% filter(Week == "Week 50")
wk1 = ggplot(USc_df2,
aes(long,lat, group=group, fill = Broiler.cut)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.1) +
scale_fill_gradientn(name = "Broiler",
colours=MyPick,
breaks=seq(0, 10, 5),
labels = c("0", "5", " 10+"),
limits=c(0,10)) +
xlab("") +
ylab("") +
facet_wrap(~Week, ncol=1, strip.position = "left") +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.position = "bottom",
#legend.position=c(0.7,0.1),
legend.direction = "horizontal",
legend.title=element_text(size=20, face="bold"),
legend.text=element_text(size=16),
strip.background = element_blank(),
strip.text = element_text(size=18, face="bold", colour="transparent"),
plot.title = element_blank(),
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank())
##Raw
Counties@data$id = rownames(Counties@data)
USc_df = fortify(Counties, region = "id")
USc_df = left_join(USc_df, Counties@data, by = "id")
USc_df$brksB = cut(USc_df$Broilers, right = F,
breaks = c(0, 100, 1000, 10000, 100000, 1000000, 10000000, 20000000, 40000000),
labels = c("[0 - 100)", "[100 - 1 k)", "[1 K - 10 k)", "[10 k - 100 k)", "[100 k - 1 m)",
"[1 m - 10 m)", "[10 m - 20 m)", "[20 m - 40 m)"))
Br.plt = ggplot(USc_df,
aes(long,lat, group=group, fill = brksB)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.5) +
scale_fill_viridis(name="", discrete=T, direction = -1) +
xlab("Longitude") +
ylab("Latitude") +
scale_x_longitude(xmin=-140, xmax=-50, step=10) +
scale_y_latitude(ymin=-10, ymax=90, step=10) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.position = "bottom",
legend.direction = "vertical",
legend.title = element_text(size=12, face="bold"),
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
plot.title = element_blank())
##Fitt
MyPick = rev(viridis(UniqueVals, option = "E"))
wk2 = ggplot(USc_df2,
aes(long,lat, group=group, fill = Fitt)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.1) +
scale_fill_gradientn(name = "Probability",
colours=MyPick,
breaks=seq(0, 1, 0.5),
#labels = c("0", "2", "4", "6", "8", " 10+"),
limits=c(0,1)) +
xlab("") +
ylab("") +
facet_wrap(~Week, ncol=1, strip.position = "left") +
#scale_x_longitude(xmin=-140, xmax=-50, step=10) +
#scale_y_latitude(ymin=-10, ymax=90, step=10) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.position = "bottom",
#legend.position=c(0.7,0.1),
#legend.direction = "horizontal",
legend.title=element_text(size=20, face="bold"),
legend.text=element_text(size=16),
strip.background = element_blank(),
strip.text = element_text(size=18, face="bold"),
plot.title = element_blank(),
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank())
grid.arrange(wk2, Br.plt, wk1, ncol=3)
Turkey Series
Mod.dat10 = Mod.dat2 %>% filter(Week <= 53)
Mod.dat10$TurkeyExp = Mod.dat10$Fit.t*Mod.dat10$Turkey
Mod.dat10$TurkeyExp = Mod.dat10$TurkeyExp/1000000
Week.levls = seq(1, 53, 3)
for(i in 1:length(Week.levls)){
tmp.wk = Mod.dat10 %>% filter(Week == Week.levls[i])
Counties$TurkeyExp = with(tmp.wk,
TurkeyExp[match(
Counties$Region,
Region)])
Counties$Week = paste("Week", Week.levls[i], sep = " ")
Counties$WeekOrd = i
Counties@data$id = rownames(Counties@data)
USc_df = fortify(Counties, region = "id")
USc_df = left_join(USc_df, Counties@data, by = "id")
if(i == 1){USc_df2 = USc_df}
else{USc_df2 = rbind(USc_df2, USc_df)}
}
pOrd = arrange(USc_df2, WeekOrd)
USc_df2$Week = factor(USc_df2$Week, levels=unique(pOrd$Week))
USc_df2 %>%
group_by(Week) %>%
summarise(Mean = mean(TurkeyExp))## # A tibble: 18 x 2
## Week Mean
## <fct> <dbl>
## 1 Week 1 0.000854
## 2 Week 4 0.000728
## 3 Week 7 0.000867
## 4 Week 10 0.00369
## 5 Week 13 0.0105
## 6 Week 16 0.0119
## 7 Week 19 0.00888
## 8 Week 22 0.00390
## 9 Week 25 0.00192
## 10 Week 28 0.00188
## 11 Week 31 0.00214
## 12 Week 34 0.00475
## 13 Week 37 0.00810
## 14 Week 40 0.00648
## 15 Week 43 0.00407
## 16 Week 46 0.000929
## 17 Week 49 0.000665
## 18 Week 52 0.000648USc_df2$TurkeyExp.cut = USc_df2$TurkeyExp
USc_df2$TurkeyExp.cut[USc_df2$TurkeyExp.cut >= 2] = 2
USc_df2$TurkeyExp.cut = ifelse(USc_df2$TurkeyExp.cut == 0, 0, ceiling(USc_df2$TurkeyExp.cut*100000)/100000)
UniqueVals = length(unique(USc_df2$TurkeyExp.cut))
MyPick = rev(viridis(UniqueVals, option = "C"))
MyPick = MyPick[25:length(MyPick)]
ggplot(USc_df2,
aes(long,lat, group=group, fill = TurkeyExp)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.1) +
scale_fill_gradientn(name = "Turkey (millions)",
colours=MyPick,
breaks=seq(0, 2, 1),
labels = c("0", "1", " 2+"),
limits=c(0,2)) +
xlab("") +
ylab("") +
facet_wrap(~Week, ncol=4) +
#scale_x_longitude(xmin=-140, xmax=-50, step=10) +
#scale_y_latitude(ymin=-10, ymax=90, step=10) +
#coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
#legend.position = "bottom",
legend.position=c(0.7,0.1),
legend.direction = "horizontal",
legend.title=element_text(size=20, face="bold"),
legend.text=element_text(size=16),
strip.background = element_blank(),
strip.text = element_text(size=12, face="bold"),
plot.title = element_blank(),
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank())
Layer Series
Mod.dat10 = Mod.dat2 %>% filter(Week <= 53)
Mod.dat10$Broiler = Mod.dat10$Fit.t*Mod.dat10$Layers
Mod.dat10$Broiler = Mod.dat10$Broiler/1000000
Week.levls = seq(1, 53, 3)
for(i in 1:length(Week.levls)){
tmp.wk = Mod.dat10 %>% filter(Week == Week.levls[i])
Counties$Broiler = with(tmp.wk,
Broiler[match(
Counties$Region,
Region)])
Counties$Week = paste("Week", Week.levls[i], sep = " ")
Counties$WeekOrd = i
Counties@data$id = rownames(Counties@data)
USc_df = fortify(Counties, region = "id")
USc_df = left_join(USc_df, Counties@data, by = "id")
if(i == 1){USc_df2 = USc_df}
else{USc_df2 = rbind(USc_df2, USc_df)}
}
pOrd = arrange(USc_df2, WeekOrd)
USc_df2$Week = factor(USc_df2$Week, levels=unique(pOrd$Week))
USc_df2$Broiler.cut = USc_df2$Broiler
USc_df2$Broiler.cut[USc_df2$Broiler.cut >= 5] = 5
USc_df2$Broiler.cut = ifelse(USc_df2$Broiler.cut == 0, 0, ceiling(USc_df2$Broiler.cut*100000)/100000)
UniqueVals = length(unique(USc_df2$Broiler.cut))
MyPick = rev(viridis(UniqueVals, option = "C"))
MyPick = MyPick[25:length(MyPick)]
ggplot(USc_df2,
aes(long,lat, group=group, fill = Broiler.cut)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.1) +
scale_fill_gradientn(name = "Layer (millions)",
colours=MyPick,
breaks=seq(0, 5, 2.5),
labels = c("0", "2.5", " 5+"),
limits=c(0,5)) +
xlab("") +
ylab("") +
facet_wrap(~Week, ncol=4) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
#legend.position = "bottom",
legend.position=c(0.7,0.1),
legend.direction = "horizontal",
legend.title=element_text(size=20, face="bold"),
legend.text=element_text(size=16),
strip.background = element_blank(),
strip.text = element_text(size=12, face="bold"),
plot.title = element_blank(),
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank())
Pullet Series
Mod.dat10 = Mod.dat2 %>% filter(Week <= 53)
Mod.dat10$Broiler = Mod.dat10$Fit.t*Mod.dat10$Pullets
Mod.dat10$Broiler = Mod.dat10$Broiler/1000000
Week.levls = seq(1, 53, 3)
for(i in 1:length(Week.levls)){
tmp.wk = Mod.dat10 %>% filter(Week == Week.levls[i])
Counties$Broiler = with(tmp.wk,
Broiler[match(
Counties$Region,
Region)])
Counties$Week = paste("Week", Week.levls[i], sep = " ")
Counties$WeekOrd = i
Counties@data$id = rownames(Counties@data)
USc_df = fortify(Counties, region = "id")
USc_df = left_join(USc_df, Counties@data, by = "id")
if(i == 1){USc_df2 = USc_df}
else{USc_df2 = rbind(USc_df2, USc_df)}
}
pOrd = arrange(USc_df2, WeekOrd)
USc_df2$Week = factor(USc_df2$Week, levels=unique(pOrd$Week))
USc_df2$Broiler.cut = USc_df2$Broiler
USc_df2$Broiler.cut[USc_df2$Broiler.cut >= 2] = 2
USc_df2$Broiler.cut = ifelse(USc_df2$Broiler.cut == 0, 0, ceiling(USc_df2$Broiler.cut*100000)/100000)
UniqueVals = length(unique(USc_df2$Broiler.cut))
MyPick = rev(viridis(UniqueVals, option = "C"))
MyPick = MyPick[25:length(MyPick)]
ggplot(USc_df2,
aes(long,lat, group=group, fill = Broiler.cut)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.1) +
scale_fill_gradientn(name = "Pullet (millions)",
colours=MyPick,
breaks=seq(0, 2, 1),
labels = c("0", "1", " 2+"),
limits=c(0,2)) +
xlab("") +
ylab("") +
facet_wrap(~Week, ncol=4) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
#legend.position = "bottom",
legend.position=c(0.7,0.1),
legend.direction = "horizontal",
legend.title=element_text(size=20, face="bold"),
legend.text=element_text(size=16),
strip.background = element_blank(),
strip.text = element_text(size=12, face="bold"),
plot.title = element_blank(),
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank())
Rooster Series
Mod.dat10 = Mod.dat2 %>% filter(Week <= 53)
Mod.dat10$Broiler = Mod.dat10$Fit.t*Mod.dat10$Roosters
Mod.dat10$Broiler = Mod.dat10$Broiler/10000
Week.levls = seq(1, 53, 3)
for(i in 1:length(Week.levls)){
tmp.wk = Mod.dat10 %>% filter(Week == Week.levls[i])
Counties$Broiler = with(tmp.wk,
Broiler[match(
Counties$Region,
Region)])
Counties$Week = paste("Week", Week.levls[i], sep = " ")
Counties$WeekOrd = i
Counties@data$id = rownames(Counties@data)
USc_df = fortify(Counties, region = "id")
USc_df = left_join(USc_df, Counties@data, by = "id")
if(i == 1){USc_df2 = USc_df}
else{USc_df2 = rbind(USc_df2, USc_df)}
}
pOrd = arrange(USc_df2, WeekOrd)
USc_df2$Week = factor(USc_df2$Week, levels=unique(pOrd$Week))
USc_df2$Broiler.cut = USc_df2$Broiler
USc_df2$Broiler.cut[USc_df2$Broiler.cut >= 8] = 8
USc_df2$Broiler.cut = ifelse(USc_df2$Broiler.cut == 0, 0, ceiling(USc_df2$Broiler.cut*100000)/100000)
UniqueVals = length(unique(USc_df2$Broiler.cut))
MyPick = rev(viridis(UniqueVals, option = "C"))
MyPick = MyPick[25:length(MyPick)]
ggplot(USc_df2,
aes(long,lat, group=group, fill = Broiler.cut)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.1) +
scale_fill_gradientn(name = "Rooster (thousands)",
colours=MyPick,
breaks=seq(0, 8, 4),
labels = c("0", "40", " 80+"),
limits=c(0,8)) +
xlab("") +
ylab("") +
facet_wrap(~Week, ncol=4) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.position=c(0.7,0.1),
legend.direction = "horizontal",
legend.title=element_text(size=20, face="bold"),
legend.text=element_text(size=16),
strip.background = element_blank(),
strip.text = element_text(size=12, face="bold"),
plot.title = element_blank(),
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank())
Duck Series
Mod.dat10 = Mod.dat2 %>% filter(Week <= 53)
Mod.dat10$Broiler = Mod.dat10$Fit.t*Mod.dat10$Ducks
Mod.dat10$Broiler = Mod.dat10$Broiler/100000
Week.levls = seq(1, 53, 3)
for(i in 1:length(Week.levls)){
tmp.wk = Mod.dat10 %>% filter(Week == Week.levls[i])
Counties$Broiler = with(tmp.wk,
Broiler[match(
Counties$Region,
Region)])
Counties$Week = paste("Week", Week.levls[i], sep = " ")
Counties$WeekOrd = i
Counties@data$id = rownames(Counties@data)
USc_df = fortify(Counties, region = "id")
USc_df = left_join(USc_df, Counties@data, by = "id")
if(i == 1){USc_df2 = USc_df}
else{USc_df2 = rbind(USc_df2, USc_df)}
}
pOrd = arrange(USc_df2, WeekOrd)
USc_df2$Week = factor(USc_df2$Week, levels=unique(pOrd$Week))
USc_df2$Broiler.cut = USc_df2$Broiler
USc_df2$Broiler.cut[USc_df2$Broiler.cut >= 0.5] = 0.5
USc_df2$Broiler.cut = ifelse(USc_df2$Broiler.cut == 0, 0, ceiling(USc_df2$Broiler.cut*100000)/100000)
UniqueVals = length(unique(USc_df2$Broiler.cut))
MyPick = rev(viridis(UniqueVals, option = "C"))
MyPick = MyPick[25:length(MyPick)]
ggplot(USc_df2,
aes(long,lat, group=group, fill = Broiler.cut)) +
geom_polygon(col="transparent") +
geom_polygon(data=USs_df,
aes(long, lat,
group=group), fill="transparent",
col = "white", size=0.1) +
scale_fill_gradientn(name = "Duck (thousands)",
colours=MyPick,
breaks=seq(0, 0.5, 0.25),
labels = c("0", "25", "50+"),
limits=c(0,0.5)) +
xlab("") +
ylab("") +
facet_wrap(~Week, ncol=4) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.position=c(0.7,0.1),
legend.direction = "horizontal",
legend.title=element_text(size=20, face="bold"),
legend.text=element_text(size=16),
strip.background = element_blank(),
strip.text = element_text(size=12, face="bold"),
plot.title = element_blank(),
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank())
4.6 Flyway Time Variation
Seasonality threshold assigned with respect to horizontal line at 0. Highlighting differences in estimating migration timing.
mic.d = as.data.frame(Model.st$summary.random$Week1)[,c(1:6)]
names(mic.d) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
mic.d$Flyway = rep(c("Atlantic Flyway", "Central Flyway", "Mississippi Flyway", "Pacific Flyway"), each = 53)
XXX = mic.d %>% filter(Flyway == "Atlantic Flyway")
Myzeros = FindZero(XXX$Mean, XXX$ID)
dateRanges <- data.frame(
from=c(1, Myzeros[1], Myzeros[2], Myzeros[3], Myzeros[4]),
to=c( Myzeros[1], Myzeros[2], Myzeros[3], Myzeros[4], 53),
Season = as.factor(c("Overwinter", "Migration", "Breeding", "Migration", "Overwinter")),
Flyway = "Atlantic Flyway")
XXX = mic.d %>% filter(Flyway == "Central Flyway")
Myzeros = FindZero(XXX$Mean, XXX$ID)
dateRanges1 <- data.frame(
from=c(1, Myzeros[2], Myzeros[3], Myzeros[4], Myzeros[5]),
to=c(Myzeros[2], Myzeros[3], Myzeros[4], Myzeros[5], 53),
Season = as.factor(c("Overwinter", "Migration", "Breeding", "Migration", "Overwinter")),
Flyway = "Central Flyway")
XXX = mic.d %>% filter(Flyway == "Mississippi Flyway")
Myzeros = FindZero(XXX$Mean, XXX$ID)
dateRanges2 <- data.frame(
from=c(1, Myzeros[1], Myzeros[2], Myzeros[3], Myzeros[4]),
to=c( Myzeros[1], Myzeros[2], Myzeros[3], Myzeros[4], 53),
Season = as.factor(c("Overwinter", "Migration", "Breeding", "Migration", "Overwinter")),
Flyway = "Mississippi Flyway")
XXX = mic.d %>% filter(Flyway == "Pacific Flyway")
Myzeros = FindZero(XXX$Mean, XXX$ID)
dateRanges3 <- data.frame(
from=c(1, Myzeros[1], Myzeros[2], Myzeros[3], Myzeros[4]),
to=c( Myzeros[1], Myzeros[2], Myzeros[3], Myzeros[4], 53),
Season = as.factor(c("Overwinter", "Migration", "Breeding", "Migration", "Overwinter")),
Flyway = "Pacific Flyway")
dateRanges = rbind(dateRanges, dateRanges1, dateRanges2, dateRanges3)
myCol = c("darkgreen", "tan", "darkred")
names(myCol) = levels(factor(dateRanges$Season))
colScale = scale_fill_manual(name = "Season", values = myCol,
labels = c("Breeding", "Migration", "Overwinter"))
Lplot = ggplot() + geom_rect(data = dateRanges,
aes(xmin = from, xmax = to,
ymin = -Inf, ymax = Inf,
fill = Season), alpha=0.1) +
colScale
Lplot + geom_smooth(data = mic.d,
aes(mic.d$ID, mic.d$Mean),
col = "black",
linetype= "solid",
method = "loess",
span = 0.15,
se = FALSE,
lwd = 1.1) +
geom_smooth(data = mic.d, aes(mic.d$ID, mic.d$Q025),
col = "darkgrey",
method = "loess",
span = 0.15,
se = FALSE,
linetype= "dashed",
lwd = 0.7) +
geom_smooth(data = mic.d, aes(mic.d$ID, mic.d$Q975),
col = "darkgrey",
method = "loess",
span = 0.15,
se = FALSE,
linetype= "dashed",
lwd = 0.7) +
geom_hline(yintercept = 0,
linetype = "solid",
col = "darkgray",
size = 0.5) +
facet_wrap(~Flyway, ncol=2) +
xlab("Week of Year ") +
ylab("Temporal Effect (logit)") +
theme_classic() +
#scale_y_continuous(breaks = seq(-0.5, 1.75, 0.5)) +
scale_x_continuous(breaks = seq(1, 53, 13),
limits =c(1,53)) +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "bottom",
legend.title=element_text(size=18),
legend.text=element_text(size=16),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
#panel.border = element_rect(colour = "black"),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
4.7 Exposure Line Graphs
Exposure by flyway, time, and poultry class.
Mod.dat$Fit = Model.st$summary.fitted.values$mean
Plt.set = Mod.dat
Plt.set$FlywayN = ifelse(Mod.dat$FlywayI == 1, "Atlantic",
ifelse(Mod.dat$FlywayI == 2, "Central",
ifelse(Mod.dat$FlywayI == 3, "Mississippi",
ifelse(Mod.dat$FlywayI == 4, "Pacific", NA))))
Plt.set2 =as.data.frame(
Plt.set %>%
group_by(Week, FlywayN, County.ANSI) %>%
summarise(Prob = mean(Fit),
Brol = Prob*sum(Broilers),
Lay = Prob*sum(Layers),
Pull = Prob*sum(Pullets),
Roo = Prob*sum(Roosters),
Turk = Prob*sum(Turkey),
Duck = Prob*sum(Ducks)) %>%
group_by(Week, FlywayN) %>%
summarise(Brol = sum(Brol)/1000000,
Lay = sum(Lay)/1000000,
Pull = sum(Pull)/1000000,
Roo = sum(Roo)/1000000,
Turk = sum(Turk)/1000000,
Duck = sum(Duck)/1000000))
Plt.set2m = melt(Plt.set2, c("Week","FlywayN"))
Tmp.frm = as.data.frame(
Plt.set2m %>%
group_by(FlywayN, variable) %>%
summarise(Mn = min(value),
Mx = max(value)))
levels(Plt.set2m$variable) = c("Broiler", "Layer", "Pullet", "Rooster", "Turkey", "Duck")
Plt.set2m$Flyway = Plt.set2m$FlywayN
ggplot(Plt.set2m, aes(Week, value,
group = Flyway,
col = Flyway)) +
geom_line(lwd=0.5) +
facet_wrap(~variable, ncol=3, scales = "free_y") +
scale_y_continuous(breaks = scales::pretty_breaks(8), limits = c(0, NA)) +
theme_classic() +
ylab("Commercial Poultry Head (millions)") +
xlab("Week of Year") +
scale_x_continuous(breaks = seq(1, 53, 13),
limits =c(1,53)) +
scale_color_manual(values=c("#332288", "#999933", "#882255", "#CC6677")) +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "bottom",
legend.title=element_text(size=18),
legend.text=element_text(size=16),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
#panel.border = element_rect(colour = "black"),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22))
4.8 AI Surveilance
Load IRD Data and then compare to estimated exposure by flyway.
IRD.df = read.csv("./IRD/IRD_BWTE_021219.csv", stringsAsFactors=FALSE)
IRD2.df = IRD.df %>%
mutate(Long = as.numeric(Longitude),
Lat = as.numeric(Latitude),
ObsDate = ymd(Collection.Date),
Year = as.integer(substr(ObsDate,0, 4)),
Month = as.integer(format(ObsDate, "%m")),
Week = as.integer(strftime(ObsDate, format = "%V")),
Day = as.integer(format(ObsDate, "%d")),
AI = Subtype,
Test = Flu.Test.Status) %>%
filter(is.na(Long) == FALSE,
is.na(Lat) == FALSE,
is.na(Collection.Date) == FALSE,
Country == "USA") %>%
select(Long, Lat, AI, ObsDate, Year, Month, Week, Day, Test)
range(IRD2.df$Year)## [1] 1986 2018Em.pnts = SpatialPointsDataFrame(IRD2.df[,c("Long", "Lat")], IRD2.df)
proj4string(Em.pnts) = proj4string(States)
Em.pnts$State = extract(States, Em.pnts)[,"atlas_name"]
Em.pnts$County = extract(Counties, Em.pnts)[,"atlas_name"]
Em.pnts$FlywayI = as.integer(factor(extract(Flyways, Em.pnts)[,"NAME"]))
Em.pnts$Flyway = ifelse(Em.pnts$FlywayI == 1, "Atlantic",
ifelse(Em.pnts$FlywayI == 2, "Central",
ifelse(Em.pnts$FlywayI == 3, "Mississippi",
ifelse(Em.pnts$FlywayI == 4, "Pacific", NA))))
Plt.set2m2 = as.data.frame(
Plt.set2m %>%
group_by(Flyway, Week) %>%
summarise(Total = sum(value)))
Yavg = as.data.frame(
Plt.set2m2 %>%
group_by(Flyway) %>%
summarise(Mean = mean(Total)))
PltPnts = as.data.frame(
Em.pnts@data %>%
filter(Test == "Positive") %>%
group_by(Flyway, Week) %>%
summarise(Pos = length(Week)))
PltPnts$Yaxis = ifelse(PltPnts$Flyway == "Atlantic", Yavg[1,2],
ifelse(PltPnts$Flyway == "Central", Yavg[2,2],
ifelse(PltPnts$Flyway == "Mississippi", Yavg[3,2],
ifelse(PltPnts$Flyway == "Pacific", Yavg[4,2], NA))))
PltPnts$Scale = as.numeric(ifelse(PltPnts$Pos == 1, 1,
ifelse(PltPnts$Pos > 0 & PltPnts$Pos <= 25, 2,
ifelse(PltPnts$Pos > 25 & PltPnts$Pos <= 50, 3,
ifelse(PltPnts$Pos > 50 & PltPnts$Pos <= 75, 4,
ifelse(PltPnts$Pos > 75 & PltPnts$Pos <= 100, 5,
ifelse(PltPnts$Pos > 100 & PltPnts$Pos <= 1000, 6,
6)))))))
PltPnts$Col = "red"
PltPnts$MyShape = 19
NPltPnts = as.data.frame(
Em.pnts@data %>%
filter(Test == "Negative") %>%
group_by(Flyway, Week) %>%
summarise(Neg = length(Week)))
NPltPnts = merge(PltPnts, NPltPnts, by=c("Flyway","Week"), all=TRUE)
NPltPnts[is.na(NPltPnts)] = 0
NPltPnts$Keep = ifelse(NPltPnts$Pos > 0, 0, 1)
NPltPnts = NPltPnts %>% filter(Keep == 1)
sum(NPltPnts$Pos)## [1] 0sum(PltPnts$Pos)## [1] 1497sum(NPltPnts$Neg)## [1] 197NPltPnts$Yaxis = ifelse(NPltPnts$Flyway == "Atlantic", Yavg[1,2],
ifelse(NPltPnts$Flyway == "Central", Yavg[2,2],
ifelse(NPltPnts$Flyway == "Mississippi", Yavg[3,2],
ifelse(NPltPnts$Flyway == "Pacific", Yavg[4,2], NA))))
dim(NPltPnts)## [1] 27 9dim(PltPnts)## [1] 35 7NPltPnts = NPltPnts %>%
mutate(Scale = 1,
Col = "black",
MyShape = 4) %>%
select(Flyway, Week, Yaxis, Scale, Col, MyShape)
PltPnts = PltPnts %>%
select(Flyway, Week, Yaxis, Scale, Col, MyShape)
All.pnts = rbind(NPltPnts, PltPnts)
All.pnts$MyShape = as.character(All.pnts$MyShape)
ShpScale = scale_shape_manual(name = "Serology", values = c(1,4),
labels = c("Positive", "Negative"))
colScale = scale_colour_manual(name = "Strain", values = c("red", "black"),
labels = c("Positive", "Negative"))
ggplot(Plt.set2m2, aes(Week, Total)) +
geom_line(lwd=0.5) +
geom_point(data = All.pnts,
aes(Week, Yaxis,
shape=factor(All.pnts$MyShape),
#col = factor(All.pnts$Col),
size = All.pnts$Scale),
position = position_jitter(w = 0.5, h = 0)) +
ShpScale +
#scale_colour_manual(values=c("red","blue")) +
facet_wrap(~Flyway, ncol=2, scales = "free_y") +
scale_size_continuous(name = "Count",
breaks = c(1:6),
labels = c("1", "2 - 25", "26 - 50", "51 - 75", "75 - 100", "101 +")) +
scale_y_continuous(breaks = scales::pretty_breaks(8), limits = c(0, NA)) +
theme_classic() +
ylab("Commercial Poultry Head (millions)") +
xlab("Week of Year") +
scale_x_continuous(breaks = seq(1, 53, 13),
limits =c(1,53)) +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "right",
legend.title=element_text(size=18),
legend.text=element_text(size=16),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
#panel.border = element_rect(colour = "black"),
title = element_text(face="bold", size=18, hjust=0.5),
axis.text=element_text(size=16),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(face="bold", size = 22),
axis.text.x = element_text(face="bold", size=14, vjust=0.5,
hjust=0.5, angle=0),
axis.title.x = element_text(face="bold", size = 22)) 