#read in data set for presence/absence of two plants
presabs.plants=read.table("http://www.public.iastate.edu/~pcaragea/S40608/Compute/plants.txt",header=TRUE))
attach(presabs.plants)
#names here are of the form S68 or H78 which indicate the
#plant type (Sorrell or Honeysuckle) and the year (1968 through 1978)
yspan=1968:1978
#plot the locations of the data, stored in xccord and ycoord
plot(xcoord,ycoord,pch=20)


#construct neighborhood using the distance
#it finds all locations farther than d1 but closer than d2. In general, set d1=0.
nb.pl=dnearneigh(x=locs.FE,d1=0,d2=13)
#this function created a neighborhood object using just the x and y coordinates of the data

#next, create the weights object based on this neighborhood
col.W <- nb2listw(nb.pl, style="W")#constructs necessary summaries for calling moran and other functions


#with these, we can now compute Moran's I for each year of the data.

Observed=S70;i=3
MI.PL=moran(Observed, col.W, length(nb.pl), Szero(col.W))
moran.plot(Observed,col.W,xlab="Obs",ylab="Neighb")
title(main=paste("Moran's I for year", yspan[i],"=",round(MI.PL$I,3)))
moran.test(Observed,col.W,rank=TRUE)
moran.mc(Observed,col.W,nsim=999)


#or we can apply another test, the join count test. 
#tests whether same-colour joins (in this case we have two colors, 0 or 1) occur more frequently than would be expected if the zones were labelled in a spatially random way. 

 
#join count statistics. Note that the data needs to be a factor.
joincount.test(factor(S70),col.W)
#could repeat this for each year, each plant species