#Load required librarries into R library(raster); library(spdep); library(ncf) ; library(usdm) setwd("C:/Users/cd833/Desktop/docs/newpapers/disease/matlabnew") library(XLConnect) # load XLConnect package wk = loadWorkbook("datz1.xls") df = readWorksheet(wk, sheet="Sheet1") lat=df[1] lon=df[2] x1=df[3] x2=df[4] x3=df[5] x4=df[6] x5=df[7] x6=df[8] x7=df[9] x8=df[10] x9=df[11] x10=df[12] x11=df[13] x12=df[14] x13=df[15] x14=df[16] x15=df[2] x16=df[17] x17=df[18] x18=df[19] y1 <-as.numeric(x1[[1]]) y2 <-as.numeric(x2[[1]]) y3 <-as.numeric(x3[[1]]) y4 <-as.numeric(x4[[1]]) y5 <-as.numeric(x5[[1]]) y6 <-as.numeric(x6[[1]]) y7 <-as.numeric(x7[[1]]) y8 <-as.numeric(x8[[1]]) y9 <-as.numeric(x9[[1]]) y10 <-as.numeric(x10[[1]]) y11 <-as.numeric(x11[[1]]) y12 <-as.numeric(x12[[1]]) y13 <-as.numeric(x13[[1]]) y14 <-as.numeric(x14[[1]]) y15 <-as.numeric(x15[[1]]) y16 <-as.numeric(x16[[1]]) y17 <-as.numeric(x17[[1]]) y18 <-as.numeric(x18[[1]]) #y1, y4, y6, y8 #Load COOR as a table of coordinates (in lat long) this is a dummy case #COOR=cbind( runif(31,120,121), runif(31,50,51)) COOR=coordinates(cbind( lat, lon)) # JID # popden # pop change # perc ag # past/current home range # past/current fecal diffusivity # past home range # past fecal diffusivity # %current home range # %current fecal diffusivity # mean rainfall # change in rainfall # current species richness # past species richness # latitude # current biomass weighted species richness (N*weight) # change in biomass weighted species richness #table S3 model FD dz = cbind(y1, y2, y3, y6) #table S3 model HR dz = cbind(y1, y2, y3, y5) #table S3 model null dz = cbind(y1, y2, y3) #Load the data into R, here a dummy variable and a set of dummy results PARAM=dz RESULTS=y16 #Create potential neighborhoods Neig5=knn2nb(knearneigh(COOR,5, longlat = T)); W5 = nb2listw(Neig5,style="S",zero.policy =T) Neig10=dnearneigh(COOR,0,10, longlat = T); W10 = nb2listw(Neig10,style="S",zero.policy =T) Neig40=dnearneigh(COOR,0,40, longlat = T); W40 = nb2listw(Neig40,style="S",zero.policy =T) Neig80=dnearneigh(COOR,0,80, longlat = T); W80 = nb2listw(Neig80,style="S",zero.policy =T) Neig100=dnearneigh(COOR,0,100, longlat = T); W100 = nb2listw(Neig100,style="S",zero.policy =T) Neig150=dnearneigh(COOR,0,150, longlat = T); W150 = nb2listw(Neig150,style="S",zero.policy =T) Neig180=dnearneigh(COOR,0,180, longlat = T); W180 = nb2listw(Neig180,style="S",zero.policy =T) Neig200=dnearneigh(COOR,0,200, longlat = T); W200 = nb2listw(Neig200,style="S",zero.policy =T) Neig220=dnearneigh(COOR,0,220, longlat = T); W220 = nb2listw(Neig220,style="S",zero.policy =T) Neig250=dnearneigh(COOR,0,250, longlat = T); W250 = nb2listw(Neig250,style="S",zero.policy =T) Neig300=dnearneigh(COOR,0,300, longlat = T); W300 = nb2listw(Neig300,style="S",zero.policy =T) NEIG250W=nb2listw(dnearneigh(COOR, 0,250, longlat = T),style="W",zero.policy =T) NEIG250B=nb2listw(dnearneigh(COOR, 0,250, longlat = T),style="B",zero.policy =T) NEIG250C=nb2listw(dnearneigh(COOR, 0,250, longlat = T),style="C",zero.policy =T) NEIG250U=nb2listw(dnearneigh(COOR, 0,250, longlat = T),style="U",zero.policy =T) NEIG250S=nb2listw(dnearneigh(COOR, 0,250, longlat = T),style="S",zero.policy =T) #Compare the models for a lm with the different neighborhoods M0=lm(RESULTS~PARAM) M100=errorsarlm(RESULTS~PARAM, ,listw=W100, tol.solve = 1e-12, zero.policy=T) M150=errorsarlm(RESULTS~PARAM, ,listw=W150, tol.solve = 1e-12, zero.policy=T) M180=errorsarlm(RESULTS~PARAM, ,listw=W180, tol.solve = 1e-12, zero.policy=T) M200=errorsarlm(RESULTS~PARAM, ,listw=W200, tol.solve = 1e-12, zero.policy=T) M220=errorsarlm(RESULTS~PARAM, ,listw=W220, tol.solve = 1e-12, zero.policy=T) M250=errorsarlm(RESULTS~PARAM, ,listw=W250, tol.solve = 1e-12, zero.policy=T) M300=errorsarlm(RESULTS~PARAM, ,listw=W300, tol.solve = 1e-12, zero.policy=T) M5=errorsarlm(RESULTS~PARAM, ,listw=W5, tol.solve = 1e-12, zero.policy=T) M10=errorsarlm(RESULTS~PARAM, ,listw=W10, tol.solve = 1e-12,zero.policy =T) M40=errorsarlm(RESULTS~PARAM, ,listw=W40, tol.solve = 1e-12, zero.policy=T) #M75=errorsarlm(RESULTS~PARAM, ,listw=W75, tol.solve = 1e-12, zero.policy=T) M80=errorsarlm(RESULTS~PARAM, ,listw=W80, tol.solve = 1e-12, zero.policy=T) #M85=errorsarlm(RESULTS~PARAM, ,listw=W85, tol.solve = 1e-12, zero.policy=T) M95=errorsarlm(RESULTS~PARAM, ,listw=W95, tol.solve = 1e-12, zero.policy=T) #Find best model #AIC(M0, M20, M40, M75, M80, M85, M90, M95, M160) AIC( M5, M10, M40, M80, M95, M100, M150, M180, M200, M220, M250, M300) #The next are based on M3 being the best. In this dummy example. None of them make any sense which is ovbious by looking at summary(errorsarlm(RESULTS~PARAM, listw=W200, tol.solve = 1e-12, zero.policy=T)) summary(lm(RESULTS~PARAM)) vif(PARAM) v1 <- vifcor(PARAM, th=0.9) # identify collinear variables that should be excluded #Plot correlelogram. The important element is a low autocorrelation for the first steps #corr = correlog(COOR[,1], COOR[,2],M200$residuals, increment =300, resamp=100,latlon=T) #corr #plot(corr) #Find the pseudo R2 of explained by parameters, neighborhood or both #The funciton predict.sarlm(M3) gives the predictions for fit which is the sum of spatial and paramters, trend which is the parameter based and signal which is the spatial. #It has an annoying format which i generally handle by exporting and importing it. write.table(predict.sarlm(M200, zero.policy=T), "TEMPORARY_SAR.txt"); read.table("TEMPORARY_SAR.txt")->PREDICT;file.remove("TEMPORARY_SAR.txt") PseudoSQ_Total=cor(RESULTS, PREDICT[,1])^2 PseudoSQ_Param=cor(RESULTS, PREDICT[,2])^2 PseudoSQ_Spatial=cor(RESULTS, PREDICT[,3])^2 v1 PseudoSQ_Param x<-M200$residuals write.csv(x,'test.csv') x1<-M200$residuals write.csv(x1,'test2.csv')