tutorial.R
In spacemakeR: Spatial modelling

### R code from vignette source 'tutorial.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: t1
###################################################
options(prompt=" ", continue=" ", width = 70)
options(digits=4)


###################################################
### code chunk number 2: ni1
###################################################
library(spdep)
#help(package="spdep")


###################################################
### code chunk number 3: ni2
###################################################
xygrid<-expand.grid(1:10,1:8)
plot(xygrid)


###################################################
### code chunk number 4: ni3
###################################################
#?cell2nb
nb1<-cell2nb(10,8,type="queen")

plot(nb1,xygrid,col="green",pch=20,cex=2)
title(main="queen neighborhood",cex.main=2)

nb1


###################################################
### code chunk number 5: ni4
###################################################
nb2<-cell2nb(10,8,type="rook")

plot(nb2,xygrid,col="blue",pch=20,cex=2)
title(main="rook neighborhood",cex.main=2)

nb2


###################################################
### code chunk number 6: ni5
###################################################
xytransect<-expand.grid(1:20,1)
nb3=cell2nb(20,1)

plot(nb3,xytransect,col="red",pch=20,cex=2)
title(main="transect of 20 sites")

summary(nb3)


###################################################
### code chunk number 7: ni6
###################################################
set.seed(3)
xyir<-matrix(runif(20),10,2)
plot(xyir,pch=20,cex=1.5,main="irregular sampling with 10 sites")


###################################################
### code chunk number 8: ni7
###################################################
#?dnearneigh
nbnear1<-dnearneigh(xyir,0,0.2)
nbnear2<-dnearneigh(xyir,0,0.3)
nbnear3<-dnearneigh(xyir,0,0.5)
nbnear4<-dnearneigh(xyir,0,1.5)

par(mfrow=c(2,2))
plot(nbnear1,xyir,col="red",pch=20,cex=2)
title(main="neighbors if 0<d<0.2")
plot(nbnear2,xyir,col="red",pch=20,cex=2)
title(main="neighbors if 0<d<0.3")
plot(nbnear3,xyir,col="red",pch=20,cex=2)
title(main="neighbors if 0<d<0.5")
plot(nbnear4,xyir,col="red",pch=20,cex=2)
title(main="neighbors if 0<d<1.5")


###################################################
### code chunk number 9: ni8
###################################################
nbnear1
nbnear2


###################################################
### code chunk number 10: ni9
###################################################
nbnear3
nbnear4


###################################################
### code chunk number 11: ni10
###################################################
#?knearneigh
knn1<-knearneigh(xyir,k=1)
nbknn1<- knn2nb(knn1,sym=T)
plot(nbknn1,xyir,col="red",pch=20,cex=2)
title(main="nearest neighbors (k=1)",cex.main=2)


###################################################
### code chunk number 12: ni11
###################################################
knn2<-knearneigh(xyir,k=2)
nbknn2<- knn2nb(knn2,sym=T)
plot(nbknn2,xyir,col="red",pch=20,cex=2)
title(main="nearest neighbors (k=2)",cex.main=2)


###################################################
### code chunk number 13: ni12
###################################################
n.comp.nb(nbknn1)
n.comp.nb(nbknn2)


###################################################
### code chunk number 14: ni13
###################################################
nbtri<-tri2nb(xyir)
nbgab<-graph2nb(gabrielneigh(xyir), sym=TRUE)
nbrel<-graph2nb(relativeneigh(xyir), sym=TRUE)
nbsoi<-graph2nb(soi.graph(nbtri,xyir), sym=TRUE)

par(mfrow=c(2,2))
plot(nbtri,xyir,col="red",pch=20,cex=2)
title(main="Delaunay triangulation")
plot(nbgab,xyir,col="red",pch=20,cex=2)
title(main="Gabriel Graph")
plot(nbrel,xyir,col="red",pch=20,cex=2)
title(main="Relative Neighbor Graph")
plot(nbsoi,xyir,col="red",pch=20,cex=2)
title(main=" Sphere of Influence Graph")


###################################################
### code chunk number 15: ni14
###################################################
diffnb(nbsoi,nbrel)


###################################################
### code chunk number 16: ni15
###################################################
str(nbsoi)
str(include.self(nbsoi))


###################################################
### code chunk number 17: ni16
###################################################
library(maptools)
columbus <- readShapePoly(system.file("etc/shapes/columbus.shp", package="spdep")[1])
xx <- poly2nb(columbus)
plot(columbus, border="grey")
plot(xx, coordinates(columbus), add=TRUE,pch=20,cex=2,col="red")
title(main="Neighborhood for polygons")


###################################################
### code chunk number 18: ni17
###################################################
#?nb2listw


###################################################
### code chunk number 19: ni18
###################################################
nbgab[[1]]


###################################################
### code chunk number 20: ni19
###################################################
round(dist(xyir),3)


###################################################
### code chunk number 21: ni20
###################################################
distgab<-nbdists(nbgab,xyir)
str(distgab)


###################################################
### code chunk number 22: ni21
###################################################
fdist<-lapply(distgab,function(x) 1-x/max(dist(xyir)))


###################################################
### code chunk number 23: ni22
###################################################
listwgab<-nb2listw(nbgab,glist=fdist,style="B")
listwgab
names(listwgab)
listwgab$neighbours[[1]]
listwgab$weights[[1]]


###################################################
### code chunk number 24: ni23
###################################################
print(listw2mat(listwgab),digits=3)


###################################################
### code chunk number 25: ni24
###################################################
library(spacemakeR)
eigengab=scores.listw(listwgab, echo=TRUE)
barplot(eigengab$values,main="Eigenvalues of spatial weighting matrix")



###################################################
### code chunk number 26: ni25
###################################################
par(mfrow=c(3,3))
for(i in 1:dim(eigengab$vectors)[2])
s.value(xyir,eigengab$vectors[,i],addaxes=F, sub=paste("Eigenvector",i,"(",round(eigengab$values[i],3),")") ,cleg=0,csub=2,neig=neig(list=nbgab))


###################################################
### code chunk number 27: ni25b
###################################################
moranI<-test.scores(eigengab,listwgab,99)
moranI


###################################################
### code chunk number 28: ni25c
###################################################
plot(eigengab$values,moranI$stat,ylab="Moran's I",xlab="Eigenvalues",pch=20,cex=2)
text(eigengab$values,moranI$stat,row.names(moranI),pos=4)
text(-1,0.5,paste("correlation =",cor(moranI$stat,eigengab$values)))
abline(a=0,b=nrow(xyir)/sum(listw2mat(listwgab)),lty=3)


###################################################
### code chunk number 29: ni26
###################################################
signi=which(moranI$p<0.05)
signi
par(mfrow=n2mfrow(length(signi)))
for(i in signi)
s.value(xyir,eigengab$vectors[,i],addaxes=F,sub=paste("ev",i),csub=2,neig=neig(list=nbgab))


###################################################
### code chunk number 30: ni27
###################################################
data(oribatid)
fau <-sqrt(oribatid$fau/outer(apply(oribatid$fau,1,sum), rep(1,ncol(oribatid$fau)),"*"))
faudt <- resid(lm(as.matrix(fau)~as.matrix(oribatid$xy)))


###################################################
### code chunk number 31: ni28
###################################################
nbtri <- tri2nb(as.matrix(oribatid$xy))
sc.tri <- scores.listw(nb2listw(nbtri,style="B"))
AIC.tri <- ortho.AIC(faudt,sc.tri$vec)
AIC.tri


###################################################
### code chunk number 32: ni29
###################################################
min(AIC.tri,na.rm=TRUE)
which.min(AIC.tri)


###################################################
### code chunk number 33: ni30
###################################################
AIC.tri <- ortho.AIC(faudt,sc.tri$vec, ord.var=TRUE)
AIC.tri$AICc[1:10]
AIC.tri$ord[1:10]


###################################################
### code chunk number 34: ni31
###################################################
tri.res <- test.W(faudt,nbtri)
names(tri.res)
names(tri.res$best)


###################################################
### code chunk number 35: ni32
###################################################
f2 <- function(x,dmax,y) {1-(x^y)/(dmax)^y}
maxi <- max(unlist(nbdists(nbtri,as.matrix(oribatid$xy))))
tri.f2 <- test.W(faudt,nbtri, f=f2,y=2:10,dmax=maxi,xy=as.matrix(oribatid$xy))


###################################################
### code chunk number 36: ni33
###################################################
names(tri.f2$best)


###################################################
### code chunk number 37: ni34
###################################################
mvspec<-variogmultiv(faudt,oribatid$xy,nclass=20)
mvspec
plot(mvspec$d,mvspec$var,ty='b',pch=20,xlab="Distance", ylab="C(distance)")


###################################################
### code chunk number 38: ni35
###################################################
dxy <- seq(give.thresh(dist(oribatid$xy)),4,le=10)
nbdnnlist<-lapply(dxy,dnearneigh, x=as.matrix(oribatid$xy),d1=0)


###################################################
### code chunk number 39: ni36
###################################################
dnn.bin <- lapply(nbdnnlist,test.W,Y=faudt)


###################################################
### code chunk number 40: ni37
###################################################
length(dnn.bin)


###################################################
### code chunk number 41: ni38
###################################################
minAIC<-sapply(dnn.bin, function(x) min(x$best$AICc,na.rm=T))


###################################################
### code chunk number 42: ni38b
###################################################
which.min(minAIC)
dxy[which.min(minAIC)]


###################################################
### code chunk number 43: ni39
###################################################
f2 <- function(x,dmax,y) {1-(x^y)/(dmax)^y}


###################################################
### code chunk number 44: ni40
###################################################
dnn.f2 <- lapply(nbdnnlist,function(x) test.W(x,Y=faudt, f=f2,y=2:10,dmax= max(unlist(nbdists(x,as.matrix(oribatid$xy)))),xy=as.matrix(oribatid$xy)))
minAIC<-sapply(dnn.f2, function(x) min(x$best$AICc,na.rm=T))
min(minAIC)
which.min(minAIC)
dnn.f2[[which.min(minAIC)]]$all


###################################################
### code chunk number 45: ni41
###################################################
par(mfrow=c(1,3))
s.value(oribatid$xy, dnn.f2[[7]]$best$vectors[,dnn.f2[[7]]$best$ord[1]],cleg=0,sub=paste("vector",dnn.f2[[7]]$best$ord[1]),csub=3)
s.value(oribatid$xy, dnn.f2[[7]]$best$vectors[,dnn.f2[[7]]$best$ord[2]],cleg=0,sub=paste("vector",dnn.f2[[7]]$best$ord[2]),csub=3)
s.value(oribatid$xy, dnn.f2[[7]]$best$vectors[,dnn.f2[[7]]$best$ord[3]],cleg=0,sub=paste("vector",dnn.f2[[7]]$best$ord[3]),csub=3)