code_05_check_daily_autocorrelation.R

library(magrittr)
library(mapview)

## Loading required package: leaflet

library(fst)
library(sp)
library(gstat)
library(automap)
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

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

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

set.seed(1)
for(i in 1:length(days)) {
  
  # Subset
  tmp = dat[dat$day == days[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.09315068

# 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")