--- #####-----Project: Vivi-Ovi-----##### title: "Spatial Analysis" author: "Liang Ma" date: "July 20, 2017" output: word_document --- #####-----Spatial Regression Analysis for Temp, Lat and Ele-----##### ```{r} setwd("D:/Liang/Projects/Vivi_vs_ovi/Coding/Spatial Regression Analysis/") memory.limit(size = 100000) data<-read.csv("data.csv") inds=sample(1:nrow(data), 5000, replace = FALSE, prob = NULL) data= data[inds,] # Establish basic linear model mod<-lm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation) summary(mod) coords= as.matrix(data[,c("Longitude","Latitude")]) # Develop neighborhood lists for different maximum distances library(spdep) m.nb20 <- dnearneigh(coords,0,20,longlat=TRUE) m.nb400 <- dnearneigh(coords,0,400,longlat=TRUE) m.nb800 <- dnearneigh(coords,0,800,longlat=TRUE) m.nb1500 <- dnearneigh(coords,0,1500,longlat=TRUE) ``` ```{r} summary(m.nb20) summary(m.nb400) summary(m.nb800) summary(m.nb1500) ``` # Plot neighborhood matrices ## NOTE: this is set up for a global analysis, so will need distance adjustments ```{r} plot(m.nb20, coords,points=TRUE) plot(m.nb400, coords,points=TRUE) plot(m.nb800, coords,points=TRUE) plot(m.nb1500, coords,points=TRUE) ``` # establish weights for the given list of neighbors, "B" = binary (all equal), "W" = row standardized ```{r} #m.sw20 <- nb2listw(m.nb20, glist=NULL, style="W", zero.policy=TRUE) m.sw400 <- nb2listw(m.nb400, glist=NULL, style="W", zero.policy=TRUE) m.sw800 <- nb2listw(m.nb800, glist=NULL, style="W", zero.policy=TRUE) m.sw1500 <- nb2listw(m.nb1500, glist=NULL, style="W", zero.policy=TRUE) ``` # Test for spatial autocorrelation in a single (e.g. response) variable, given the neighborhod links and weights ```{r} #moran.test(data$Prop_vivi, m.sw20, zero.policy = TRUE) #moran.plot(data$Prop_vivi, m.sw20, zero.policy = TRUE) moran.test(data$Prop_vivi, m.sw400, zero.policy = TRUE) moran.plot(data$Prop_vivi, m.sw400, zero.policy = TRUE) moran.test(data$Prop_vivi, m.sw800, zero.policy = TRUE) moran.plot(data$Prop_vivi, m.sw800, zero.policy = TRUE) moran.test(data$Prop_vivi, m.sw1500, zero.policy = TRUE) moran.plot(data$Prop_vivi, m.sw1500, zero.policy = TRUE) ``` # Test for spatial autocorrelation of the linear model residuals, given the neighborhod links and weights ```{r} #lm.morantest(mod, m.sw20, zero.policy=TRUE) lm.morantest(mod, m.sw400, zero.policy=TRUE) #lm.morantest(mod, m.sw800, zero.policy=TRUE) #lm.morantest(mod, m.sw1500, zero.policy=TRUE) ``` # Calculate correlograms of nonspatial model residuals giving Moran's I for each distance lag # Spline correlogram gives smooth fit with 95% ci ```{r} library(ncf) cor.lmresid <- correlog(data$Longitude,data$Latitude,z=residuals(mod),increment=100,resamp=1,latlon=TRUE) corspline.lmresid <- spline.correlog(data$Longitude,data$Latitude,z=residuals(mod),resamp = 50, latlon=TRUE) ``` ```{r, echo=FALSE} # Plot correlograms par(mfrow=c(2,2)) plot(cor.lmresid$mean.of.class[1:150], cor.lmresid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="linear model",ylim=c(0,0.6),xlim=c(0,2000)) plot(corspline.lmresid) ``` # Lagrange Multiplier diagnostics for spatial dependence in linear models #The statistics are the simple LM test for error dependence (LMerr), the simple LM test for a missing spatially lagged dependent variable (LMlag), variants of these robust #to the presence of the other (RLMerr, RLMlag - RLMerr tests for error dependence in the possible presence of a missing lagged dependent variable, RLMlag the other way round), and a portmanteau #test (SARMA, in fact LMerr + RLMlag). #In a specification search, you would first assess which of the LMerr and LMlag are significant. In this case, both are, which necessitates the comparison of the robust forms. ```{r} #lm.LMtests(mod, m.sw20, zero.policy=TRUE, test=c("LMerr", "LMlag", "RLMerr","RLMlag", "SARMA")) lm.LMtests(mod, m.sw400, zero.policy=TRUE, test=c("LMerr", "LMlag", "RLMerr","RLMlag", "SARMA")) lm.LMtests(mod, m.sw800, zero.policy=TRUE, test=c("LMerr", "LMlag", "RLMerr","RLMlag", "SARMA")) lm.LMtests(mod, m.sw1500, zero.policy=TRUE, test=c("LMerr", "LMlag", "RLMerr","RLMlag", "SARMA")) ``` # Perform spatial simultaneous autoregressive lag model estimation ```{r} #m_x.lagSAR20 <- lagsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw20, method="eigen", quiet=FALSE, zero.policy=TRUE, interval=c(-0.5,0.5)) m_x.lagSAR400 <- lagsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw400, method="eigen", quiet=FALSE, zero.policy=TRUE, tol.solve = 1e-15,interval=c(-0.5,0.5)) #m_x.lagSAR800 <- lagsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw800, method="eigen", quiet=FALSE, zero.policy=TRUE, interval=c(-0.5,0.5)) #m_x.lagSAR1500 <- lagsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw1500, method="eigen", quiet=FALSE, zero.policy=TRUE, interval=c(-0.5,0.5)) #summary(m_x.lagSAR20) summary(m_x.lagSAR400) #summary(m_x.lagSAR800) #summary(m_x.lagSAR1500) ``` # Perform spatial simultaneous autoregressive error model estimation ```{r} #m_x.errorSAR20 <- errorsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw20, method="eigen", quiet=FALSE, zero.policy=TRUE, interval=c(-0.5,0.5)) m_x.errorSAR400 <- errorsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw400, method="eigen", quiet=FALSE, zero.policy=TRUE, tol.solve = 1e-15,interval=c(-0.5,0.5)) #m_x.errorSAR800 <- errorsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw800, method="eigen", quiet=FALSE, zero.policy=TRUE, interval=c(-0.5,0.5)) #m_x.errorSAR1500 <- errorsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw1500, method="eigen", quiet=FALSE, zero.policy=TRUE, interval=c(-0.5,0.5)) #summary(m_x.errorSAR20) summary(m_x.errorSAR400) #summary(m_x.errorSAR800) #summary(m_x.errorSAR1500) ``` #Perform spatial simultaneous autoregressive lag + error (mixed) model (i.e. Durbin model) ```{r} #m_x.mixedSAR20 <- lagsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw20, type="mixed",method="eigen", quiet=FALSE, zero.policy=TRUE, interval=c(-0.5,0.5)) m_x.mixedSAR400 <- lagsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw400, type="mixed",method="eigen", quiet=FALSE, zero.policy=TRUE, tol.solve = 1e-15,interval=c(-0.5,0.5)) #m_x.mixedSAR800 <- lagsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw800, type="mixed",method="eigen", quiet=FALSE, zero.policy=TRUE, interval=c(-0.5,0.5)) #m_x.mixedSAR1500 <- lagsarlm(data$Prop_vivi~data$Dev_temp*data$Latitude*data$Elevation, data=data, m.sw1500, type="mixed",method="eigen", quiet=FALSE, zero.policy=TRUE, interval=c(-0.5,0.5)) #summary(m_x.mixedSAR20) summary(m_x.mixedSAR400) #summary(m_x.mixedSAR800) #summary(m_x.mixedSAR1500) ``` ```{r} #Calculate correlograms of spatial model residuals #cor.lagSAR20resid <- correlog(data$Longitude,data$Latitude,z=residuals(m_x.lagSAR20),increment=100,resamp=1,latlon=TRUE) cor.lagSAR400resid <- correlog(data$Longitude,data$Latitude,z=residuals(m_x.lagSAR400),increment=100,resamp=1,latlon=TRUE) #cor.lagSAR800resid <- correlog(data$Longitude,data$Latitude,z=residuals(m_x.lagSAR800),increment=100,resamp=1,latlon=TRUE) #cor.lagSAR1500resid <- correlog(data$Longitude,data$Latitude,z=residuals(m_x.lagSAR1500),increment=100,resamp=1,latlon=TRUE) #cor.errorSAR20resid <- correlog(data$Longitude,data$Latitude,z=residuals(m_x.errorSAR20),increment=100,resamp=1,latlon=TRUE) cor.errorSAR400resid <- correlog(data$Longitude,data$Latitude,z=residuals(m_x.errorSAR400),increment=100,resamp=1,latlon=TRUE) #cor.errorSAR800resid <- correlog(data$Longitude,data$Latitude,z=residuals(m_x.errorSAR800),increment=100,resamp=1,latlon=TRUE) #cor.errorSAR1500resid <-correlog(data$Longitude,data$Latitude,z=residuals(m_x.errorSAR1500),increment=100,resamp=1,latlon=TRUE) #cor.mixedSAR20resid <- correlog(data$Longitude,data$Latitude,z=residuals(m_x.mixedSAR20),increment=100,resamp=1,latlon=TRUE) cor.mixedSAR400resid <- correlog(data$Longitude,data$Latitude,z=residuals(m_x.mixedSAR400),increment=100,resamp=1,latlon=TRUE) #cor.mixedSAR800resid <- correlog(data$Longitude,data$Latitude,z=residuals(m_x.mixedSAR800),increment=100,resamp=1,latlon=TRUE) #cor.mixedSAR1500resid <-correlog(data$Longitude,data$Latitude,z=residuals(m_x.mixedSAR1500),increment=100,resamp=1,latlon=TRUE) ``` ```{r} # Plot correlograms par(mfrow=c(2,2)) #plot(cor.lagSAR20resid$mean.of.class[1:150], cor.lagSAR20resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",ylim=c(-0.4,0.6)) plot(cor.lagSAR400resid$mean.of.class[1:150], cor.lagSAR400resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",xlim=c(0,2000),ylim=c(-0.4,0.6)) #plot(cor.lagSAR800resid$mean.of.class[1:150], cor.lagSAR800resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",ylim=c(-0.4,0.6)) #plot(cor.lagSAR1500resid$mean.of.class[1:150], cor.lagSAR1500resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",ylim=c(-0.4,0.6)) #plot(cor.errorSAR20resid$mean.of.class[1:150], cor.errorSAR20resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",ylim=c(-0.4,0.6)) plot(cor.errorSAR400resid$mean.of.class[1:150], cor.errorSAR400resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",xlim=c(0,2000),ylim=c(-0.4,0.6)) #plot(cor.errorSAR800resid$mean.of.class[1:150], cor.errorSAR800resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",ylim=c(-0.4,0.6)) #plot(cor.errorSAR1500resid$mean.of.class[1:150], cor.errorSAR1500resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",ylim=c(-0.4,0.6)) #plot(cor.mixedSAR20resid$mean.of.class[1:150], cor.mixedSAR20resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",ylim=c(-0.4,0.6)) plot(cor.mixedSAR400resid$mean.of.class[1:150], cor.mixedSAR400resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",xlim=c(0,2000),ylim=c(-0.4,0.6)) #plot(cor.mixedSAR800resid$mean.of.class[1:150], cor.mixedSAR800resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",ylim=c(-0.4,0.6)) #plot(cor.mixedSAR1500resid$mean.of.class[1:150], cor.mixedSAR1500resid$correlation[1:150],xlab="Distance class",ylab="Moran's I",type="b",main="spatial model",ylim=c(-0.4,0.6)) #plot(corspline.lagSAR800resid) ``` ```{r} # Perform global Moran's test on residuals. Not really correct way of doing it! moran.test(residuals(m_x.lagSAR400), m.sw400, zero.policy = TRUE) # use normalization test to account for residuals, Anselin P11 moran.test(residuals(m_x.lagSAR400), m.sw400, randomisation=FALSE,alternative="two.sided", zero.policy = TRUE) moran.test(residuals(m_x.errorSAR400), m.sw400, zero.policy = TRUE) # use normalization test to account for residuals, Anselin P11 moran.test(residuals(m_x.errorSAR400), m.sw400, randomisation=FALSE,alternative="two.sided", zero.policy = TRUE) moran.test(residuals(m_x.mixedSAR400), m.sw400, zero.policy = TRUE) # use normalization test to account for residuals, Anselin P11 moran.test(residuals(m_x.mixedSAR400), m.sw400, randomisation=FALSE,alternative="two.sided", zero.policy = TRUE) # Model comparison AIC(mod) #AIC(m_x.lagSAR20) AIC(m_x.lagSAR400) #AIC(m_x.lagSAR800) #AIC(m_x.lagSAR1500) #AIC(m_x.errorSAR20) AIC(m_x.errorSAR400) #AIC(m_x.errorSAR800) #AIC(m_x.errorSAR1500) #AIC(m_x.mixedSAR20) AIC(m_x.mixedSAR400) #AIC(m_x.mixedSAR800) #AIC(m_x.mixedSAR1500) ```