code_07_check_daily_autocorrelation_in_bimon_error.R

library(magrittr)
library(mapview)

## Loading required package: leaflet

library(fst)
library(sp)
library(gstat)
library(automap)
library(plyr)
library(spdep)

## Loading required package: Matrix

setwd("/home/michael/Dropbox/BGU/Itai_Kloog/p_46_FRK_Kriging_test/")

# Read data
dat = read.fst("01_raw/modX2009.fst")

# Rename
dat$lon = dat$Longitude
dat$lat = dat$Latitude
dat$Longitude = NULL
dat$Latitude = NULL

# Classify 'bimon'
dat$month = dat$day %>% format("%m") %>% as.numeric
dat$bimon[dat$month %in% c(1:2)] = "01-02"
dat$bimon[dat$month %in% c(3:4)] = "03-04"
dat$bimon[dat$month %in% c(5:6)] = "05-06"
dat$bimon[dat$month %in% c(7:8)] = "07-08"
dat$bimon[dat$month %in% c(9:10)] = "09-10"
dat$bimon[dat$month %in% c(11:12)] = "11-12"

# Aggregate
dat = ddply(
  dat,
  c("station", "bimon"),
  summarize,
  lon = mean(lon),
  lat = mean(lat),
  error = mean(error)
)

# To Spatial
coordinates(dat) = ~ lon + lat
proj4string(dat) = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

# Subset bimon
bimon = dat$bimon %>% unique %>% sort
s = rep(NA, length(bimon))
p = rep(NA, length(bimon))

set.seed(1)
for(i in 1:length(bimon)) {
  
  # Subset
  tmp = dat[dat$bimon == bimon[i], ]
  
  # Moran's I test
  w = 1/as.matrix(dist(coordinates(tmp)))
  diag(w) = 0
  fit = moran.test(tmp$error, mat2listw(w))
  s[i] = fit$statistic
  p[i] = fit$p.value
  
}

# Proportion of days with significant positive autocorrelation
mean(s > 0 & p < 0.05)

## [1] 0

# Plot - Estimates
hist(s, breaks = 50, main = "Moran's I estimates, per day")
abline(v = 0, col = "red")
text(0.05, 40, "I = 0", col = "red")

# Plot - Estimates
hist(p, breaks = 50, main = "Moran's I p-values, per day")
abline(v = 0.05, col = "red")
text(0.05, 23, "p = 0.05", col = "red")