library(tidyverse) library(raster) date="2-15-2021" run=3 V= "G" # Do you want to use the good and marginal EVI observations (GM) or just the good (G) data.grouped=read.table("/home/wh350/climate_variability/evi_sddata2-15-20213.grouped.txt") data.grouped=data.grouped%>%filter(temp_sd.1950!=0&temp_sd.2000!=0 & pre_sd.1950!=0&pre_sd.2000!=0 & vpd_sd.1950!=0&vpd_sd.2000!=0) data.grouped=data.grouped %>% filter(real_change.temp>0 | real_change.pre>0 |real_change.vpd>0) spatial.weights = function(data.x, data.y){ myco=cbind(data.x, data.y) s=sp::SpatialPoints(myco,proj4string = crs("+proj=longlat +ellps=clrk66 +no_defs")) myco.knn=spdep::knearneigh(s,k=24,longlat = T) myco.nb=spdep::knn2nb(myco.knn) myco.l2=spdep::nb2listw(myco.nb) return(myco.l2) } spatial.weights=spatial.weights(data.x = data.grouped$x, data.y = data.grouped$y) bc <- MASS::boxcox(data.grouped$evi_sd ~ scale(poly(data.grouped$temp_sd.2000,2)) + scale(poly(data.grouped$pre_sd.2000,2)) +scale(poly(data.grouped$vpd_sd.2000,2))+data.grouped$leaf.habit+data.grouped$leaf.shape+data.grouped$biome) (lambda <- bc$x[which.max(bc$y)]) lm.modern = lm(scale(((evi_sd^lambda - 1)/lambda)) ~ biome + leaf.shape + scale(poly(temp_sd.2000, 2)) + scale(poly(vpd_sd.2000, 2)),data=data.grouped) print(summary(lm.modern)) sp.error.modern <- spatialreg::errorsarlm(scale(((evi_sd^lambda - 1)/lambda)) ~ biome + leaf.shape + scale(poly(temp_sd.2000, 2)) + scale(poly(vpd_sd.2000, 2)), data=data.grouped, spatial.weights) print(summary(sp.error.modern)) sp.lag.modern <- spatialreg::lagsarlm(scale(((evi_sd^lambda - 1)/lambda)) ~ biome + leaf.shape + scale(poly(temp_sd.2000, 2)) + scale(poly(vpd_sd.2000, 2)),data=data.grouped, spatial.weights) print(summary(sp.lag.modern)) print(spdep::moran.mc(resid(lm.modern), spatial.weights, nsim=999, zero.policy=TRUE)) print(spdep::moran.mc(resid(sp.error.modern), spatial.weights, nsim=999, zero.policy=TRUE)) print(spdep::moran.mc(resid(sp.lag.modern), spatial.weights, nsim=999, zero.policy=TRUE)) print("Finished modern analysis")