##Landscape Structure Influences Urban Vegetation Vertical Structure ##Matthew Mitchell, Dan Wu, Kasper Johansen, Martine Maron, Clive McAlpine & Jonathan Rhodes ##University of Queensland ##School of Geography, Planning and Environmental Management ##June 24, 2016 ##set the working directory to 1km scale data setwd("~/Desktop/UQ/ARC ES Project/Brisbane Veg Structure Analysis/Current Data") ##Load required libraries for model selection analysis library(car) library(MuMIn) library(gstat) library(spdep) library(arm) source("SAR_Script.R") #import 1ha resolution data hdata <- read.csv("data_1ha_50_20150815.csv") hdata.1 <- hdata hdata.1$X <- NULL ##Gets rid of extra variable in the data hdata.1$road_conn[is.na(hdata.1$road_conn)] <- 0 #Turns NAs in road connectivity to zeros hdata.1<-hdata.1[complete.cases(hdata.1),] #removes 1ha cells that have missing data (could not be calculated) hdata.1 <- hdata.1[(hdata.1$clumpy>-0.5),] #removes outliers with very small clumpiness values hdata.1 <- hdata.1[(hdata.1$mb_dwel_dens<400),] #removes outlier where dwelling density is very large hdata.1 <- hdata.1[(hdata.1$lot_size!=0),] #removes 1ha cells where lot size is zero ##Identify any 1ha cells without neighbours and remove from analysis coords <- as.matrix(cbind(hdata.1$x, hdata.1$y)) #creates coordinate dataset of 1ha cells neigh <- dnearneigh(coords, 0, 150) #calculated neighbouring cells for dataset based on 150m distance between cell centers neigh #shows which 1ha cells are missing neighbours and must be removed ##Remove these 1ha cells hdata.1 <- hdata.1[c(1:19133,19135:28123,28125:39723,39725:43316,43318:52420,52422:52623),] ##################################################################### #####Model selection for vegetation density between 0.15-1m ##################################################################### ####################################### #### Physical variables - vegetation density 0.15-1m #### ##Rescale variables hdata.1$aspect_cos.s <- rescale(hdata.1$aspect_cos) hdata.1$aspect_sin.s <- rescale(hdata.1$aspect_sin) hdata.1$elev.s <- rescale(hdata.1$elev) hdata.1$slope.s <- rescale(hdata.1$slope) ##Create null, linear, and polynomial models h.phys.0151.lm <- lm(log(dens_015_1+0.01)~1, data=hdata.1) h.phys.0151.lm.1 <- lm(log(dens_015_1+0.01)~aspect_cos.s+aspect_sin.s+elev.s+slope.s, data=hdata.1) h.phys.0151.lm.2 <- lm(log(dens_015_1+0.01)~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2), data=hdata.1) ##Comparison three models from above to choose best based on AICc and look at VIFs comp.models(h.phys.0151.lm.1, h.phys.0151.lm.2) ##Create error, lag, and mixed SAR models, compare using AIC and ANOVA, and calculate Moran's I for each h.sar.phys.0151 <- mm.spatial.sparse(hdata.1, 150, h.phys.0151.lm.2, "MC") ##View results mm.sum.sp(h.sar.phys.0151) ##Plot residuals against fitted values, residuals spatially, and semivariograms spatial.plot(h.phys.0151.lm.2, hdata.1, "0.15-1m Veg Density - Polynomial physical non-spatial model") spatial.plot(h.sar.phys.0151$mixed, hdata.1, "0.15-1m Veg Density - Mixed physical SAR model") ##Create best model and run diagnostics h.best.phys.0151 <- lagsarlm(log(dens_015_1+0.01) ~ poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2), listw=h.sar.phys.0151$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") #AICc value AICc(h.best.phys.0151) ##summary of best model and Nagelkerke R2 summary(h.best.phys.0151, Nagelkerke=TRUE) ##Plots of residuals vs. fitted values, residuals spatially, and semivariogram diag.best(best.phys.0151, data.1, sar.phys.0151$weight, "0.15-1m Veg Density - Best SAR model") ##Creation of null spatial model (accounts for spatial autoregression only) h.best.phys.0151.null <- lagsarlm(log(dens_015_1+0.01) ~ 1, listw=h.sar.phys.0151$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") ##Calculation of Nagelkerke R2 values spatial.diag(h.phys.0151.lm, h.phys.0151.lm.2, h.best.phys.0151, h.best.phys.0151.null, hdata.1) ##Raw SAR model for plotting h.best.phys.0151.raw <- lagsarlm(log(dens_015_1+0.01) ~ poly(aspect_cos.s,2,raw=T)+poly(aspect_sin.s,2,raw=T)+ poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T), listw=h.sar.phys.0151$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") ##Partial residual plots plot.partial(h.best.phys.0151, h.best.phys.0151.raw, hdata.1, 1, "aspect_cos", "aspect_cos.s", "Northness", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.phys.0151, h.best.phys.0151.raw, hdata.1, 2, "aspect_sin", "aspect_sin.s", "Eastness", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.phys.0151, h.best.phys.0151.raw, hdata.1, 3, "elev", "elev.s", "Elevation (m)", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.phys.0151, h.best.phys.0151.raw, hdata.1, 4, "slope", "slope.s", "Slope (degrees)", "0.15-1m Veg. Density Partial Residuals") ####################################### ####Soil variables - vegetation density 0.15-1m #### hdata.1$water_cap.s <- rescale(hdata.1$water_cap) hdata.1$tot_p.s <- rescale(hdata.1$tot_p) hdata.1$soc.s <- rescale(hdata.1$soc) hdata.1$tot_n.s <- rescale(hdata.1$tot_n) h.soil.0151.lm <- lm(log(dens_015_1+0.01)~1, data=hdata.1) h.soil.0151.lm.1 <- lm(log(dens_015_1+0.01)~water_cap.s+tot_p.s+soc.s+tot_n.s, data=hdata.1) h.soil.0151.lm.2 <- lm(log(dens_015_1+0.01)~poly(water_cap.s,2)+poly(tot_p.s,2)+poly(soc.s,2)+poly(tot_n.s,2), data=hdata.1) comp.models(h.soil.0151.lm.1, h.soil.0151.lm.2) h.sar.soil.0151 <- mm.spatial.sparse(hdata.1, 150, h.soil.0151.lm.2, "MC") mm.sum.sp(h.sar.soil.0151) spatial.plot(soil.0151.lm.2, data.1, "0.15-1m Veg Density - Soil polynomial non-spatial model") spatial.plot(sar.soil.0151$mixed, data.1, "0.15-1m Veg Density - Soil mixed SAR model") h.best.soil.0151 <- lagsarlm(log(dens_015_1+0.01) ~ poly(soc.s,2)+poly(tot_n.s,2)+poly(water_cap.s,2)+poly(tot_p.s,2), listw=h.sar.soil.0151$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.soil.0151) summary(h.best.soil.0151, Nagelkerke=TRUE) diag.best(best.soil.0151, data.1, sar.soil.0151$weight, "0.15-1m Veg Density - Best Soil SAR model") h.best.soil.0151.null <- lagsarlm(log(dens_015_1+0.01) ~ 1, listw=h.sar.soil.0151$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.soil.0151.lm, h.soil.0151.lm.2, h.best.soil.0151, h.best.soil.0151.null, hdata.1) h.best.soil.0151.raw <- lagsarlm(log(dens_015_1+0.01) ~ poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+ poly(water_cap.s,2,raw=T)+poly(tot_p.s,2,raw=T), listw=h.sar.soil.0151$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") plot.partial(h.best.soil.0151, h.best.soil.0151.raw, hdata.1, 1, "soc", "soc.s", "Soil organic carbon (%)", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.soil.0151, h.best.soil.0151.raw, hdata.1, 2, "tot_n", "tot_n.s", "Total soil N (%)", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.soil.0151, h.best.soil.0151.raw, hdata.1, 3, "water_cap", "water_cap.s", "Soil water capacity (%)", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.soil.0151, h.best.soil.0151.raw, hdata.1, 4, "tot_p", "tot_p.s", "Total soil P (%)", "0.15-1m Veg. Density Partial Residuals") ####################################### ####Census variables - vegetation density 0.15-1m #### hdata.1$prop_noncauc.s <- rescale(hdata.1$prop_noncauc) hdata.1$sa1_medage.s <- rescale(hdata.1$sa1_medage) hdata.1$sa1_medtothinc.s <- rescale(hdata.1$sa1_medtothinc) hdata.1$sa1_avghouse.s <- rescale(hdata.1$sa1_avghouse) h.cens.0151.lm <- lm(log(dens_015_1+0.01)~1, data=hdata.1) h.cens.0151.lm.1 <- lm(log(dens_015_1+0.01)~prop_noncauc.s+sa1_medage.s+sa1_medtothinc.s+sa1_avghouse.s, data=hdata.1) h.cens.0151.lm.2 <- lm(log(dens_015_1+0.01)~poly(prop_noncauc.s,2)+poly(sa1_medage.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2), data=hdata.1) comp.models(h.cens.0151.lm.1, h.cens.0151.lm.2) h.sar.cens.0151 <- mm.spatial.sparse(hdata.1, 150, h.cens.0151.lm.2, "MC") mm.sum.sp(h.sar.cens.0151) spatial.plot(cens.0151.lm.2, data.1, "0.15-1m Veg Density - Census polynomial non-spatial model") spatial.plot(sar.cens.0151$mixed, data.1, "0.15-1m Veg Density - Census mixed SAR model") h.best.cens.0151 <- lagsarlm(log(dens_015_1+0.01) ~ poly(prop_noncauc.s,2)+poly(sa1_medage.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2), listw=h.sar.cens.0151$weight, type="lag", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.cens.0151) summary(h.best.cens.0151, Nagelkerke=TRUE) diag.best(best.cens.0151, data.1, sar.cens.0151$weight, "0.15-1m Veg Density - Best Soil SAR model") h.best.cens.0151.null <- lagsarlm(log(dens_015_1+0.01) ~ 1, listw=h.sar.cens.0151$weight, type="lag", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.cens.0151.lm, h.cens.0151.lm.2, h.best.cens.0151, h.best.cens.0151.null, hdata.1) h.best.cens.0151.raw <- lagsarlm(log(dens_015_1+0.01)~poly(prop_noncauc.s,2,raw=T)+poly(sa1_medage.s,2,raw=T)+ poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_avghouse.s,2,raw=T), listw=h.sar.cens.0151$weight, type="lag", method="MC", data=hdata.1, na.action="na.fail") plot.partial(h.best.cens.0151, h.best.cens.0151.raw, hdata.1, 2, "sa1_medage", "sa1_medage.s", "Median pop. age", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.cens.0151, h.best.cens.0151.raw, hdata.1, 3, "sa1_medtothinc", "sa1_medtothinc.s", "Median weekly household income ($)", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.cens.0151, h.best.cens.0151.raw, hdata.1, 1, "prop_noncauc", "prop_noncauc.s", "Proportion Non-European Descent", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.cens.0151, h.best.cens.0151.raw, hdata.1, 4, "sa1_avghouse", "sa1_avghouse.s", "Avg. household size", "0.15-1m Veg. Density Partial Residuals") ####################################### ####Urban form variables - vegetation density 0.15-1m #### hdata.1$road_density <- hdata.1$road_length/1000 hdata.1$road_density.s <- rescale(hdata.1$road_density) hdata.1$mb_dweldens.s <- rescale(hdata.1$mb_dwel_dens) hdata.1$park_prop.s <- rescale(hdata.1$park_prop) hdata.1$lot_size.s <- rescale(hdata.1$lot_size) hdata.1$log_lot_size <- log10(hdata.1$lot_size) hdata.1$log_lot_size.s <- rescale(hdata.1$log_lot_size) hdata.1$park_prop_per <- hdata.1$park_prop/100 h.urb.0151.lm <- lm(log(dens_015_1+0.01)~1, data=hdata.1) h.urb.0151.lm.1 <- lm(log(dens_015_1+0.01)~mb_dweldens.s+road_density.s+park_prop.s+lot_size.s, data=hdata.1) h.urb.0151.lm.2 <- lm(log(dens_015_1+0.01)~poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), data=hdata.1) comp.models(h.urb.0151.lm.1, h.urb.0151.lm.2) h.sar.urb.0151 <- mm.spatial.sparse(hdata.1, 150, h.urb.0151.lm.2, "MC") mm.sum.sp(h.sar.urb.0151) spatial.plot(urb.0151.lm.2, data.1, "0.15-1m Veg Density - Urban form polynomial non-spatial model") spatial.plot(sar.cens.0151$error, data.1, "0.15-1m Veg Density - Urban form mixed SAR model") h.best.urb.0151 <- lagsarlm(log(dens_015_1+0.01) ~ poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), listw=h.sar.urb.0151$weight, method="MC", type="mixed", data=hdata.1, na.action="na.fail") AICc(h.best.urb.0151) summary(h.best.urb.0151, Nagelkerke=TRUE) diag.best(best.urb.0151, data.1, sar.urb.0151$weight, "0.15-1m Veg Density - Best Urban Form SAR model") h.best.urb.0151.null <- lagsarlm(log(dens_015_1+0.01) ~ 1, listw=h.sar.urb.0151$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.urb.0151.lm, h.urb.0151.lm.2, h.best.urb.0151, h.best.urb.0151.null, hdata.1) h.best.urb.0151.raw <- lagsarlm(log(dens_015_1+0.01)~poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+ poly(park_prop.s,2,raw=T)+poly(lot_size.s,2,raw=T), listw=h.sar.urb.0151$weight, data=hdata.1, method="MC", type="mixed", na.action="na.fail") plot.partial(h.best.urb.0151, h.best.urb.0151.raw, hdata.1, 2, "road_density", "road_density.s", "Street density (km/km2)", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.urb.0151, h.best.urb.0151.raw, hdata.1, 3, "park_prop_per", "park_prop.s", "Proportion parkland (%)", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.urb.0151, h.best.urb.0151.raw, hdata.1, 1, "mb_dwel_dens", "mb_dweldens.s", "Dwelling density (dwellings/ha)", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.urb.0151, h.best.urb.0151.raw, hdata.1, 4, "log_lot_size", "log_lot_size.s", "Log average lot size (m2)", "0.15-1m Veg. Density Partial Residuals") ####################################### ####Landscape variables - vegetation density 0.15-1m #### hdata.1$clumpy.s <- rescale(hdata.1$clumpy) hdata.1$patch_area_mn.s <- rescale(hdata.1$patch_area_mn) hdata.1$perimeter_area_mn.s <- rescale(hdata.1$perimeter_area_mn) hdata.1$perimeter_area_div <- hdata.1$perimeter_area_mn/10000 hdata.1$tree_area_per <- hdata.1$tree_area*100 hdata.1$tree_area.s <- rescale(hdata.1$tree_area*100) hdata.1$number_patches.s <- rescale(hdata.1$number_patches) h.land.0151.lm <- lm(log(dens_015_1+0.01)~1, data=hdata.1) h.land.0151.lm.1 <- lm(log(dens_015_1+0.01)~clumpy.s+number_patches.s+perimeter_area_mn.s+tree_area.s, data=hdata.1) h.land.0151.lm.2 <- lm(log(dens_015_1+0.01)~poly(clumpy.s,2)+poly(number_patches.s,2)+poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), data=hdata.1) comp.models(h.land.0151.lm.1, h.land.0151.lm.2) h.sar.land.0151 <- mm.spatial.sparse(hdata.1, 150, h.land.0151.lm.2, "MC") mm.sum.sp(h.sar.land.0151) spatial.plot(land.0151.lm.2, data.1, "0.15-1m Veg Density - Landscape polynomial non-spatial model") spatial.plot(sar.land.0151$error, data.1, "0.15-1m Veg Density - Landscape mixed SAR model") h.best.land.0151 <- lagsarlm(log(dens_015_1+0.01) ~ poly(clumpy.s,2)+poly(number_patches.s,2)+ poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), listw=h.sar.land.0151$weight, data=hdata.1, method="MC", type="mixed", na.action="na.fail") AICc(h.best.land.0151) summary(h.best.land.0151, Nagelkerke=TRUE) diag.best(best.land.0151, data.1, sar.land.0151$weight, "0.15-1m Veg Density - Best Landscape SAR model") h.best.land.0151.null <- lagsarlm(log(dens_015_1+0.01) ~ 1, listw=h.sar.land.0151$weight, data=hdata.1, method="MC", type="mixed", na.action="na.fail") spatial.diag(h.land.0151.lm, h.land.0151.lm.2, h.best.land.0151, h.best.land.0151.null, hdata.1) h.best.land.0151.raw <- lagsarlm(log(dens_015_1+0.01)~poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+ poly(perimeter_area_mn.s,2,raw=T)+poly(tree_area.s,2,raw=T), listw=h.sar.land.0151$weight, data=hdata.1, method="MC", type="mixed", na.action="na.fail") plot.partial(h.best.land.0151, h.best.land.0151.raw, hdata.1, 1, "clumpy", "clumpy.s", "Clumpiness", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.land.0151, h.best.land.0151.raw, hdata.1, 2, "number_patches", "number_patches.s", "Number forest fragments", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.land.0151, h.best.land.0151.raw, hdata.1, 3, "perimeter_area_div", "perimeter_area_mn.s", "Forest perimeter:area ratio", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.land.0151, h.best.land.0151.raw, hdata.1, 4, "tree_area", "tree_area.s", "Proportion forest (%)", "0.15-1m Veg. Density Partial Residuals") ####################################### ####Combined model - vegetation density 0.15-1m #### h.comb.0151.lm <- lm(log(dens_015_1+0.01)~1, data=hdata.1) h.comb.0151.lm.1 <- lm(log(dens_015_1+0.01)~aspect_cos.s+aspect_sin.s+elev.s+slope.s+ tot_p.s+water_cap.s+soc.s+tot_n.s+sa1_medtothinc.s+sa1_avghouse.s+ sa1_medage.s+prop_noncauc.s+mb_dweldens.s+road_density.s+lot_size.s+park_prop.s+ clumpy.s+number_patches.s+tree_area.s+perimeter_area_mn.s, data=hdata.1) h.comb.0151.lm.2 <- lm(log(dens_015_1+0.01)~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), data=hdata.1) comp.models(h.comb.0151.lm.1, h.comb.0151.lm.2) h.sar.comb.0151 <- mm.spatial.sparse(hdata.1, 150, h.comb.0151.lm.2, "MC") mm.sum.sp(h.sar.comb.0151) spatial.plot(comb.0151.lm.2, data.1, "0.15-1m Veg Density - Combined polynomial non-spatial model") spatial.plot(sar.comb.0151$mixed, data.1, "0.15-1m Veg Density - Comnbined mixed SAR model") h.best.comb.0151 <- lagsarlm(log(dens_015_1+0.01) ~ poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), type="mixed", listw=h.sar.comb.0151$weight, data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.comb.0151) summary(h.best.comb.0151) diag.best(best.comb.0151, data.1, sar.comb.0151$weight, "0.15-1m Veg Density - Best Combined SAR model") h.best.comb.0151.null <- lagsarlm(log(dens_015_1+0.01) ~ 1, type="mixed", listw=h.sar.comb.0151$weight, data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.comb.0151.lm, h.comb.0151.lm.2, h.best.comb.0151, h.best.comb.0151.null, hdata.1) ##Creation of combined model with all possible variables h.best.comb.0151.raw <- lagsarlm(log(dens_015_1+0.01)~poly(aspect_cos.s,2,raw=T)+poly(aspect_sin.s,2,raw=T)+poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(tot_p.s,2,raw=T)+poly(water_cap.s,2,raw=T)+poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_avghouse.s,2,raw=T)+ poly(sa1_medage.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T)+poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+poly(lot_size.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+poly(tree_area.s,2,raw=T)+poly(perimeter_area_mn.s,2,raw=T), type="mixed", listw=h.sar.comb.0151$weight, data=hdata.1, method="MC", na.action="na.fail") summary(best.comb.0151.raw) ##Creation of all possible 4-variable models for comparison h.best.comb.0151.raw.avg <- dredge(h.best.comb.0151.raw, m.max=4) h.best.comb.0151.raw.avg2 <- model.avg(h.best.comb.0151.raw.avg) summary(h.best.comb.0151.raw.avg2) ##Selection and summary of best combined model h.comb.0151.avg.best <- get.models(h.best.comb.0151.raw.avg,1)[[1]] summary(get.models(h.best.comb.0151.raw.avg,1)[[1]]) ##best model includes tree area, aspect cos, aspect sin, slope ##Creation of best combined standard and raw models and best combined null model h.best.simple.0151 <- lagsarlm(log(dens_015_1+0.01)~poly(tree_area.s,2)+poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+ poly(slope.s,2), listw=h.sar.comb.0151$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.simple.0151.raw <- lagsarlm(log(dens_015_1+0.01)~poly(tree_area.s,2,raw=T)+poly(aspect_cos.s,2,raw=T)+poly(aspect_sin.s,2,raw=T)+ poly(slope.s,2,raw=T), listw=h.sar.comb.0151$weight,data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.null <- lagsarlm(log(dens_015_1+0.01)~1, listw=h.sar.comb.0151$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") AICc(h.best.simple.0151) summary(h.best.simple.0151) summary(h.best.simple.0151.raw) spatial.diag(h.comb.0151.lm, h.comb.0151.lm.2, h.best.simple.0151, h.best.null, hdata.1) ##Calculation of SAR and LM Moran's I values h.comb.0151.best.lm.2 <- lm(log(dens_015_1+0.01)~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(slope.s,2)+ poly(tree_area.s,2), data=hdata.1) h.sar.best.0151 <- mm.spatial.sparse(hdata.1, 150, h.comb.0151.best.lm.2, "MC") mm.sum.sp(h.sar.best.0151) plot.partial(h.best.simple.0151, h.best.simple.0151.raw, hdata.1, 1, "tree_area", "tree_area.s", "Proportion tree cover (%)", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.simple.0151, h.best.simple.0151.raw, hdata.1, 2, "aspect_cos", "aspect_cos.s", "North aspect", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.simple.0151, h.best.simple.0151.raw, hdata.1, 3, "aspect_sin", "aspect_sin.s", "East aspect", "0.15-1m Veg. Density Partial Residuals") plot.partial(h.best.simple.0151, h.best.simple.0151.raw, hdata.1, 4, "slope", "slope.s", "Slope (degrees)", "0.15-1m Veg. Density Partial Residuals") ##################################################################### #####Model selection for vegetation density between 1-2m ##################################################################### ####################################### ####Physical variables - vegetation density 1-2m #### h.phys.12.lm <- lm(log(dens_1_2+0.01)~1, data=hdata.1) h.phys.12.lm.1 <- lm(log(dens_1_2+0.01)~aspect_cos.s+aspect_sin.s+elev.s+slope.s, data=hdata.1) h.phys.12.lm.2 <- lm(log(dens_1_2+0.01)~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2), data=hdata.1) comp.models(h.phys.12.lm.1, h.phys.12.lm.2) h.sar.phys.12 <- mm.spatial.sparse(hdata.1, 150, h.phys.12.lm.2, "MC") mm.sum.sp(h.sar.phys.12) spatial.plot(phys.12.lm.2, data.1, "1-2m Veg Density - Polynomial non-spatial model") spatial.plot(sar.phys.12$mixed, data.1, "1-2m Veg Density - Mixed SAR model") h.best.phys.12 <- lagsarlm(log(dens_1_2+0.01) ~ poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(slope.s,2)+poly(elev.s,2), listw=h.sar.phys.12$weight, data=hdata.1, method="MC", type="mixed", na.action="na.fail") AICc(h.best.phys.12) summary(h.best.phys.12) diag.best(best.phys.12, data.1, sar.phys.12$weight, "1-2m Veg Density - Best SAR model") h.best.phys.12.null <- lagsarlm(log(dens_1_2+0.01) ~ 1, listw=h.sar.phys.12$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") spatial.diag(h.phys.12.lm, h.phys.12.lm.2, h.best.phys.12, h.best.phys.12.null, hdata.1) h.best.phys.12.raw <- lagsarlm(log(dens_1_2+0.01) ~ poly(aspect_cos.s,2,raw=T)+poly(aspect_sin.s,2,raw=T)+ poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T), listw=h.sar.phys.12$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.phys.12, h.best.phys.12.raw, hdata.1, 1, "aspect_cos", "aspect_cos.s", "Northness", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.phys.12, h.best.phys.12.raw, hdata.1, 2, "aspect_sin", "aspect_sin.s", "Eastness", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.phys.12, h.best.phys.12.raw, hdata.1, 3, "elev", "elev.s", "Elevation (m)", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.phys.12, h.best.phys.12.raw, hdata.1, 4, "slope", "slope.s", "Slope (degrees)", "1-2m Veg. Density Partial Residuals") ####################################### ####Soil variables - vegetation density 1-2m #### h.soil.12.lm <- lm(log(dens_1_2+0.01)~1, data=hdata.1) h.soil.12.lm.1 <- lm(log(dens_1_2+0.01)~water_cap.s+tot_p.s+soc.s+tot_n.s, data=hdata.1) h.soil.12.lm.2 <- lm(log(dens_1_2+0.01)~poly(water_cap.s,2)+poly(tot_p.s,2)+poly(soc.s,2)+poly(tot_n.s,2), data=hdata.1) comp.models(h.soil.12.lm.1, h.soil.12.lm.2) h.sar.soil.12 <- mm.spatial.sparse(hdata.1, 150, h.soil.12.lm.2, "MC") mm.sum.sp(h.sar.soil.12) spatial.plot(soil.12.lm.2, data.1, "1-2m Veg Density - Soil polynomial non-spatial model") spatial.plot(sar.soil.12$mixed, data.1, "1-2m Veg Density - Soil mixed SAR model") h.best.soil.12 <- lagsarlm(log(dens_1_2+0.01) ~ poly(soc.s,2)+poly(tot_n.s,2)+poly(water_cap.s,2)+poly(tot_p.s,2), listw=h.sar.soil.12$weight, data=hdata.1, type="lag", method="MC", na.action="na.fail") AICc(h.best.soil.12) summary(h.best.soil.12) diag.best(best.soil.12, data.1, sar.soil.12$weight, "1-2m Veg Density - Best Soil SAR model") h.best.soil.12.null <- lagsarlm(log(dens_1_2+0.01) ~ 1, listw=h.sar.soil.12$weight, type="lag", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.soil.12.lm, h.soil.12.lm.2, h.best.soil.12, h.best.soil.12.null, hdata.1) h.best.soil.12.raw <- lagsarlm(log(dens_1_2+0.01) ~ poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+ poly(water_cap.s,2,raw=T)+poly(tot_p.s,2,raw=T), listw=h.sar.soil.12$weight, data=hdata.1, type="lag", method="MC", na.action="na.fail") plot.partial(h.best.soil.12, h.best.soil.12.raw, hdata.1, 1, "soc", "soc.s", "Soil organic carbon (%)", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.soil.12, h.best.soil.12.raw, hdata.1, 2, "tot_n", "tot_n.s", "Total soil N (%)", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.soil.12, h.best.soil.12.raw, hdata.1, 3, "water_cap", "water_cap.s", "Soil water capacity (%)", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.soil.12, h.best.soil.12.raw, hdata.1, 4, "tot_p", "tot_p.s", "Total soil P (%)", "1-2m Veg. Density Partial Residuals") ####################################### ####Census variables - vegetation density 1-2m #### h.cens.12.lm <- lm(log(dens_1_2+0.01)~1, data=hdata.1) h.cens.12.lm.1 <- lm(log(dens_1_2+0.01)~prop_noncauc.s+sa1_medtothinc.s+sa1_medage.s+sa1_avghouse.s, data=hdata.1) h.cens.12.lm.2 <- lm(log(dens_1_2+0.01)~poly(prop_noncauc.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2), data=hdata.1) comp.models(h.cens.12.lm.1, h.cens.12.lm.2) h.sar.cens.12 <- mm.spatial.sparse(hdata.1, 150, h.cens.12.lm.2, "MC") mm.sum.sp(h.sar.cens.12) spatial.plot(cens.12.lm.2, data.1, "1-2m Veg Density - Census polynomial non-spatial model") spatial.plot(sar.cens.12$mixed, data.1, "0.15-1m Veg Density - Census mixed SAR model") h.best.cens.12 <- lagsarlm(log(dens_1_2+0.01) ~ poly(prop_noncauc.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2), listw=h.sar.cens.12$weight, method="MC", type="lag", data=hdata.1, na.action="na.fail") AICc(h.best.cens.12) summary(h.best.cens.12) diag.best(best.cens.12, data.2, sar.cens.12$weight, "1-2m Veg Density - Best Soil SAR model") h.best.cens.12.null <- lagsarlm(log(dens_1_2+0.01) ~ 1, listw=h.sar.cens.12$weight, type="lag", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.cens.12.lm, h.cens.12.lm.2, h.best.cens.12, h.best.cens.12.null, hdata.1) h.best.cens.12.raw <- lagsarlm(log(dens_1_2+0.01)~poly(prop_noncauc.s,2,raw=T)+poly(sa1_medtothinc.s,2,raw=T)+ poly(sa1_avghouse.s,2,raw=T)+poly(sa1_medage.s,2,raw=T), listw=h.sar.cens.12$weight, type="lag", method="MC", data=hdata.1, na.action="na.fail") plot.partial(h.best.cens.12, h.best.cens.12.raw, hdata.1, 4, "sa1_medage", "sa1_medage.s", "Median pop. age", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.cens.12, h.best.cens.12.raw, hdata.1, 2, "sa1_medtothinc", "sa1_medtothinc.s", "Median weekly household income ($)", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.cens.12, h.best.cens.12.raw, hdata.1, 1, "prop_noncauc", "prop_noncauc.s", "Proportion Non-European Descent", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.cens.12, h.best.cens.12.raw, hdata.1, 3, "sa1_avghouse", "sa1_avghouse.s", "Avg. household size", "1-2m Veg. Density Partial Residuals") ####################################### ####Urban form variables - vegetation density 1-2m #### h.urb.12.lm <- lm(log(dens_1_2+0.01)~1, data=hdata.1) h.urb.12.lm.1 <- lm(log(dens_1_2+0.01)~mb_dweldens.s+road_density.s+park_prop.s+lot_size.s, data=hdata.1) h.urb.12.lm.2 <- lm(log(dens_1_2+0.01)~poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), data=hdata.1) comp.models(h.urb.12.lm.1, h.urb.12.lm.2) h.sar.urb.12 <- mm.spatial.sparse(hdata.1, 150, h.urb.12.lm.2, "MC") mm.sum.sp(h.sar.urb.12) spatial.plot(urb.12.lm.2, data.1, "1-2m Veg Density - Urban form polynomial non-spatial model") spatial.plot(sar.cens.12$error, data.1, "1-2m Veg Density - Urban form mixed SAR model") h.best.urb.12 <- lagsarlm(log(dens_1_2+0.01) ~ poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), listw=h.sar.urb.12$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.urb.12) summary(h.best.urb.12) diag.best(best.urb.12, data.1, sar.urb.12$weight, "1-2m Veg Density - Best Urban Form SAR model") h.best.urb.12.null <- lagsarlm(log(dens_1_2+0.01) ~ 1, listw=h.sar.urb.12$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.urb.12.lm, h.urb.12.lm.2, h.best.urb.12, h.best.urb.12.null, hdata.1) h.best.urb.12.raw <- lagsarlm(log(dens_1_2+0.01)~poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+ poly(park_prop.s,2,raw=T)+poly(lot_size.s,2,raw=T), listw=h.sar.urb.12$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.urb.12, h.best.urb.12.raw, hdata.1, 2, "road_density", "road_density.s", "Street density (km/km2)", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.urb.12, h.best.urb.12.raw, hdata.1, 3, "park_prop_per", "park_prop.s", "Proportion parkland (%)", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.urb.12, h.best.urb.12.raw, hdata.1, 1, "mb_dwel_dens", "mb_dweldens.s", "Dwelling density (dwellings/ha)", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.urb.12, h.best.urb.12.raw, hdata.1, 4, "log_lot_size", "log_lot_size.s", "Mean lot size", "1-2m Veg. Density Partial Residuals") ####################################### ####Landscape variables - vegetation density 1-2m #### h.land.12.lm <- lm(log(dens_1_2+0.01)~1, data=hdata.1) h.land.12.lm.1 <- lm(log(dens_1_2+0.01)~clumpy.s+number_patches.s+perimeter_area_mn.s+tree_area.s, data=hdata.1) h.land.12.lm.2 <- lm(log(dens_1_2+0.01)~poly(clumpy.s,2)+poly(number_patches.s,2)+poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), data=hdata.1) comp.models(h.land.12.lm.1, h.land.12.lm.2) h.sar.land.12 <- mm.spatial.sparse(hdata.1, 150, h.land.12.lm.2, "MC") mm.sum.sp(h.sar.land.12) spatial.plot(land.12.lm.2, data.1, "1-2m Veg Density - Landscape polynomial non-spatial model") spatial.plot(sar.land.12$error, data.1, "1-2m Veg Density - Landscape error SAR model") h.best.land.12 <- lagsarlm(log(dens_1_2+0.01) ~ poly(clumpy.s,2)+poly(number_patches.s,2)+ poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), listw=h.sar.land.12$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") AICc(h.best.land.12) summary(h.best.land.12) diag.best(best.land.12, data.1, sar.land.12$weight, "1-2m Veg Density - Best Landscape SAR model") h.best.land.12.null <- lagsarlm(log(dens_1_2+0.01) ~ 1, listw=h.sar.land.12$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") spatial.diag(h.land.12.lm, h.land.12.lm.2, h.best.land.12, h.best.land.12.null, hdata.1) h.best.land.12.raw <- lagsarlm(log(dens_1_2+0.01)~poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+ poly(perimeter_area_mn.s,2,raw=T)+poly(tree_area.s,2,raw=T), listw=h.sar.land.12$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.land.12, h.best.land.12.raw, hdata.1, 1, "clumpy", "clumpy.s", "Clumpiness", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.land.12, h.best.land.12.raw, hdata.1, 2, "number_patches", "number_patches.s", "Number forest patches", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.land.12, h.best.land.12.raw, hdata.1, 3, "perimeter_area_mn", "perimeter_area_mn.s", "Forest perimeter:area ratio", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.land.12, h.best.land.12.raw, hdata.1, 4, "tree_area", "tree_area.s", "Proportion forest (%)", "1-2m Veg. Density Partial Residuals") ####################################### ####Combined model - vegetation density 1-2m #### h.comb.12.lm <- lm(log(dens_1_2+0.01)~1, data=hdata.1) h.comb.12.lm.1 <- lm(log(dens_1_2+0.01)~aspect_cos.s+aspect_sin.s+elev.s+slope.s+ tot_p.s+water_cap.s+soc.s+tot_n.s+sa1_medtothinc.s+sa1_avghouse.s+ sa1_medage.s+prop_noncauc.s+mb_dweldens.s+road_density.s+lot_size.s+park_prop.s+ clumpy.s+number_patches.s+tree_area.s+perimeter_area_mn.s, data=hdata.1) h.comb.12.lm.2 <- lm(log(dens_1_2+0.01)~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), data=hdata.1) comp.models(h.comb.12.lm.1, h.comb.12.lm.2) h.sar.comb.12 <- mm.spatial.sparse(hdata.1, 150, h.comb.12.lm.2, "MC") mm.sum.sp(h.sar.comb.12) spatial.plot(comb.12.lm.2, data.1, "1-2m Veg Density - Combined polynomial non-spatial model") spatial.plot(sar.comb.12$mixed, data.1, "1-2m Veg Density - Comnbined mixed SAR model") h.best.comb.12 <- lagsarlm(log(dens_1_2+0.01) ~ poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), type="mixed", listw=h.sar.comb.12$weight, data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.comb.12) summary(h.best.comb.12) h.best.comb.12.raw.avg <- dredge(h.best.comb.12.raw, m.max=4, trace=2) h.best.comb.12.raw.avg2 <- model.avg(h.best.comb.12.raw.avg) summary(h.best.comb.12.raw.avg2) h.comb.12.avg.best <- get.models(h.best.comb.12.raw.avg,1)[[1]] summary(get.models(h.best.comb.12.raw.avg,1)[[1]]) ##best model includes slope, clumpy, park prop, tree area h.best.simple.12 <- lagsarlm(log(dens_1_2+0.01)~poly(tree_area.s,2)+poly(clumpy.s,2)+poly(park_prop.s,2)+ poly(slope.s,2), listw=h.sar.comb.12$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.simple.12.raw <- lagsarlm(log(dens_1_2+0.01)~poly(tree_area.s,2,raw=T)+poly(clumpy.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(slope.s,2,raw=T), listw=h.sar.comb.12$weight,data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.null <- lagsarlm(log(dens_1_2+0.01)~1, listw=h.sar.comb.12$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") AICc(h.best.simple.12) summary(h.best.simple.12) summary(h.best.simple.12.raw) spatial.diag(h.comb.12.lm, h.comb.12.lm.2, h.best.simple.12, h.best.null, hdata.1) h.comb.12.best.lm.2 <- lm(log(dens_1_2+0.01)~poly(tree_area.s,2)+poly(clumpy.s,2)+poly(park_prop.s,2)+ poly(slope.s,2), data=hdata.1) h.sar.best.12 <- mm.spatial.sparse(hdata.1, 150, h.comb.12.best.lm.2, "MC") mm.sum.sp(h.sar.best.12) plot.partial(h.best.simple.12, h.best.simple.12.raw, hdata.1, 2, "clumpy", "clumpy.s", "Tree cover clumpiness", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.simple.12, h.best.simple.12.raw, hdata.1, 3, "park_prop", "park_prop.s", "Proportion parkland (%)", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.simple.12, h.best.simple.12.raw, hdata.1, 4, "slope", "slope.s", "Slope (degrees)", "1-2m Veg. Density Partial Residuals") plot.partial(h.best.simple.12, h.best.simple.12.raw, hdata.1, 1, "tree_area", "tree_area.s", "Proportion tree cover (%)", "1-2m Veg. Density Partial Residuals") ##################################################################### #####Model selection for vegetation density between 2-5m ##################################################################### ####################################### ####Physical variables - vegetation density 2-5m #### h.phys.25.lm <- lm(dens_2_5~1, data=hdata.1) h.phys.25.lm.1 <- lm(dens_2_5~aspect_cos.s+aspect_sin.s+elev.s+slope.s, data=hdata.1) h.phys.25.lm.2 <- lm(dens_2_5~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2), data=hdata.1) comp.models(h.phys.25.lm.1, h.phys.25.lm.2) h.sar.phys.25 <- mm.spatial.sparse(hdata.1, 150, h.phys.25.lm.2, "MC") mm.sum.sp(h.sar.phys.25) spatial.plot(phys.25.lm.2, data.1, "2-5m Veg Density - Polynomial non-spatial model") spatial.plot(sar.phys.25$error, data.1, "2-5m Veg Density - Lag SAR model") h.best.phys.25 <- lagsarlm(dens_2_5 ~ poly(elev.s,2)+poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+ poly(slope.s,2), listw=h.sar.phys.25$weight, type="lag", data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.phys.25) summary(h.best.phys.25) diag.best(best.phys.25, data.1, sar.phys.25$weight, "2-5m Veg Density - Best SAR model") h.best.phys.25.null <- lagsarlm(dens_2_5 ~ 1, listw=h.sar.phys.25$weight, data=hdata.1, type="lag", method="MC", na.action="na.fail") spatial.diag(h.phys.25.lm, h.phys.25.lm.2, h.best.phys.25, h.best.phys.25.null, hdata.1) h.best.phys.25.raw <- lagsarlm(dens_2_5 ~ poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(aspect_cos.s,2,raw=T)+poly(aspect_sin,2,raw=T), listw=h.sar.phys.25$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.phys.25, h.best.phys.25.raw, hdata.1, 3, "aspect_cos", "aspect_cos.s", "Northness", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.phys.25, h.best.phys.25.raw, hdata.1, 4, "aspect_sin", "aspect_sin.s", "Eastness", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.phys.25, h.best.phys.25.raw, hdata.1, 1, "elev", "elev.s", "Elevation (m)", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.phys.25, h.best.phys.25.raw, hdata.1, 2, "slope", "slope.s", "Slope (degrees)", "2-5m Veg. Density Partial Residuals") ####################################### ####Soil variables - vegetation density 2-5m #### h.soil.25.lm <- lm(dens_2_5~1, data=hdata.1) h.soil.25.lm.1 <- lm(dens_2_5~water_cap.s+tot_p.s+soc.s+tot_n.s, data=hdata.1) h.soil.25.lm.2 <- lm(dens_2_5~poly(water_cap.s,2)+poly(tot_p.s,2)+poly(soc.s,2)+poly(tot_n.s,2), data=hdata.1) comp.models(h.soil.25.lm.1, h.soil.25.lm.2) h.sar.soil.25 <- mm.spatial.sparse(hdata.1, 150, h.soil.25.lm.2, "MC") mm.sum.sp(h.sar.soil.25) spatial.plot(soil.25.lm.2, data.1, "2-5m Veg Density - Soil polynomial non-spatial model") spatial.plot(sar.soil.25$mixed, data.1, "2-5m Veg Density - Soil mixed SAR model") h.best.soil.25 <- lagsarlm(dens_2_5 ~ poly(soc.s,2)+poly(tot_n.s,2)+poly(water_cap.s,2)+poly(tot_p.s,2), listw=h.sar.soil.25$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.soil.25) summary(h.best.soil.25) diag.best(best.soil.25, data.1, sar.soil.25$weight, "2-5m Veg Density - Best Soil SAR model") h.best.soil.25.null <- lagsarlm(dens_2_5 ~ 1, listw=h.sar.soil.25$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.soil.25.lm, h.soil.25.lm.2, h.best.soil.25, h.best.soil.25.null, hdata.1) h.best.soil.25.raw <- lagsarlm(dens_2_5 ~ poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+ poly(water_cap.s,2,raw=T)+poly(tot_p.s,2,raw=T), listw=h.sar.soil.25$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") plot.partial(h.best.soil.25, h.best.soil.25.raw, hdata.1, 1, "soc", "soc.s", "Soil organic carbon (%)", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.soil.25, h.best.soil.25.raw, hdata.1, 2, "tot_n", "tot_n.s", "Total soil N (%)", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.soil.25, h.best.soil.25.raw, hdata.1, 3, "water_cap", "water_cap.s", "Soil water capacity (%)", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.soil.25, h.best.soil.25.raw, hdata.1, 7, "tot_p", "tot_p.s", "Total soil P (%)", "2-5m Veg. Density Partial Residuals") ####################################### ####Census variables - vegetation density 2-5m #### h.cens.25.lm <- lm(dens_2_5~1, data=hdata.1) h.cens.25.lm.1 <- lm(dens_2_5~prop_noncauc.s+sa1_medtothinc.s+sa1_medage.s+sa1_avghouse, data=hdata.1) h.cens.25.lm.2 <- lm(dens_2_5~poly(prop_noncauc.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2), data=hdata.1) comp.models(h.cens.25.lm.1, h.cens.25.lm.2) h.sar.cens.25 <- mm.spatial.sparse(hdata.1, 150, h.cens.25.lm.2, "MC") mm.sum.sp(h.sar.cens.25) spatial.plot(cens.25.lm.2, data.1, "2-5m Veg Density - Census polynomial non-spatial model") spatial.plot(sar.cens.25$mixed, data.1, "2-5m Veg Density - Census mixed SAR model") h.best.cens.25 <- lagsarlm(dens_2_5 ~ poly(prop_noncauc.s,2)+poly(sa1_medtothinc.s,2)+ poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2), listw=h.sar.cens.25$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") AICc(h.best.cens.25) summary(h.best.cens.25) diag.best(best.cens.25, data.1, sar.cens.25$weight, "2-5m Veg Density - Best Soil SAR model") h.best.cens.25.null <- lagsarlm(dens_2_5 ~ 1, listw=h.sar.cens.25$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") spatial.diag(h.cens.25.lm, h.cens.25.lm.2, h.best.cens.25, h.best.cens.25.null, hdata.1) h.best.cens.25.raw <- lagsarlm(dens_2_5~poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_medage.s,2,raw=T)+ poly(sa1_avghouse.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T), listw=h.sar.cens.25$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.cens.25, h.best.cens.25.raw, hdata.1, 2, "sa1_medage", "sa1_medage.s", "Median pop. age", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.cens.25, h.best.cens.25.raw, hdata.1, 1, "sa1_medtothinc", "sa1_medtothinc.s", "Median weekly household income ($)", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.cens.25, h.best.cens.25.raw, hdata.1, 4, "prop_noncauc", "prop_noncauc.s", "Proportion Non-European Descent", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.cens.25, h.best.cens.25.raw, hdata.1, 3, "sa1_avghouse", "sa1_avghouse.s", "Avg. household size", "2-5m Veg. Density Partial Residuals") ####################################### ####Urban form variables - vegetation density 2-5m #### h.urb.25.lm <- lm(dens_2_5~1, data=hdata.1) h.urb.25.lm.1 <- lm(dens_2_5~mb_dweldens.s+road_density.s+park_prop.s+lot_size.s, data=hdata.1) h.urb.25.lm.2 <- lm(dens_2_5~poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), data=hdata.1) comp.models(h.urb.25.lm.1, h.urb.25.lm.2) h.sar.urb.25 <- mm.spatial.sparse(hdata.1, 150, h.urb.25.lm.2, "MC") mm.sum.sp(h.sar.urb.25) spatial.plot(urb.25.lm.2, data.1, "2-5m Veg Density - Urban form polynomial non-spatial model") spatial.plot(sar.cens.25$error, data.1, "2-5m Veg Density - Urban form mixed SAR model") h.best.urb.25 <- lagsarlm(dens_2_5 ~ poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), listw=h.sar.urb.25$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.urb.25) summary(h.best.urb.25) diag.best(best.urb.25, data.1, sar.urb.25$weight, "2-5m Veg Density - Best Urban Form SAR model") h.best.urb.25.null <- lagsarlm(dens_2_5 ~ 1, listw=h.sar.urb.25$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.urb.25.lm, h.urb.25.lm.2, h.best.urb.25, h.best.urb.25.null, hdata.1) h.best.urb.25.raw <- lagsarlm(dens_2_5~poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+ poly(park_prop.s,2,raw=T)+poly(lot_size.s,2,raw=T), listw=h.sar.urb.25$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.urb.25, h.best.urb.25.raw, hdata.1, 2, "road_density", "road_density.s", "Street density (km/km2)", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.urb.25, h.best.urb.25.raw, hdata.1, 3, "park_prop_per", "park_prop.s", "Proportion parkland (%)", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.urb.25, h.best.urb.25.raw, hdata.1, 1, "mb_dwel_dens", "mb_dweldens.s", "Dwelling density (dwellings/ha)", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.urb.25, h.best.urb.25.raw, hdata.1, 4, "log_lot_size", "log_lot_size.s", "Mean lot size (m2)", "2-5m Veg. Density Partial Residuals") ####################################### ####Landscape variables - vegetation density 2-5mm #### h.land.25.lm <- lm(dens_2_5~1, data=hdata.1) h.land.25.lm.1 <- lm(dens_2_5~clumpy.s+number_patches.s+perimeter_area_mn.s+tree_area.s, data=hdata.1) h.land.25.lm.2 <- lm(dens_2_5~poly(clumpy.s,2)+poly(number_patches.s,2)+poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), data=hdata.1) comp.models(h.land.25.lm.1, h.land.25.lm.2) h.sar.land.25 <- mm.spatial.sparse(hdata.1, 150, h.land.25.lm.2, "MC") mm.sum.sp(h.sar.land.25) spatial.plot(land.25.lm.2, data.1, "2-5m Veg Density - Landscape polynomial non-spatial model") spatial.plot(sar.land.25$lag, data.1, "2-5m Veg Density - Landscape error SAR model") h.best.land.25 <- lagsarlm(dens_2_5 ~ poly(number_patches.s,2)+poly(clumpy.s,2)+ poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), listw=h.sar.land.25$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.land.25) summary(h.best.land.25) diag.best(best.land.25, data.1, sar.land.25$weight, "2-5m Veg Density - Best Landscape SAR model") h.best.land.25.null <- lagsarlm(dens_2_5 ~ 1, listw=h.sar.land.25$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.land.25.lm, h.land.25.lm.2, h.best.land.25, h.best.land.25.null, hdata.1) h.best.land.25.raw <- lagsarlm(dens_2_5~poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+ poly(perimeter_area_mn.s,2,raw=T)+poly(tree_area.s,2,raw=T), listw=h.sar.land.25$weight, type="lag", data=hdata.1, method="MC", na.action="na.fail") plot.partial(h.best.land.25, h.best.land.25.raw, hdata.1, 4, "clumpy", "clumpy.s", "Clumpiness", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.land.25, h.best.land.25.raw, hdata.1, 1, "number_patches", "number_patches.s", "Number forest fragments", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.land.25, h.best.land.25.raw, hdata.1, 2, "perimeter_area_mn", "perimeter_area_mn.s", "Forest perimeter:area ratio", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.land.25, h.best.land.25.raw, hdata.1, 3, "tree_area", "tree_area.s", "Proportion forest (%)", "2-5m Veg. Density Partial Residuals") ####################################### ####Combined model - vegetation density 2-5m #### h.comb.25.lm <- lm(dens_2_5~1, data=hdata.1) h.comb.25.lm.1 <- lm(dens_2_5~aspect_cos.s+aspect_sin.s+elev.s+slope.s+ tot_p.s+water_cap.s+soc.s+tot_n.s+sa1_medtothinc.s+sa1_avghouse.s+ sa1_medage.s+prop_noncauc.s+mb_dweldens.s+road_density.s+lot_size.s+park_prop.s+ clumpy.s+number_patches.s+tree_area.s+perimeter_area_mn.s, data=hdata.1) h.comb.25.lm.2 <- lm(dens_2_5~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), data=hdata.1) comp.models(h.comb.25.lm.1, h.comb.25.lm.2) h.sar.comb.25 <- mm.spatial.sparse(hdata.1, 150, h.comb.25.lm.2, "MC") mm.sum.sp(h.sar.comb.25) spatial.plot(comb.25.lm.2, data.1, "2-5m Veg Density - Combined polynomial non-spatial model") spatial.plot(sar.comb.25$mixed, data.1, "2-5m Veg Density - Comnbined mixed SAR model") h.best.comb.25 <- lagsarlm(dens_2_5 ~ poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), type="mixed", listw=h.sar.comb.25$weight, data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.comb.25) summary(h.best.comb.25) diag.best(best.comb.25, data.1, sar.comb.25$weight, "2-5m Veg Density - Best Combined SAR model") h.best.comb.25.null <- lagsarlm(dens_2_5 ~ 1, type="mixed", listw=h.sar.comb.25$weight, data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.comb.25.lm, h.comb.25.lm.2, h.best.comb.25, h.best.comb.25.null, hdata.1) h.best.comb.25.raw <- lagsarlm(dens_2_5~poly(aspect_cos.s,2,raw=T)+poly(aspect_sin.s,2,raw=T)+poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(tot_p.s,2,raw=T)+poly(water_cap.s,2,raw=T)+poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_avghouse.s,2,raw=T)+ poly(sa1_medage.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T)+poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+poly(lot_size.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+poly(tree_area.s,2,raw=T)+poly(perimeter_area_mn.s,2,raw=T), type="mixed", listw=h.sar.comb.25$weight, data=hdata.1, method="MC", na.action="na.fail") summary(best.comb.25.raw) h.best.comb.25.raw.avg <- dredge(h.best.comb.25.raw, m.max=4, trace=2) h.best.comb.25.raw.avg2 <- model.avg(h.best.comb.25.raw.avg) summary(h.best.comb.25.raw.avg2) h.comb.25.avg.best <- get.models(h.best.comb.25.raw.avg,1)[[1]] summary(get.models(h.best.comb.25.raw.avg,1)[[1]]) ##includes dwelling density, tree area, park prop, number patches h.best.simple.25 <- lagsarlm(dens_2_5~poly(tree_area.s,2)+poly(mb_dweldens.s,2)+poly(park_prop.s,2)+ poly(number_patches.s,2), listw=h.sar.comb.25$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.simple.25.raw <- lagsarlm(dens_2_5~poly(tree_area.s,2,raw=T)+poly(mb_dweldens.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(number_patches.s,2,raw=T), listw=h.sar.comb.25$weight,data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.null <- lagsarlm(dens_2_5~1, listw=h.sar.comb.25$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") AICc(h.best.simple.25) summary(h.best.simple.25) summary(h.best.simple.25.raw) spatial.diag(h.comb.25.lm, h.comb.25.lm.2, h.best.simple.25, h.best.null, hdata.1) h.comb.25.best.lm.2 <- lm(dens_2_5~poly(tree_area.s,2)+poly(mb_dweldens.s,2)+poly(park_prop.s,2)+ poly(number_patches.s,2), data=hdata.1) h.sar.best.25 <- mm.spatial.sparse(hdata.1, 150, h.comb.25.best.lm.2, "MC") mm.sum.sp(h.sar.best.25) plot.partial(h.best.simple.25, h.best.simple.25.raw, hdata.1, 2, "mb_dwel_dens", "mb_dweldens.s", "Dwelling density (dwellings/ha)", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.simple.25, h.best.simple.25.raw, hdata.1, 3, "park_prop", "park_prop.s", "Proportion parkland (%)", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.simple.25, h.best.simple.25.raw, hdata.1, 4, "number_patches", "number_patches.s", "Number of tree cover patches", "2-5m Veg. Density Partial Residuals") plot.partial(h.best.simple.25, h.best.simple.25.raw, hdata.1, 1, "tree_area", "tree_area.s", "Proportion tree cover (%)", "2-5m Veg. Density Partial Residuals") ##################################################################### #####Model selection for vegetation density between 5-10m ##################################################################### ####################################### ####Physical variables - vegetation density 5-10m #### h.phys.510.lm <- lm(dens_5_10~1, data=hdata.1) h.phys.510.lm.1 <- lm(dens_5_10~aspect_cos.s+aspect_sin.s+elev.s+slope.s, data=hdata.1) h.phys.510.lm.2 <- lm(dens_5_10~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2), data=hdata.1) comp.models(h.phys.510.lm.1, h.phys.510.lm.2) h.sar.phys.510 <- mm.spatial.sparse(hdata.1, 150, h.phys.510.lm.2, "MC") mm.sum.sp(h.sar.phys.510) spatial.plot(phys.510.lm.2, data.1, "5-10m Veg Density - Polynomial non-spatial model") spatial.plot(sar.phys.510$error, data.1, "5-10m Veg Density - Lag SAR model") h.best.phys.510 <- lagsarlm(dens_5_10 ~ poly(elev.s,2)+poly(aspect_cos.s,2)+poly(slope.s,2)+poly(aspect_sin.s,2), listw=h.sar.phys.510$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") AICc(h.best.phys.510) summary(h.best.phys.510) diag.best(best.phys.510, data.1, sar.phys.510$weight, "5-10m Veg Density - Best SAR model") h.best.phys.510.null <- lagsarlm(dens_5_10 ~ 1, listw=h.sar.phys.510$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") spatial.diag(h.phys.510.lm, h.phys.510.lm.2, h.best.phys.510, h.best.phys.510.null, hdata.1) h.best.phys.510.raw <- lagsarlm(dens_5_10 ~ poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(aspect_cos.s,2,raw=T)+poly(aspect_sin,2,raw=T), listw=h.sar.phys.510$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") ##Partial residual plotting to look at trends plot.partial(h.best.phys.510, h.best.phys.510.raw, hdata.1, 3, "aspect_cos", "aspect_cos.s", "Northness", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.phys.510, h.best.phys.510.raw, hdata.1, 4, "aspect_sin", "aspect_sin.s", "Eastness", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.phys.510, h.best.phys.510.raw, hdata.1, 1, "elev", "elev.s", "Elevation (m)", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.phys.510, h.best.phys.510.raw, hdata.1, 2, "slope", "slope.s", "Slope (degrees)", "5-10m Veg. Density Partial Residuals") ####################################### ####Soil variables - vegetation density 5-10m #### h.soil.510.lm <- lm(dens_5_10~1, data=hdata.1) h.soil.510.lm.1 <- lm(dens_5_10~water_cap.s+tot_p.s+soc.s+tot_n.s, data=hdata.1) h.soil.510.lm.2 <- lm(dens_5_10~poly(water_cap.s,2)+poly(tot_p.s,2)+poly(soc.s,2)+poly(tot_n.s,2), data=hdata.1) comp.models(h.soil.510.lm.1, h.soil.510.lm.2) h.sar.soil.510 <- mm.spatial.sparse(hdata.1, 150, h.soil.510.lm.2, "MC") mm.sum.sp(h.sar.soil.510) spatial.plot(soil.510.lm.2, data.1, "5-10m Veg Density - Soil polynomial non-spatial model") spatial.plot(sar.soil.510$mixed, data.1, "5-10m Veg Density - Soil mixed SAR model") h.best.soil.510 <- lagsarlm(dens_5_10 ~ poly(soc.s,2)+poly(tot_n.s,2)+poly(water_cap.s,2)+poly(tot_p.s,2), listw=h.sar.soil.510$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.soil.510) summary(h.best.soil.510) diag.best(best.soil.510, data.1, sar.soil.510$weight, "5-10m Veg Density - Best Soil SAR model") h.best.soil.510.null <- lagsarlm(dens_5_10 ~ 1, listw=h.sar.soil.510$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.soil.510.lm, h.soil.510.lm.2, h.best.soil.510, h.best.soil.510.null, hdata.1) h.best.soil.510.raw <- lagsarlm(dens_5_10 ~ poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+ poly(water_cap.s,2,raw=T)+poly(tot_p.s,2,raw=T), listw=h.sar.soil.510$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") plot.partial(h.best.soil.510, h.best.soil.510.raw, hdata.1, 1, "soc", "soc.s", "Soil organic carbon (%)", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.soil.510, h.best.soil.510.raw, hdata.1, 2, "tot_n", "tot_n.s", "Total soil N (%)", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.soil.510, h.best.soil.510.raw, hdata.1, 3, "water_cap", "water_cap.s", "Soil water capacity (%)", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.soil.510, h.best.soil.510.raw, hdata.1, 4, "tot_p", "tot_p.s", "Total soil P (%)", "5-10m Veg. Density Partial Residuals") ####################################### ####Census variables - vegetation density 5-10m #### h.cens.510.lm <- lm(dens_5_10~1, data=hdata.1) h.cens.510.lm.1 <- lm(dens_5_10~prop_noncauc.s+sa1_medtothinc.s+sa1_medage.s+sa1_avghouse.s, data=hdata.1) h.cens.510.lm.2 <- lm(dens_5_10~poly(prop_noncauc.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2), data=hdata.1) comp.models(h.cens.510.lm.1, h.cens.510.lm.2) h.sar.cens.510 <- mm.spatial.sparse(hdata.1, 150, h.cens.510.lm.2, "MC") mm.sum.sp(h.sar.cens.510) spatial.plot(cens.510.lm.2, data.1, "5-10m Veg Density - Census polynomial non-spatial model") spatial.plot(sar.cens.510$lag, data.1, "5-10m Veg Density - Census mixed SAR model") h.best.cens.510 <- lagsarlm(dens_5_10 ~ poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2)+poly(prop_noncauc.s,2), listw=h.sar.cens.510$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.cens.510) summary(h.best.cens.510) diag.best(best.cens.510, data.1, sar.cens.510$weight, "5-10m Veg Density - Best Soil SAR model") h.best.cens.510.null <- lagsarlm(dens_5_10 ~ 1, listw=h.sar.cens.510$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.cens.510.lm, h.cens.510.lm.2, h.best.cens.510, h.best.cens.510.null, hdata.1) h.best.cens.510.raw <- lagsarlm(dens_5_10~poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_medage.s,2,raw=T)+ poly(sa1_avghouse.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T), listw=h.sar.cens.510$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") plot.partial(h.best.cens.510, h.best.cens.510.raw, hdata.1, 2, "sa1_medage", "sa1_medage.s", "Median pop. age", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.cens.510, h.best.cens.510.raw, hdata.1, 1, "sa1_medtothinc", "sa1_medtothinc.s", "Median weekly household income ($)", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.cens.510, h.best.cens.510.raw, hdata.1, 4, "prop_noncauc", "prop_noncauc.s", "Proportion Non-European Descent", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.cens.510, h.best.cens.510.raw, hdata.1, 3, "sa1_avghouse", "sa1_avghouse.s", "Avg. household size", "5-10m Veg. Density Partial Residuals") ####################################### ####Urban form variables - vegetation density 5-10m #### h.urb.510.lm <- lm(dens_5_10~1, data=hdata.1) h.urb.510.lm.1 <- lm(dens_5_10~mb_dweldens.s+road_density.s+park_prop.s+lot_size.s, data=hdata.1) h.urb.510.lm.2 <- lm(dens_5_10~poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), data=hdata.1) comp.models(h.urb.510.lm.1, h.urb.510.lm.2) h.sar.urb.510 <- mm.spatial.sparse(hdata.1, 150, h.urb.510.lm.2, "MC") mm.sum.sp(h.sar.urb.510) spatial.plot(urb.510.lm.2, data.1, "5-10m Veg Density - Urban form polynomial non-spatial model") spatial.plot(sar.cens.510$mixed, data.1, "5-10m Veg Density - Urban form mixed SAR model") h.best.urb.510 <- lagsarlm(dens_5_10 ~ poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), listw=h.sar.urb.510$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.urb.510) summary(h.best.urb.510) diag.best(best.urb.510, data.1, sar.urb.510$weight, "5-10m Veg Density - Best Urban Form SAR model") h.best.urb.510.null <- lagsarlm(dens_5_10 ~ 1, listw=h.sar.urb.510$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.urb.510.lm, h.urb.510.lm.2, h.best.urb.510, h.best.urb.510.null, hdata.1) h.best.urb.510.raw <- lagsarlm(dens_5_10~poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+ poly(park_prop.s,2,raw=T)+poly(lot_size.s,2,raw=T), listw=h.sar.urb.510$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.urb.510, h.best.urb.510.raw, hdata.1, 2, "road_density", "road_density.s", "Street density (km/km2)", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.urb.510, h.best.urb.510.raw, hdata.1, 3, "park_prop_per", "park_prop.s", "Proportion parkland (%)", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.urb.510, h.best.urb.510.raw, hdata.1, 1, "mb_dwel_dens", "mb_dweldens.s", "Dwelling density (dwellings/ha)", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.urb.510, h.best.urb.510.raw, hdata.1, 4, "log_lot_size", "log_lot_size.s", "Mean lot size (m2)", "5-10m Veg. Density Partial Residuals") ####################################### ####Landscape variables - vegetation density 5-10m #### h.land.510.lm <- lm(dens_5_10~1, data=hdata.1) h.land.510.lm.1 <- lm(dens_5_10~clumpy.s+number_patches.s+perimeter_area_mn.s+tree_area.s, data=hdata.1) h.land.510.lm.2 <- lm(dens_5_10~poly(clumpy.s,2)+poly(number_patches.s,2)+poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), data=hdata.1) comp.models(h.land.510.lm.1, h.land.510.lm.2) h.sar.land.510 <- mm.spatial.sparse(hdata.1, 150, h.land.510.lm.2, "MC") mm.sum.sp(h.sar.land.510) spatial.plot(land.510.lm.2, data.1, "5-10m Veg Density - Landscape polynomial non-spatial model") spatial.plot(sar.land.510$mixed, data.1, "5-10m Veg Density - Landscape error SAR model") h.best.land.510 <- lagsarlm(dens_5_10 ~ poly(clumpy.s,2)+poly(number_patches.s,2)+ poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), listw=h.sar.land.510$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.land.510) summary(h.best.land.510) diag.best(best.land.510, data.1, sar.land.510$weight, "5-10m Veg Density - Best Landscape SAR model") h.best.land.510.null <- lagsarlm(dens_5_10 ~ 1, listw=h.sar.land.510$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.land.510.lm, h.land.510.lm.2, h.best.land.510, h.best.land.510.null, hdata.1) h.best.land.510.raw <- lagsarlm(dens_5_10~poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+ poly(perimeter_area_mn.s,2,raw=T)+poly(tree_area.s,2,raw=T), listw=h.sar.land.510$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") plot.partial(h.best.land.510, h.best.land.510.raw, hdata.1, 1, "clumpy", "clumpy.s", "Clumpiness", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.land.510, h.best.land.510.raw, hdata.1, 2, "number_patches", "number_patches.s", "Avg. forest patch size (ha)", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.land.510, h.best.land.510.raw, hdata.1, 3, "perimeter_area_mn", "perimeter_area_mn.s", "Forest perimeter:area ratio", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.land.510, h.best.land.510.raw, hdata.1, 4, "tree_area", "tree_area.s", "Proportion forest (%)", "5-10m Veg. Density Partial Residuals") ####################################### ####Combined model - vegetation density 5-10m #### h.comb.510.lm <- lm(dens_5_10~1, data=hdata.1) h.comb.510.lm.1 <- lm(dens_5_10~aspect_cos.s+aspect_sin.s+elev.s+slope.s+ tot_p.s+water_cap.s+soc.s+tot_n.s+sa1_medtothinc.s+sa1_avghouse.s+ sa1_medage.s+prop_noncauc.s+mb_dweldens.s+road_density.s+lot_size.s+park_prop.s+ clumpy.s+number_patches.s+tree_area.s+perimeter_area_mn.s, data=hdata.1) h.comb.510.lm.2 <- lm(dens_5_10~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), data=hdata.1) comp.models(h.comb.510.lm.1, h.comb.510.lm.2) h.sar.comb.510 <- mm.spatial.sparse(hdata.1, 150, h.comb.510.lm.2, "MC") mm.sum.sp(h.sar.comb.510) spatial.plot(comb.510.lm.2, data.1, "5-10m Veg Density - Combined polynomial non-spatial model") spatial.plot(sar.comb.510$mixed, data.1, "5-10m Veg Density - Comnbined mixed SAR model") h.best.comb.510 <- lagsarlm(dens_5_10 ~ poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), type="mixed", listw=h.sar.comb.510$weight, data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.comb.510) summary(h.best.comb.510) diag.best(best.comb.510, data.1, sar.comb.510$weight, "5-10m Veg Density - Best Combined SAR model") h.best.comb.510.null <- lagsarlm(dens_5_10 ~ 1, type="mixed", listw=h.sar.comb.510$weight, data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.comb.510.lm, h.comb.510.lm.2, h.best.comb.510, h.best.comb.510.null, hdata.1) h.best.comb.510.raw <- lagsarlm(dens_5_10~poly(aspect_cos.s,2,raw=T)+poly(aspect_sin.s,2,raw=T)+poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(tot_p.s,2,raw=T)+poly(water_cap.s,2,raw=T)+poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_avghouse.s,2,raw=T)+ poly(sa1_medage.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T)+poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+poly(lot_size.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+poly(tree_area.s,2,raw=T)+poly(perimeter_area_mn.s,2,raw=T), type="mixed", listw=h.sar.comb.510$weight, data=hdata.1, method="MC", na.action="na.fail") h.best.comb.510.raw.avg <- dredge(h.best.comb.510.raw, m.max=4, trace=2) h.best.comb.510.raw.avg2 <- model.avg(h.best.comb.510.raw.avg) summary(h.best.comb.510.raw.avg2) h.comb.510.avg.best <- get.models(h.best.comb.510.raw.avg,1)[[1]] summary(get.models(h.best.comb.510.raw.avg,1)[[1]]) ##includes clumpy, dwel dens, avg house, tree area h.best.simple.510 <- lagsarlm(dens_5_10~poly(tree_area.s,2)+poly(sa1_avghouse.s,2)+poly(mb_dweldens.s,2)+ poly(clumpy.s,2), listw=h.sar.comb.510$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.simple.510.raw <- lagsarlm(dens_5_10~poly(tree_area.s,2,raw=T)+poly(sa1_avghouse.s,2,raw=T)+poly(mb_dweldens.s,2,raw=T)+ poly(clumpy.s,2,raw=T), listw=h.sar.comb.510$weight,data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.510.null <- lagsarlm(dens_5_10~1, listw=h.sar.comb.510$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") AICc(h.best.simple.510) summary(h.best.simple.510) summary(h.best.simple.510.raw) spatial.diag(h.comb.510.lm, h.comb.510.lm.2, h.best.simple.510, h.best.510.null, hdata.1) h.comb.510.best.lm.2 <- lm(dens_5_10~poly(tree_area.s,2)+poly(sa1_avghouse.s,2)+poly(mb_dweldens.s,2)+ poly(clumpy.s,2), data=hdata.1) h.sar.best.510 <- mm.spatial.sparse(hdata.1, 150, h.comb.510.best.lm.2, "MC") mm.sum.sp(h.sar.best.510) plot.partial(h.best.simple.510, h.best.simple.510.raw, hdata.1, 2, "sa1_avghouse", "sa1_avghouse.s", "Household size (number persons)", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.simple.510, h.best.simple.510.raw, hdata.1, 3, "mb_dwel_dens", "mb_dweldens.s", "Dwelling density (dwellings/ha)", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.simple.510, h.best.simple.510.raw, hdata.1, 4, "clumpy", "clumpy.s", "Tree cover clumpiness", "5-10m Veg. Density Partial Residuals") plot.partial(h.best.simple.510, h.best.simple.510.raw, hdata.1, 1, "tree_area", "tree_area.s", "Proportion tree cover (%)", "5-10m Veg. Density Partial Residuals") ##################################################################### #####Model selection for vegetation density >10m ##################################################################### ####################################### ####Physical variables - vegetation density >10m #### h.phys.10.lm <- lm(log(dens_above10+0.01)~1, data=hdata.1) h.phys.10.lm.1 <- lm(log(dens_above10+0.01)~aspect_cos.s+aspect_sin.s+elev.s+slope.s, data=hdata.1) h.phys.10.lm.2 <- lm(log(dens_above10+0.01)~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2), data=hdata.1) comp.models(h.phys.10.lm.1, h.phys.10.lm.2) h.sar.phys.10 <- mm.spatial.sparse(hdata.1, 150, h.phys.10.lm.2, "MC") mm.sum.sp(h.sar.phys.10) spatial.plot(phys.10.lm.2, data.1, ">10m Veg Density - Polynomial non-spatial model") spatial.plot(sar.phys.10$mixed, data.1, ">10m Veg Density - Mixed SAR model") h.best.phys.10 <- lagsarlm(log(dens_above10+0.01) ~ poly(elev.s,2)+poly(aspect_sin.s,2)+poly(slope.s,2)+poly(aspect_cos.s,2), listw=h.sar.phys.10$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.phys.10) summary(h.best.phys.10) diag.best(best.phys.10, data.1, sar.phys.10$weight, ">10m Veg Density - Best SAR model") h.best.phys.10.null <- lagsarlm(log(dens_above10+0.01) ~ 1, listw=h.sar.phys.10$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.phys.10.lm, h.phys.10.lm.2, h.best.phys.10, h.best.phys.10.null, hdata.1) h.best.phys.10.raw <- lagsarlm(log(dens_above10+0.01) ~ poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(aspect_sin,2,raw=T)+poly(aspect_cos.s,2,raw=T), listw=h.sar.phys.10$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") plot.partial(h.best.phys.10, h.best.phys.10.raw, hdata.1, 4, "aspect_cos", "aspect_cos.s", "Northness", ">10m Veg. Density Partial Residuals") plot.partial(h.best.phys.10, h.best.phys.10.raw, hdata.1, 3, "aspect_sin", "aspect_sin.s", "Eastness", ">10m Veg. Density Partial Residuals") plot.partial(h.best.phys.10, h.best.phys.10.raw, hdata.1, 1, "elev", "elev.s", "Elevation (m)", ">10m Veg. Density Partial Residuals") plot.partial(h.best.phys.10, h.best.phys.10.raw, hdata.1, 2, "slope", "slope.s", "Slope (degrees)", ">10m Veg. Density Partial Residuals") ####################################### ####Soil variables - vegetation density >10m #### h.soil.10.lm <- lm(log(dens_above10+0.01)~1, data=hdata.1) h.soil.10.lm.1 <- lm(log(dens_above10+0.01)~water_cap.s+tot_p.s+soc.s+tot_n.s, data=hdata.1) h.soil.10.lm.2 <- lm(log(dens_above10+0.01)~poly(water_cap.s,2)+poly(tot_p.s,2)+poly(soc.s,2)+poly(tot_n.s,2), data=hdata.1) comp.models(h.soil.10.lm.1, h.soil.10.lm.2) h.sar.soil.10 <- mm.spatial.sparse(hdata.1, 150, h.soil.10.lm.2, "MC") mm.sum.sp(h.sar.soil.10) spatial.plot(soil.10.lm.2, data.1, ">10m Veg Density - Soil polynomial non-spatial model") spatial.plot(sar.soil.10$error, data.1, ">10m Veg Density - Soil mixed SAR model") h.best.soil.10 <- lagsarlm(log(dens_above10+0.01) ~ poly(soc.s,2)+poly(tot_n.s,2)+poly(water_cap.s,2)+poly(tot_p.s,2), listw=h.sar.soil.10$weight, data=hdata.1, type="lag", method="MC", na.action="na.fail") AICc(h.best.soil.10) summary(h.best.soil.10) diag.best(best.soil.10, data.1, sar.soil.10$weight, ">10m Veg Density - Best Soil SAR model") h.best.soil.10.null <- lagsarlm(log(dens_above10+0.01) ~ 1, listw=h.sar.soil.10$weight, data=hdata.1, type="lag", method="MC", na.action="na.fail") spatial.diag(h.soil.10.lm, h.soil.10.lm.2, h.best.soil.10, h.best.soil.10.null, hdata.1) h.best.soil.10.raw <- lagsarlm(log(dens_above10+0.01) ~ poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+ poly(water_cap.s,2,raw=T)+poly(tot_p.s,2,raw=T), listw=h.sar.soil.510$weight, data=hdata.1, type="lag", method="MC", na.action="na.fail") plot.partial(h.best.soil.10, h.best.soil.10.raw, hdata.1, 1, "soc", "soc.s", "Soil organic carbon (%)", ">10m Veg. Density Partial Residuals") plot.partial(h.best.soil.10, h.best.soil.10.raw, hdata.1, 2, "tot_n", "tot_n.s", "Total soil N (%)", ">10m Veg. Density Partial Residuals") plot.partial(h.best.soil.10, h.best.soil.10.raw, hdata.1, 3, "water_cap", "water_cap.s", "Soil water capacity (%)", ">10m Veg. Density Partial Residuals") plot.partial(h.best.soil.10, h.best.soil.10.raw, hdata.1, 4, "tot_p", "tot_p.s", "Total soil P (%)", ">10m Veg. Density Partial Residuals") ####################################### ####Census variables - vegetation density >10m #### h.cens.10.lm <- lm(log(dens_above10+0.01)~1, data=hdata.1) h.cens.10.lm.1 <- lm(log(dens_above10+0.01)~prop_noncauc.s+sa1_medtothinc.s+sa1_medage.s+sa1_avghouse.s, data=hdata.1) h.cens.10.lm.2 <- lm(log(dens_above10+0.01)~poly(prop_noncauc.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2), data=hdata.1) comp.models(h.cens.10.lm.1, h.cens.10.lm.2) h.sar.cens.10 <- mm.spatial.sparse(hdata.1, 150, h.cens.10.lm.2, "MC") mm.sum.sp(h.sar.cens.10) spatial.plot(cens.10.lm.2, data.1, ">10m Veg Density - Census polynomial non-spatial model") spatial.plot(sar.cens.10$mixed, data.1, ">10m Veg Density - Census error SAR model") h.best.cens.10 <- lagsarlm(log(dens_above10+0.01) ~ poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(sa1_avghouse.s,2), listw=h.sar.cens.10$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.cens.10) summary(h.best.cens.10) diag.best(best.cens.10, data.1, sar.cens.10$weight, ">10m Veg Density - Best Soil SAR model") h.best.cens.10.null <- lagsarlm(log(dens_above10+0.01) ~ 1, listw=h.sar.cens.10$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.cens.10.lm, h.cens.10.lm.2, h.best.cens.10, h.best.cens.10.null, hdata.1) h.best.cens.10.raw <- lagsarlm(log(dens_above10+0.01)~poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_medage.s,2,raw=T)+ poly(sa1_avghouse.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T), listw=h.sar.cens.10$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") plot.partial(h.best.cens.10, h.best.cens.10.raw, hdata.1, 2, "sa1_medage", "sa1_medage.s", "Median pop. age", ">10m Veg. Density Partial Residuals") plot.partial(h.best.cens.10, h.best.cens.10.raw, hdata.1, 1, "sa1_medtothinc", "sa1_medtothinc.s", "Median weekly household income ($)", ">10m Veg. Density Partial Residuals") plot.partial(h.best.cens.10, h.best.cens.10.raw, hdata.1, 4, "prop_noncauc", "prop_noncauc.s", "Proportion Non-European Descent", ">10m Veg. Density Partial Residuals") plot.partial(h.best.cens.10, h.best.cens.10.raw, hdata.1, 3, "sa1_avghouse", "sa1_avghouse.s", "Avg. household size", ">10m Veg. Density Partial Residuals") ####################################### ####Urban form variables - vegetation density >10m #### h.urb.10.lm <- lm(log(dens_above10+0.01)~1, data=hdata.1) h.urb.10.lm.1 <- lm(log(dens_above10+0.01)~mb_dweldens.s+road_density.s+park_prop.s+lot_size.s, data=hdata.1) h.urb.10.lm.2 <- lm(log(dens_above10+0.01)~poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), data=hdata.1) comp.models(h.urb.10.lm.1, h.urb.10.lm.2) h.sar.urb.10 <- mm.spatial.sparse(hdata.1, 150, h.urb.10.lm.2, "MC") mm.sum.sp(h.sar.urb.10) spatial.plot(urb.10.lm.2, data.1, ">10m Veg Density - Urban form polynomial non-spatial model") spatial.plot(sar.cens.10$mixed, data.1, ">10m Veg Density - Urban form mixed SAR model") h.best.urb.10 <- lagsarlm(log(dens_above10+0.01) ~ poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), listw=h.sar.urb.10$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.urb.10) summary(h.best.urb.10) diag.best(best.urb.10, data.1, sar.urb.10$weight, ">10m Veg Density - Best Urban Form SAR model") h.best.urb.10.null <- lagsarlm(log(dens_above10+0.01) ~ 1, listw=h.sar.urb.10$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.urb.10.lm, h.urb.10.lm.2, h.best.urb.10, h.best.urb.10.null, hdata.1) h.best.urb.10.raw <- lagsarlm(log(dens_above10+0.01)~poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+ poly(park_prop.s,2,raw=T)+poly(lot_size.s,2,raw=T), listw=h.sar.urb.10$weight, data=hdata.1, method="MC", type="mixed", na.action="na.fail") plot.partial(h.best.urb.10, h.best.urb.10.raw, hdata.1, 2, "road_density", "road_density.s", "Street density (km/km2)", ">10m Veg. Density Partial Residuals") plot.partial(h.best.urb.10, h.best.urb.10.raw, hdata.1, 3, "park_prop_per", "park_prop.s", "Proportion parkland (%)", ">10m Veg. Density Partial Residuals") plot.partial(h.best.urb.10, h.best.urb.10.raw, hdata.1, 1, "mb_dwel_dens", "mb_dweldens.s", "Dwelling density (dwellings/ha)", ">10m Veg. Density Partial Residuals") plot.partial(h.best.urb.10, h.best.urb.10.raw, hdata.1, 4, "log_lot_size", "log_lot_size.s", "Mean lot size (m2)", ">10m Veg. Density Partial Residuals") ####################################### ####Landscape variables - vegetation density >10m #### h.land.10.lm <- lm(log(dens_above10+0.01)~1, data=hdata.1) h.land.10.lm.1 <- lm(log(dens_above10+0.01)~clumpy.s+number_patches.s+perimeter_area_mn.s+tree_area.s, data=hdata.1) h.land.10.lm.2 <- lm(log(dens_above10+0.01)~poly(clumpy.s,2)+poly(number_patches.s,2)+poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), data=hdata.1) comp.models(h.land.10.lm.1, h.land.10.lm.2) h.sar.land.10 <- mm.spatial.sparse(hdata.1, 150, h.land.10.lm.2, "MC") mm.sum.sp(h.sar.land.10) spatial.plot(land.10.lm.2, data.1, ">10m Veg Density - Landscape polynomial non-spatial model") spatial.plot(sar.land.10$mixed, data.1, ">10m Veg Density - Landscape mixed SAR model") h.best.land.10 <- lagsarlm(log(dens_above10+0.01) ~ poly(clumpy.s,2)+poly(number_patches.s,2)+ poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), listw=h.sar.land.10$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.land.10) summary(h.best.land.10) diag.best(best.land.10, data.1, sar.land.10$weight, ">10m Veg Density - Best Landscape SAR model") h.best.land.10.null <- lagsarlm(log(dens_above10+0.01) ~ 1, listw=h.sar.land.10$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.land.10.lm, h.land.10.lm.2, h.best.land.10, h.best.land.10.null, hdata.1) h.best.land.10.raw <- lagsarlm(log(dens_above10+0.01)~poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+ poly(perimeter_area_mn.s,2,raw=T)+poly(tree_area.s,2,raw=T), listw=h.sar.land.10$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") plot.partial(h.best.land.10, h.best.land.10.raw, hdata.1, 1, "clumpy", "clumpy.s", "Clumpiness", ">10m Veg. Density Partial Residuals") plot.partial(h.best.land.10, h.best.land.10.raw, hdata.1, 2, "number_patches", "number_patches.s", "Number forest fragments", ">10m Veg. Density Partial Residuals") plot.partial(h.best.land.10, h.best.land.10.raw, hdata.1, 3, "perimeter_area_mn", "perimeter_area_mn.s", "Forest perimeter:area ratio", ">10m Veg. Density Partial Residuals") plot.partial(h.best.land.10, h.best.land.10.raw, hdata.1, 4, "tree_area", "tree_area.s", "Proportion forest (%)", ">10m Veg. Density Partial Residuals") ####################################### ####Combined model - vegetation density >10m #### h.comb.10.lm <- lm(log(dens_above10+0.01)~1, data=hdata.1) h.comb.10.lm.1 <- lm(log(dens_above10+0.01)~aspect_cos.s+aspect_sin.s+elev.s+slope.s+ tot_p.s+water_cap.s+soc.s+tot_n.s+sa1_medtothinc.s+sa1_avghouse.s+ sa1_medage.s+prop_noncauc.s+mb_dweldens.s+road_density.s+lot_size.s+park_prop.s+ clumpy.s+number_patches.s+tree_area.s+perimeter_area_mn.s, data=hdata.1) h.comb.10.lm.2 <- lm(log(dens_above10+0.01)~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), data=hdata.1) comp.models(h.comb.10.lm.1, h.comb.10.lm.2) h.sar.comb.10 <- mm.spatial.sparse(hdata.1, 150, h.comb.10.lm.2, "MC") mm.sum.sp(h.sar.comb.10) spatial.plot(comb.10.lm.2, data.1, ">10m Veg Density - Combined polynomial non-spatial model") spatial.plot(sar.comb.10$mixed, data.1, ">10m Veg Density - Comnbined mixed SAR model") h.best.comb.10 <- lagsarlm(log(dens_above10+0.01) ~ poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), type="mixed", listw=h.sar.comb.510$weight, data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.comb.10) summary(h.best.comb.10) diag.best(best.comb.10, data.1, sar.comb.10$weight, ">10m Veg Density - Best Combined SAR model") h.best.comb.10.null <- lagsarlm(log(dens_above10+0.01) ~ 1, type="mixed", listw=h.sar.comb.10$weight, data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.comb.10.lm, h.comb.10.lm.2, h.best.comb.10, h.best.comb.10.null, hdata.1) h.best.comb.10.raw <- lagsarlm(log(dens_above10+0.01)~poly(aspect_cos.s,2,raw=T)+poly(aspect_sin.s,2,raw=T)+poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(tot_p.s,2,raw=T)+poly(water_cap.s,2,raw=T)+poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_avghouse.s,2,raw=T)+ poly(sa1_medage.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T)+poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+poly(lot_size.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+poly(tree_area.s,2,raw=T)+poly(perimeter_area_mn.s,2,raw=T), type="mixed", listw=h.sar.comb.510$weight, data=hdata.1, method="MC", na.action="na.fail") h.best.comb.10.raw.avg <- dredge(h.best.comb.10.raw, m.max=4, trace=2) h.best.comb.10.raw.avg2 <- model.avg(h.best.comb.10.raw.avg) summary(h.best.comb.10.raw.avg2) h.comb.10.avg.best <- get.models(h.best.comb.10.raw.avg,1)[[1]] summary(get.models(h.best.comb.10.raw.avg,1)[[1]]) ##includes tree area, clumpy, number patches, park prop h.best.simple.10 <- lagsarlm(log(dens_above10+0.01)~poly(tree_area.s,2)+poly(number_patches.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2), listw=h.sar.comb.10$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.simple.10.raw <- lagsarlm(log(dens_above10+0.01)~poly(tree_area.s,2,raw=T)+poly(number_patches.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(clumpy.s,2,raw=T), listw=h.sar.comb.10$weight,data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.10.null <- lagsarlm(log(dens_above10+0.01)~1, listw=h.sar.comb.10$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") AICc(h.best.simple.10) summary(h.best.simple.10) summary(h.best.simple.10.raw) spatial.diag(h.comb.10.lm, h.comb.10.lm.2, h.best.simple.10, h.best.10.null, hdata.1) h.comb.10.best.lm.2 <- lm(log(dens_above10+0.01)~poly(tree_area.s,2)+poly(number_patches.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2), data=hdata.1) h.sar.best.10 <- mm.spatial.sparse(hdata.1, 150, h.comb.10.best.lm.2, "MC") mm.sum.sp(h.sar.best.10) plot.partial(h.best.simple.10, h.best.simple.10.raw, hdata.1, 4, "clumpy", "clumpy.s", "Tree cover clumpiness", ">10m Veg. Density Partial Residuals") plot.partial(h.best.simple.10, h.best.simple.10.raw, hdata.1, 3, "park_prop", "park_prop.s", "Proportion parkland (%)", ">10m Veg. Density Partial Residuals") plot.partial(h.best.simple.10, h.best.simple.10.raw, hdata.1, 2, "number_patches", "number_patches.s", "Number of tree cover patches", ">10m Veg. Density Partial Residuals") plot.partial(h.best.simple.10, h.best.simple.10.raw, hdata.1, 1, "tree_area", "tree_area.s", "Proportion tree cover (%)", ">10m Veg. Density Partial Residuals") ##################################################################### #####Model selection for Foliage Projective Cover ##################################################################### ####################################### ####Physical variables - FPC #### h.phys.fpc.lm <- lm(fpc~1, data=hdata.1) h.phys.fpc.lm.1 <- lm(fpc~aspect_cos.s+aspect_sin.s+elev.s+slope.s, data=hdata.1) h.phys.fpc.lm.2 <- lm(fpc~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2), data=hdata.1) comp.models(h.phys.fpc.lm.1, h.phys.fpc.lm.2) h.sar.phys.fpc <- mm.spatial.sparse(hdata.1, 150, h.phys.fpc.lm.2, "MC") mm.sum.sp(h.sar.phys.fpc) spatial.plot(phys.fpc.lm.2, data.1, "FPC - Polynomial non-spatial model") spatial.plot(sar.phys.fpc$error, data.1, "FPC - Mixed SAR model") h.best.phys.fpc <- lagsarlm(fpc ~ poly(elev.s,2)+poly(aspect_cos.s,2)+poly(slope.s,2)+poly(aspect_sin.s,2), listw=h.sar.phys.fpc$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") AICc(h.best.phys.fpc) summary(h.best.phys.fpc) diag.best(best.phys.fpc, data.1, sar.phys.fpc$weight, "FPC - Best SAR model") h.best.phys.fpc.null <- lagsarlm(fpc ~ 1, listw=h.sar.phys.fpc$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") spatial.diag(h.phys.fpc.lm, h.phys.fpc.lm.2, h.best.phys.fpc, h.best.phys.fpc.null, hdata.1) h.best.phys.fpc.raw <- lagsarlm(fpc ~ poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(aspect_cos.s,2,raw=T)+poly(aspect_sin,2,raw=T), listw=h.sar.phys.fpc$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.phys.fpc, h.best.phys.fpc.raw, hdata.1, 3, "aspect_cos", "aspect_cos.s", "Northness", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.phys.fpc, h.best.phys.fpc.raw, hdata.1, 4, "aspect_sin", "aspect_sin.s", "Eastness", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.phys.fpc, h.best.phys.fpc.raw, hdata.1, 1, "elev", "elev.s", "Elevation (m)", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.phys.fpc, h.best.phys.fpc.raw, hdata.1, 2, "slope", "slope.s", "Slope (degrees)", "Foliage Projective Cover Partial Residuals") ####################################### ####Soil variables - FPC #### h.soil.fpc.lm <- lm(fpc~1, data=hdata.1) h.soil.fpc.lm.1 <- lm(fpc~water_cap.s+tot_p.s+soc.s+tot_n.s, data=hdata.1) h.soil.fpc.lm.2 <- lm(fpc~poly(water_cap.s,2)+poly(tot_p.s,2)+poly(soc.s,2)+poly(tot_n.s,2), data=hdata.1) comp.models(h.soil.fpc.lm.1, h.soil.fpc.lm.2) h.sar.soil.fpc <- mm.spatial.sparse(hdata.1, 150, h.soil.fpc.lm.2, "MC") mm.sum.sp(h.sar.soil.fpc) spatial.plot(soil.fpc.lm.2, data.1, "FPC - Soil polynomial non-spatial model") spatial.plot(sar.soil.fpc$mixed, data.1, "FPC - Soil mixed SAR model") h.best.soil.fpc <- lagsarlm(fpc ~ poly(soc.s,2)+poly(tot_p.s,2)+poly(water_cap.s,2)+poly(tot_n.s,2), listw=h.sar.soil.fpc$weight, type="lag", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.soil.fpc) summary(h.best.soil.fpc) diag.best(best.soil.fpc, data.1, sar.soil.fpc$weight, "FPC - Best Soil SAR model") h.best.soil.fpc.null <- lagsarlm(fpc ~ 1, listw=h.sar.soil.fpc$weight, type="lag", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.soil.fpc.lm, h.soil.fpc.lm.2, h.best.soil.fpc, h.best.soil.fpc.null, hdata.1) h.best.soil.fpc.raw <- lagsarlm(fpc ~ poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+ poly(water_cap.s,2,raw=T)+poly(tot_p.s,2,raw=T), listw=h.sar.soil.fpc$weight, type="lag", data=hdata.1, method="MC", na.action="na.fail") plot.partial(h.best.soil.fpc, h.best.soil.fpc.raw, hdata.1, 1, "soc", "soc.s", "Soil organic carbon (%)", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.soil.fpc, h.best.soil.fpc.raw, hdata.1, 2, "tot_n", "tot_n.s", "Total soil N (%)", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.soil.fpc, h.best.soil.fpc.raw, hdata.1, 3, "water_cap", "water_cap.s", "Soil water capacity (%)", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.soil.fpc, h.best.soil.fpc.raw, hdata.1, 4, "tot_p", "tot_p.s", "Total soil P (%)", "Foliage Projective Cover Partial Residuals") ####################################### ####Census variables - FPC #### h.cens.fpc.lm <- lm(fpc~1, data=hdata.1) h.cens.fpc.lm.1 <- lm(fpc~prop_noncauc.s+sa1_medtothinc.s+sa1_medage.s+sa1_avghouse.s, data=hdata.1) h.cens.fpc.lm.2 <- lm(fpc~poly(prop_noncauc.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2), data=hdata.1) comp.models(h.cens.fpc.lm.1, h.cens.fpc.lm.2) h.sar.cens.fpc <- mm.spatial.sparse(hdata.1, 150, h.cens.fpc.lm.2, "MC") mm.sum.sp(h.sar.cens.fpc) spatial.plot(cens.fpc.lm.2, data.1, "FPC - Census polynomial non-spatial model") spatial.plot(sar.cens.fpc$mixed, data.1, "FPC - Census mixed SAR model") h.best.cens.fpc <- lagsarlm(fpc ~ poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2)+poly(prop_noncauc.s,2), listw=h.sar.cens.fpc$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.cens.fpc) summary(h.best.cens.fpc) diag.best(best.cens.fpc, data.1, sar.cens.fpc$weight, "FPC - Best Soil SAR model") h.best.cens.fpc.null <- lagsarlm(fpc ~ 1, listw=h.sar.cens.fpc$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.cens.fpc.lm, h.cens.fpc.lm.2, h.best.cens.fpc, h.best.cens.fpc.null, hdata.1) h.best.cens.fpc.raw <- lagsarlm(fpc~poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_medage.s,2,raw=T)+ poly(sa1_avghouse.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T), listw=h.sar.cens.fpc$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") plot.partial(h.best.cens.fpc, h.best.cens.fpc.raw, hdata.1, 2, "sa1_medage", "sa1_medage.s", "Median pop. age", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.cens.fpc, h.best.cens.fpc.raw, hdata.1, 1, "sa1_medtothinc", "sa1_medtothinc.s", "Median weekly household income ($)", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.cens.fpc, h.best.cens.fpc.raw, hdata.1, 4, "prop_noncauc", "prop_noncauc.s", "Proportion Non-European Descent", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.cens.fpc, h.best.cens.fpc.raw, hdata.1, 3, "sa1_avghouse", "sa1_avghouse.s", "Avg. household size", "Foliage Projective Cover Partial Residuals") ####################################### ####Urban form variables - vegetation density FPC #### h.urb.fpc.lm <- lm(fpc~1, data=hdata.1) h.urb.fpc.lm.1 <- lm(fpc~mb_dweldens.s+road_density.s+park_prop.s+lot_size.s, data=hdata.1) h.urb.fpc.lm.2 <- lm(fpc~poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), data=hdata.1) comp.models(h.urb.fpc.lm.1, h.urb.fpc.lm.2) h.sar.urb.fpc <- mm.spatial.sparse(hdata.1, 150, h.urb.fpc.lm.2, "MC") mm.sum.sp(h.sar.urb.fpc) spatial.plot(urb.fpc.lm.2, data.1, "Foliage Projective Cover - Urban form polynomial non-spatial model") spatial.plot(sar.cens.fpc$mixed, data.1, "Foliage Projective Cover - Urban form mixed SAR model") h.best.urb.fpc <- lagsarlm(fpc ~ poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), listw=h.sar.urb.fpc$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.urb.fpc) summary(h.best.urb.fpc) diag.best(best.urb.fpc, data.1, sar.urb.fpc$weight, "Foliage Projective Cover - Best Urban Form SAR model") h.best.urb.fpc.null <- lagsarlm(fpc ~ 1, listw=h.sar.urb.fpc$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.urb.fpc.lm, h.urb.fpc.lm.2, h.best.urb.fpc, h.best.urb.fpc.null, hdata.1) h.best.urb.fpc.raw <- lagsarlm(fpc~poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+ poly(park_prop.s,2,raw=T)+poly(lot_size.s,2,raw=T), listw=h.sar.urb.fpc$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.urb.fpc, h.best.urb.fpc.raw, hdata.1, 2, "road_density", "road_density.s", "Street density (km/km2)", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.urb.fpc, h.best.urb.fpc.raw, hdata.1, 3, "park_prop_per", "park_prop.s", "Proportion parkland (%)", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.urb.fpc, h.best.urb.fpc.raw, hdata.1, 1, "mb_dwel_dens", "mb_dweldens.s", "Dwelling density (dwellings/ha)", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.urb.fpc, h.best.urb.fpc.raw, hdata.1, 4, "log_lot_size", "log_lot_size.s", "Mean lot size (m2)", "Foliage Projective Cover Partial Residuals") ####################################### ####Landscape variables - FPC #### h.land.fpc.lm <- lm(fpc~1, data=hdata.1) h.land.fpc.lm.1 <- lm(fpc~clumpy.s+number_patches.s+perimeter_area_mn.s+tree_area.s, data=hdata.1) h.land.fpc.lm.2 <- lm(fpc~poly(clumpy.s,2)+poly(number_patches.s,2)+poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), data=hdata.1) comp.models(h.land.fpc.lm.1, h.land.fpc.lm.2) h.sar.land.fpc <- mm.spatial.sparse(hdata.1, 150, h.land.fpc.lm.2, "MC") mm.sum.sp(h.sar.land.fpc) spatial.plot(land.fpc.lm.2, data.1, "FPC - Landscape polynomial non-spatial model") spatial.plot(sar.land.fpc$mixed, data.1, "FPC - Landscape mixed SAR model") h.best.land.fpc <- lagsarlm(fpc ~ poly(clumpy.s,2)+ poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), listw=h.sar.land.fpc$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.land.fpc) summary(h.best.land.fpc) diag.best(best.land.fpc, data.1, sar.land.fpc$weight, "FPC - Best Landscape SAR model") h.best.land.fpc.null <- lagsarlm(fpc ~ 1, listw=h.sar.land.fpc$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.land.fpc.lm, h.land.fpc.lm.2, h.best.land.fpc, h.best.land.fpc.null, hdata.1) h.best.land.fpc.raw <- lagsarlm(fpc~poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+ poly(perimeter_area_mn.s,2,raw=T)+poly(tree_area.s,2,raw=T), listw=h.sar.land.fpc$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") plot.partial(h.best.land.fpc, h.best.land.fpc.raw, hdata.1, 1, "clumpy", "clumpy.s", "Clumpiness", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.land.fpc, h.best.land.fpc.raw, hdata.1, 2, "number_patches", "number_patches.s", "Number forest fragments", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.land.fpc, h.best.land.fpc.raw, hdata.1, 3, "perimeter_area_mn", "perimeter_area_mn.s", "Forest perimeter:area ratio", "Foliage Projective Cover Partial Residuals") plot.partial(h.best.land.fpc, h.best.land.fpc.raw, hdata.1, 4, "tree_area", "tree_area.s", "Proportion forest (%)", "Foliage Projective Cover Partial Residuals") ####################################### ####Combined model - vegetation density FPC #### h.comb.fpc.lm <- lm(fpc~1, data=hdata.1) h.comb.fpc.lm.1 <- lm(fpc~aspect_cos.s+aspect_sin.s+elev.s+slope.s+ tot_p.s+water_cap.s+soc.s+tot_n.s+sa1_medtothinc.s+sa1_avghouse.s+ sa1_medage.s+prop_noncauc.s+mb_dweldens.s+road_density.s+lot_size.s+park_prop.s+ clumpy.s+number_patches.s+tree_area.s+perimeter_area_mn.s, data=hdata.1) h.comb.fpc.lm.2 <- lm(fpc~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), data=hdata.1) comp.models(h.comb.fpc.lm.1, h.comb.fpc.lm.2) h.sar.comb.fpc <- mm.spatial.sparse(hdata.1, 150, h.comb.fpc.lm.2, "MC") mm.sum.sp(h.sar.comb.fpc) spatial.plot(comb.fpc.lm.2, data.1, "FPC - Combined polynomial non-spatial model") spatial.plot(sar.comb.fpc$mixed, data.1, "FPC - Comnbined mixed SAR model") h.best.comb.fpc <- lagsarlm(fpc ~ poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), type="mixed", listw=h.sar.comb.fpc$weight, data=hdata.1, method="MC", na.action="na.fail") AICc(h.best.comb.fpc) summary(h.best.comb.fpc) diag.best(best.comb.fpc, data.1, sar.comb.fpc$weight, "FPC - Best Combined SAR model") h.best.comb.fpc.null <- lagsarlm(fpc ~ 1, type="mixed", listw=h.sar.comb.fpc$weight, data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.comb.fpc.lm, h.comb.fpc.lm.2, h.best.comb.fpc, h.best.comb.fpc.null, hdata.1) h.best.comb.fpc.raw <- lagsarlm(fpc~poly(aspect_cos.s,2,raw=T)+poly(aspect_sin.s,2,raw=T)+poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(tot_p.s,2,raw=T)+poly(water_cap.s,2,raw=T)+poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_avghouse.s,2,raw=T)+ poly(sa1_medage.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T)+poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+poly(lot_size.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+poly(tree_area.s,2,raw=T)+poly(perimeter_area_mn.s,2,raw=T), type="mixed", listw=h.sar.comb.fpc$weight, data=hdata.1, method="MC", na.action="na.fail") h.best.comb.fpc.raw.avg <- dredge(h.best.comb.fpc.raw, m.max=4, trace=2) h.best.comb.fpc.raw.avg2 <- model.avg(h.best.comb.fpc.raw.avg) summary(h.best.comb.fpc.raw.avg2) h.comb.fpc.avg.best <- get.models(h.best.comb.fpc.raw.avg,1)[[1]] summary(get.models(h.best.comb.fpc.raw.avg,1)[[1]]) ##includes tree area, clumpy, park prop, avg house h.best.simple.fpc <- lagsarlm(fpc~poly(tree_area.s,2)+poly(clumpy.s,2)+poly(park_prop.s,2)+ poly(sa1_avghouse.s,2), listw=h.sar.comb.fpc$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.simple.fpc.raw <- lagsarlm(fpc~poly(tree_area.s,2,raw=T)+poly(clumpy.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(sa1_avghouse.s,2,raw=T), listw=h.sar.comb.fpc$weight,data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.fpc.null <- lagsarlm(fpc~1, listw=h.sar.comb.fpc$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") AICc(h.best.simple.fpc) summary(h.best.simple.fpc) summary(h.best.simple.fpc.raw) spatial.diag(h.comb.fpc.lm, h.comb.fpc.lm.2, h.best.simple.fpc, h.best.fpc.null, hdata.1) h.comb.fpc.best.lm.2 <- lm(fpc~poly(tree_area.s,2)+poly(clumpy.s,2)+poly(park_prop.s,2)+ poly(sa1_avghouse.s,2), data=hdata.1) h.sar.best.fpc <- mm.spatial.sparse(hdata.1, 150, h.comb.fpc.best.lm.2, "MC") mm.sum.sp(h.sar.best.fpc) plot.partial(h.best.simple.fpc, h.best.simple.fpc.raw, hdata.1, 2, "clumpy", "clumpy.s", "Clumpiness", "FPC Partial Residuals") plot.partial(h.best.simple.fpc, h.best.simple.fpc.raw, hdata.1, 4, "sa1_avghouse", "sa1_avghouse.s", "Average household size (persons)", "FPC Partial Residuals") plot.partial(h.best.simple.fpc, h.best.simple.fpc.raw, hdata.1, 3, "park_prop", "park_prop.s", "Proportion parkland (%)", "FPC Partial Residuals") plot.partial(h.best.simple.fpc, h.best.simple.fpc.raw, hdata.1, 1, "tree_area_per", "tree_area.s", "Proportion tree cover (%)", "FPC Partial Residuals") ##################################################################### #####Model selection for FHD ##################################################################### ##Calculate Shannon Index based on vegetation densities of different layers library(vegan) data.shannon <- hdata.1[,c(12:16)] sv <- diversity(data.shannon, index="shannon") hdata.9 <- cbind(hdata.1,sv) ####################################### ####Physical variables - FHD #### hdata.9$aspect_cos.s <- rescale(hdata.9$aspect_cos) hdata.9$aspect_sin.s <- rescale(hdata.9$aspect_sin) hdata.9$elev.s <- rescale(hdata.9$elev) hdata.9$slope.s <- rescale(hdata.9$slope) h.phys.sv.lm <- lm(sv~1, data=hdata.9) h.phys.sv.lm.1 <- lm(sv~aspect_cos.s+aspect_sin.s+elev.s+slope.s, data=hdata.9) h.phys.sv.lm.2 <- lm(sv~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2), data=hdata.9) comp.models(h.phys.sv.lm.1, h.phys.sv.lm.2) h.sar.phys.sv <- mm.spatial.sparse(hdata.9, 150, h.phys.sv.lm.2, "MC") mm.sum.sp(h.sar.phys.sv) spatial.plot(phys.sv.lm.2, data.9, "FHD - Polynomial non-spatial model") spatial.plot(sar.phys.sv$mixed, data.9, "FHD - Mixed SAR model") h.best.phys.sv <- lagsarlm(sv ~ poly(elev.s,2)+poly(slope.s,2)+poly(aspect_sin.s,2)+poly(aspect_cos.s,2), listw=h.sar.phys.sv$weight, type="mixed", method="MC", data=hdata.9, na.action="na.fail") AICc(h.best.phys.sv) summary(h.best.phys.sv) diag.best(best.phys.sv, data.9, sar.phys.sv$weight, "FHD - Best SAR model") h.best.phys.sv.null <- lagsarlm(sv ~ 1, listw=h.sar.phys.sv$weight, type="mixed", method="MC", data=hdata.9, na.action="na.fail") spatial.diag(h.phys.sv.lm, h.phys.sv.lm.2, h.best.phys.sv, h.best.phys.sv.null, hdata.9) h.best.phys.sv.raw <- lagsarlm(sv ~ poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(aspect_cos.s,2,raw=T)+poly(aspect_sin,2,raw=T), listw=h.sar.phys.sv$weight, type="mixed", method="MC", data=hdata.9, na.action="na.fail") ##Partial residual plotting to look at trends plot.partial(h.best.phys.sv, h.best.phys.sv.raw, hdata.9, 3, "aspect_cos", "aspect_cos.s", "Northness", "Foliage Height Diversity Partial Residuals") plot.partial(h.best.phys.sv, h.best.phys.sv.raw, hdata.9, 4, "aspect_sin", "aspect_sin.s", "Eastness", "Foliage Height Diversity Partial Residuals") plot.partial(h.best.phys.sv, h.best.phys.sv.raw, hdata.9, 1, "elev", "elev.s", "Elevation (m)", "Foliage Height Diversity Partial Residuals") plot.partial(h.best.phys.sv, h.best.phys.sv.raw, hdata.9, 2, "slope", "slope.s", "Slope (degrees)", "Foliage Height Diversity Partial Residuals") ####################################### ####Soil variables - FHD #### hdata.9$water_cap.s <- rescale(hdata.9$water_cap) hdata.9$tot_p.s <- rescale(hdata.9$tot_p) hdata.9$soc.s <- rescale(hdata.9$soc) hdata.9$tot_n.s <- rescale(hdata.9$tot_n) h.soil.sv.lm <- lm(sv~1, data=hdata.9) h.soil.sv.lm.1 <- lm(sv~water_cap.s+tot_p.s+soc.s+tot_n.s, data=hdata.9) h.soil.sv.lm.2 <- lm(sv~poly(water_cap.s,2)+poly(tot_p.s,2)+poly(soc.s,2)+poly(tot_n.s,2), data=hdata.9) comp.models(h.soil.sv.lm.1, h.soil.sv.lm.2) h.sar.soil.sv <- mm.spatial.sparse(hdata.9, 150, h.soil.sv.lm.2, "MC") mm.sum.sp(h.sar.soil.sv) spatial.plot(soil.sv.lm.2, data.9, "FHD - Soil polynomial non-spatial model") spatial.plot(sar.soil.sv$mixed, data.9, "FHD - Soil mixed SAR model") h.best.soil.sv <- lagsarlm(sv ~ poly(soc.s,2)+poly(tot_p.s,2)+poly(tot_n.s,2)+poly(water_cap.s,2), listw=h.sar.soil.sv$weight, type="mixed", method="MC", data=hdata.9, na.action="na.fail") AICc(h.best.soil.sv) summary(h.best.soil.sv) diag.best(best.soil.sv, data.9, sar.soil.sv$weight, "FHD - Best Soil SAR model") h.best.soil.sv.null <- lagsarlm(sv ~ 1, listw=h.sar.soil.sv$weight, type="mixed", method="MC", data=hdata.9, na.action="na.fail") spatial.diag(h.soil.sv.lm, h.soil.sv.lm.2, h.best.soil.sv, h.best.soil.sv.null, hdata.9) h.best.soil.sv.raw <- lagsarlm(sv ~ poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+ poly(water_cap.s,2,raw=T)+poly(tot_p.s,2,raw=T), listw=h.sar.soil.sv$weight, type="mixed", method="MC", data=hdata.9, na.action="na.fail") plot.partial(h.best.soil.sv, h.best.soil.sv.raw, hdata.9, 1, "soc", "soc.s", "Soil organic carbon (%)", "Foliage Height Diversity Partial Residuals") plot.partial(h.best.soil.sv, h.best.soil.sv.raw, hdata.9, 2, "tot_n", "tot_n.s", "Total soil N (%)", "Foliage Height Diversity Partial Residuals") plot.partial(h.best.soil.sv, h.best.soil.sv.raw, hdata.9, 3, "water_cap", "water_cap.s", "Soil water capacity (%)", "Foliage Height Diversity Partial Residuals") plot.partial(h.best.soil.sv, h.best.soil.sv.raw, hdata.9, 4, "tot_p", "tot_p.s", "Total soil P (%)", "Foliage Height Diversity Partial Residuals") ####################################### ####Census variables - FHD #### hdata.9$prop_noncauc.s <- rescale(hdata.9$prop_noncauc) hdata.9$sa1_avghouse.s <- rescale(hdata.9$sa1_avghouse) hdata.9$sa1_medage.s <- rescale(hdata.9$sa1_medage) hdata.9$sa1_medtothinc.s <- rescale(hdata.9$sa1_medtothinc) h.cens.sv.lm <- lm(sv~1, data=hdata.9) h.cens.sv.lm.1 <- lm(sv~prop_noncauc.s+sa1_medtothinc.s+sa1_medage.s+sa1_avghouse.s, data=hdata.9) h.cens.sv.lm.2 <- lm(sv~poly(prop_noncauc.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2), data=hdata.9) comp.models(h.cens.sv.lm.1, h.cens.sv.lm.2) h.sar.cens.sv <- mm.spatial.sparse(hdata.9, 150, h.cens.sv.lm.2, "MC") mm.sum.sp(h.sar.cens.sv) spatial.plot(cens.sv.lm.2, data.9, "FHD - Census polynomial non-spatial model") spatial.plot(sar.cens.sv$lag, data.9, "FHD - Census mixed SAR model") h.best.cens.sv <- lagsarlm(sv ~ poly(prop_noncauc.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2), listw=h.sar.cens.sv$weight, type="mixed", method="MC", data=hdata.9, na.action="na.fail") AICc(h.best.cens.sv) summary(h.best.cens.sv) diag.best(best.cens.sv, data.9, sar.cens.sv$weight, "FHD - Best Soil SAR model") h.best.cens.sv.null <- lagsarlm(sv ~ 1, listw=h.sar.cens.sv$weight, type="mixed", method="MC", data=hdata.9, na.action="na.fail") spatial.diag(h.cens.sv.lm, h.cens.sv.lm.2, h.best.cens.sv, h.best.cens.sv.null, hdata.9) h.best.cens.sv.raw <- lagsarlm(sv~poly(sa1_avghouse.s,2,raw=T)+poly(sa1_medage.s,2,raw=T)+ poly(sa1_medtothinc.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T), listw=h.sar.cens.sv$weight, type="mixed", data=hdata.9, method="MC", na.action="na.fail") plot.partial(h.best.cens.sv, h.best.cens.sv.raw, hdata.9, 2, "sa1_medage", "sa1_medage.s", "Median pop. age", "Foliage Height Diversity Partial Residuals") plot.partial(h.best.cens.sv, h.best.cens.sv.raw, hdata.9, 3, "sa1_medtothinc", "sa1_medtothinc.s", "Median weekly household income ($)", "Foliage Height Diversity Partial Residuals") plot.partial(h.best.cens.sv, h.best.cens.sv.raw, hdata.9, 4, "prop_noncauc", "prop_noncauc.s", "Proportion Non-European Descent", "Foliage Height Diversity Partial Residuals") plot.partial(h.best.cens.sv, h.best.cens.sv.raw, hdata.9, 1, "sa1_avghouse", "sa1_avghouse.s", "Avg. household size", "Foliage Height Diversity Partial Residuals") ####################################### ####Urban form variables - vegetation density FHD #### hdata.9$road_density <- hdata.9$road_length/1000 hdata.9$road_density.s <- rescale(hdata.9$road_density) hdata.9$mb_dweldens.s <- rescale(hdata.9$mb_dwel_dens) hdata.9$park_prop.s <- rescale(hdata.9$park_prop) hdata.9$lot_size.s <- rescale(hdata.9$lot_size) hdata.9$log_lot_size <- log10(hdata.9$lot_size) hdata.9$log_lot_size.s <- rescale(hdata.9$log_lot_size) hdata.9$park_prop_per <- hdata.9$park_prop/10000 h.urb.sv.lm <- lm(sv~1, data=hdata.9) h.urb.sv.lm.1 <- lm(sv~mb_dweldens.s+road_density.s+park_prop.s+lot_size.s, data=hdata.9) h.urb.sv.lm.2 <- lm(sv~poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), data=hdata.9) comp.models(h.urb.sv.lm.1, h.urb.sv.lm.2) h.sar.urb.sv <- mm.spatial.sparse(hdata.9, 150, h.urb.sv.lm.2, "MC") mm.sum.sp(h.sar.urb.sv) spatial.plot(urb.sv.lm.2, data.9, "FHD - Urban form polynomial non-spatial model") spatial.plot(sar.cens.sv$mixed, data.9, "FHD - Urban form mixed SAR model") h.best.urb.sv <- lagsarlm(sv ~ poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), listw=h.sar.urb.sv$weight, type="mixed", method="MC", data=hdata.9, na.action="na.fail") AICc(h.best.urb.sv) summary(h.best.urb.sv) diag.best(best.urb.sv, data.9, sar.urb.sv$weight, "FHD - Best Urban Form SAR model") h.best.urb.sv.null <- lagsarlm(sv ~ 1, listw=h.sar.urb.sv$weight, type="mixed", data=hdata.9, method="MC", na.action="na.fail") spatial.diag(h.urb.sv.lm, h.urb.sv.lm.2, h.best.urb.sv, h.best.urb.sv.null, hdata.9) h.best.urb.sv.raw <- lagsarlm(sv~poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+ poly(park_prop.s,2,raw=T)+poly(lot_size.s,2,raw=T), listw=h.sar.urb.sv$weight, data=hdata.9, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.urb.sv, h.best.urb.sv.raw, hdata.9, 2, "road_density", "road_density.s", "Street density (km/km2)", "FHD Partial Residuals") plot.partial(h.best.urb.sv, h.best.urb.sv.raw, hdata.9, 3, "park_prop_per", "park_prop.s", "Proportion parkland (%)", "FHD Partial Residuals") plot.partial(h.best.urb.sv, h.best.urb.sv.raw, hdata.9, 1, "mb_dwel_dens", "mb_dweldens.s", "Dwelling density (dwellings/ha)", "FHD Partial Residuals") plot.partial(h.best.urb.sv, h.best.urb.sv.raw, hdata.9, 4, "log_lot_size", "log_lot_size.s", "Mean lot size (m2)", "FHD Partial Residuals") ####################################### ####Landscape variables - FHD #### hdata.9$clumpy.s <- rescale(hdata.9$clumpy) hdata.9$patch_area_mn.s <- rescale(hdata.9$patch_area_mn) hdata.9$perimeter_area_mn.s <- rescale(hdata.9$perimeter_area_mn) hdata.9$tree_area.s <- rescale(hdata.9$tree_area) hdata.9$number_patches.s <- rescale(hdata.9$number_patches) h.land.sv.lm <- lm(sv~1, data=hdata.9) h.land.sv.lm.1 <- lm(sv~clumpy.s+number_patches.s+perimeter_area_mn.s+tree_area.s, data=hdata.9) h.land.sv.lm.2 <- lm(sv~poly(clumpy.s,2)+poly(number_patches.s,2)+poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), data=hdata.9) comp.models(h.land.sv.lm.1, h.land.sv.lm.2) h.sar.land.sv <- mm.spatial.sparse(hdata.9, 150, h.land.sv.lm.2, "MC") mm.sum.sp(h.sar.land.sv) spatial.plot(land.sv.lm.2, data.9, "FHD - Landscape polynomial non-spatial model") spatial.plot(sar.land.sv$error, data.9, "FHD - Landscape error SAR model") h.best.land.sv <- lagsarlm(sv ~ poly(clumpy.s,2)+poly(number_patches.s,2)+ poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), listw=h.sar.land.sv$weight, data=hdata.9, type="mixed", method="MC", na.action="na.fail") AICc(h.best.land.sv) summary(h.best.land.sv) diag.best(best.land.sv, data.9, sar.land.sv$weight, "FHD - Best Landscape SAR model") h.best.land.sv.null <- lagsarlm(sv ~ 1, listw=h.sar.land.sv$weight, data=hdata.9, type="mixed", method="MC", na.action="na.fail") spatial.diag(h.land.sv.lm, h.land.sv.lm.2, h.best.land.sv, h.best.land.sv.null, hdata.9) h.best.land.sv.raw <- lagsarlm(sv~poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+ poly(perimeter_area_mn.s,2,raw=T)+poly(tree_area.s,2,raw=T), listw=h.sar.land.sv$weight, data=hdata.9, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.land.sv, h.best.land.sv.raw, hdata.9, 1, "clumpy", "clumpy.s", "Clumpiness", "FHD Partial Residuals") plot.partial(h.best.land.sv, h.best.land.sv.raw, hdata.9, 2, "number_patches", "number_patches.s", "Number forest fragments", "FHD Partial Residuals") plot.partial(h.best.land.sv, h.best.land.sv.raw, hdata.9, 3, "perimeter_area_mn", "perimeter_area_mn.s", "Forest perimeter:area ratio", "FHD Partial Residuals") plot.partial(h.best.land.sv, h.best.land.sv.raw, hdata.9, 4, "tree_area", "tree_area.s", "Proportion forest (%)", "FHD Partial Residuals") ####################################### ####Combined model - vegetation density FHD #### h.comb.sv.lm <- lm(sv~1, data=hdata.9) h.comb.sv.lm.1 <- lm(sv~aspect_cos.s+aspect_sin.s+elev.s+slope.s+ tot_p.s+water_cap.s+soc.s+tot_n.s+sa1_medtothinc.s+sa1_avghouse.s+ sa1_medage.s+prop_noncauc.s+mb_dweldens.s+road_density.s+lot_size.s+park_prop.s+ clumpy.s+number_patches.s+tree_area.s+perimeter_area_mn.s, data=hdata.9) h.comb.sv.lm.2 <- lm(sv~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), data=hdata.9) comp.models(h.comb.sv.lm.1, h.comb.sv.lm.2) h.sar.comb.sv <- mm.spatial.sparse(hdata.9, 150, h.comb.sv.lm.2, "MC") mm.sum.sp(h.sar.comb.sv) spatial.plot(comb.sv.lm.2, data.9, "FHD - Combined polynomial non-spatial model") spatial.plot(sar.comb.sv$mixed, data.9, "FHD - Comnbined mixed SAR model") h.best.comb.sv <- lagsarlm(sv ~ poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), listw=h.sar.comb.sv$weight, data=hdata.9, type="mixed", method="MC", na.action="na.fail") AICc(h.best.comb.sv) summary(h.best.comb.sv) diag.best(best.comb.sv, data.9, sar.comb.sv$weight, "FHD - Best Combined SAR model") h.best.comb.sv.null <- lagsarlm(sv ~ 1, listw=h.sar.comb.sv$weight, data=hdata.9, type="mixed", method="MC", na.action="na.fail") spatial.diag(h.comb.sv.lm, h.comb.sv.lm.2, h.best.comb.sv, h.best.comb.sv.null, hdata.9) h.best.comb.sv.raw <- lagsarlm(sv~poly(aspect_cos.s,2,raw=T)+poly(aspect_sin.s,2,raw=T)+poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(tot_p.s,2,raw=T)+poly(water_cap.s,2,raw=T)+poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_avghouse.s,2,raw=T)+ poly(sa1_medage.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T)+poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+poly(lot_size.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+poly(tree_area.s,2,raw=T)+poly(perimeter_area_mn.s,2,raw=T), listw=h.sar.comb.sv$weight, data=hdata.9, type="mixed", method="MC", na.action="na.fail") h.best.comb.sv.raw.avg <- dredge(h.best.comb.sv.raw, m.max=4, trace=2) h.best.comb.sv.raw.avg2 <- model.avg(h.best.comb.sv.raw.avg) summary(h.best.comb.sv.raw.avg2) h.comb.sv.avg.best <- get.models(h.best.comb.sv.raw.avg,1)[[1]] summary(get.models(h.best.comb.sv.raw.avg,1)[[1]]) ##includes tree area, number patches, park prop, slope h.best.simple.sv <- lagsarlm(sv~poly(tree_area.s,2)+poly(number_patches.s,2)+poly(park_prop.s,2)+ poly(slope.s,2), listw=h.sar.comb.sv$weight, data=hdata.9, type="mixed", na.action="na.fail", method="MC") h.best.simple.sv.raw <- lagsarlm(sv~poly(tree_area.s,2,raw=T)+poly(number_patches.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(slope.s,2,raw=T), listw=h.sar.comb.sv$weight,data=hdata.9, type="mixed", na.action="na.fail", method="MC") h.best.sv.null <- lagsarlm(sv~1, listw=h.sar.comb.sv$weight, data=hdata.9, type="mixed", na.action="na.fail", method="MC") AICc(h.best.simple.sv) summary(h.best.simple.sv) summary(h.best.simple.sv.raw) spatial.diag(h.comb.sv.lm, h.comb.sv.lm.2, h.best.simple.sv, h.best.sv.null, hdata.9) h.comb.sv.best.lm.2 <- lm(sv~poly(tree_area.s,2)+poly(number_patches.s,2)+poly(park_prop.s,2)+ poly(slope.s,2), data=hdata.9) h.sar.best.sv <- mm.spatial.sparse(hdata.1, 150, h.comb.sv.best.lm.2, "MC") mm.sum.sp(h.sar.best.sv) hdata.9$tree_area_per <- hdata.9$tree_area*100 plot.partial(h.best.simple.sv, h.best.simple.sv.raw, hdata.9, 1, "tree_area_per", "tree_area.s", "Proportion tree cover (%)", "FHD Partial Residuals") plot.partial(h.best.simple.sv, h.best.simple.sv.raw, hdata.9, 2, "number_patches", "number_patches.s", "Number of tree cover patches", "FHD Partial Residuals") plot.partial(h.best.simple.sv, h.best.simple.sv.raw, hdata.9, 3, "park_prop_per", "park_prop.s", "Proportion parkland (%)", "FHD Partial Residuals") plot.partial(h.best.simple.sv, h.best.simple.sv.raw, hdata.9, 4, "slope", "slope.s", "Slope (degrees)", "FHD Partial Residuals") ##################################################################### #####Model selection for Canopy Height ##################################################################### ####################################### ####Physical variables - Canopy Height #### h.phys.h95.lm <- lm(log(ht_p95+0.01)~1, data=hdata.1) h.phys.h95.lm.1 <- lm(log(ht_p95+0.01)~aspect_cos.s+aspect_sin.s+elev.s+slope.s, data=hdata.1) h.phys.h95.lm.2 <- lm(log(ht_p95+0.01)~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2), data=hdata.1) comp.models(h.phys.h95.lm.1, h.phys.h95.lm.2) h.sar.phys.h95 <- mm.spatial.sparse(hdata.1, 150, h.phys.h95.lm.2, "MC") mm.sum.sp(h.sar.phys.h95) spatial.plot(phys.h95.lm.2, data.1, "Canopy Height - Polynomial non-spatial model") spatial.plot(sar.phys.h95$mixed, data.1, "Canopy Height - Mixed SAR model") h.best.phys.h95 <- lagsarlm(log(ht_p95+0.01) ~ poly(elev.s,2)+poly(slope.s,2)+poly(aspect_cos.s,2)+poly(aspect_sin.s,2), listw=h.sar.phys.h95$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.phys.h95) summary(h.best.phys.h95) diag.best(best.phys.h95, data.1, sar.phys.h95$weight, "Canopy Height - Best SAR model") h.best.phys.h95.null <- lagsarlm(log(ht_p95+0.01) ~ 1, listw=h.sar.phys.h95$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.phys.h95.lm, h.phys.h95.lm.2, h.best.phys.h95, h.best.phys.h95.null, hdata.1) h.best.phys.h95.raw <- lagsarlm(log(ht_p95+0.01) ~ poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(aspect_cos.s,2,raw=T)+poly(aspect_sin,2,raw=T), listw=h.sar.phys.h95$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") ##Partial residual plotting to look at trends plot.partial(h.best.phys.h95, h.best.phys.h95.raw, hdata.1, 3, "aspect_cos", "aspect_cos.s", "Northness", "Canopy Height Partial Residuals") plot.partial(h.best.phys.h95, h.best.phys.h95.raw, hdata.1, 4, "aspect_sin", "aspect_sin.s", "Eastness", "Canopy Height Partial Residuals") plot.partial(h.best.phys.h95, h.best.phys.h95.raw, hdata.1, 1, "elev", "elev.s", "Elevation (m)", "Canopy Height Partial Residuals") plot.partial(h.best.phys.h95, h.best.phys.h95.raw, hdata.1, 2, "slope", "slope.s", "Slope (degrees)", "Canopy Height Partial Residuals") ####################################### ####Soil variables - Canopy Height #### h.soil.h95.lm <- lm(log(ht_p95+0.01)~1, data=hdata.1) h.soil.h95.lm.1 <- lm(log(ht_p95+0.01)~water_cap.s+tot_p.s+soc.s+tot_n.s, data=hdata.1) h.soil.h95.lm.2 <- lm(log(ht_p95+0.01)~poly(water_cap.s,2)+poly(tot_p.s,2)+poly(soc.s,2)+poly(tot_n.s,2), data=hdata.1) comp.models(h.soil.h95.lm.1, h.soil.h95.lm.2) h.sar.soil.h95 <- mm.spatial.sparse(hdata.1, 150, h.soil.h95.lm.2, "MC") mm.sum.sp(h.sar.soil.h95) spatial.plot(soil.h95.lm.2, data.1, "Canopy height - Soil polynomial non-spatial model") spatial.plot(sar.soil.h95$error, data.1, "Canopy height - Soil mixed SAR model") h.best.soil.h95 <- lagsarlm(log(ht_p95+0.01) ~ poly(soc.s,2)+poly(tot_n.s,2)+poly(tot_p.s,2)+poly(water_cap.s,2), listw=h.sar.soil.h95$weight, data=hdata.1, type="lag", method="MC", na.action="na.fail") AICc(h.best.soil.h95) summary(h.best.soil.h95) diag.best(best.soil.h95, data.1, sar.soil.h95$weight, "Canopy Height - Best Soil SAR model") h.best.soil.h95.null <- lagsarlm(log(ht_p95+0.01) ~ 1, listw=h.sar.soil.h95$weight, data=hdata.1, type="lag", method="MC", na.action="na.fail") spatial.diag(h.soil.h95.lm, h.soil.h95.lm.2, h.best.soil.h95, h.best.soil.h95.null, hdata.1) h.best.soil.h95.raw <- lagsarlm(log(ht_p95+0.01) ~ poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+ poly(water_cap.s,2,raw=T)+poly(tot_p.s,2,raw=T), listw=h.sar.soil.h95$weight, data=hdata.1, type="lag", method="MC", na.action="na.fail") ##Partial residual plotting to look at trends plot.partial(h.best.soil.h95, h.best.soil.h95.raw, hdata.1, 1, "soc", "soc.s", "Soil organic carbon (%)", "Canopy Height Partial Residuals") plot.partial(h.best.soil.h95, h.best.soil.h95.raw, hdata.1, 2, "tot_n", "tot_n.s", "Total soil N (%)", "Canopy Height Partial Residuals") plot.partial(h.best.soil.h95, h.best.soil.h95.raw, hdata.1, 3, "water_cap", "water_cap.s", "Soil water capacity (%)", "Canopy Height Partial Residuals") plot.partial(h.best.soil.h95, h.best.soil.h95.raw, hdata.1, 4, "tot_p", "tot_p.s", "Total soil P (%)", "Canopy Height Partial Residuals") ####################################### ####Census variables - Canopy Height #### h.cens.h95.lm <- lm(log(ht_p95+0.01)~1, data=hdata.1) h.cens.h95.lm.1 <- lm(log(ht_p95+0.01)~prop_noncauc.s+sa1_medtothinc.s+sa1_medage.s+sa1_avghouse.s, data=hdata.1) h.cens.h95.lm.2 <- lm(log(ht_p95+0.01)~poly(prop_noncauc.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2), data=hdata.1) comp.models(h.cens.h95.lm.1, h.cens.h95.lm.2) h.sar.cens.h95 <- mm.spatial.sparse(hdata.1, 150, h.cens.h95.lm.2, "MC") mm.sum.sp(h.sar.cens.h95) spatial.plot(cens.h95.lm.2, data.1, "Canopy Height - Census polynomial non-spatial model") spatial.plot(sar.cens.h95$mixed, data.1, "Canopy Height - Census mixed SAR model") h.best.cens.h95 <- lagsarlm(log(ht_p95+0.01) ~ poly(prop_noncauc.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_medage.s,2)+poly(sa1_avghouse.s,2), listw=h.sar.cens.h95$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.cens.h95) summary(h.best.cens.h95) diag.best(best.cens.h95, data.1, sar.cens.h95$weight, "Canopy Height - Best Soil SAR model") h.best.cens.h95.null <- lagsarlm(log(ht_p95+0.01) ~ 1, listw=h.sar.cens.h95$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") spatial.diag(h.cens.h95.lm, h.cens.h95.lm.2, h.best.cens.h95, h.best.cens.h95.null, hdata.1) h.best.cens.h95.raw <- lagsarlm(log(ht_p95+0.01)~poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_medage.s,2,raw=T)+ poly(sa1_avghouse.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T), listw=h.sar.cens.h95$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") plot.partial(h.best.cens.h95, h.best.cens.h95.raw, hdata.1, 2, "sa1_medage", "sa1_medage.s", "Median pop. age", "Canopy Height Partial Residuals") plot.partial(h.best.cens.h95, h.best.cens.h95.raw, hdata.1, 1, "sa1_medtothinc", "sa1_medtothinc.s", "Median weekly household income ($)", "Canopy Height Partial Residuals") plot.partial(h.best.cens.h95, h.best.cens.h95.raw, hdata.1, 4, "prop_noncauc", "prop_noncauc.s", "Proportion Non-European Descent", "Canopy Height Partial Residuals") plot.partial(h.best.cens.h95, h.best.cens.h95.raw, hdata.1, 3, "sa1_avghouse", "sa1_avghouse.s", "Avg. household size", "Canopy Height Partial Residuals") ####################################### ####Urban form variables - Canopy Height #### h.urb.h95.lm <- lm(log(ht_p95+0.01)~1, data=hdata.1) h.urb.h95.lm.1 <- lm(log(ht_p95+0.01)~mb_dweldens.s+road_density.s+park_prop.s+lot_size.s, data=hdata.1) h.urb.h95.lm.2 <- lm(log(ht_p95+0.01)~poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), data=hdata.1) comp.models(h.urb.h95.lm.1, h.urb.h95.lm.2) h.sar.urb.h95 <- mm.spatial.sparse(hdata.1, 150, h.urb.h95.lm.2, "MC") mm.sum.sp(h.sar.urb.h95) spatial.plot(urb.h95.lm.2, data.1, "FHD - Urban form polynomial non-spatial model") spatial.plot(sar.cens.h95$mixed, data.1, "FHD - Urban form mixed SAR model") h.best.urb.h95 <- lagsarlm(log(ht_p95+0.01) ~ poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(park_prop.s,2)+poly(lot_size.s,2), listw=h.sar.urb.h95$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.urb.h95) summary(h.best.urb.h95) diag.best(best.urb.h95, data.1, sar.urb.h95$weight, "Canopy Height - Best Urban Form SAR model") h.best.urb.h95.null <- lagsarlm(log(ht_p95+0.01) ~ 1, listw=h.sar.urb.h95$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.urb.h95.lm, h.urb.h95.lm.2, h.best.urb.h95, h.best.urb.h95.null, hdata.1) h.best.urb.h95.raw <- lagsarlm(log(ht_p95+0.01)~poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+ poly(park_prop.s,2,raw=T)+poly(lot_size.s,2,raw=T), listw=h.sar.urb.h95$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") plot.partial(h.best.urb.h95, h.best.urb.h95.raw, hdata.1, 2, "road_density", "road_density.s", "Street density (km/km2)", "Canopy Height Partial Residuals") plot.partial(h.best.urb.h95, h.best.urb.h95.raw, hdata.1, 3, "park_prop_per", "park_prop.s", "Proportion parkland (%)", "Canopy Height Partial Residuals") plot.partial(h.best.urb.h95, h.best.urb.h95.raw, hdata.1, 1, "mb_dwel_dens", "mb_dweldens.s", "Dwelling density (dwellings/ha)", "Canopy Height Partial Residuals") plot.partial(h.best.urb.h95, h.best.urb.h95.raw, hdata.1, 4, "log_lot_size", "log_lot_size.s", "Mean lot size (m2)", "Canopy Height Partial Residuals") ####################################### ####Landscape variables - Canopy Height #### h.land.h95.lm <- lm(log(ht_p95+0.01)~1, data=hdata.1) h.land.h95.lm.1 <- lm(log(ht_p95+0.01)~clumpy.s+number_patches.s+perimeter_area_mn.s+tree_area.s, data=hdata.1) h.land.h95.lm.2 <- lm(log(ht_p95+0.01)~poly(clumpy.s,2)+poly(number_patches.s,2)+poly(perimeter_area_mn.s,2)+poly(tree_area.s,2), data=hdata.1) comp.models(h.land.h95.lm.1, h.land.h95.lm.2) h.sar.land.h95 <- mm.spatial.sparse(hdata.1, 150, h.land.h95.lm.2, "MC") mm.sum.sp(h.sar.land.h95) spatial.plot(land.h95.lm.2, data.1, "Canopy Height - Landscape polynomial non-spatial model") spatial.plot(sar.land.h95$mixed, data.1, "Canopy Height - Landscape error SAR model") h.best.land.h95 <- lagsarlm(log(ht_p95+0.01) ~ poly(clumpy.s,2)+poly(number_patches.s,2)+ poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), listw=h.sar.land.h95$weight, type="mixed", method="MC", data=hdata.1, na.action="na.fail") AICc(h.best.land.h95) summary(h.best.land.h95) diag.best(best.land.h95, data.1, sar.land.h95$weight, "Canopy Height - Best Landscape SAR model") h.best.land.h95.null <- lagsarlm(log(ht_p95+0.01) ~ 1, listw=h.sar.land.h95$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") spatial.diag(h.land.h95.lm, h.land.h95.lm.2, h.best.land.h95, h.best.land.h95.null, hdata.1) h.best.land.h95.raw <- lagsarlm(log(ht_p95+0.01)~poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+ poly(perimeter_area_mn.s,2,raw=T)+poly(tree_area.s,2,raw=T), listw=h.sar.land.h95$weight, type="mixed", data=hdata.1, method="MC", na.action="na.fail") plot.partial(h.best.land.h95, h.best.land.h95.raw, hdata.1, 1, "clumpy", "clumpy.s", "Clumpiness", "Canopy Height Partial Residuals") plot.partial(h.best.land.h95, h.best.land.h95.raw, hdata.1, 2, "number_patches", "number_patches.s", "Number forest fragments", "Canopy Height Partial Residuals") plot.partial(h.best.land.h95, h.best.land.h95.raw, hdata.1, 3, "perimeter_area_mn", "perimeter_area_mn.s", "Forest perimeter:area ratio", "Canopy Height Partial Residuals") plot.partial(h.best.land.h95, h.best.land.h95.raw, hdata.1, 4, "tree_area", "tree_area.s", "Proportion forest (%)", "Canopy Height Partial Residuals") ####################################### ####Combined model - Canopy Height #### h.comb.h95.lm <- lm(log(ht_p95+0.01)~1, data=hdata.1) h.comb.h95.lm.1 <- lm(log(ht_p95+0.01)~aspect_cos.s+aspect_sin.s+elev.s+slope.s+ tot_p.s+water_cap.s+soc.s+tot_n.s+sa1_medtothinc.s+sa1_avghouse.s+ sa1_medage.s+prop_noncauc.s+mb_dweldens.s+road_density.s+lot_size.s+park_prop.s+ clumpy.s+number_patches.s+tree_area.s+perimeter_area_mn.s, data=hdata.1) h.comb.h95.lm.2 <- lm(log(ht_p95+0.01)~poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), data=hdata.1) comp.models(h.comb.h95.lm.1, h.comb.h95.lm.2) h.sar.comb.h95 <- mm.spatial.sparse(hdata.1, 150, h.comb.h95.lm.2, "MC") mm.sum.sp(h.sar.comb.h95) spatial.plot(comb.h95.lm.2, data.1, "Canopy Height - Combined polynomial non-spatial model") spatial.plot(sar.comb.h95$mixed, data.1, "Canopy Height - Comnbined mixed SAR model") h.best.comb.h95 <- lagsarlm(log(ht_p95+0.01) ~ poly(aspect_cos.s,2)+poly(aspect_sin.s,2)+poly(elev.s,2)+poly(slope.s,2)+ poly(tot_p.s,2)+poly(water_cap.s,2)+poly(soc.s,2)+poly(tot_n.s,2)+poly(sa1_medtothinc.s,2)+poly(sa1_avghouse.s,2)+ poly(sa1_medage.s,2)+poly(prop_noncauc.s,2)+poly(mb_dweldens.s,2)+poly(road_density.s,2)+poly(lot_size.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2)+poly(number_patches.s,2)+poly(tree_area.s,2)+poly(perimeter_area_mn.s,2), listw=h.sar.comb.h95$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") AICc(h.best.comb.h95) summary(h.best.comb.h95) diag.best(best.comb.h95, data.1, sar.comb.h95$weight, "FHD - Best Combined SAR model") h.best.comb.h95.null <- lagsarlm(log(ht_p95+0.01) ~ 1, listw=h.sar.comb.h95$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") spatial.diag(h.comb.h95.lm, h.comb.h95.lm.2, h.best.comb.h95, h.best.comb.h95.null, hdata.1) h.best.comb.h95.raw <- lagsarlm(log(ht_p95+0.01)~poly(aspect_cos.s,2,raw=T)+poly(aspect_sin.s,2,raw=T)+poly(elev.s,2,raw=T)+poly(slope.s,2,raw=T)+ poly(tot_p.s,2,raw=T)+poly(water_cap.s,2,raw=T)+poly(soc.s,2,raw=T)+poly(tot_n.s,2,raw=T)+poly(sa1_medtothinc.s,2,raw=T)+poly(sa1_avghouse.s,2,raw=T)+ poly(sa1_medage.s,2,raw=T)+poly(prop_noncauc.s,2,raw=T)+poly(mb_dweldens.s,2,raw=T)+poly(road_density.s,2,raw=T)+poly(lot_size.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(clumpy.s,2,raw=T)+poly(number_patches.s,2,raw=T)+poly(tree_area.s,2,raw=T)+poly(perimeter_area_mn.s,2,raw=T), listw=h.sar.comb.h95$weight, data=hdata.1, type="mixed", method="MC", na.action="na.fail") h.best.comb.h95.raw.avg <- dredge(h.best.comb.h95.raw, m.max=4, trace=2) h.best.comb.h95.raw.avg2 <- model.avg(h.best.comb.h95.raw.avg) summary(h.best.comb.h95.raw.avg2) h.comb.h95.avg.best <- get.models(h.best.comb.h95.raw.avg,1)[[1]] summary(get.models(h.best.comb.h95.raw.avg,1)[[1]]) ##includes tree area, clumpy, number patches, park prop h.best.simple.h95 <- lagsarlm(log(ht_p95+0.01)~poly(tree_area.s,2)+poly(number_patches.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2), listw=h.sar.comb.h95$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.simple.h95.raw <- lagsarlm(log(ht_p95+0.01)~poly(tree_area.s,2,raw=T)+poly(number_patches.s,2,raw=T)+poly(park_prop.s,2,raw=T)+ poly(clumpy.s,2,raw=T), listw=h.sar.comb.h95$weight,data=hdata.1, type="mixed", na.action="na.fail", method="MC") h.best.h95.null <- lagsarlm(log(ht_p95+0.01)~1, listw=h.sar.comb.h95$weight, data=hdata.1, type="mixed", na.action="na.fail", method="MC") AICc(h.best.simple.h95) summary(h.best.simple.h95) summary(h.best.simple.h95.raw) spatial.diag(h.comb.h95.lm, h.comb.h95.lm.2, h.best.simple.h95, h.best.h95.null, hdata.1) h.comb.h95.best.lm.2 <- lm(log(ht_p95+0.01)~poly(tree_area.s,2)+poly(number_patches.s,2)+poly(park_prop.s,2)+ poly(clumpy.s,2), data=hdata.1) h.sar.best.h95 <- mm.spatial.sparse(hdata.1, 150, h.comb.h95.best.lm.2, "MC") mm.sum.sp(h.sar.best.h95) plot.partial(h.best.simple.h95, h.best.simple.h95.raw, hdata.1, 1, "tree_area_per", "tree_area.s", "Proportion tree cover (%)", "Canopy Height Partial Residuals") plot.partial(h.best.simple.h95, h.best.simple.h95.raw, hdata.1, 2, "number_patches", "number_patches.s", "Number of tree cover patches", "Canopy Height Partial Residuals") plot.partial(h.best.simple.h95, h.best.simple.h95.raw, hdata.1, 3, "park_prop_per", "park_prop.s", "Proportion parkland (%)", "Canopy Height Partial Residuals") plot.partial(h.best.simple.h95, h.best.simple.h95.raw, hdata.1, 4, "clumpy", "clumpy.s", "Tree cover clumpiness", "Canopy Height Partial Residuals")