# load package spdep – need to do from GUI drop-down menu# ### getting the data in ### lsg <- read.csv(file.choose(),header=TRUE,skip=1) attach (lsg) summary(lsg) ### neighbor matrix and inverse dist Wij ### #create a neighbor list, calling all pairs neighbors that are within 0km and the extent of the study area - 300km lsg.nb <- dnearneigh (cbind (X, Y), 0, 300000) #calculate neighbor distances between all neighbors identified by the neighbor list (i.e. all pairs as defined above) nbd <- nbdists(lsg.nb, cbind (X,Y)) #create a geographic weight list that weights each connection in the neighbor list according to a function - here the inverse of the distance calculated in the step above gl <- lapply(nbd, function(x) 1/x) #create a neighbor weighting list - here includes the list of neighbors (defining all pairs within the study extent to be neighbors) and applying a wieght list - here based on the inverse distance weighting calculated above based on distances between all pairs lsg.listw <- nb2listw(lsg.nb, glist=gl) ### run spatial models ### # construct spatial lag regression as with normal regression, but here the weights matrix is part of the data - the lag model command is what determines how the weights go into the model # in summary Rho gives the spatial autoregressive parameter # and LR gives results of a likelihood ratio test of the model with and without the Rho parameter # use these to determine the significance of the contribution of autocorrelation to the model # Rho can be estimated with more precision and greater significance than if the variable is manually lagged and included directly in a standard regression (i.e. as with directly including a parameter to modify prevalence by distance from core) # model for sPCA1 #extremes of the model - null and full lsg1.p <- lagsarlm ( IDaxis1 ~ 1 , data = lsg, lsg.listw) lsg1.full <- lagsarlm ( IDaxis1 ~ ECO + CAN + ADJ + RIV + INT + HWY , data = lsg, lsg.listw) #compare AIC within set AIC (lsg1.p) AIC (lsg1.full) #PCA1 models with HABITAT variables lsg1.X2 <- lagsarlm ( IDaxis1 ~ ECO + CAN + ADJ , data = lsg, lsg.listw) lsg1.X3 <- lagsarlm ( IDaxis1 ~ ECO + CAN , data = lsg, lsg.listw) lsg1.X4 <- lagsarlm ( IDaxis1 ~ ECO + ADJ , data = lsg, lsg.listw) lsg1.X5 <- lagsarlm ( IDaxis1 ~ CAN + ADJ , data = lsg, lsg.listw) lsg1.X6 <- lagsarlm ( IDaxis1 ~ ECO , data = lsg, lsg.listw) lsg1.X7 <- lagsarlm ( IDaxis1 ~ CAN , data = lsg, lsg.listw) lsg1.X8 <- lagsarlm ( IDaxis1 ~ ADJ , data = lsg, lsg.listw) #compare AIC within HABITAT set AIC (lsg1.X2) AIC (lsg1.X3) AIC (lsg1.X4) AIC (lsg1.X5) AIC (lsg1.X6) AIC (lsg1.X7) AIC (lsg1.X8) #sPCA1 models with BARRIER variables lsg1.X9 <- lagsarlm ( IDaxis1 ~ INT + HWY + RIV , data = lsg, lsg.listw) lsg1.X10 <- lagsarlm ( IDaxis1 ~ INT + HWY , data = lsg, lsg.listw) lsg1.X11 <- lagsarlm ( IDaxis1 ~ INT + RIV , data = lsg, lsg.listw) lsg1.X12 <- lagsarlm ( IDaxis1 ~ HWY + RIV , data = lsg, lsg.listw) lsg1.X13 <- lagsarlm ( IDaxis1 ~ INT , data = lsg, lsg.listw) lsg1.X14 <- lagsarlm ( IDaxis1 ~ HWY , data = lsg, lsg.listw) lsg1.X15 <- lagsarlm ( IDaxis1 ~ RIV , data = lsg, lsg.listw) #compare AIC within BARRIER set AIC (lsg1.X9) AIC (lsg1.X10) AIC (lsg1.X11) AIC (lsg1.X12) AIC (lsg1.X13) AIC (lsg1.X14) AIC (lsg1.X15) #best COMBO's after initial screening lsg1.C1 <- lagsarlm ( IDaxis1 ~ + INT + HWY + RIV + ECO + CAN , data = lsg, lsg.listw) lsg1.C2 <- lagsarlm ( IDaxis1 ~ + INT + HWY + RIV + ECO + ADJ , data = lsg, lsg.listw) lsg1.C3 <- lagsarlm ( IDaxis1 ~ + INT + HWY + RIV + ECO , data = lsg, lsg.listw) #compare AIC within COMBO set AIC (lsg1.C1) AIC (lsg1.C2) AIC (lsg1.C3) #summary and parameter output for top contenders (within 2AIC of top model) summary(lsg1.X9 ) summary(lsg1.C1 ) summary(lsg1.C3 ) # model for sPCA2 #extremes of the model - null and full lsg2.p <- lagsarlm ( IDaxis2 ~ 1 , data = lsg, lsg.listw) lsg2.full <- lagsarlm ( IDaxis2 ~ ECO + CAN + ADJ + RIV + INT + HWY , data = lsg, lsg.listw) # compare AIC within set AIC (lsg2.p) AIC (lsg2.full) #PCA1 models with HABITAT variables lsg2.X2 <- lagsarlm ( IDaxis2 ~ ECO + CAN + ADJ , data = lsg, lsg.listw) lsg2.X3 <- lagsarlm ( IDaxis2 ~ ECO + CAN , data = lsg, lsg.listw) lsg2.X4 <- lagsarlm ( IDaxis2 ~ ECO + ADJ , data = lsg, lsg.listw) lsg2.X5 <- lagsarlm ( IDaxis2 ~ CAN + ADJ , data = lsg, lsg.listw) lsg2.X6 <- lagsarlm ( IDaxis2 ~ ECO , data = lsg, lsg.listw) lsg2.X7 <- lagsarlm ( IDaxis2 ~ CAN , data = lsg, lsg.listw) lsg2.X8 <- lagsarlm ( IDaxis2 ~ ADJ , data = lsg, lsg.listw) #compare AIC within HABITAT set AIC (lsg2.X2) AIC (lsg2.X3) AIC (lsg2.X4) AIC (lsg2.X5) AIC (lsg2.X6) AIC (lsg2.X7) AIC (lsg2.X8) #sPCA1 models with BARRIER variables lsg2.X9 <- lagsarlm ( IDaxis2 ~ INT + HWY + RIV , data = lsg, lsg.listw) lsg2.X10 <- lagsarlm ( IDaxis2 ~ INT + HWY , data = lsg, lsg.listw) lsg2.X11 <- lagsarlm ( IDaxis2 ~ INT + RIV , data = lsg, lsg.listw) lsg2.X12 <- lagsarlm ( IDaxis2 ~ HWY + RIV , data = lsg, lsg.listw) lsg2.X13 <- lagsarlm ( IDaxis2 ~ INT , data = lsg, lsg.listw) lsg2.X14 <- lagsarlm ( IDaxis2 ~ HWY , data = lsg, lsg.listw) lsg2.X15 <- lagsarlm ( IDaxis2 ~ RIV , data = lsg, lsg.listw) #sPCA1 models with BARRIER variables AIC (lsg2.X9) AIC (lsg2.X10) AIC (lsg2.X11) AIC (lsg2.X12) AIC (lsg2.X13) AIC (lsg2.X14) AIC (lsg2.X15) #best COMBO's after initial screening lsg2.C1 <- lagsarlm ( IDaxis2 ~ + INT + HWY + RIV + ECO + CAN + ADJ , data = lsg, lsg.listw) lsg2.C2 <- lagsarlm ( IDaxis2 ~ + INT + HWY + RIV + ECO + ADJ , data = lsg, lsg.listw) lsg2.C3 <- lagsarlm ( IDaxis2 ~ , data = lsg, lsg.listw) #compare AIC within COMBO set AIC (lsg2.C1) AIC (lsg2.C2) AIC (lsg2.C3) #summary and parameter output for top contenders (within 2AIC of top model) summary(lsg2.X9 ) summary(lsg2.C1 ) summary(lsg2.C2 ) ### test autocorrelation of residuals ### # requires regression object and listw object # use regression object and weights object for moran test lsg1.moran <- lm.morantest (lsg1.C1, lsg.listw) print (lsg1.moran) lsg2.moran <- lm.morantest (lsg2.C1, lsg.listw) print (lsg2.moran) # or for 2 sided test name.moran <- lm.morantest (name.lm, lsg.listw, alternative = "two.sided") print (name.moran)