##### CODE FOR ADJUSTING FOR SPATIAL AUTOCORRELATION IN HURRICANE SANDY RESIDUALS ##### Note just illustrating the principle with a submodel (response + parents) rather than a lavaan model # load example data file sandy.dat <- read.csv("SEM_11_4_Spatial_Autocorr_data.csv") names(sandy.dat) dim(sandy.dat) #shows there are 107 rows of data ### Model for adjustment (selected by AICc model comparisons in prior analyses) # mod1 mod1 <-lm(change ~ dist.code + I(dist.code^2) + pos + rate.c + surge, data = sandy.dat) summary(mod1) ### Capture residuals change.res1 <- resid(mod1) ### Determine nearest neighbors from x, y coordinates ## Create x,y coordinate matrix x <- sandy.dat$latitude y <- sandy.dat$longitude xydat <- cbind(x, y) ### Load needed library library(spdep) vignette(package = "spdep") # identify k= nearest neighbors xy.knn <- knearneigh(xydat,k=2) # returns neighbors list xy.nb <- knn2nb(xy.knn) ### Moran's test of model residuals (for mod1) {requires residuals and neighbors list} MoransResults <- moran.test(change.res1, nb2listw(xy.nb, style="W")) print(MoransResults) ### Computation of effective sample size n = 107 eff_ss = round(n/(1+abs(MoransResults$estimate[1]))) print(eff_ss) ##### Adjust standard errors of model (mod1) n = 107; n.adj = 91 # Formula for adjustment se.adj = sqrt(n/n.adj) # Model mod1 <-lm(change ~ dist.code + I(dist.code^2) + pos + rate.c + surge, data = sandy.dat) # Extract standard errors std.errors <- coef(summary(mod1))[, 2] # Adjust standard errors se1.adj = std.errors[1]*se.adj; se2.adj = std.errors[2]*se.adj; se3.adj = std.errors[3]*se.adj; se4.adj = std.errors[4]*se.adj; se5.adj = std.errors[5]*se.adj; se6.adj = std.errors[6]*se.adj; # Print list adjusted standard errors std.errors.adjusted <- data.frame(se1.adj, se2.adj, se3.adj, se4.adj, se5.adj, se6.adj) std.errors.adj.t <- t(std.errors.adjusted) print(std.errors.adj.t) ### Calculate adjusted p-values for mod1 # get parameter estimates mod1.ests <- coef(summary(mod1))[, 1] ; mod1.ests mod1.est1 <- mod1.ests[1] # compute adjusted p-value for intercept se1.adj = 9.8787912 #from previous step n.adj = 91 t1.adj = mod1.est1/se1.adj; p1.adj <- 2*pt(-abs(t1.adj),df=n.adj-1); round(p1.adj, 4) # compare to unadjusted (just for information sake) round(coef(summary(mod1))[, 4], 4) # proportional increase in p-value 0.9869/0.9858 #very very slight affect ##### Adjusting AICc Comparisons ### Competing Models # Model 1 mod1 <-lm(change ~ dist.code + I(dist.code^2) + pos + rate.c + surge, data = sandy.dat) # Model 2 - drop position variable mod2 <-lm(change ~ dist.code + I(dist.code^2) + rate.c + surge, data = sandy.dat) # AICc table library(AICcmodavg) # compare models using AIC aictab(list(mod1, mod2), c("Model 1", "Model 2")) ### Adjustment # Extract AIC values (these are not influenced by sample size) mod1.AIC <- AIC(mod1) #AIC = 843.588 mod2.AIC <- AIC(mod2) #AIC = 851.877 # Bring in adjusted sample size estimate n.adj = 91 # Extract number of parameters (K) in model from AICc table mod1.num.params <- aictab(list(mod1), "Model 1")$K #K1 = 7 mod2.num.params <- aictab(list(mod2), "Model 2")$K #k2 = 6 # Compute adjusted AICc mod1.AICc.adj <- mod1.AIC +((2*mod1.num.params*(mod1.num.params+1)) /(n.adj - mod1.num.params-1)) mod2.AICc.adj <- mod2.AIC +((2*mod2.num.params*(mod2.num.params+1)) /(n.adj - mod2.num.params-1)) # Compute AICc difference (model 1 minus model 2) mod1.AICc.adj - mod2.AICc.adj