library(spdep) library(akima)#need for interp library(fields)#need for image.plot drumlin=read.table("http://www.public.iastate.edu/~pcaragea/S40608/Compute/drumlin_example.txt",header=TRUE) attach(drumlin) #look at each region's counts summary(cnt1) summary(cnt2) summary(cnt3) #plot these data cnt1.interp=interp(U,V,cnt1) image(cnt1.interp,xlab="U",ylab="V")#this graph is not really showing the areas (i.e. observational units). #To fix this, change the interpolation grid to match the one on which data was collected (here, 11 by 11 lattice) #and since I will be using this output grid over and over, I would like to save it grid.u=seq(min(U),max(U),length=11) grid.v=seq(min(V),max(V),length=11) cnt1.interp=interp(U,V,cnt1,xo=grid.u,yo=grid.v) image(cnt1.interp,xlab="U",ylab="V") #could add the counts on each cell text(U,V,cnt1) #add a title title(main="Region 1, Original counts") #could do this in a gray scale: image(cnt1.interp,xlab="U",ylab="V",col=gray(seq(1,0,length=6))) #note here that white corresponds to counts of zero and black corresponds to 5. If my counts were more diverse than that I could have used a finer grey scale #even better, I can use the command image.plot in the package "fields" and get a legend as well image.plot(cnt1.interp,xlab="U",ylab="V",col=gray(seq(1,0,length=6))) title(main="Region 1, Original counts") #could add the the counts text(U,V,cnt1)#problem is you can't see them on black. Could change the grey scale a little to stop it from going all the way to black (BTW, black is grey(0) and grey(1) is white) image.plot(cnt1.interp,xlab="U",ylab="V",col=gray(seq(1,.3,length=6))) title(main="Region 1, Original counts") #could add the the counts text(U,V,cnt1)#p #to compute any of the indeces of spatial depdendence we need to have the neighborhood---->weights matrix #For a regular lattice we can do this using cell2nb function in spdep #compute neighborhoods; 4 nearest neighbors : rook nb.drum <- cell2nb(11,11)#for an 11 by 11 regular lattice summary(nb.drum) xyc <- attr(nb.drum, "region.id") xy <- matrix(as.integer(unlist(strsplit(xyc, ":"))), ncol=2, byrow=TRUE) plot(nb.drum, xy,cex=2)#neighbors are connected by lines # 8 nearest neighbors: Queen nb.drum.Q <- cell2nb(11,11,type="queen")#for an 11 by 11 regular lattice summary(nb.drum.Q) xyc <- attr(nb.drum.Q, "region.id") xy <- matrix(as.integer(unlist(strsplit(xyc, ":"))), ncol=2, byrow=TRUE) plot(nb.drum.Q, xy,cex=2)#neighbors are connected #compute Moran's I. In spdep, we construct the weights matrix and other important summaries by using the function nb2listw col.W <- nb2listw(nb.drum, style="W")#constructs necessary summaries for calling moran and other functions names(col.W)#it's good to know what it returns #the weights are stored in col.W$weights #by default, it standardizes by row (divides 1 by the number of neighbors) #therefore, the sum of all weights is going to be n*1. Here n is 121 (because we have an 11 by 11 lattice) #you have several options here: W,B,C,U,S #W is the default, row standarized #B is basic binary coding (0 or 1) #C is globally standardized #U is C/nr. of neighbors #S is a variance stabilizing coding scheme #We are now ready to compute Moran's I using the function moran() MI.1=moran(cnt1, col.W, length(nb.drum), Szero(col.W)) MI.2=moran(cnt2, col.W, length(nb.drum), Szero(col.W)) MI.3=moran(cnt3, col.W, length(nb.drum), Szero(col.W)) #plot par(mfrow=c(2,2))#plots a matrix of 2 by 2 graphs moran.plot(cnt1,col.W,main=paste("Moran's I for Region 1:",round(MI.1$I,3))) moran.plot(cnt2,col.W,main=paste("Moran's I for Region 2:",round(MI.2$I,3))) moran.plot(cnt3,col.W,main=paste("Moran's I for Region 3:",round(MI.3$I,3))) # 8 nearest neighbors #compute Moran's I col.WQ <- nb2listw(nb.drum.Q, style="W") MI.1Q=moran(cnt1, col.WQ, length(nb.drum.Q), Szero(col.WQ)) MI.2Q=moran(cnt2, col.WQ, length(nb.drum.Q), Szero(col.WQ)) MI.3Q=moran(cnt3, col.WQ, length(nb.drum.Q), Szero(col.WQ)) #plot par(mfrow=c(2,2)) moran.plot(cnt1,col.WQ,main=paste("Moran's I for Region 1:",round(MI.1Q$I,3))) moran.plot(cnt2,col.WQ,main=paste("Moran's I for Region 2:",round(MI.2Q$I,3))) moran.plot(cnt3,col.WQ,main=paste("Moran's I for Region 3:",round(MI.3Q$I,3))) #try several weighting schemes (B and C, for example) #alternatively, Geary's C: GI.1=geary(cnt1, col.W, n=length(nb.drum),n1=length(nb.drum)-1,Szero(col.W)) GI.2=geary(cnt2, col.W, n=length(nb.drum),n1=length(nb.drum)-1,Szero(col.W)) GI.3=geary(cnt3, col.W, n=length(nb.drum),n1=length(nb.drum)-1,Szero(col.W)) #tests of significance for Moran's I #Null hypothesis: Moran's I is 0 (alternative default >0) T1=moran.test(cnt1,col.W,rank=TRUE)#for region 1 #the p-value is stored in T1$p.value and can be interpreted in the usual way #note here that "rank" is FALSE by default (for continuous data), and should be changed to TRUE for data that is not continuous. It uses an adaptation of Moran's I T1$p.value #in theory, a small p-value is an indication of spatial dependence #you could use a Monte Carlo based test, in which you need to specify the number of permutations of the data you'd like to look at T1.mc=moran.mc(cnt1,col.W,nsim=1000) #again, the p-value is stored under p.value. This p-value is based on the rank of the observed Moran's I in the series of permutations you have performed. Therefore, you need to have a large number of permutations to ensure T1.mc$p.value #now for regions 2 and 3 T2=moran.test(cnt2,col.W,rank=TRUE)#for region 2 T2$p.value T3=moran.test(cnt3,col.W,rank=TRUE)#for region 3 T3$p.value #and monte carlo tests T2.mc=moran.mc(cnt2,col.W,nsim=1000) T2.mc$p.value T3.mc=moran.mc(cnt3,col.W,nsim=1000) T3.mc$p.value #you could do the same type of testing for Geary's C #the functions would be geary.test and geary.mc #local moran's I (LISA) T1.LISA=localmoran(cnt1,col.W) T2.LISA=localmoran(cnt2,col.W) T3.LISA=localmoran(cnt3,col.W) #stores statistics on the first column, the z-score on the fourth and p-value on the fifth #when the local Moran I statistics are positive we have an indication of high-high or low-low patterns (i.e. positive spatial dependence), while negative values indicate high-low or low-high patterns (i.e. negative spatial dependence) #Moran's I is the average of the local Moran's I statistics #plot the surface of z-scores. Interpolate first T1.z=interp(U,V,T1.LISA[,4],xo=grid.u,yo=grid.v) image.plot(T1.z,ylab="V",xlab="U")#uses the standard color scheme title(main="Region 1, map of LISA statistics") image.plot(T1.z,col=cm.colors(30),ylab="V",xlab="U")#another color scheme #we can do the same thing for p-values and obtain what is often called a LISA significance map T1.p=interp(U,V,T1.LISA[,5],xo=grid.u,yo=grid.v) image.plot(T1.p,col=grey(seq(1,0,l=10)),breaks=seq(0,1,l=11),ylab="V",xlab="U") title(main="Region 1, Significance map") #it may be helpful to look at the p-values directly on this map text(U,V,round(T1.LISA[,5],2)) #a cluster map is a combination of the two previous maps:significant locations are color coded by type of spatial autocorrelation: high-high, low-low, high-low and low-high #but I haven't figured out how to get this in R in the simplest way. I'll look into this a little more. #could look at the boxplot of local Moran statistics (to see if any locations show very different local autocorrelation patterns) boxplot(T1.LISA[,1],T2.LISA[,1],T3.LISA[,1])#side-by-side boxplots for the three regions of the study #since Moran's I is the average of these local statistics, we should be weary of skewed distributions or a large number of outliers.