# R code for models, tables and figures of the manuscript: 
# "Climate differentiates forest structure across a residential macrosystem" 
# by Alessandro Ossola & Matthew Hopton, 
# US-EPA, NRMRL-ORD-STD, 26 W. Martin Luther King Dr., Cincinnati, OH, 45268, USA.
# Corresponding Author, Email: hopton.matthew@epa.gov, T: +1 513 569 7718
# version 11/12/2017

.libPaths("C:/users/aossola/R/library")
library(ggplot2)
library(gtable)   #for dual axes plots
library(grid)     #for dual axes plots
library(reshape)  #for pivot tables

# uploading data
data.ori<-read.csv("M:/Net MyDocuments/Products/Pubs/EcolAppl/Data.csv", sep = ",", header=T, na.string = "NA")
attach(data.ori)
names(data.ori)
dim(data.ori) 

# organizing data
# excluding census tracts with less than 100 parcels/km2 and tracts with residential cover <10%
data<-subset(data.ori, ParcDens>100 & PercResid>10 & PercResid<100)
data$Rent <- as.numeric(as.character(data$Rent))
attach(data)
#sapply(data,class)

# retrieving some general stats
dim(data)[1]           #tot n. tracts
table(data$City)       #n. tracts/city 
sum(ParcelNo)          #tot n. parcels
aggregate(ParcelNo~City, data=data, na.action = na.omit, FUN=sum)
sum(ALAND10/1000000)   #tot area tracts (land) in km2
aggregate(ALAND10/1000000~City, data=data, na.action = na.omit, FUN=sum)
sum(Area_m2/1000000)   #tot area residential parcels in km2
aggregate(Area_m2/1000000~City, data=data, na.action = na.omit, FUN=sum)
mean(ParcDens)         #average parcel density
aggregate(ParcDens~City, data=data, na.action = na.omit, FUN=mean)
mean(PercResid)        #average percent residential land cover 
aggregate(PercResid~City, data=data, na.action = na.omit, FUN=mean)
mean(MeanParArea)      #average parcel area in m2
aggregate(MeanParArea~City, data=data, na.action = na.omit, FUN=mean)
sum(Population)        #tot population
aggregate(Population~City, data=data, na.action = na.omit, FUN=sum)
mean(Population/(ALAND10/1000000))        # population density
aggregate(Population/(ALAND10/1000000)~City, data=data, na.action = na.omit, FUN=mean)

####################################################################################################################
################# GLM modeling with VIF threshold=2, R2=0.50 #######################################################

#selecting variables for multicollinearity checking
sub.vif<-na.omit(subset(data, select=c(ET, MeanParArea, ParcDens, PercResid, Decade, VacantUnits., Age, RaceDiv, 
                        RenterOcc., HighSchool, Degree, BelowPoverty, Gini, Income, Unemployed, Rent)))

# checking variable multicollinearity using Variance Inflation Factors (VIF) and excluding variables which have
# a correlation higher than ~0.5 (VIF threshold =2, Zuur et al. (2010). A protocol for data exploration to avoid 
# common statistical problems. Methods in Ecology and Evolution 1: 3-14.). 
library(usdm)
sub.vif = data.frame(sub.vif)
vif(sub.vif)
vifstep(sub.vif, th=2)

cor.test(MeanParArea, ParcDens)
cor.test(MeanParArea, PercResid)
cor.test(MeanParArea, Decade)
cor.test(Income, VacantUnits.)
cor.test(Income, Age)
cor.test(Income, RaceDiv)
cor.test(Income, Gini)
cor.test(Income, Unemployed)
cor.test(Income, Rent)

#selecting variables for modeling
sub<-na.omit(subset(data, select=c(City, TvolAREA, TcanAREA, ET, ParcDens, PercResid, Decade, VacantUnits., Age, RaceDiv,  
                   Gini, Unemployed, Rent)))

# Prior to modeling GLMs, preliminary GAMs have been fitted to check the likely type of response between outcome 
# and predictor variables (e.g. TvolAREA and TcanAREA). These analyses are not presented in this code.

library(caret)    # for variable importance varImp()

#######################modeling TvolAREA########################
glm.TvolAREA.global<-glm(log(TvolAREA+1)~poly(log(ET+1),2) + poly(log(ParcDens+1),2) + poly(PercResid,2) + 
                    I(Decade) + log(VacantUnits.+1) +  Age + RaceDiv + 
                    Gini + poly(log(Unemployed+1),2) + poly(log(Rent+1),2), data=sub, 
                    family=gaussian(link=log))
summary(glm.TvolAREA.global)
par(mfrow=c(2,2))
plot(glm.TvolAREA.global)
mod_fit.TvolAREA.global <- train(log(TvolAREA+1)~poly(log(ET+1),2) + poly(log(ParcDens+1),2) + poly(PercResid,2) + 
                    I(Decade) + log(VacantUnits.+1) +  Age + RaceDiv + 
                    Gini + poly(log(Unemployed+1),2) + poly(log(Rent+1),2), data=sub, method="glm", family=gaussian(link=log))
varImp(mod_fit.TvolAREA.global)

glm.TvolAREA.clim.morpho<-glm(log(TvolAREA+1)~poly(log(ET+1),2) + poly(log(ParcDens+1),2) + poly(PercResid,2) + 
                    I(Decade), data=sub, family=gaussian(link=log))
summary(glm.TvolAREA.clim.morpho)
par(mfrow=c(2,2))
plot(glm.TvolAREA.clim.morpho)

glm.TvolAREA.clim.socio<-glm(log(TvolAREA+1)~poly(log(ET+1),2) + log(VacantUnits.+1) +  Age + RaceDiv + 
                    Gini + poly(log(Unemployed+1),2) + poly(log(Rent+1),2), 
                    data=sub, family=gaussian(link=log))
summary(glm.TvolAREA.clim.socio)
par(mfrow=c(2,2))
plot(glm.TvolAREA.clim.socio)

glm.TvolAREA.morpho.socio<-glm(log(TvolAREA+1)~ poly(log(ParcDens+1),2) + poly(PercResid,2) + 
                    I(Decade) + log(VacantUnits.+1) +  Age + RaceDiv + 
                    Gini + poly(log(Unemployed+1),2) + poly(log(Rent+1),2), data=sub, 
                    family=gaussian(link=log))
summary(glm.TvolAREA.morpho.socio)
par(mfrow=c(2,2))
plot(glm.TvolAREA.morpho.socio)

glm.TvolAREA.clim<-glm(log(TvolAREA+1)~poly(log(ET+1),2), data=sub, 
                    family=gaussian(link=log))
summary(glm.TvolAREA.clim)
par(mfrow=c(2,2))
plot(glm.TvolAREA.clim)

glm.TvolAREA.morpho<-glm(log(TvolAREA+1)~ poly(log(ParcDens+1),2) + poly(PercResid,2) + 
                    I(Decade), data=sub, family=gaussian(link=log))
summary(glm.TvolAREA.morpho)
par(mfrow=c(2,2))
plot(glm.TvolAREA.morpho)

glm.TvolAREA.socio<-glm(log(TvolAREA+1)~log(VacantUnits.+1) +  Age + RaceDiv + 
                   Gini + poly(log(Unemployed+1),2) + poly(log(Rent+1),2), data=sub, 
                    family=gaussian(link=log))
summary(glm.TvolAREA.socio)
par(mfrow=c(2,2))
plot(glm.TvolAREA.socio)

library(AICcmodavg)
Can.mod<-list()
Can.mod[[1]]<-glm.TvolAREA.global
Can.mod[[2]]<-glm.TvolAREA.clim.morpho
Can.mod[[3]]<-glm.TvolAREA.clim.socio
Can.mod[[4]]<-glm.TvolAREA.morpho.socio
Can.mod[[5]]<-glm.TvolAREA.clim
Can.mod[[6]]<-glm.TvolAREA.morpho
Can.mod[[7]]<-glm.TvolAREA.socio
aictab(cand.set = Can.mod, sort = TRUE, second.ord=FALSE)


#######################modeling TcanAREA########################
#GLMs binomial logit
# for tree percent canopy cover logistic regression rather than linear has more predictive power (Bigsby et al., 2014),
# thus the binomial family (logit link) has been selected for fitting GLMs. 

glm.TcanAREA.global<-glm(log(TcanAREA+1)~poly(log(ET+1),2) + poly(log(ParcDens+1),2) + poly(PercResid,2) + 
                    I(Decade) + log(VacantUnits.+1) +  Age + RaceDiv + 
                    Gini + poly(log(Unemployed+1),2) + log(Rent+1), data=sub, 
                    family=binomial(link=logit))
summary(glm.TcanAREA.global)
par(mfrow=c(2,2))
plot(glm.TcanAREA.global)

glm.TcanAREA.clim.morpho<-glm(log(TcanAREA+1)~poly(log(ET+1),2) + poly(log(ParcDens+1),2) + poly(PercResid,2) + 
                    I(Decade), data=sub, family=binomial(link=logit))
summary(glm.TcanAREA.clim.morpho)
par(mfrow=c(2,2))
plot(glm.TcanAREA.clim.morpho)
mod_fit.TcanAREA.clim.morpho <- train(log(TcanAREA+1)~poly(log(ET+1),2) + poly(log(ParcDens+1),2) + poly(PercResid,2) + 
                    I(Decade), data=sub, method="glm", family=binomial(link=logit))
varImp(mod_fit.TcanAREA.clim.morpho)

glm.TcanAREA.clim.socio<-glm(log(TcanAREA+1)~poly(log(ET+1),2) + + log(VacantUnits.+1) +  Age + RaceDiv + 
                    Gini + poly(log(Unemployed+1),2) + log(Rent+1), 
                    data=sub, family=binomial(link=logit))
summary(glm.TcanAREA.clim.socio)
par(mfrow=c(2,2))
plot(glm.TcanAREA.clim.socio)

glm.TcanAREA.morpho.socio<-glm(log(TcanAREA+1)~ poly(log(ParcDens+1),2) + poly(PercResid,2) + 
                    I(Decade) + log(VacantUnits.+1) +  Age + RaceDiv + Gini + poly(log(Unemployed+1),2) + 
                    log(Rent+1), data=sub, 
                    family=binomial(link=logit))
summary(glm.TcanAREA.morpho.socio)
par(mfrow=c(2,2))
plot(glm.TcanAREA.morpho.socio)

glm.TcanAREA.clim<-glm(log(TcanAREA+1)~poly(log(ET+1),2), data=sub, 
                    family=binomial(link=logit))
summary(glm.TcanAREA.clim)
par(mfrow=c(2,2))
plot(glm.TcanAREA.clim)

glm.TcanAREA.morpho<-glm(log(TcanAREA+1)~ poly(log(ParcDens+1),2) + poly(PercResid,2) + 
                    I(Decade), data=sub, family=binomial(link=logit))
summary(glm.TcanAREA.morpho)
par(mfrow=c(2,2))
plot(glm.TcanAREA.morpho)

glm.TcanAREA.socio<-glm(log(TcanAREA+1)~log(VacantUnits.+1) +  Age + RaceDiv + Gini + poly(log(Unemployed+1),2) + 
                    log(Rent+1), data=sub, 
                    family=binomial(link=logit))
summary(glm.TcanAREA.socio)
par(mfrow=c(2,2))
plot(glm.TcanAREA.socio)

library(AICcmodavg)
Can.mod<-list()
Can.mod[[1]]<-glm.TcanAREA.global
Can.mod[[2]]<-glm.TcanAREA.clim.morpho
Can.mod[[3]]<-glm.TcanAREA.clim.socio
Can.mod[[4]]<-glm.TcanAREA.morpho.socio
Can.mod[[5]]<-glm.TcanAREA.clim
Can.mod[[6]]<-glm.TcanAREA.morpho
Can.mod[[7]]<-glm.TcanAREA.socio
aictab(cand.set = Can.mod, sort = TRUE, second.ord=FALSE)

########################### codes for figures  #################################################################
####### Figure 2 #####

dev.new(width=17, height=17)
fig2 <- ggplot(data = data, aes(y = TvolAREA, x = TcanAREA))
fig2<- fig2 + geom_point(aes(shape = factor(City), col=factor(City))) +  theme_bw() + 
  scale_shape_manual(name="     City", values = c(0,1,2,3,4,5,6,8,10))  + scale_colour_discrete(name  ="     City") +
  theme(legend.position = c(0.15, 0.7)) + #theme(legend.position="none") +
  ylab(bquote('Residential forest volume ('*~m^3~m^-2*')')) + xlim(0,1) + ylim(0,12) +
  xlab(bquote('Residential forest cover ('*m^2~m^-2*')')) + ggtitle("") +
  theme(legend.key.size = unit(0.8, "cm")) +
  stat_function(fun=function(x)0.5911*exp(3.6348*x), geom="line", size=0.8) #+  #R2 is 0.7959
  #annotate(geom="text", x=0.8, y=1, label="y = 0.5911 · exp (3.6348 · x)", color="black")
fig2

#ggsave(filename="C:/Users/aossola/Desktop/fig2.pdf", plot=fig2, dpi=300)
#ggsave(filename="C:/Users/aossola/Desktop/fig2.png", plot=fig2, dpi=300)

####### Figure 3 #####

dev.new(width=17, height=17)
fig3A <- ggplot(data = data, aes(x = ET, y = TcanAREA)) 
fig3A<- fig3A + geom_point(aes(shape = factor(City), col=factor(City)), position = position_jitter(width = 7)) + 
     scale_shape_manual(name="     City", values = c(0,1,2,3,4,5,6,8,10))  + scale_colour_discrete(name  ="     City") +
     theme_bw() + #theme(legend.position="none") +
     xlab(bquote('Mean annual ET ('*mm*')')) + ylim(0,1) + 
     ylab(bquote('Residential forest cover ('*~m^2~m^-2*')')) + ggtitle("") + 
     geom_boxplot(aes(fill = factor(City)), outlier.shape=NA, lwd=0.7, fatten=1.5) + 
     scale_fill_manual(values = alpha(c("lightgrey","lightgrey","lightgrey","lightgrey","lightgrey","lightgrey",
     "lightgrey","lightgrey","lightgrey"), 0.002), name="     City") 
fig3A

dev.new(width=17, height=17)
fig3B <- ggplot(data = data, aes(x = ET, y = TvolAREA)) 
fig3B<- fig3B + geom_point(aes(shape = factor(City), col=factor(City)), position = position_jitter(width = 7)) + 
     scale_shape_manual(name="     City", values = c(0,1,2,3,4,5,6,8,10))  + scale_colour_discrete(name  ="     City") +
     theme_bw() + #theme(legend.position="none") +
     xlab(bquote('Mean annual ET ('*mm*')')) + ylim(0,12) + 
     ylab(bquote('Residential forest volume ('*~m^3~m^-2*')')) + ggtitle("") + 
     geom_boxplot(aes(fill = factor(City)), outlier.shape=NA, lwd=0.7, fatten=1.5) + 
     scale_fill_manual(values = alpha(c("lightgrey","lightgrey","lightgrey","lightgrey","lightgrey","lightgrey",
     "lightgrey","lightgrey","lightgrey"), 0.002), name="     City") 
fig3B

library(ggplot2)
library(gridExtra)
library(grid)

grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {

  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x + theme(legend.position="none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)

  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl),
                                            legend,
                                            ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend,
                                           ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))

  grid.newpage()
  grid.draw(combined)

  # return gtable invisibly
  invisible(combined)

}

dev.new(width=20, height=23)
fig3<- grid_arrange_shared_legend(fig3A, fig3B, ncol = 1, nrow = 2, position="right")
fig3

ggsave(filename="C:/Users/MQ20174608/Desktop/fig3.pdf", width=8, height=10, plot=fig3, dpi=300)
ggsave(filename="C:/Users/MQ20174608/Desktop/fig3.png", width=8, height=10, plot=fig3, dpi=300)

# Dunn test post Kruskal to assign letters to similar mean values
library(PMCMR)
can<-posthoc.kruskal.dunn.test(x=TcanAREA, g=ET, p.adjust.method="bonferroni")
can = can$p.value
library(rcompanion)
can1 = fullPTable(can)
can1
library(multcompView)
multcompLetters(can1, 
                compare="<", 
                threshold=0.05,
                Letters=letters,
                reversed = FALSE)

vol<-posthoc.kruskal.dunn.test(x=TvolAREA, g=ET, p.adjust.method="bonferroni")
vol = vol$p.value
library(rcompanion)
vol1 = fullPTable(vol)
vol1
library(multcompView)
multcompLetters(vol1, 
                compare="<", 
                threshold=0.05,
                Letters=letters,
                reversed = FALSE)


####### Figure 4 #####
library(ggplot2)
library(gridExtra)
library(grid)

dev.new(width=13, height=13)
fig4A <- ggplot(data = data, aes(x = ParcDens, y = TcanAREA))
fig4A<- fig4A + geom_point(aes(shape = factor(City), col=factor(City))) +  theme_bw() + 
     scale_shape_manual(name="     City", values = c(0,1,2,3,4,5,6,8,10))  + scale_colour_discrete(name  ="     City") +
     theme(legend.position = c(0.7, 0.7)) + #theme(legend.position="none") +
     xlab(bquote('Residential parcel density ('*units ~km^-2*')')) + 
     ylim(0,1) + xlim(0,3600) + theme(legend.key.size = unit(0.4, "cm")) +
     ylab(bquote('Residential forest cover ('*~m^2~m^-2*')')) + 
     ggtitle("") + 
     theme(axis.title.x=element_blank())
fig4A

dev.new(width=13, height=13)
fig4B <- ggplot(data = data, aes(x = PercResid, y = TcanAREA))
fig4B<- fig4B + geom_point(aes(shape = factor(City), col=factor(City))) +  theme_bw() + 
     scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme(legend.position="none") +
     #xlab(bquote('Residential land cover (%)')) + 
     ylim(0,1) + xlim(0,100) +
     #ylab(bquote('Residential forest cover ('*~m^2~m^-2*')')) + 
     ggtitle("") + 
     theme(axis.title.x=element_blank(), axis.title.y=element_blank())
fig4B

dev.new(width=13, height=13)
fig4C <- ggplot(data = data, aes(x = Decade, y = TcanAREA)) 
fig4C<- fig4C + geom_point(aes(shape = factor(City), col=factor(City)), position = position_jitter(width = 2)) + 
     scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme_bw() + #theme(legend.position="none") +
     #xlab(bquote('Maximum housing development (decade)')) + 
     ylim(0,1) + xlim(1925,2010) +
     #ylab(bquote('Residential forest cover ('*~m^2~m^-2*')')) + 
     geom_boxplot(aes(fill = factor(Decade)), outlier.shape=NA, lwd=0.7, fatten=1.5) + 
     scale_fill_manual(values = alpha(c("lightgrey","lightgrey","lightgrey","lightgrey","lightgrey","lightgrey",
     "lightgrey","lightgrey","lightgrey"), 0.002)) + theme(legend.position="none") +
     ggtitle("") + 
     theme(axis.title.x=element_blank(), axis.title.y=element_blank())
fig4C

dev.new(width=13, height=13)
fig4D <- ggplot(data = data, aes(x = ParcDens, y = TvolAREA))
fig4D<- fig4D + geom_point(aes(shape = factor(City), col=factor(City))) +  theme_bw() + 
     scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme(legend.position="none") +
     xlab(bquote('Residential parcel density ('*units ~km^-2*')')) + ylim(0,12) + xlim(0,3600) +
     ylab(bquote('Residential forest volume ('*~m^3~m^-2*')')) + ggtitle("") 
fig4D

dev.new(width=13, height=13)
fig4E <- ggplot(data = data, aes(x = PercResid, y = TvolAREA))
fig4E<- fig4E + geom_point(aes(shape = factor(City), col=factor(City))) +  theme_bw() + 
     scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme(legend.position="none") +
     xlab(bquote('Residential land cover (%)')) + 
     ylim(0,12) + xlim(0,100) +
     #ylab(bquote('Residential forest volume ('*~m^3~m^-2*')')) + 
     ggtitle("") + 
     theme(axis.title.y=element_blank())
fig4E

dev.new(width=13, height=13)
fig4F <- ggplot(data = data, aes(x = Decade, y = TvolAREA)) 
fig4F<- fig4F + geom_point(aes(shape = factor(City), col=factor(City)), position = position_jitter(width = 2)) + 
     scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme_bw() + #theme(legend.position="none") +
     xlab(bquote('Maximum housing development (decade)')) + ylim(0,12) + xlim(1925,2010) +
     #ylab(bquote('Residential forest volume ('*~m^3~m^-2*')')) + 
     geom_boxplot(aes(fill = factor(Decade)), outlier.shape=NA, lwd=0.7, fatten=1.5) + 
     scale_fill_manual(values = alpha(c("lightgrey","lightgrey","lightgrey","lightgrey","lightgrey","lightgrey",
     "lightgrey","lightgrey","lightgrey"), 0.002)) + theme(legend.position="none") +
     ggtitle("") + 
     theme(axis.title.y=element_blank())
fig4F

#compose fig4
dev.new(width=150, height=120)
fig4<- grid.arrange(fig4A, fig4B, fig4C, fig4D, fig4E, fig4F, ncol=3)
fig4

ggsave(filename="C:/Users/aossola/Desktop/fig4.pdf", width=10, height=8, plot=fig4, dpi=300)
ggsave(filename="C:/Users/aossola/Desktop/fig4.png", width=10, height=8, plot=fig4, dpi=300)

#get color palette
#ggplot_build(fig4F)$data


####### Figure 5 and Appendix E #####
library(devtools)
source_gist("524eade46135f6348140")
# from https://gist.github.com/kdauria/524eade46135f6348140
stat_smooth_func <- function(mapping = NULL, data = NULL,
                        geom = "smooth", position = "identity",
                        ...,
                        method = "auto",
                        formula = y ~ x,
                        se = TRUE,
                        n = 80,
                        span = 0.75,
                        fullrange = FALSE,
                        level = 0.95,
                        method.args = list(),
                        na.rm = FALSE,
                        show.legend = NA,
                        inherit.aes = TRUE,
                        xpos = NULL,
                        ypos = NULL) {
  layer(
    data = data,
    mapping = mapping,
    stat = StatSmoothFunc,
    geom = geom,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list(
      method = method,
      formula = formula,
      se = se,
      n = n,
      fullrange = fullrange,
      level = level,
      na.rm = na.rm,
      method.args = method.args,
      span = span,
      xpos = xpos,
      ypos = ypos,
      ...
    )
  )
}


StatSmoothFunc <- ggproto("StatSmooth", Stat,
                      
                      setup_params = function(data, params) {
                        # Figure out what type of smoothing to do: loess for small datasets,
                        # gam with a cubic regression basis for large data
                        # This is based on the size of the _largest_ group.
                        if (identical(params$method, "auto")) {
                          max_group <- max(table(data$group))
                          
                          if (max_group < 1000) {
                            params$method <- "loess"
                          } else {
                            params$method <- "gam"
                            params$formula <- y ~ s(x, bs = "cs")
                          }
                        }
                        if (identical(params$method, "gam")) {
                          params$method <- mgcv::gam
                        }
                        
                        params
                      },
                      
                      compute_group = function(data, scales, method = "auto", formula = y~x,
                                               se = TRUE, n = 80, span = 0.75, fullrange = FALSE,
                                               xseq = NULL, level = 0.95, method.args = list(),
                                               na.rm = FALSE, xpos=NULL, ypos=NULL) {
                        if (length(unique(data$x)) < 2) {
                          # Not enough data to perform fit
                          return(data.frame())
                        }
                        
                        if (is.null(data$weight)) data$weight <- 1
                        
                        if (is.null(xseq)) {
                          if (is.integer(data$x)) {
                            if (fullrange) {
                              xseq <- scales$x$dimension()
                            } else {
                              xseq <- sort(unique(data$x))
                            }
                          } else {
                            if (fullrange) {
                              range <- scales$x$dimension()
                            } else {
                              range <- range(data$x, na.rm = TRUE)
                            }
                            xseq <- seq(range[1], range[2], length.out = n)
                          }
                        }
                        # Special case span because it's the most commonly used model argument
                        if (identical(method, "loess")) {
                          method.args$span <- span
                        }
                        
                        if (is.character(method)) method <- match.fun(method)
                        
                        base.args <- list(quote(formula), data = quote(data), weights = quote(weight))
                        model <- do.call(method, c(base.args, method.args))
                        
                        m = model
                        eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
                                         list(a = format(coef(m)[1], digits = 3), 
                                              b = format(coef(m)[2], digits = 3), 
                                              r2 = format(summary(m)$r.squared, digits = 3)))
                        func_string = as.character(as.expression(eq))
                        
                        if(is.null(xpos)) xpos = min(data$x)*0.9
                        if(is.null(ypos)) ypos = max(data$y)*0.9
                        data.frame(x=xpos, y=ypos, label=func_string)
                        
                      },
                      
                      required_aes = c("x", "y")
)

#order cities by increasing ET
data$city = factor(data$City, levels=c('Fargo-Moorhead','Milwaukee','Newark','Boston', 'Washington', 'Norfolk', 
                   'Raleigh-Durham', 'Birmingham', 'Tallahassee'))

dev.new(width=23, height=20)
fig5<- ggplot(data = data, aes(x = Income/1000, y = TvolAREA)) +
  stat_smooth_func(geom="text",method="lm",size=3, xpos=120, ypos=-1, parse=TRUE) +
  geom_smooth(method="lm",se=FALSE) +
  geom_point(aes(shape = factor(city), col=factor(city))) +  theme_bw() + 
  scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme(legend.position="none") +
  xlab(bquote('Median annual income (k$)')) + ylim(-2,12) + xlim(0,250) +
  ylab(bquote('Residential forest volume ('*~m^3~m^-2*')')) + ggtitle("") + facet_wrap( ~ city, ncol=3)
fig5
ggsave(filename="C:/Users/aossola/Desktop/fig5.pdf", plot=fig5, dpi=300)
ggsave(filename="C:/Users/aossola/Desktop/fig5.png", plot=fig5, dpi=300)

dev.new(width=23, height=20)
AppendixE<- ggplot(data = data, aes(x = Income/1000, y = TcanAREA)) +
  stat_smooth_func(geom="text",method="lm",size=3, xpos=120, ypos=-0.05, parse=TRUE) +
  geom_smooth(method="lm",se=FALSE) +
  geom_point(aes(shape = factor(city), col=factor(city))) +  theme_bw() + 
  scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme(legend.position="none") +
  xlab(bquote('Median annual income (k$)')) + ylim(-0.1,1) + xlim(0,250) +
  ylab(bquote('Residential forest cover ('*~m^2~m^-2*')')) + ggtitle("") + facet_wrap( ~ city, ncol=3)
AppendixE
ggsave(filename="C:/Users/aossola/Desktop/AppendixE.pdf", plot=AppendixE, dpi=300)
ggsave(filename="C:/Users/aossola/Desktop/AppendixE.png", plot=AppendixE, dpi=300)


# extra plot TvolAREA~Rent
dev.new(width=23, height=20)
TvolAREARent<- ggplot(data = data, aes(x = Rent, y = TvolAREA)) +
  stat_smooth_func(geom="text",method="lm",size=3, xpos=1000, ypos=-1, parse=TRUE) +
  geom_smooth(method="lm",se=FALSE) +
  geom_point(aes(shape = factor(city), col=factor(city))) +  theme_bw() + 
  scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme(legend.position="none") +
  xlab(bquote('Median monthly rent ($)')) + ylim(-2,12) + xlim(0,2000) +
  ylab(bquote('Residential forest volume ('*~m^3~m^-2*')')) + ggtitle("") + facet_wrap( ~ city, ncol=3)
TvolAREARent
ggsave(filename="C:/Users/aossola/Desktop/TvolAREARent.pdf", plot=TvolAREARent, dpi=300)
ggsave(filename="C:/Users/aossola/Desktop/TvolAREARent.png", plot=TvolAREARent, dpi=300)

# extra plot TcanAREA~Rent
dev.new(width=23, height=20)
TcanAREARent<- ggplot(data = data, aes(x = Rent, y = TcanAREA)) +
  stat_smooth_func(geom="text",method="lm",size=3, xpos=1000, ypos=-0.05, parse=TRUE) +
  geom_smooth(method="lm",se=FALSE) +
  geom_point(aes(shape = factor(city), col=factor(city))) +  theme_bw() + 
  scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme(legend.position="none") +
  xlab(bquote('Median monthly rent ($)')) + ylim(-0.1,1) + xlim(0,2000) +
  ylab(bquote('Residential forest cover ('*~m^2~m^-2*')')) + ggtitle("") + facet_wrap( ~ city, ncol=3)
TcanAREARent
ggsave(filename="C:/Users/aossola/Desktop/TcanAREARent.pdf", plot=TcanAREARent, dpi=300)
ggsave(filename="C:/Users/aossola/Desktop/TcanAREARent.png", plot=TcanAREARent, dpi=300)
 

####
IncomeNorm<-vector("numeric", length = dim(data)[1])
cbind(data, IncomeNorm)
data$IncomeNorm[City=="Fargo-Moorhead"] <- c(Income[City=="Fargo-Moorhead"] /median(Income[City=="Fargo-Moorhead"]))
data$IncomeNorm[City=="Milwaukee"] <- c(Income[City=="Milwaukee"] /median(Income[City=="Milwaukee"]))
data$IncomeNorm[City=="Boston"] <- c(Income[City=="Boston"] /median(Income[City=="Boston"]))
data$IncomeNorm[City=="Newark"] <- c(Income[City=="Newark"] /median(Income[City=="Newark"]))
data$IncomeNorm[City=="Washington"] <- c(Income[City=="Washington"] /median(Income[City=="Washington"]))
data$IncomeNorm[City=="Norfolk"] <- c(Income[City=="Norfolk"] /median(Income[City=="Norfolk"]))
data$IncomeNorm[City=="Raleigh-Durham"] <- c(Income[City=="Raleigh-Durham"] /median(Income[City=="Raleigh-Durham"]))
data$IncomeNorm[City=="Birmingham"] <- c(Income[City=="Birmingham"] /median(Income[City=="Birmingham"]))
data$IncomeNorm[City=="Tallahassee"] <- c(Income[City=="Tallahassee"] /median(Income[City=="Tallahassee"]))

dev.new(width=23, height=20)
AppendixE<- ggplot(data = data, aes(x = IncomeNorm, y = TcanAREA)) +
  stat_smooth_func(geom="text",method="lm",size=3, xpos=120, ypos=-0.05, parse=TRUE) +
  geom_smooth(method="lm",se=FALSE) +
  geom_point(aes(shape = factor(city), col=factor(city))) +  theme_bw() + 
  scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme(legend.position="none") +
  xlab(bquote('Median annual income normalized by median')) + ylim(-0.1,1) + xlim(0,5) +
  ylab(bquote('Residential forest cover ('*~m^2~m^-2*')')) + ggtitle("") + facet_wrap( ~ city, ncol=3)
AppendixE

####
IncomeScaled<-vector("numeric", length = dim(data)[1])
cbind(data, IncomeScaled)
data$IncomeScaled[City=="Fargo-Moorhead"] <- c(scale(Income[City=="Fargo-Moorhead"]))
data$IncomeScaled[City=="Milwaukee"] <- c(scale(Income[City=="Milwaukee"]))
data$IncomeScaled[City=="Boston"] <- c(scale(Income[City=="Boston"]))
data$IncomeScaled[City=="Newark"] <- c(scale(Income[City=="Newark"]))
data$IncomeScaled[City=="Washington"] <- c(scale(Income[City=="Washington"]))
data$IncomeScaled[City=="Norfolk"] <- c(scale(Income[City=="Norfolk"]))
data$IncomeScaled[City=="Raleigh-Durham"] <- c(scale(Income[City=="Raleigh-Durham"]))
data$IncomeScaled[City=="Birmingham"] <- c(scale(Income[City=="Birmingham"]))
data$IncomeScaled[City=="Tallahassee"] <- c(scale(Income[City=="Tallahassee"]))

dev.new(width=23, height=20)
AppendixE<- ggplot(data = data, aes(x = IncomeScaled, y = TcanAREA)) +
  stat_smooth_func(geom="text",method="lm",size=3, xpos=120, ypos=-0.05, parse=TRUE) +
  geom_smooth(method="lm",se=FALSE) +
  geom_point(aes(shape = factor(city), col=factor(city))) +  theme_bw() + 
  scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme(legend.position="none") +
  xlab(bquote('Median annual income (scaled)')) + ylim(-0.1,1) + xlim(-1,1) +
  ylab(bquote('Residential forest cover ('*~m^2~m^-2*')')) + ggtitle("") + facet_wrap( ~ city, ncol=3)
AppendixE


####### Figure 6A #####
fit <- nls(data=data, formula = TcanAREA ~ p1*log(MeanParArea) + p2, start = list(p1 = 0.5, p2 = 3000))
summary(fit) #pseudo-r2=0.5876

dev.new(width=17, height=17)
fig6A <- ggplot(data = data, aes(y = TcanAREA, x = MeanParArea))
fig6A<- fig6A + geom_point(aes(shape = factor(City), col=factor(City))) +  theme_bw() + 
  scale_shape_manual(name="     City", values = c(0,1,2,3,4,5,6,8,10))  + scale_colour_discrete(name  ="     City") +
  theme(legend.position = c(0.87, 0.3)) + #theme(legend.position="none") +
  xlab(bquote('Mean residential parcel size ('*m^2*')')) + ylim(0,1) + xlim(0,5000) +
  ylab(bquote('Residential forest cover ('*m^2~m^-2*')')) + ggtitle("") +
  theme(legend.key.size = unit(0.7, "cm")) +
  stat_function(fun=function(x)0.207126*log(x)-0.959923, geom="line", size=0.8, na.rm = TRUE) #+  
  #annotate(geom="text", x=0.8, y=1, label="y = 0.5911 · exp (3.6348 · x)", color="black")
fig6A

ggsave(filename="C:/Users/aossola/Desktop/fig6A.pdf", plot=fig6A, dpi=300)
ggsave(filename="C:/Users/aossola/Desktop/fig6A.png", plot=fig6A, dpi=300)

####### Figure 6B #####

dev.new(width=13, height=13)
fig6B <- ggplot(data = data, aes(x = Decade, y = MeanParArea)) 
fig6B<- fig6B + geom_point(aes(shape = factor(City), col=factor(City)), position = position_jitter(width = 2)) + 
     scale_shape_manual(values = c(0,1,2,3,4,5,6,8,10)) + theme_bw() + #theme(legend.position="none") +
     xlab(bquote('Maximum housing development (decade)')) + ylim(0,5000) + xlim(1925,2010) +
     ylab(bquote('Mean residential parcel size ('*m^2*')')) + 
     geom_boxplot(aes(fill = factor(Decade)), outlier.shape=NA, lwd=0.7, fatten=1.5) + 
     scale_fill_manual(values = alpha(c("lightgrey","lightgrey","lightgrey","lightgrey","lightgrey","lightgrey",
     "lightgrey","lightgrey","lightgrey"), 0.002)) + theme(legend.position="none") +
     ggtitle("") 
fig6B

ggsave(filename="C:/Users/aossola/Desktop/fig6B.pdf", plot=fig6B, dpi=300)
ggsave(filename="C:/Users/aossola/Desktop/fig6B.png", plot=fig6B, dpi=300)