23 Day 23 (April 16)

23.1 Announcements

23.2 Spatio-temporal models for disease data

  • Example

    library(sf)
library(sp)
library(raster)

url <- "https://www.dropbox.com/scl/fi/9ymxt900s77uq50ca6dgc/Enders-et-al.-2018-data.csv?rlkey=0rxjwleenhgu0gvzow5p0x9xf&dl=1"
df1 <- read.csv(url)
df1 <- df1[,c(1,4,5,8,9,10)] # Keep only the data on bird cherry-oat aphid

head(df1)
    ##   BCOA BYDV.totalpos.BCOA BCOA.totaltested year      long      lat
## 1   12                  2               10 2014 -95.16269 37.86238
## 2    1                  0                1 2014 -95.28463 38.29669
## 3    2                  0                2 2014 -95.33038 39.59482
## 4    0                 NA                0 2014 -95.32098 39.50696
## 5    8                  0                8 2014 -98.55469 38.48455
## 6    1                  0                1 2014 -98.84792 38.32772
    # Download shapefile of Kansas from census.gov
download.file("http://www2.census.gov/geo/tiger/GENZ2015/shp/cb_2015_us_state_20m.zip", destfile = "states.zip")
unzip("states.zip")
sf.us <- st_read("cb_2015_us_state_20m.shp",quiet = TRUE)
sf.kansas <- sf.us[48,6]
sf.kansas <- as(sf.kansas, 'Spatial')
#plot(sf.kansas,main="",col="white")

# Make SpatialPoints data frame
pts.sample <- data.frame(long = df1$long,lat = df1$lat, 
                     count = df1$BCOA,
                     BYDV.pos = df1$BYDV.totalpos.BCOA,
                     BYDV.tot = df1$BCOA.totaltested,
                     BYDV.prop = df1$BYDV.totalpos.BCOA/df1$BCOA.totaltested)
coordinates(pts.sample) =~ long + lat
proj4string(pts.sample) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")


# Plot counts of Bird cherry-oat aphid
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(sf.kansas,main="Abundance of Bird Cherry-oat Aphid")
points(pts.sample[,1],col=rgb(0.4,0.8,0.5,0.9),pch=ifelse(pts.sample$count>0,20,4),cex=pts.sample$count/50+0.5)
legend("right",inset=c(-0.25,0),legend = c(0,1,10,20,40,60), bty = "n", text.col = "black", 
   pch=c(4,20,20,20,20,20), cex=1.3,pt.cex=c(0,1,10,20,40,60)/50+0.5,col=rgb(0.4,0.8,0.5,0.9))

    # Plot proportion of number of Bird cherry-oat aphid infected with BYDV
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(sf.kansas, main="Proportion infected")
points(pts.sample[,4],col=rgb(0.8,0.5,0.4,0.9),pch=ifelse(is.na(pts.sample$BYDV.prop)==FALSE,20,4),cex=ifelse(is.na(pts.sample$BYDV.prop)==FALSE,pts.sample$BYDV.prop,0)/0.5+0.5)
legend("right",inset=c(-0.25,0),legend = c("NA",0,0.25,0.50,0.75,1.00), bty = "n", text.col = "black", pch=c(4,20,20,20,20,20), cex=1.3,pt.cex=c(0,0,0.25,0.50,0.75,1.00)/0.5+0.5,col=rgb(0.8,0.5,0.4,0.9))

  • Study goals

    • Make accurate predictions of vector abundance and probability of BYDV infection at times and locations where data was not collected.
    • Understand the environmental factors (e.g., temperature) that influence the vector abundance and probability of BYDV infection.

  • Auxiliary data

    • Weather/climate data
    • Land cover data
    • Crop data

  • Model choice: dynamic vs. descriptive spatio-temporal model?

    • Class discussion

  • Next steps

    • Write out the statistical model(s)
    • Determine how we will evaluate how model(s)
    • Determine how we will make inference

23.3 Earthquake data

  • New York Times article

  • New York Times articles

  • Earthquake data from Kansas

    library(sf)
library(sp)
library(raster)
library(lubridate)
library(animation)    
library(gifski)


# Download shapefile of Kansas from census.gov
download.file("http://www2.census.gov/geo/tiger/GENZ2015/shp/cb_2015_us_state_20m.zip", destfile = "states.zip")
unzip("states.zip")
sf.us <- st_read("cb_2015_us_state_20m.shp",quiet = TRUE)
sf.kansas <- sf.us[48,6]
sf.kansas <- as(sf.kansas, 'Spatial')

# Load earthquake data
url <- "https://www.dropbox.com/scl/fi/uhicc7qo4zcxeq79y1eh9/ks_earthquake_data.csv?rlkey=lbq8kxzx9g067domp46ffh7jw&dl=1"
df.eq <- read.csv(url)
df.eq$Date.time <- ymd_hms(df.eq$Date.time)
pts.eq <- SpatialPointsDataFrame(coords = df.eq[,c(3,2)], data = df.eq[,c(1,4)], proj4string = crs(sf.kansas))

# Plot spatial map earthquake data
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(sf.kansas,main="")
points(pts.eq,col=rgb(0.4,0.8,0.5,0.3),pch=21,cex=pts.eq$Magnitude/3)
legend("right",inset=c(-0.25,0),legend = c(1,2,3,4,5), bty = "n", text.col = "black", 
   pch=21, cex=1.3,pt.cex=c(1,2,3,4,5)/3,col=rgb(0.4,0.8,0.5,0.6))

    # Plot timeseries of earthquake data
plot(as.numeric(names(table(year(df.eq$Date.time)))),c(table(year(df.eq$Date.time))),
 xlab="Year",ylab="Number of Earthquakes in KS",pch=20)

    plot(df.eq$Date.time,df.eq$Magnitude,pch=20,cex=.2,xlab="Year",ylab="Earthquake magnitude")

    • Animation of Kansas Earthquake data
    # Make animation of earthquake data
date <- seq(as.Date("1977-1-1"), as.Date("2020-10-31"), by = "year")
t <- round(time_length(interval(min(date),pts.eq$Date.time),"year"))

for(i in 1:length(date)){
  par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
  plot(sf.kansas,main=year(date[i]))
  legend("right",inset=c(-0.25,0),legend = c(1,2,3,4,5), bty = "n", text.col = "black", 
     pch=20, cex=1.3,pt.cex=c(1,2,3,4,5)/3,col=rgb(0.4,0.8,0.5,1))
  if(length(which(t==i))>0){
  points(pts.eq[which(t==i),],col=rgb(0.4,0.8,0.5,1),pch=20,cex=pts.eq[which(t==i),]$Magnitude/3)
                        }
                    }

23.4 Spatio-temporal models for earthquake data

  • What are the goals of our study?
    • Prediction
      • Make point and areal level predictions
    • Inference
      • Understand the spatio-temporal covariates that may increase or decrease the risk of an earthquake
  • What auxiliary data do we need?
  • Kansas oil and gas well data
  • Descriptive vs. dynamic approach
  • Write out statistical model for earthquake data
  • Semi-live example (Download R code)