setwd("E:/spring16")
library(sp)
## Warning: package 'sp' was built under R version 3.2.5
load("wUSPrec_Bnd.Rdata")
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.2.5
library(ggmap)
## Warning: package 'ggmap' was built under R version 3.2.5
sp_map <- fortify(wUSBoundary)
## Regions defined for each Polygons
wUSPrec.df <- data.frame(wUSPrec)
gg <- ggplot()
gg <- gg + geom_map(data=sp_map,map=sp_map,
                    aes(x=long, y=lat, map_id=id),
                    fill= "white",color="black")
gg <- gg + geom_point(data=wUSPrec.df, aes(x=long, y=lat, size=prec), 
                      color="orange")
gg <- gg + scale_size(breaks=seq(0,300,by=50))
gg <- gg + labs(x=NULL, y=NULL)
gg <- gg + theme(panel.grid=element_blank(), panel.border=element_blank())
gg <- gg + theme(axis.ticks=element_blank(), axis.text=element_blank())
gg <- gg + theme(legend.position="top")
gg

Look at that map! Right away we can see that large dots are near large dots and small dots are near small dots, there must me some spatial autocorrelation going on.

Scatterplot matrix

library(gstat)
## Warning: package 'gstat' was built under R version 3.2.5
library(sp)
library(nlme)
library(spatstat)
## Warning: package 'spatstat' was built under R version 3.2.5
## 
## spatstat 1.45-0       (nickname: 'One Lakh') 
## For an introduction to spatstat, type 'beginner'
## 
## Attaching package: 'spatstat'
## The following object is masked from 'package:gstat':
## 
##     idw
library(ncf)
## Warning: package 'ncf' was built under R version 3.2.4
library(Matrix)
library(spdep)
## Warning: package 'spdep' was built under R version 3.2.5
library(raster)
## Warning: package 'raster' was built under R version 3.2.5
## 
## Attaching package: 'raster'
## The following objects are masked from 'package:spatstat':
## 
##     area, rotate, shift
## The following object is masked from 'package:nlme':
## 
##     getData
hscat(prec~1,data=wUSPrec,breaks = seq(0,1e5,by=2e4))

This is telling us what we saw on the map, just with some r-squared values. We can see that near neighbors are more similar than distant neighbors.Seems like after 80km there is not much correlation. Correlations are great and all but these do not tell us anything about direction, that coming up.

Variograms

preVar <- variogram(prec~1, data=wUSPrec, cloud = FALSE)
  plot(preVar,pch=20,cex=1.5,col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Precipitation")

anivar=variogram(prec~1, data=wUSPrec,cloud=FALSE, alpha=c(0,90,180,270))
plot(anivar,pch=20,cex=1.25,main="Precipitation by direction")

I created two variograms. The first is the standard one which shows us squared variance by distance. The second is the also squared variance by distance but incorporating direction as well (at least I am petty sure that is what this code did). It looks like there is strong autocorrelation to the north and south that gradually decreases with distance. However, autocorrelation to the east and west is strong to a point, maybe around 100km or less, then are no longer correlated.This is Nifty! because it tells us that direction is important to the structure of these data.

Moran’s I (That sounds like it should be a chapter in some science fiction novel)

p.1 <- 1/as.matrix(dist(coordinates(wUSPrec)))
diag(p.1) <- 0
moran.test(wUSPrec$prec,mat2listw(p.1))
## 
##  Moran I test under randomisation
## 
## data:  wUSPrec$prec  
## weights: mat2listw(p.1)  
## 
## Moran I statistic standard deviate = 11.886, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1657032288     -0.0089285714      0.0002158469
n <- 6
res <- data.frame(k=2^(1:n),I=rep(NA,n))
for(i in 1:n){
  w <- knn2nb(knearneigh(wUSPrec,k=2^i))
  res$I[i] <- moran.test(wUSPrec$prec,nb2listw(w))$estimate[1]
}
plot(res,type="b",main="Moran's I by number of neighbors",pch=20,cex=1.5)

Here we have I as a function of number of neighbors. We see that as I decreases with increasing neighbor number. This is cool and all, but not exactly helpful for environmental data sets.

Now by distance

precip1 <- spline.correlog(x=coordinates(wUSPrec)[,1], y=coordinates(wUSPrec)[,2],
                         z=wUSPrec$prec, resamp=100, quiet=TRUE)
plot(precip1)

Looks like there is positive correlation up until about 175km then there may or may not be some correlation between there and about 700km. Then there seems to be some negative correlation after 700km. This is also pretty nifty. It would be interesting to see this plot by direction since we know N,S,E,w are correlated to different degrees at different distances.