localregression.Rmd

---
title: "Local Regression analysis"
author: "Rosie"
date: "7/22/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(lubridate)
library(tidyverse)
library(vegan)
library(sf)
if (!require("rspatial")) devtools::install_github('rspatial/rspatial')
library(rspatial)
library( spgwr )
library(nlme)
library(boot)
library(itsadug)
library(mgcv)
library(raster)

#load the data
temps = readRDS("Temp_filtered (1).rds")


#first calculate the daily means
tempmean = temps %>%
  filter(Date > as.Date("2014-1-1"), Station != "RYF") %>%
  mutate(julian = yday(Date)) %>%
  group_by(Station, julian) %>%
  summarize(Tempave = mean(Temp, na.rm = T), TempMax = max(Temp, na.rm = T), 
            TempMin = min(Temp, na.rm = T))


tempmeanx = temps %>%
  filter(Date > as.Date("2014-1-1"), Station != "RYF") %>%
  mutate(julian = yday(Date), Year = year(Date)) %>%
  group_by(Station, julian, Year, Date) %>%
  summarize(Tempave = mean(Temp, na.rm = T), TempMax = max(Temp, na.rm = T), 
            TempMin = min(Temp, na.rm = T), Temprange = TempMax-TempMin)

save(tempmeanx, file = "tempmeanx.RData")

```

## Local regression

One of the cool ways you can address differences across space is run a regression that calculates different coefficients over space. I ran the model over the continuous stations, then extrapolated the results to the rest of the Detla.

First we have to do some data manipulation to get everything into a spatial format. 

```{r}
#read in shapefile of the delta
delta = read_sf("DeltaShapefile/hydro_delta_marsh.shp")

#add lat/longs for the stations
stas = read.csv("StationLatLongs.csv")


#attached lat/longs to mean temperature
tempmean2 = merge(tempmean, stas)

#Specify a coordinate reference system and turn it into a spatial object
alb <- CRS("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")

#Filter the Delta shapefile so it doesn't include areas where we dont have a lot of data
delta2 = sf::as_Spatial(dplyr::filter(delta, HNAME != "SAN FRANCISCO BAY", HNAME != "SAN PABLO BAY"))

#I'm not positive we have to do this transformation, but it was in teh example so I'm doing it.
ctst <- spTransform(delta2, alb)
```

Now I'll filter out just the data from January-June and find the optimal badwidth for the interpolation. I think this basically tells you how much "smoothing" to put between your sampling points.
```{r}
tempmean2.1 = dplyr::filter(tempmean2, julian <190)
tempmeansub = sample_n(tempmean2.1, 1000)
sp = tempmeansub
coordinates(sp) = ~ Longitude + Latitude
crs(sp) <- "+proj=longlat +datum=NAD83"
spt <- spTransform(sp, alb)

#find optimal bandwith
#THIS TAKES FOREVER!
bw <- gwr.sel( Temp~ julian, data= spt)


```

To make the pretty map, I need to create a blank raster from the Delta shapefile. Once I calculate the predicted values for the regression, I'll fill it in.


```{r}
r <- raster(ctst,  res=100)
rdelta2 <- rasterize(ctst, r)
newpts <- rasterToPoints(rdelta2)

#Turn Newpoints into a spatial grid
newpts2 = as.data.frame(newpts)
coordinates(newpts2) = ~ x  + y
crs(newpts2) = alb

save(newpts2, rdelta2, r, delta2, file = "spatialdata.RData")

```

Now I run the GWR function to calculate local regressions

```{r warning=F}
#

g <- gwr(Temp ~ julian, data= spt,
         bandwidth=bw, fit.points=newpts[, 1:2])
g

resultsg = as.data.frame(g$SDF)

#sigTest = abs(g$SDF$julian) -2 * g$SDF$julian_se 

#Note: i'd like to run it with "se.fit = T" so I can calculate the 
#standard error, but it doens't seem to be working. I had it running 
#all weekend, and noting. 

```

Put the fitted points back on my raster and plot it!

```{r fig.width=15}
slope <- rdelta2
intercept <- rdelta2
slope[!is.na(slope)] <- g$SDF$julian
intercept[!is.na(intercept)] <- g$SDF$'(Intercept)'
s <- stack(intercept, slope)
names(s) <- c('intercept', 'slope')
plot(s, xlim = c(-200000, -100000), ylim = c(-50000,50000))


```

##Now the cool part - how things have changed over the historical record

```{r}


#first calculate the daily means
temps$Year = year(temps$Date)
tempmax = temps %>%
  mutate(julian = yday(Date)) %>%
  group_by(Station, julian, Year, Date) %>%
  summarize(Tempave = mean(Temp, na.rm = T),
            maxTemp = max(Temp, na.rm = T),
            minTemp = min(Temp, na.rm = T))

tempmax2 = merge(tempmax, stas)

#tempmax2.1 = dplyr::filter(tempmax2, julian <250)
#spx = sample_n(tempmax2.1, 10000)
spx = tempmax2

#turn it into a spatial object
coordinates(spx) = ~ Longitude + Latitude
crs(spx) <- "+proj=longlat +datum=NAD83"
sptx <- spTransform(spx, alb)

#find badwidth
bwx <- gwr.sel( maxTemp~ julian + Year, data= sptx)


gx <- gwr(maxTemp~ julian + Year, data= sptx,
         bandwidth=bwx, fit.points=newpts[, 1:2])
gx


```

Put the fitted points back on my raster and plot it!

```{r fig.width=15}
slopeDay <- rdelta2
slopeYear <- rdelta2
intercept <- rdelta2
slopeDay[!is.na(slopeDay)] <- gx$SDF$julian
slopeYear[!is.na(slopeYear)] <- gx$SDF$Year
intercept[!is.na(intercept)] <- gx$SDF$'(Intercept)'
s <- stack(intercept, slopeDay, slopeYear)
names(s) <- c('intercept', 'Seasonal Change', "Inter-annual Change")
plot(s, xlim = c(-200000, -100000), ylim = c(-50000,50000))


```

## Temporal autocorrelation

I'm not sure I need to worry about temporal autocorrelation if time is part of my model, but I'll do some testing.

```{r}

testdat = filter(tempmax2.1, Station == first(Station)) %>%
arrange(Date)

acf(testdat$maxTemp)

pacf(testdat$maxTemp)

lmtest = lm(maxTemp~ julian + Year, data = testdat)

summary(lmtest)

acf(lmtest$residuals)

pacf(lmtest$residuals)

#find the start value of rho
acf(resid(lmtest), plot=FALSE)$acf[2]

# I'm confused and don't know what i'm diong

lmtest2 = gls(maxTemp~ julian + Year, data = testdat, correlation = corAR1(form = ~ Date))

acf(lmtest2$residuals)
acf_resid(lmtest2)


#try a gam
m1 = bam(maxTemp ~ s(julian) + s(Year), data=testdat)
acf(m1$residuals)

summary(m1)

#calculate start value of Rho
r1 = acf(resid(m1), plot=FALSE)$acf[2]

#identify the start of the time series
testdat = start_event(testdat, column="Date", event=c("Station", "Year"), label.event="Event")

m1AR1 = bam(maxTemp ~ s(julian) + s(Year), data=testdat, rho=r1, AR.start=testdat$start.event)

summary(m1AR1)
acf_resid(m1AR1)


```
So, I may want to try and use a GAM or something. Not sure.


## Package GWmodel

Before I investigate spatial GAMs, I want to try a different package for the geographically weighted regression.

First, compute the distances between the points on the grid where the model will be estimated, and the station locations. This takes a while, but if it is pre-computed, successive applications of GWR functions are speeded up.
```{r}
library(GWmodel)

DM <- gw.dist(dp.locat=coordinates(spt),rp.locat=coordinates(newpts2))


#now the model
gwr.res <- gwr.basic(Temp~ julian, data= spt, regression.points=newpts2, bw=bw, dMat=DM,kernel='gaussian')

gwr.res

```

Put the fitted points back on my raster and plot it!

```{r fig.width=15}
slopeDay <- rdelta2
intercept <- rdelta2
slopeDay[!is.na(slopeDay)] <- gwr.res$SDF$julian
#slopeYear[!is.na(slopeYear)] <- gx$SDF$Year
intercept[!is.na(intercept)] <- gwr.res$SDF$Intercept
s <- stack(intercept, slopeDay)
names(s) <- c('intercept', 'Seasonal Change')
plot(s, xlim = c(-200000, -100000), ylim = c(-50000,50000))

```

Estimate standard errors using bootstrapping

```{r}

bsm.res1 <- gwr.bootstrap(Temp~julian, kernel = "gaussian", sptx, R=99)
bsm.res1


set.seed(4676)

gwrcoef <- function(hpdf,i) gwr.basic, data=spt[i,], regression.points=newpts2, bw=bw, dMat=DM[i,],kernel='gaussian')$SDF$PROF

bootres <- boot(spt,gwrcoef,100)
gwr.res$SDF$bsePROF <- sqrt(apply(bootres$t,2,var))
image(gwr.res$SDF,'bsePROF')
plot(londonborough,add=TRUE)

###UGH!! NONE OF THIS IS WORKING!!!

```
Maybe I should just use seasonal means so I don' tneed to worry about the temporal autocorrelation

```{r}

#filter out just the summer months and calculate mean, min max for the season
summer = dplyr::filter(tempmax2, julian <244, julian >153) %>%
  group_by(Year, Station, Latitude, Longitude) %>%
  summarize(maxTemp = mean(maxTemp, na.rm = T), minTemp = mean(minTemp, na.rm = T), Temp = mean(Temp, na.rm = T))
#spx = sample_n(tempmax2.1, 10000)
sppt = summer

#turn it into a spatial object
coordinates(sppt) = ~ Longitude + Latitude
crs(sppt) <- "+proj=longlat +datum=NAD83"
sppt <- spTransform(sppt, alb)

#find badwidth
bwx <- gwr.sel( maxTemp~ Year, data= sppt)


DM <- gw.dist(dp.locat=coordinates(sppt),rp.locat=coordinates(newpts2))


#now the model
gwr.res <- gwr.basic(Temp~ Year, data= sppt, regression.points=newpts2, bw=bwx, dMat=DM,kernel='gaussian')

gwr.res

#see what the other package gives us
gx <- gwr(maxTemp~ Year, data= sppt,
         bandwidth=bwx, fit.points=newpts[, 1:2])
gx

slopeYear <- rdelta2
intercept <- rdelta2
slopeYear[!is.na(slopeYear)] <- gwr.res$SDF$Year
intercept[!is.na(intercept)] <- gwr.res$SDF$Intercept
s <- stack(intercept, slopeYear)
names(s) <- c('intercept', 'long-term change')
plot(s, xlim = c(-200000, -100000), ylim = c(-50000,50000))



bsm.res1 <- gwr.bootstrap(Temp~Year, kernel = "gaussian", sppt, R=99)
bsm.res1
#UGH still not working


```
## GAMs

```{r}
tempmax2 = arrange(tempmax2, Station, Date) %>%
  filter(Year >1995)

g1 = bam(Temp ~ s(julian, bs="cc" ) + s(Year) + s(Latitude, Longitude), 
data = tempmax2, method = "REML")

summary(g1)
plot(g1)
vis.gam(g1, view = c("Longitude", "Latitude"), plot.type = "contour",
        too.far = 0.1)

vis.gam(g1, view = c("Longitude", "Latitude"), plot.type = "persp",
        too.far = 0.1, se = 2)

#use a tensor product smooth for an interaction between 
#space and time

g2 = bam(Temp ~ s(julian, bs="cc" ) + s(Year) + 
           te(Latitude, Longitude, Year), 
data = tempmax2, method = "REML")

summary(g2)

#now the interaction version
g3 = bam(Tempave ~ s(julian, bs="cc" ) + s(Year) + 
           ti(Latitude, Longitude, Year) +
           ti(Latitude, Longitude, julian), 
data = tempmax2, method = "REML")

summary(g3)
plot(g3)
```

try adding the temporal correlation term

```{r}
#calculate start value of Rho
r1 = acf(resid(g3), plot=FALSE)$acf[2]

#identify the start of the time series
testdat = start_event(tempmax2, column="Date", event=c("Station", "Year"), label.event="Event")

g4 = bam(Tempave ~ s(julian, bs = "cc") + s(Year) + 
           ti(Latitude, Longitude, Year)+ 
           ti(Latitude, Longitude, julian), 
data =testdat, method = "REML",  rho=r1, AR.start=testdat$start.event)

summary(g4)
plot(g4)
acf(resid(g4))
acf_resid(g4, split_pred = c("Year"))
```
now do it with maximum temperatures

```{r}
g3x = bam(maxTemp ~ s(julian) + s(Year) + 
           ti(Latitude, Longitude, Year) +
           ti(Latitude, Longitude, julian), 

data = testdat, method = "REML")
#calculate start value of Rho
r1x = acf(resid(g3x), plot=FALSE)$acf[2]


g4x = bam(maxTemp ~ s(julian) + s(Year) + 
           ti(Latitude, Longitude, Year)+ 
           ti(Latitude, Longitude, julian), 
data =testdat, method = "REML",  rho=r1x, AR.start=testdat$start.event)

summary(g4x)
plot(g4x)
acf_resid(g4x)

vis.gam(g4x, view = c("Longitude", "Latitude"), plot.type = "contour",
        too.far = 0.1)

#gamtabs(g4x, type = 'html')


g4x2 = bam(maxTemp ~ 
           te(Latitude, Longitude, Year)+ 
           te(Latitude, Longitude, julian), 
data =testdat, method = "REML",  rho=r1x, AR.start=testdat$start.event)

summary(g4x2)
plot(g4x2)
acf_resid(g4x2)

vis.gam(g4x, view = c("Longitude", "Latitude"), plot.type = "contour",
        too.far = 0.1)


#gamtabs(g4x, type = 'html')
```

Just the stations with long-term data

```{r}
load("Tempmax.RData")
tempmax_sum = group_by(tempmax, Station) %>%
  summarize(minyear = min(Year), maxyear = max(Year))

tempmax_long = left_join(testdat, tempmax_sum) %>%
  filter(minyear < 2001)

g3x2 = bam(maxTemp ~ s(julian, bs = "cc") + s(Year), 

data = tempmax_long, method = "REML")
#calculate start value of Rho
r1x2 = acf(resid(g3x2), plot=FALSE)$acf[2]


g4x2 = bam(maxTemp ~ 
           s(Year)+ 
          s(julian, bs = "cc"), 
data =tempmax_long, method = "REML",  rho=r1x2, AR.start=tempmax_long$start.event)

summary(g4x2)
plot(g4x2)


#now try the Minimum temperatures, like peggy suggested

g3x2m = bam(minTemp ~ s(julian, bs = "cc") + s(Year), 

data = tempmax_long, method = "REML")
#calculate start value of Rho
r1x2m = acf(resid(g3x2m), plot=FALSE)$acf[2]


g4x2m = bam(minTemp ~ 
           s(Year)+ 
          s(julian, bs = "cc"), 
data =tempmax_long, method = "REML",  rho=r1x2m, AR.start=tempmax_long$start.event)

summary(g4x2m)
plot(g4x2m)
#hm. Minimum temperatures were not very interesting either
```

Just winter
```{r}
#let's try just the winter, since peggy and sam both thought that 
#might be useful

winter = mutate(tempmax_long, Month = month(Date)) %>%
  filter(Month %in% c(11))

gwin = bam(minTemp ~ 
           s(Year)+ 
          s(julian), 
data =winter, method = "REML",  rho=r1x2m, AR.start=winter$start.event)

summary(gwin)
plot(gwin)


gwin2 = bam(Tempave ~ 
           s(Year)+ 
          s(julian), 
data =winter, method = "REML",  rho=r1x2m, AR.start=winter$start.event)

summary(gwin2)
plot(gwin2)

```
If I don't have a lot of stations with the long time series, ditch spatial structure and move to random effects. 

```{r}
save(tempmax, file = "Tempmax.RData")
save(tempmax2.1, file = "tempmax2.1.RData")
save(tempmean, file = "tempmean.RData")
save(tempmean2.1, file = "tempmean.RData")
save(testdat, file = "testdat.RData")

load("testdat.RData")
g4x2 = gamm(maxTemp ~ 
           s(Year)+ 
           s(julian, bs = "cc"), 
data =tempmax_long, method = "REML", correlation = corAR1(form = ~1|Station))
#Error: cannot allocate vector of size 10.0 Gb
#meh
```

We wouldn't expect to see a ton of change-over-time anyway. Maybe concentrate on spatial differences
```{r}
g5 =  bam(maxTemp ~ s(julian, bs = "cc") +  
           te(Latitude, Longitude), 
data =testdat, method = "REML")

r5 = acf(resid(g5), plot=FALSE)$acf[2]



g5 =  bam(maxTemp ~ s(julian, bs = "cc") +
           te(Latitude, Longitude), 
data =testdat, method = "REML",  rho=r5, AR.start=testdat$start.event)

plot(g5)
```

I gotta figure out how to put these results on a map

```{r}
load("spatialdata.RData")
regions = read_sf("RosieRegions/shpExport.shp")
regions = st_transform(regions, crs = 4326)

stas = read.csv("StationLatLongs.csv")
coordinates(stas) = ~ Longitude + Latitude
crs(stas) <- "+proj=longlat +datum=NAD83"


delta2x = spTransform(delta2, CRS("+init=epsg:4326"))
r <- raster(delta2x,  res=100)
rdelta3 <- rasterize(delta2, r)
newpts3 <- rasterToPoints(rdelta3)

predict(newpts, g5)


delta = st_transform(delta,crs=4326)
stas = st_as_sf(stas) %>%
  st_transform(stas, crs=4326)

WQ_pred<-function(Full_data=Data,
                  Delta_subregions = regions,
                  Delta_water=delta,
                  Stations = stas,
                  n=100, 
                  Years=round(seq(min(Full_data$Year)+2, max(Full_data$Year)-2, length.out=9)),
                  Julian_days=yday(ymd(paste("2001", c(1,4,7,10), "15", sep="-"))) #Jan, Apr, Jul, and Oct 15 for a non-leap year

){
  
  # Create point locations on a grid for predictions
  Points<-st_make_grid(Delta_subregions, n=n)%>%
    st_as_sf(crs=st_crs(Delta_subregions))%>%
    st_join(Delta_water)%>% 
    
    # Joining a map of delta waterways (from my spacetools package) to ensure all these points are over water.
        st_coordinates()%>%
    as_tibble()%>%
    mutate(Location=1:nrow(.))%>%
    dplyr::select(Longitude=X, Latitude=Y, Location)
  
  
  # Create full dataset for predictions
  newdata<-expand.grid(Year= Years,
                       Location=1:nrow(Points),
                       Julian_day=Julian_days)%>% # Create all combinations of predictor variables
    left_join(Points, by="Location")%>% #Add Lat/Longs to each location
    mutate(Latitude_s=(Latitude-mean(Full_data$Latitude, na.rm=T))/sd(Full_data$Latitude, na.rm=T), # Standardize each variable based on full dataset for model
           Longitude_s=(Longitude-mean(Full_data$Longitude, na.rm=T))/sd(Full_data$Longitude, na.rm=T),
           Year_s=(Year-mean(Full_data$Year, na.rm=T))/sd(Full_data$Year, na.rm=T),
           Julian_day_s = (Julian_day-mean(Full_data$julian, na.rm=T))/sd(Full_data$julian, na.rm=T),
           Year_fac=factor(Year)) %>%
    st_as_sf(coords=c("Longitude", "Latitude"), crs=4326, remove=FALSE)%>% # Turn into sf object
      
    st_transform(crs=st_crs(Delta_subregions))%>% # transform to crs of Delta shapefile
      
    st_join(Delta_subregions, join = st_intersects)

  return(newdata)
}


newdata_year <- WQ_pred(Full_data=tempmax_long, 
                        Julian_days = yday(ymd(paste("2001", 1:12, "15", sep="-"))),
                        Years=round(min(tempmax_long$Year):max(tempmax_long$Year)))

newdata_year = rename(newdata_year, julian = Julian_day)

modellc4_predictions<-predict(g5, newdata=newdata_year, type="response", se.fit=TRUE, discrete=T) # Create predictions

# Predictions stored as "modellc4_predictions.Rds"



newdata<-newdata_year%>%
  mutate(Prediction=modellc4_predictions$fit)%>%
  mutate(SE=modellc4_predictions$se.fit,
         L95=Prediction-SE*1.96,
         U95=Prediction+SE*1.96)%>%
  mutate(Date=as.Date(julian, origin=as.Date(paste(Year, "01", "01", sep="-")))) # Create Date variable from Julian Day and Year

library(stars)
# Function to rasterize all dates. Creates a 3D raster Latitude x Longitude x Date 
Rasterize_all <- function(data, var, out_crs=4326, n=100){
  var<-rlang::enquo(var)
  rlang::as_name(var)
  preds<-map(unique(data$Date), function(x) st_rasterize(data%>%
                                                           filter(Date==x)%>%
                                                           dplyr::select(!!var), 
                                                         template=st_as_stars(st_bbox(delta2), dx=diff(st_bbox(delta2)[c(1, 3)])/n, dy=diff(st_bbox(delta2)[c(2, 4)])/n, values = NA_real_))%>%
               st_warp(crs=out_crs))
  
  # Then bind all dates together into 1 raster
  out <- exec(c, !!!preds, along=list(Date=unique(data$Date)))
  return(out)
}

# Create full rasterization of all predictions for interactive visualizations
rastered_preds<-Rasterize_all(newdata, Prediction)
# Same for SE
rastered_SE<-Rasterize_all(newdata, SE)
# Bind SE and predictions together
rastered_predsSE<-c(rastered_preds, rastered_SE)

newdata_rast_season <- newdata%>%
  mutate(Month=month(Date))%>%
  filter(Year%in%round(seq(min(newdata$Year)+2, max(newdata$Year)-2, length.out=9)) & Month%in%c(1,4,7,10))



raster_plot<-function(data, Years=unique(newdata$Year), labels="All"){
  ggplot()+
    geom_stars(data=data)+
   # facet_wrap(julian~.)+
    scale_fill_viridis_c(name="Temperature", na.value="white", breaks=seq(6,26,by=1), labels= function(x) ifelse((x/2)==as.integer(x/2), as.character(x), ""),
                         guide = guide_colorbar(direction="horizontal", title.position = "top", barwidth = 4, ticks.linewidth = 2,
                                                barheight=0.4, title.hjust=0.5, label.position="bottom", label.theme=element_text(size=8), 
                                                title.theme=element_text(size=10)))+
    coord_sf()+
    ylab("Latitude")+
    xlab("Longitude")+
    theme_bw()
  }

raster_plot(rastered_preds)



#########################################
rp2 = as.data.frame(rastered_preds)
 ggplot()+
    geom_stars(data=rastered_preds)+
    scale_fill_viridis_c(name="Temperature", na.value="white", breaks=seq(6,26,by=1), labels= function(x) ifelse((x/2)==as.integer(x/2), as.character(x), ""),
                         guide = guide_colorbar(direction="horizontal", title.position = "top", barwidth = 4, ticks.linewidth = 2,
                                                barheight=0.4, title.hjust=0.5, label.position="bottom", label.theme=element_text(size=8), 
                                                title.theme=element_text(size=10)))+
    coord_sf()+
    ylab("Latitude")+
    xlab("Longitude")+
    theme_bw()

```