index.Rmd

---
title: Vegetation communities on commercial developments are heterogenous and determined by development and landscaping decisions, not neighborhood socioeconomics
bibliography: ../../../bibliography.bib
output: 
    word_document:
        reference_docx: PeerJ_Preprint.docx
always_allow_html: yes
---

```{r VEG-SETUP, include=FALSE, echo = FALSE}

# Written by Karen Dyson: @_kdyson; karenldyson@gmail.com
# 
# Note that re-running this code to produce the manuscript will result in slightly different p-values for PERMANOVA each time as it is estimating the underlying distribution of F to estimate p.


knitr::opts_chunk$set(
    echo = FALSE,
    warning = FALSE,
    message = FALSE
    )
# Previous titles:
# Heterogeneity within and between land uses: commercial office vegetation is determined by development and landscaping decisions
# Woody vegetation communities on office developments are heterogenous and meaningfully influenced by development and landscaping choices
# Development and landscaping choices differentiate heterogenous tree and shrub communities on office developments

# This manuscript .Rmd code relies on two sourced R codes:
# 'vegetation_data_processing_Combined.R' and 'vegetation_data_analysis_Combine.R.' These do the
# heavy lifting before importing results here for ggplot-ifying and presenting
# in Results and Discussion.
# 
# This .Rmd does not create the title page required for PLOS ONE

## Data processing --> analysis --> Rmarkdown (this doc)


## Run sourced R codes

  source("vegetation_data_analysis_Combine.R")
  


## Load required libraries
  library(ggplot2)
  library(rgl)
  library(png)
  library(grid)
  library(gridExtra)
  library(ggrepel)
  library(tibble)
  library(pander)
  library(vegan)
  library(dplyr)

  options(digits = 3)

    levels(vegetation.clusters$tree.cluster.name) <- c("Native", "Ornamental")
    levels(vegetation.clusters$shrub.cluster.name) <- c("Ornamental", "Native")
    
    
    population.covariates.subset <-
    population.covariates[population.covariates$vegetation_class %in%
    c("HH", "MC", "MD", "MM", "LL"),]
    
    
    PA.tree.sp.only <- matrify.tree.sp.only
        PA.tree.sp.only[PA.tree.sp.only > 0] <- 1
    PA.shrub.site <- matrify.shrub.site
        PA.shrub.site[PA.shrub.site > 0] <- 1
        
  
```


# Abstract

In urban ecosystems, woody vegetation communities and the ecosystem functions and habitat they provide are largely controlled by humans. 
These communities are assembled during development, landscaping, and maintenance processes according to decisions made by human actors.
While vegetation communities on residential land uses are increasingly well studied, these efforts generally have not extended to other land uses, including commercial property.
To fill this gap, I surveyed tree and shrub communities on office developments located in Redmond and Bellevue, Washington, USA, and explored whether aggregated neighborhood and parcel scale socio-economic variables or variables describing the outcome of development and landscaping actions better explained variation in vegetation communities.
I found that both tree and shrub communities on office developments are heterogenous, with sites characterized by native or ornamental vegetation.
The heterogeneity I observed in vegetation communities within one land use suggests that different ecosystem functions, habitat quality, and habitat quantities are provided on office developments.
Greater provision of e.g. native conifer habitat is possible using currently existing developments as models.
Additionally, the outcome of development and landscaping decisions explained more variation in community composition than the socio-economic factors found significant on residential property.
Together with previous research showing that residential property owner attitudes and actions are more important than socio-economic descriptors, my results suggest that individual motivators, including intended audience, may be the primary determinant of urban vegetation communities. 
Future urban ecology research should consider sampling the vegetation gradient within land uses, better understanding individual motivation for vegetation management, and creating models of the urban ecosystems that account for alternate decision pathways on different land uses.


# Introduction

Woody vegetation community composition, structure, and distribution are largely controlled by human decisions and actions in urban ecosystems [@faeth2011; @avolio2018; @pickett2008; @gibb2002; @mullaney2015; @peters2011; @sharpe1986]. 
Development, landscaping, and ongoing maintenance are important milestones for management decisions that determine vegetation community characteristics. 
Changes to the vegetation community alter ecosystem service provision and habitat quality and quantity [@faeth2011; @byrne2007; @lehmann2014]. 

Development has replaced fire as the primary disturbance driver and precursor to new forest stands in the coastal Puget Sound region of Washington [@gibb2002; @walcott1899; @sharpe1986; @halpern1995].
The mechanisms of disturbance when clearing and grading land for development include removing vegetation, removing topsoil, and compacting soil with heavy equipment [Fig \@ref(fig:FIG-development-pic); @andres2004; @dorney1984; @mckinney2002; @grimm2017; @turner2005].
Decisions made by developers and landowners at the time of development determine the extent of disturbance and influence future site conditions (Fig 1). 
For example, choosing to preserve existing trees determines legacy vegetation and influences stand characteristics like age and size [@dorney1984]. 


```{r FIG-development-pic, fig.cap='**Commercial development project located in Redmond, Washington**. Depicted: a. clearing the site of vegetation and b. grading the site and digging the foundation. Photo credit: K. Dyson.'}

	knitr::include_graphics("Fig1.tif")


```


Vegetation succession in urban ecosystems is determined through ecological processes such as dispersal and regeneration from seed banks and through landscaping and ongoing maintenance decisions made by developers and landowners [@zipperer2010process].
The latter has become the dominant process [@grimm2017; @faeth2011; @widrlechner1990; @dorney1984; @vanheezik2014; @goodness2018; @kendal2012].
Ornamental introduced shrubs, trees, or grasses are often chosen for landscaping, though using native species is becoming more common [@vanheezik2014; @faeth2011; @germaine1998; @burghardt2009; @blair1996; @mckinney2002]. 
Once planted, these require significant ongoing maintenance inputs to arrest succession and maintain the desired aesthetic [@zipperer2010process; @faeth2011; @lebauer2008; @lepczyk2004]. 
Along with trees retained through tree preservation policies, landscape plantings represent a significant portion of the vegetation on site and of the habitat quality and quantity available to other organisms [@faeth2011; @avolio2018].

Drivers determining vegetation management decisions, actions, and outcomes are multi-scalar, and include policy, community social pressures, aggregated neighborhood socio-economic status, and the motivations and preferences of individual landowners [@cook2012]. 
Relevant public policies include clearing and grading permitting processes, impervious surface maximums and minimums via parking space requirements, tree protection policies, canopy cover goals, and vegetation planting policies [@young2011; @EPA2011; @delaria2008]. These policies are frequently enacted to protect ecosystem services, including carbon sequestration and aesthetic preferences [@wolf2005; @collins2011; @elmendorf2008; @goodness2018; @conway2016; @dorney1984].

Community drivers include the social norms and customs that influence individual behavior [@cook2012]. 
On residential properties, homeowners alter preferences for their own yards in response to the choices of nearby neighbor's yards [@nassauer2009], though their assumptions about neighbor preference are not always accurate [@peterson2012]. 
On commercial properties, owners may alter preferences to appeal to prospective and existing tenants [@laverne2003; @levy2013].

Neighborhood socio-economic status is often identified an important predictor of vegetation communities in studies of residential property.
These variables are aggregated to the neighborhood scale, though they reflect group membership of the individual.
Group membership often serves as a proxy for commonly held attitudes and ability to manipulate their enviroment [@heynen2006], however they are also inextricably linked with systematic forces of inequality influencing the spatial distribution of wealth in a city [@schell2019].
Socio-economic variables correlated with canopy cover and other vegetation metrics include: current and historic household income [@leong2018; @hope2003; @larsen2006; @boone2010; @krafft2016; @sierra2017; @clarke2013; @luck2009; @avolio2015; @heynen2006; @avolio2018], education level [@krafft2016; @luck2009; @kendal2012], ethnic composition [@luck2009; @grove2006data; @heynen2006; @leong2018], home value [@mills2016], home ownership [@heynen2006], and housing age [@jim1993; @sierra2017; @boone2010; @clarke2013; @avolio2018].
However, researchers that disaggregate socio-economic characteristics find that individual attitudes may be more important than these aggregated measures that serve as a proxy  [@shakeel2014; @kendal2012].

In municipal parks, education level and park age were only occasionally important  [@martin2004; @kendal2012]. These aggregated measures are thought to influence vegetation through neighborhood investment, advocacy, and legacy effects [@boone2010; @rigolon2018; @schell2019]. 
Individual scale drivers on other land uses are poorly studied.

Developers for all land uses are often motivated by cost and investment decisions [@almagor2017]. 
Bulk construction paired with removing existing vegetation is purportedly cheaper, though preserving vegetation may be less expensive in the long run [@mckinney2002].

These management decisions which create vegetation communities and patterns in cities also impact ecosystem function, food webs, and biodiversity [@faeth2011; @avolio2018; @faeth2005; @wittig2010; @mckinney2002; @dorney1984]. 
Different tree and shrub species have different capacity for carbon sequestration [@tenneson2013; @tang2016].
Introduced ornamentals generally do not same insect species, or the same biomass or diversity of fauna as native habitat [@burghardt2009; @mckinney2002; @crisp1998; @rebele1994; @mach2018].
These changes to habitat quality and quantity also impact higher trophic levels [@burghardt2009; @marzluff2001; @faeth2011; @alberti2003; @alberti2005; @polasky2005; @rosenzweig2003]. 
For the urban matrix to support conservation, decision makers across land uses need to take actions that support locally important vegetation habitat [@goddard2010; @miller2002]. 

While the drivers and outcomes of decision making are increasingly well studied on residential private property, other land uses have not been given the same attention [@snep2011; @bourne2013]. 
For example, commercial and industrial land uses are generally included only as independent variables in remote sensing studies of factors influencing percent canopy cover [e.g. @fan2019; @mills2016]. 
Additionally, research where the unit of analysis is defined by the area of influence of specific decision makers is also needed.
Aggregated measures, such as vegetation transects through neighborhoods or canopy cover of a census block, cannot examine specific decision outcomes as they conflate different actors and their motivations and actions, and previous research shows that motivations differ between actors [@kendal2012; @tenneson2013].

To fill this gap, I examined woody vegetation community composition on office developments in Bellevue and Redmond, Washington, USA. 
Specifically, I examined 1) tree and shrub communities present on office developments and 2) whether aggregated and parcel specific socio-economic variables or development and landscaping outcomes better explained observed variation in vegetation communities.

I hypothesized that vegetation communities on office developments would be heterogeneous.
I also hypothesized that aggregated socio-economic variables found significant in explaining vegetation patterns on residential property would not be significant on office developments [@avolio2018; @hope2003; @conway2016],
but that parcel scale variables would. 
Finally, I hypothesized that the outcome of development and landscaping actions would better explain variation in tree and shrub community structure.
I found that woody vegetation communities on office developments are heterogenous with distinct community types, and that in contrast with residential property, development and landscaping actions explain this variability better than socio-economic variables. 



# Materials and methods

This study was approved by the University of Washington Human Subjects Division under Determination of Exemption #48246. Field surveys were approved in writing by private property owners or managers of office developments. I did not perform any animal research or collect plant, animal, or other materials from a natural setting.

## Study area and site selection

Redmond (2017 population 64,000) and Bellevue (population 144,000) are located east of Seattle in King County, Washington [@Census2017]. Both cities share a similar ecological history, a similar disturbance timeline for logging and agriculture, and have grown considerably since the opening of the Evergreen Point Floating Bridge (SR 520) in 1963. They are at similar elevations (< 160 m) and experience the same climate and similar weather. 

The sampling frame was limited to Redmond and Bellevue north of I-90, excluded developments in Bellevue’s central business district, and contained parcels defined as office use by the King County Assessor’s Office (Fig \@ref(fig:FIG-study-area-map)). 
I grouped adjacent parcels built within three years of one another and with the same owner to create a unit of analysis based on human action not cadastral boundaries.
This initial population size was `r length(population.covariates$vegetation_class) ` developments.

```{r FIG-study-area-map, fig.cap='(ref:studyareacap)'}

  knitr::include_graphics("Fig2.tiff")


```

(ref:studyareacap) **Map of office development study sites in Redmond and Bellevue, Washington.** The population of office developments with High (HH), Medium Canopy (MC), Medium Diverse (MD), Medium (MM), and Low (LL) vegetation types are represented with colored circles; excluded sites (no vegetation/LP, wetlands/WW, and under construction/XX) are represented with gray triangles. Sampled sites are shown with colored diamonds.

I used disproportionate stratified random sampling to ensure that my sample included sites across the entire vegetation gradient.
I classified the vegetation at each potential study site into categories using visual estimation during site visits in winter 2014 (Fig \@ref(fig:FIG-vegetation-types), Table \@ref(tab:TAB-vegclass)). Sites with no vegetation, with wetlands, or those that were currently under construction or undergoing landscape replanting were excluded from the analysis (total `r sum(population.covariates$vegetation_class %in% c("WW", "XX", "LP")) ` sites). The remaining pool of `r sum(!(population.covariates$vegetation_class %in% c("WW", "XX", "LP")))` potential sites had no notable hydrological features on site.


```{r FIG-vegetation-types, fig.cap='(ref:vegetationtypes)', dpi=300}

  knitr::include_graphics("Fig3.tif")


```

(ref:vegetationtypes) Examples of each vegetation type. From top left to bottom right: High (HH); Medium Canopy (MC); Medium Diverse (MD); Medium (MM); Low (LL); no vegetation (LP; excluded); wetlands (WW; excluded).


```{r TAB-vegclass}

vegclasssum <- tibble(
`Vegetation Type` = c(
"High",
"Medium Canopy",
"Medium Diverse",
"Medium",
"Low",
"No Vegetation",
"Wetlands",
"Under Construction"
),
`Tree Cover` = c(
"30% native tree cover",
"30% native tree cover",
"15% tree cover",
"15% tree cover",
"< 10% tree cover",
"No trees",
"No requirement",
"No requirement"
),
`Shrub Richness` = c(
"> 5 native shrub genera",
"No requirement",
"> 5 native shrub genera",
"> 5 shrub genera",
"< 5 shrub genera",
"No shrubs",
"No requirement",
"No requirement"
),
`Strata Size` = c(
sum(population.covariates$vegetation_class == "HH"),
sum(population.covariates$vegetation_class == "MC"),
sum(population.covariates$vegetation_class == "MD"),
sum(population.covariates$vegetation_class == "MM"),
sum(population.covariates$vegetation_class == "LL"),
sum(population.covariates$vegetation_class == "LP"),
sum(population.covariates$vegetation_class == "WW"),
sum(population.covariates$vegetation_class == "XX")

),
`Sampled (n)` = c(5, 3, 4, 3, 5, 0, 0, 0),
Notes = c(
" ",
" ",
" ",
" ",
" ",
"Excluded from further analysis",
"Excluded from further analysis",
"Excluded from further analysis"
)
)

knitr::kable(vegclasssum, caption = "**Vegetation type assignment criteria and strata size.**",
             format = "pandoc", escape = TRUE)


```

I conducted stratified random sampling on sites with High, Medium Canopy, Medium Diverse, Medium, and Low vegetation types.
Site selection was restricted based on site area and surrounding impervious surfaces for concurrent studies [@dysonPhD].
Limiting sampling of these extremes reduced my ability to detect community differences along these gradients, though socio-economic variables are not covariate.

I requested property access through three mailings sent to the property owner or manager on file in the King County Assessor's database [@dyson2019].
I targeted vegetation categories underrepresented in my sample in the second and third mailings. Of 46 mailed requests, 20 (`r round(20/46 * 100, 1)`%) received no response or were not deliverable. 
Of the 26 (`r round(26/46 * 100, 1)`%) responses received, 6 (`r round(6/26*100, 1)`%) of were rejected and 20 (`r round(20/26*100, 1)`%) were accepted in writing by an individual with authority to do so [@dyson2019].

Commercial use of sample sites included light industrial, white collar office space, and medical/dental offices. 
Some sites were fully leased to tenants, while others were either partly or fully owner-occupied. 
Company size ranged from less than 10 to many thousand employees.


## Independent variables

I examined three categories of independent variables: 1. aggregated and parcel scale socio-economic variables, 2. development and landscaping outcome variables, and 3. ground cover material and maintenance regime variables (Table 2). 
Aggregated socio-economic variables, sometimes called neighborhood socio-economic variables, reflect group membership of individual landowners aggregated to a common level (in the United States, often block groups). 
Parcel scale socio-economic variables reflect individual pieces of land, including building age and value. 
Socio-economic variables were derived from existing databases [Table \@ref(tab:TAB-covariate-definitions); @ACS5yr; @taxDB; @KCDB; @xian2011; @homer2015].
Variables were chosen based on previous research and analyzed in QGIS 3.2 [@hope2003; @martin2004; @walker2009; @dana2002; @grove2014; @QGIS].


```{r TAB-covariate-definitions}

covariate.definitions <- tibble(
`Variable Name` = c(
"1. AGGREGATED AND PARCEL SCALE SOCIO-ECONOMIC VARIABLES",

"Area (acre)",
"Town",
"Building Age (in 2017)",
"Building Quality",
"Appraised Land Value per Acre",
"Impervious w/in 500 m (%)",
"Median Income",
"Percent Foreign-Born",

"2. DEVELOPMENT AND LANDSCAPING OUTCOME VARIABLES",

"Stands Predate Development",
"Median Height of Dominant Conifer",
"Density of Native Conifers",

"3. GROUND COVER MATERIAL AND MAINTENANCE REGIME VARIABLES",

"Ground Cover Types (%)",
"Dead Wood (count)",
"Irrigation",
"Mulch, Herbicide, and/or Fertilizer Application"
),

Definition = c(
" ",

"Site area, in acres.",
"Location, Bellevue or Redmond.",
"Age of building on site (or mean age for multiple buildings) in 2017.",
"Categorical 'quality class' assigned to buildings on the site",
"Appraised land value divided by site area. One missing assessed land values were replaced with population median land value.",
"Percent impervious surface within 500 m of the site's perimeter.",

"The median income of residents for the site's block group",
"The percent of residents born outside of the United States for the site's block group.",

" ",

"Binary variable indicating presence of a cluster of three+ trees that predate development. (Development)",
"Median height (m) of five dominant native conifer trees; age proxy. (Development)",
"Total density of Douglas-fir, western redcedar, and western hemlock (Development/Landscaping).",

" ",

"Ground cover types on site including lawn, mulch, and impervious surface.",
"Total abundance of stumps, logs, and snags on site.",
"Binary variable indicating whether irrigation is used during the summer months.",
"Binary variables (3) indicating whether landscaping crew applies mulch, herbicides, or fertilizers to a site."

),

`Data Source` = c(
" ",
"King County Assessor",
"King County Assessor",
"King County Assessor",
"King County Assessor",
"King County Assessor",
"National Land Cover Database 2011 Percent Developed Imperviousness dataset updated in 2014",

"American Community Survey 2014 5-year block group",
"American Community Survey 2014 5-year block group",

" ",
"Site survey",
"Site survey",
"Site survey",

" ",
"Site survey",
"Site survey",
"Interviews and site survey",
"Interviews and site survey"


),

`Population` = c(
" ",

paste0("Range: ", paste0(range(population.covariates.subset$acres)%>%round(2), collapse = "-"),
       "; Mean(SD): ", mean(population.covariates.subset$acres)%>%round(2), "(",
       sd(population.covariates.subset$acres)%>%round(2), ")"),
paste("Bellevue:", sum(population.covariates.subset$DistrictName=="BELLEVUE", na.rm = T),
      "Redmond:", sum(population.covariates.subset$DistrictName=="REDMOND", na.rm=T)),
paste0("Range: ", paste0(range(population.covariates.subset$age.2017, na.rm = T)%>%round(1),
                         collapse = "-"),
       "; Mean(SD): ", mean(population.covariates.subset$age.2017, na.rm = T)%>%round(1), "(",
       sd(population.covariates.subset$age.2017, na.rm = T)%>%round(2), ")"),

paste("Below Average:", sum(population.covariates.subset$BldgQuality_lookup %in%
                                c("LOW COST", "LOW/AVERAGE"), na.rm=T),
      "Average:", sum(population.covariates.subset$BldgQuality_lookup=="AVERAGE", na.rm = T),
      "Average/Good:", sum(population.covariates.subset$BldgQuality_lookup=="AVERAGE/GOOD", na.rm=T),
      "Good:", sum(population.covariates.subset$BldgQuality_lookup=="GOOD", na.rm=T),
      "Good/Excellent:", sum(population.covariates.subset$BldgQuality_lookup=="GOOD/EXCELLENT", na.rm=T)
      ),

paste0("Range: ", paste0(range(population.covariates.subset$landrentperacre)%>%round(0), collapse = "-"),
       "; Mean(SD): ", mean(population.covariates.subset$landrentperacre)%>%round(0), "(",
       sd(population.covariates.subset$landrentperacre)%>%round(0), ")"),
paste0("Range: ", paste0(range(population.covariates.subset$dissolved_parcel_500m_buffer_impervious_500m_mean)%>%round(1), collapse = "-"),
       "; Mean(SD): ", mean(population.covariates.subset$dissolved_parcel_500m_buffer_impervious_500m_mean)%>%round(1), "(",
       sd(population.covariates.subset$dissolved_parcel_500m_buffer_impervious_500m_mean)%>%round(1), ")"),
paste0("Range: ", paste0(range(population.covariates.subset$MedianHouseholdIncome_B19013e1)%>%round(0), collapse = "-"),
       "; Mean(SD): ", mean(population.covariates.subset$MedianHouseholdIncome_B19013e1)%>%round(0), "(",
       sd(population.covariates.subset$MedianHouseholdIncome_B19013e1)%>%round(0), ")"),
paste0("Range: ", paste0(range(population.covariates.subset$Proportion_ForeignBorn_B99051e5*100)%>%round(1), collapse = "-"),
       "; Mean(SD): ", mean(population.covariates.subset$Proportion_ForeignBorn_B99051e5*100)%>%round(1), "(",
       sd(population.covariates.subset$Proportion_ForeignBorn_B99051e5*100)%>%round(1), ")"),

" ",

"NA",
"NA",
"NA",

" ",

"NA",
"NA",
"NA",
"NA"
    
),

`Sample` = c(
" ",

paste0("Range: ", paste0(range(sample.covariates$acres)%>%round(2), collapse = "-"),
       "; Mean(SD): ", mean(sample.covariates$acres)%>%round(2), "(",
       sd(sample.covariates$acres)%>%round(2), ")"),
paste("Bellevue:", sum(sample.covariates$DistrictName=="BELLEVUE", na.rm = T),
      "Redmond:", sum(sample.covariates$DistrictName=="REDMOND", na.rm=T)),
paste0("Range: ", paste0(range(sample.covariates$age.2017), collapse = "-"),
       "; Mean(SD): ", mean(sample.covariates$age.2017)%>%round(1), "(",
       sd(sample.covariates$age.2017)%>%round(1), ")"),

paste("Below Average:", sum(sample.covariates$BldgQuality_lookup %in%
                                c("LOW COST", "LOW/AVERAGE"), na.rm=T),
      "Average:", sum(sample.covariates$BldgQuality_lookup=="AVERAGE", na.rm = T),
      "Average/Good:", sum(sample.covariates$BldgQuality_lookup=="AVERAGE/GOOD", na.rm=T),
      "Good:", sum(sample.covariates$BldgQuality_lookup=="GOOD", na.rm=T),
      "Good/Excellent:", sum(sample.covariates$BldgQuality_lookup=="GOOD/EXCELLENT", na.rm=T)
      ),

paste0("Range: ", paste0(range(sample.covariates$landrentperacre)%>%round(0), collapse = "-"),
       "; Mean(SD): ", mean(sample.covariates$landrentperacre)%>%round(0), "(",
       sd(sample.covariates$landrentperacre)%>%round(0), ")"),
paste0("Range: ", paste0(range(sample.covariates$dissolved_parcel_500m_buffer_impervious_500m_mean)%>%round(1), collapse = "-"),
       "; Mean(SD): ", mean(sample.covariates$dissolved_parcel_500m_buffer_impervious_500m_mean)%>%round(1), "(",
       sd(sample.covariates$dissolved_parcel_500m_buffer_impervious_500m_mean)%>%round(1), ")"),
paste0("Range: ", paste0(range(sample.covariates$MedianHouseholdIncome_B19013e1)%>%round(0), collapse = "-"),
       "; Mean(SD): ", mean(sample.covariates$MedianHouseholdIncome_B19013e1)%>%round(0), "(",
       sd(sample.covariates$MedianHouseholdIncome_B19013e1)%>%round(0), ")"),
paste0("Range: ", paste0(range(sample.covariates$Proportion_ForeignBorn_B99051e5*100)%>%round(1), collapse = "-"),
       "; Mean(SD): ", mean(sample.covariates$Proportion_ForeignBorn_B99051e5*100)%>%round(1), "(",
       sd(sample.covariates$Proportion_ForeignBorn_B99051e5*100)%>%round(1), ")"),

" s",

paste("Yes:", sum(management.landscaping$stand.predate.development=="yes", na.rm = T),
      "No:", sum(management.landscaping$stand.predate.development=="no", na.rm=T)),

paste0("Range: ", paste0(range(management.landscaping$height.m.median)%>%round(1), collapse = "-"),
       "; Mean(SD): ", mean(management.landscaping$height.m.median)%>%round(1), "(",
       sd(management.landscaping$height.m.median)%>%round(1), ")"),
paste0("Range: ", paste0(range(management.landscaping$tree.conifer.nat.dens)%>%round(1), collapse = "-"),
       "; Mean(SD): ", mean(management.landscaping$tree.conifer.nat.dens)%>%round(1), "(",
       sd(management.landscaping$tree.conifer.nat.dens)%>%round(1), ")"),

" s",

paste("Mean (SD) Grass:", mean(dens.matrify.gc$grass*100)%>%round(1),
            "(",
            sd(dens.matrify.gc$grass*100)%>%round(1), ");",
      "Impervious:", mean(dens.matrify.gc$impervious.sqft*100)%>%round(1),
            "(",
            sd(dens.matrify.gc$impervious.sqft*100)%>%round(1), ");",
      "Dirt/Litter:", mean(dens.matrify.gc$dirt.litter*100)%>%round(1),
            "(",
            sd(dens.matrify.gc$dirt.litter*100)%>%round(1), ")"),


paste0("Range: ", paste0(range(management.landscaping$height.m.median)%>%round(1), collapse = "-"),
       "; Mean(SD): ", mean(management.landscaping$height.m.median)%>%round(1), "(",
       sd(management.landscaping$height.m.median)%>%round(1), ")"),

paste("Yes:", sum(management.landscaping$Irrigation=="yes", na.rm = T),
      "No:", sum(management.landscaping$Irrigation=="no", na.rm=T)),

paste("Mulch Y/N:", paste(sum(management.landscaping$Mulch=="yes", na.rm = T),
                           sum(management.landscaping$Mulch=="no", na.rm=T), 
                           sep = "/"),
      "Herbicide:", paste(sum(management.landscaping$Herbicide=="yes", na.rm = T),
                           sum(management.landscaping$Herbicide=="no", na.rm=T), 
                           sep = "/"),
      "Fertilizer:", paste(sum(management.landscaping$Fertilizer=="yes", na.rm = T),
                           sum(management.landscaping$Fertilizer=="no", na.rm=T), sep = "/")
    )

)
)


knitr::kable(covariate.definitions,
       caption = '(ref:definition)',
       format = "pandoc", escape = TRUE, format.args = list(decimal.mark = '.', big.mark = ","))
 


```


(ref:definition) **Definition of independent variables used in PERMANOVA and correlation analysis [@taxDB; @KCDB; @ACS5yr; @xian2011; @homer2015].** Summary statistics for independent variables for both the population of office developments in Redmond and Bellevue and the sample of sites studied (`r sum(!(population.covariates$vegetation_class %in% c("WW", "XX", "LP"))) ` and 20 sites, respectively). Median income ($) and proportion foreign born are included to compare patterns in commercial developments with patterns found significant in residential research.


I measured the height of dominant native conifers using a Nikon Forestry Pro Laser Rangefinder.
This is a proxy measure for age as I did not collect tree cores due to liability concerns [@dyson2019].
I used historical records and site construction plans to determine whether each site retained a stand of three adjacent tree predating site development.
I used _Pseudotsuga menziesii_ (Mirb.) Franco, _Thuja plicata_ Donn ex D. Don, and _Tsuga heterophylla_ (Raf.) Sarg. counts to calculate native conifer density. 

I digitized broad ground cover material classes in QGIS to calculate area [@QGIS].
Pervious cover types recorded include dense vegetation, dirt/litter, lawn (turf grass including moss and forb species), gravel, dense ivy, mulch, and water. 
I used semi-structured interviews of property owners, managers, and landscaping services along with site visits to obtain maintenance regime variables [University of Washington Human Subjects Division Determination of Exemption #48246; @dexter1970; @harvey2011].
Irrigation, mulching, herbicide, and fertilizer application had only three "no" responses and thus could not be used to draw any well supported conclusions.


## Vegetation data collection

I censused vegetation communities during the summer of 2015, excluding saplings with DBH < 3''.
Each tree or shrub was identified to species or genus in consultation with experts at the Center for Urban Horticulture at University of Washington [@sibley2009; @dirr1997hardy; @dirr2009manual].
Some tree and shrub species were grouped at the genus level due to the abundance of very similar cultivars in the landscaping trade, including _Malus_ Mill. [@sierra2017]. 
Following previous studies, I grouped conifers under 2 m into a broad class of dwarf conifer species [@daniels2006yds]. 
`r sum(matrify.tree.sp.only[ grepl("cedar.1|cedar.2|unknown", colnames(matrify.tree.sp.only))]) ` individual trees (`r 100*sum(matrify.tree.sp.only[ grepl("cedar.1|cedar.2|unknown", colnames(matrify.tree.sp.only))])/sum(matrify.tree.sp.only) `% of total trees) and `r sum(matrify.shrub.site[grepl("Unknown", colnames(matrify.shrub.site))]) ` shrubs (`r 100*(sum(matrify.shrub.site[grepl("Unknown", colnames(matrify.shrub.site))])/sum(matrify.shrub.site) ) `% of total shrubs) could not be identified; these were given a unique identifier code for multivariate community analysis.

I assigned tree and shrub genera to one of three provenance categories---native, non-native, or ambiguous [@USGSplants; @USDAplants].
The ambiguous category was used for genera including both native and non-native cultivated species that are difficult to distinguish, and/or frequently interbred and sold as crosses.
For example, some _Mahonia_ Nutt. species are native (_M. aquifolium_ Pursh Nutt. and _M. nervosa_ Pursh Nutt.), while others originate in Asia (_Mahonia japonica_ Thumb. DC.) and many hybrids are bred and sold by nurseries (e.g. _Mahonia x media_ "Charity" Brickell). 


## Identifying and describing vegetation clusters on office developments

I standardized tree and shrub abundance data and ground cover area by total site area in acres.
This transformation preserves parcel boundaries as the unit of analysis and reflects developer and landowner actions during and following development that determine the amount of impervious surface and pervious area, the number of trees preserved, and the number of trees and shrubs planted.
All analyses were performed in R [@citeR].

To delineate vegetation community clusters on office developments, I used a flexible agglomerative nesting function (agnes {vegan}) with a beta of -0.5 to produce an ecologically interpretable dendrogram with minimal chaining [@mccune2002; @dufrene1997; @milligan1989; @breckenridge2000; @vegan]. 
Using the resulting groups, I performed indicator species analysis, which assesses the predictive values of species as indicators of the conditions at site groups, using multipatt {indicspecies} [@decaceres2010; @indicspecies; @decaceres2013] and a custom wrapper for repetition of the permutation-based function [@dysoncode]. 
To examine citywide patterns, I extrapolated cluster membership to the entire population of office developments in the study area using proportions.

I used simple univariate permutational analysis of variance (PERMANOVA) models to test if continuous variables differed between vegetation community cluster groups and Pearson’s Chi-squared test to test if categorical variables differed [adonis2 {vegan} and chisq.test {stats}; @anderson2001; @vegan]. 
PERMANOVA is a permutation-based implementation of analysis of variance (ANOVA) that avoids assumptions about underlying distributions of community structure and can be used with non-Euclidian distance matrices [@anderson2001].
Bartlett tests of homogeneity found no difference between group variances (bartlett.test {stats}).


## Explaining variation in tree and shrub community structure

I analyzed tree and shrub communities separately to determine if the two communities responded differently, and as development and landscaping outcome variables were derived from measurements of the tree community.

I used non-metric multidimensional scaling (NMDS) to evaluate relationships between development and landscaping variables and the tree community [metaMDS {vegan}; @mccune2002; @vegan].
NMDS is a rank-based ordination technique that is robust to data without identifiable distributions, can be used with any distance or dissimilarity measure [@mccune2002].
To determine the relationship between development and landscaping outcome variables and the tree community, I used convex hull plots and fitted environmental vectors [ordiplot and envfit {vegan}; @vegan].

I used PERMANOVA to test relationships between other variables and the tree community and all variables and the shrub community. 
I used a multi-step approach to avoid transforming independent variables or using ordination to collapse related variables, as these actions make results less interpretable for urban planners and other professionals.
I first tested each independent variable with community matrices in a simple multivariate PERMANOVA model.
To ensure differences in categorical variables were due to location and not dispersion, I used ANOVA to test for significant differences in dispersion [anova {stats} and betadisper {vegan}; @vegan].
I then constructed models using all variables with significant pseudo-_F_ values in all possible single and multiple multivariate model combinations. Significance was assessed at the $\alpha \leq$ 0.05 level following Holm-Bonferroni correction for multiple comparisons.
I used a custom Akaike information criterion with correction (AICc) function based on residual sums of squares to compare models and identify those with the best support [@dysoncode].


# Results and discussion

## Observed woody vegetation communities

I recorded a total of `r sum(matrify.tree.sp.only)` trees and  `r sum(matrify.shrub.site)` shrubs from `r length(unique(tree.native$tree.species)) ` and `r length(unique(shrub.native$shrub.species)) ` taxonomic groups respectively (S1 and S2 Tables). 
Only _Rhododendron_ L. were found on all 20 sites surveyed. Four tree species and nine shrub species were found on more than half of all office developments, with `r sum(colSums(PA.tree.sp.only) == 1) ` tree species and `r sum(colSums(PA.shrub.site) == 1) ` shrub taxa found only on only one development. 

Native tree species accounted for `r (100 * (sum(matrify.tree.sp.only.native) / sum(matrify.tree.sp.only)))%>%round(1) `% of total individuals observed, and three of the top five most abundant species. 
On average, native species accounted for `r (100 * mean(rowSums(matrify.tree.sp.only.native) / rowSums(matrify.tree.sp.only)))%>%round(1) `% of the trees found on each office development, though sites varied widely with 0%--99% native trees.
_Pseudotsuga menziesii_ was by far the most abundant tree species, with `r (100*( 746 / sum(matrify.tree.sp.only)))%>%round(1)`% of observed individuals.
_Thuja plicata_ (`r (100*( 245 / sum(matrify.tree.sp.only)))%>%round(1)`%),
_Acer macrophyllum_ Pursh (`r (100*( 218 / sum(matrify.tree.sp.only)))%>%round(1)`%),
_Acer rubrum_ L. (`r (100*( 132 / sum(matrify.tree.sp.only)))%>%round(1)`%), and
_Acer platanoides_ L. (`r (100*( 100 / sum(matrify.tree.sp.only)))%>%round(1)`%) complete the top five. _Prunus_ L. and _Alnus rubra_ Bong. were both widespread taxa (found on 12 and 9 sites, respectively) but were never abundant on any one site.

In contrast, native shrub species accounted for only `r (100 * (sum(matrify.shrub.site.native) / sum(matrify.shrub.site)))%>%round(1)`% individual shrubs observed.
On average, native shrubs accounted for `r (100 * mean(rowSums(matrify.shrub.site.native) / rowSums(matrify.shrub.site)))%>%round(1) `% of the shrubs observed at each office development, and never more than `r (100 * max(rowSums(matrify.shrub.site.native) / rowSums(matrify.shrub.site)))%>%round(1) `% of individual shrubs. 
The two most abundant shrub species were the native
_Gaultheria shallon_ Pursh (`r (100 * 1267 / sum(matrify.shrub.site))%>%round(1)`%), which frequently occurs in low, dense mats, and 
_Berberis Mahonia gp._ Nutt. (`r (100 * 1007 / sum(matrify.shrub.site))%>%round(1)`%) which is comprised of native, introduced, and hybrid species.
The rest of the top five most abundant shrub species were all non-native, including 
_Prunus laurocerasus_ L. (`r (100 * 682 / sum(matrify.shrub.site))%>%round(1)`%),
_Rhododendron_ (`r (100 * 612 / sum(matrify.shrub.site))%>%round(1)`%), and 
_Cornus sericea_ L. (`r (100 * 419 / sum(matrify.shrub.site))%>%round(1)`%). 

```{r beta-calc}

tree.sorensen <- vegdist(matrify.tree.sp.only, binary = TRUE)%>%round(3) # take mean for sorensen index
shrub.sorensen <- vegdist(matrify.shrub.site, binary = TRUE)%>%round(3) #ditto

tree.z.beta <- betadiver(matrify.tree.sp.only, "z")%>%round(3)
shrub.z.beta <- betadiver(matrify.shrub.site, "z")%>%round(3)


beta.tibble <- tibble(Metric = c("Tree Sorensen", "Tree Arrhenius Model z", "Shrub Sorensen", "Shrub Arrhenius Model z"),
                      Minimum = c(min(tree.sorensen), min(tree.z.beta), min(shrub.sorensen), min(shrub.z.beta)),
                      Maximum = c(max(tree.sorensen), max(tree.z.beta), max(shrub.sorensen), max(shrub.z.beta)),
                      Mean = c(mean(tree.sorensen)%>%round(3),
                               mean(tree.z.beta)%>%round(3),
                               mean(shrub.sorensen)%>%round(3),
                               mean(shrub.z.beta)%>%round(3)),
                      `S.D.` = c(sd(tree.sorensen)%>%round(3),
                                 sd(tree.z.beta)%>%round(3),
                                 sd(shrub.sorensen)%>%round(3),
                                 sd(shrub.z.beta)%>%round(3)),
                      Median = c(median(tree.sorensen), median(tree.z.beta),
                                 median(shrub.sorensen), median(shrub.z.beta)))

```

Measures of tree and shrub abundance, density, and diversity varied substantially between sites (Table \@ref(tab:TAB-site-summary-tables-veg)).
In general, total species richness and native species richness were positively correlated (Pearson's Correlation for Tree:  `r cor(descriptive.stats.site$tree.native.richness, descriptive.stats.site$tree.sp.richness)`; Shrub: `r cor(descriptive.stats.site$shrub.native.richness, descriptive.stats.site$shrub.sp.richness)`), though four sites with above average species richness had three or fewer native species planted.
Remnant large native conifer abundance, primarily _Pseudotsuga menziesii_, greatly contributed to sites with greater tree abundance (Pearson's: `r cor(descriptive.stats.site$tree.conifer.L.native, descriptive.stats.site$tree.site.abundance)`); consequently, Shannon diversity was generally lower on sites with more native trees (Pearson's: `r cor(descriptive.stats.site$tree.conifer.native.abund, descriptive.stats.site$tree.Shannon)`).


```{r TAB-site-summary-tables-veg}

diver.metrics <- data.frame(matrix(unlist(pretty.descriptive.stats.site),
                          nrow = nrow(pretty.descriptive.stats.site),
                          byrow = FALSE),
                   stringsAsFactors = FALSE)
colnames(diver.metrics) <- colnames(pretty.descriptive.stats.site)
diver.metrics <- rbind(diver.metrics,beta.tibble)

diver.metrics.display <- c("Tree Abundance",
                           "Tree Density",
                           "Tree Species Richness",
                           "Tree Shannon Diversity",
                           "Tree Effective Species Richness",
                           "Tree Sorensen",
                           "Tree Arrhenius Model z",
                           "Native Tree Abundance",
                           "Native Tree Density",
                           "Native Tree Species Richness",
                           "Native Tree Shannon Diversity",
                           "Native Tree Effective Species Richness",
                           "Native Conifer Abundance",    
                           "Native Conifer Density",
                           "Shrub Abundance",
                           "Shrub Density",
                           "Shrub Species Richness",
                           "Shrub Shannon Diversity",
                           "Shrub Effective Species Richness",
                           "Shrub Sorensen",
                           "Shrub Arrhenius Model z",
                           "Native Shrub Abundance",
                           "Native Shrub Desnity",
                           "Native Shrub Species Richness",
                           "Native Shrub Shannon Diversity",
                           "Native Shrub Effective Species Richness")
    

knitr::kable(diver.metrics[ diver.metrics$Metric %in% diver.metrics.display , c(1:3,6,4,5)],
               row.names = FALSE,
               caption = "Metrics for tree and shrub communities on sampled office developments. H' is Shannon's diversity index [@shannon1949], effective species richness = exp(H') [@jost2006], density = individuals per acre.",
               format = "pandoc", escape = TRUE,
               format.args = list(decimal.mark = '.', big.mark = ","))

  
  

```


Overall, these measures are within the ranges reported by other urban ecology studies, though differences in methodology and particularly the use of small plots [e.g. @clarke2013] and remote sensing [e.g. @luck2009] in other studies and stratified sampling along a vegetation gradient in this study make comparison more difficult.
The most abundant tree species on office developments are similar to those on residential properties in western Washington [@tenneson2013; @mills2016]. 
Measures of diversity were generally lower than residential property [@clarke2013; @martin2004]. 
Species richness was comparable to other commercial land uses and city parks [@clarke2013; @martin2004] though lower than residential land uses [@clarke2013; @sierra2017; @martin2004; @jim1993]. 
Measures of beta diversity, suggesting low similarity between locations, were also comparable [@sierra2017].


## Divergent vegetation groups found on office developments

I identified two groups of tree and shrub vegetation (flexible beta = -0.5; agglomerative coefficients of `r beta.5.tree$ac%>%round(2) ` and `r beta.5.shrub$ac ` respectively; Table \@ref(tab:TAB-groupabund)). 
Using indicator species analysis, I identified the Native Tree group (11 sites) as characterized by _Thuja plicata_, _Acer macrophyllum_ Pursh, _Arbutus menziesii_ Pursh, and _Alnus rubra_ Bong, while the Ornamental Tree group (9 sites) is characterized by _Acer rubrum_ L. The Native Shrub group (11 sites) is characterized by _Gaultheria shallon_ Pursh, _Mahonia gp._ Nutt., _Symphoricarpos_ Duham., and _Ribes sanguineum_ Pursh, and the Ornamental Shrub group (9 sites) by _Thuja occidentalis_ L. 


```{r TAB-groupabund}

tree.native$tree.star <- paste0(tree.native$tree.scientific.update,
                                ifelse(
                                    tree.native$tree.species %in%
                                        colnames(matrify.tree.sp.only.native),
                                    "*",
                                    "" )
                                )


names(tree.cluster.native) <- tree.native[ match(names(tree.cluster.native),
                                                 tree.native$tree.species),
                                           "tree.star"] 

names(tree.cluster.orndec) <- tree.native[ match(names(tree.cluster.orndec),
                                                 tree.native$tree.species),
                                           "tree.star"]


names(shrub.cluster.native) <- paste0(names(shrub.cluster.native),
                                      ifelse(
                                      names(shrub.cluster.native) %in%
                                      colnames(matrify.shrub.site.native),
                                      "*",
                                      ""
                                      ))

names(shrub.cluster.cultiv) <- paste0(names(shrub.cluster.cultiv),
                                      ifelse(
                                      names(shrub.cluster.cultiv) %in%
                                      colnames(matrify.shrub.site.native),
                                      "*",
                                      ""
                                      ))

temp <- tibble(
    `Rank` = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"),
    `Native Tree Group` = paste(names(head(
        tree.cluster.native, 10)),
        " (",
        (head(tree.cluster.native, 10) / 11) %>% round(1),
        ")" , sep = ""),
    `Ornamental Tree Group` = paste(names(head(
        tree.cluster.orndec, 10)),  " (",
        (head(tree.cluster.orndec, 10) / 9) %>% round(1),
        ")" , sep = ""),
    `Native Shrub Group` = paste(names(head(
        shrub.cluster.native, 10)), " (",
        (head(shrub.cluster.native, 10) / 11) %>% round(1),
        ")" , sep = ""),
    `Ornamental Shrub Group` = paste(names(head(
        shrub.cluster.cultiv, 10)), " (",
        (head(shrub.cluster.cultiv, 10) / 9) %>% round(1),
        ")" , sep = "")
    )

knitr::kable(temp,
             caption = "Rank abundance of tree and shrub taxa for each community group identified by flexible-beta cluster analysis. Asterisk indicates native tree and shrub species. Number in parenthesis is mean abundance of the species in the community group.", 
             row.names = FALSE,
             format = "pandoc", escape = TRUE, format.args = list(decimal.mark = '.', big.mark = ","))


  
```


```{r CALC-mean-calcs}

shrub.ac.mean <- merge(sample.covariates, vegetation.clusters,
                     by = "SiteName",
                     sort = FALSE) %>% 
    group_by(shrub.cluster.name) %>%
    summarise_if(is.numeric, .funs = mean)

tree.ac.mean <- merge(sample.covariates, vegetation.clusters,
                     by = "SiteName",
                     sort = FALSE) %>% 
    group_by(tree.cluster.name) %>%
    summarise_if(is.numeric, .funs = mean)
 
    
shrub.ml.mean <- merge(vegetation.clusters, management.landscaping,
                     by.x = "SiteName", by.y= "Site",
                     sort = FALSE) %>%
    mutate(tot.abund = rowSums(matrify.shrub.site), 
           tot.density = rowSums(dens.matrify.shrub.site)) %>%
    group_by(shrub.cluster.name) %>%
    summarise_if(is.numeric, .funs = mean)

tree.ml.mean <- mutate(management.landscaping,
                       tot.abund = rowSums(matrify.tree.sp.only),
                       tot.density = rowSums(dens.matrify.tree.sp.only)) %>%
    group_by(tree.cluster.name) %>%
    summarise_if(is.numeric, .funs = mean)    


shrub.gc.mean <- merge(dens.matrify.gc, vegetation.clusters,
                     by.x = 0, by.y = "SiteName",
                     sort = FALSE) %>% 
    group_by(shrub.cluster.name) %>%
    summarise_at(.vars = c("dirt.litter",
                            "grass", "mulch", "impervious.sqft"
                            ), .funs = mean)

tree.gc.mean <- merge(dens.matrify.gc, vegetation.clusters,
                     by.x = 0, by.y = "SiteName",
                     sort = FALSE) %>% 
    group_by(tree.cluster.name) %>%
    summarise_at(.vars = c("dirt.litter",
                            "grass", "mulch", "impervious.sqft"
                            ), .funs = mean)


```


The two groups are distinct in the average density of trees and shrubs per site 
(Native Tree mean = `r tree.ml.mean[1,'tot.density']%>%round(1) `,
Ornamental Tree mean = `r tree.ml.mean[2,"tot.density"]%>%round(1) `
with Pr(>_F_) = `r dens.tree[1,5]%>%round(3) `; 
Native Shrub mean = `r shrub.ml.mean[2,"tot.density"]%>%round(1) `,
Ornamental Shrub mean = `r shrub.ml.mean[1,"tot.density"]%>%round(1) `
with Pr(>_F_) = `r dens.shrub[1,5]%>%round(4) `
). 
The mean median height of dominant native conifers was also significantly different between clusters for trees and shrubs (Native Tree mean  = `r tree.ml.mean[1,"height.m.median"]%>%round(1) ` m,
and Ornamental mean = `r tree.ml.mean[2,"height.m.median"]%>%round(1) ` m,
with Pr(>_F_) = `r mheight.tree[1,5]%>%round(4) `;
Native Shrub mean = `r shrub.ml.mean[2,"height.m.median"]%>%round(1) ` m,
and Ornamental mean = `r shrub.ml.mean[1,"height.m.median"]%>%round(1)` m,
with Pr(>_F_) = `r mheight.shrub[1,5]%>%round(3)`).
However, there was no difference in area between Native and Ornamental clusters for either trees or shrubs (
tree Pr(>_F_) = `r area.tree[1,5]%>%round(3)`;
shrub Pr(>_F_) = `r area.shrub[1,5]%>%round(3)`). 
Dead wood abundance was significantly greater on Native Tree sites than Ornamental Tree sites 
(Native Tree mean = `r tree.ml.mean[1,"dead.wood"]%>%round(1)`,
and Ornamental mean = `r tree.ml.mean[2,"dead.wood"]%>%round(1)`,
with Pr(>_F_) = `r dw.tree[1,5]%>%round(3)`),
but not between shrub groups. 

Only impervious surface cover between Native and Ornamental Tree sites differed significantly 
(Native Tree mean = `r 100*tree.gc.mean[1, "impervious.sqft"]%>%round(1)`,
and Ornamental mean = `r 100*tree.gc.mean[2, "impervious.sqft"] %>%round(1)`,
with Pr(>_F_) = `r imp.tree[1,5]%>%round(4)`).
No other ground covers differed.

There was also substantial co-occurence between Native and Ornamental groups. 
Of the 20 office developments surveyed, nine sites belong to both Native Tree and Shrub community groups, and seven sites belong to both Ornamental Tree and Shrub community groups.
This suggests that the sequential decisions made concerning tree preservation, tree plantings, and shrub plantings are related. 
The observed differences in species composition, between group differences, and  high turnover (beta diversity) support the conclusion that woody vegetation communities on office developments are heterogenous.

Native Tree and Shrub communities are more rare than Ornamental Tree and Shrub communities. Extrapolation suggests there are approximately 70 Native Tree and 335 Ornamental Tree developments (`r 7000/405 `%), and 152 Native Shrub and 253 Ornamental Shrub developments (`r 15200/405`%). 
The accuracy of these estimates is influenced by the Medium vegetation type, as it is large and proportionally under sampled, and the relatively small sample size.


## Socio-economic variables poorly explain variation in tree or shrub community composition

Neither aggregated neighborhood scale measures of residential socio-economic status nor parcel scale measures of economic value and the built environment explained variation in tree and shrub community composition on office developments following Holm-Bonferroni correction (S3 Table). 
For the tree community, median household income was significant only before correction.
This largely supports my hypothesis that aggregated socio-economic variables describing neighborhoods are not important for adjacent commercial properties as well [@hope2003; @leong2018], though additional research is needed.

Theoretically, for residential socio-economic variables to drive vegetation on office developments, the adjacent residential context must influence developer and commercial landowner vegetation choices.
Generally, zoning code in Bellevue and Redmond seeks to screen land uses from one another. 
Owners of office developments are likely signaling to prospective and existing tenants  [@laverne2003; @levy2013], in contrast with owners of residential properties, who use vegetation choices to signal to their neighbors of similar socio-economic status [@nassauer2009; @cook2012; @peterson2012]. 
Other explanations include proximity to desirable amenities–that is, both residential and commercial properties near amenities are appealing to more wealthy neighbors/tenants, the influence of city design review boards, or neighborhood overlay districts.

Studies examining why decision makers on commercial property make planting decisions are fewer in number than residential homeowners, though existing studies provide important early insight. 
For example, for landscape architects in Toronto factors like site aspect, appearance, and available space rated more highly in species selection than whether species are native or nearby canopy composition [@conway2016].

Parcel scale measures of economic value and the built environment were not significant, which fails to support my hypothesis. Previous research found that property value explained variation in the woody vegetation community [@mills2016]. 
Similarly, site age was suggested as a determinant of woody vegetation community composition by studies on residential properties [@boone2010; @avolio2018], landscaping professionals I interviewed, and my examination of contemporaneous landscaping plans filed with the cities of Bellevue and Redmond. Landscaping professionals mentioned trends in plant popularity, including _Pieris japonica_ (Thunb.) D. Don ex G. Don in the late 1980s and increasing use of native plants like _Ribes sanguineum_ since 2000. 
Alternative explanations for this finding include that building age is a poor measure for landscaping age due to replanting; an interaction between age and landscaping budget; or that a subset of office developments are planted with in vogue landscape plants, such as the common but under sampled Medium vegetation type. 

Differences in study design may also be responsible for these divergent results. 
Other studies use index response variables and univariate regression [@hope2003; @martin2004], measures dependent on effort [@karlik2001; @martin2004], and plot or transect designs which confound different actors and outcomes [@clarke2013; @bourne2013].


## Development and landscaping outcomes are related to tree and shrub community composition

Multiple variables describing the outcome of development and landscaping actions explain variation in tree and shrub community composition. For the tree community, NMDS with convex hulls and fitted environmental vectors found strong relationships with median dominant native conifer height, native conifer density, the presence of stands predating development, and dead wood abundance (particularly stump abundance, Fig \@ref(fig:FIG-tree-management-NMDS)). 
These variables were also included in the best supported PERMANOVA models for the shrub community (Table \@ref(tab:TAB-shrub-model-comparison)). 




```{r FIG-tree-management-NMDS, fig.cap = '(ref:NMDS)', dpi=300, fig.width=6.5, fig.height=5.5}

  spd.col <- ifelse(management.landscaping$stand.predate.development == "yes", "black", "gray70")
      


      par(mfrow = c(2,2), mar = c(2,3,2,1))
      
      temp <- ordisurf(NMDS.tree.2d ~ management.landscaping$height.m.median, plot = F)
        plot(temp, ylab = "NMDS2", xlab = "")
        points(NMDS.tree.2d,
               display = "sites", 
               cex = 1.2, pch = 16)
        ordihull(NMDS.tree.2d, 
                 groups = vegetation.clusters$tree.cluster.name,
                 col = "gray",
                 draw = "polygon", 
                 label = TRUE,
                 alpha = .2)
        #plot(envfit(NMDS.tree.2d, management.landscaping$height.m.median),  labels = "")
      title(main = "Median Domin. Native Conif. Ht. (m)")
      
      temp <- ordisurf(NMDS.tree.2d ~ management.landscaping$tree.conifer.nat.dens,
               plot = F)
        plot(temp, ylab = "NMDS2", xlab = "")
        points(NMDS.tree.2d,
               display = "sites",
               cex = 1.2,
               pch = 16)
        ordihull(NMDS.tree.2d,
                 groups = vegetation.clusters$tree.cluster.name,
                 col = "gray", draw = "polygon",
                 label = TRUE, alpha = .2)
        #plot(envfit(NMDS.tree.2d, management.landscaping$tree.conifer.nat.dens),  labels = "")
      title(main = "Native Conifer Density")

      temp <- plot(NMDS.tree.2d, ylab = "NMDS2",
           xlab = "NMDS1", display = "sites")
          points(NMDS.tree.2d, display = "sites",
                 col = spd.col, cex = 1.2, pch = 16)
          ordihull(NMDS.tree.2d, groups = vegetation.clusters$tree.cluster.name,
                   col = "gray", draw = "polygon",
                   label = TRUE, alpha = .2)
          #plot(envfit(NMDS.tree.2d, data.frame(as.factor(management.landscaping$stand.predate.development)),
          #           na.rm = T), labels = c("No SPD", "Yes SPD"))
      title(main = "Stands Predating Development")
      
      temp <-
          ordisurf(NMDS.tree.2d ~ management.landscaping$dead.wood,
                   plot = F)
          plot(temp, ylab = "NMDS2", xlab = "NMDS1")
          points(NMDS.tree.2d,
                 display = "sites",
                 cex = 1.2,
                 pch = 16)
          ordihull(NMDS.tree.2d,
                   groups = vegetation.clusters$tree.cluster.name,
                   col = "gray",
                   draw = "polygon",
                   label = TRUE,
                   alpha = .2
                  )
          # plot(envfit(NMDS.tree.2d, data.frame(management.landscaping$dead.wood, 
          #                                      management.landscaping$Stump.count,
          #                                      management.landscaping$Snag.count)),  
          #      labels = c("Total", "Stumps", "Snags"))
          title(main = "Dead Wood Abundance")


```

(ref:NMDS) Two dimensional NMDS representation of tree community composition. Median dominant native conifer height, native conifer density, and the presence of stands predating development are associated with the first NMDS axis. Dead wood is associated with both axes. Black dots represent sites with stands predating development, gray dots sites without. Ordination has not been rotated prior to plotting.


```{r TAB-shrub-model-comparison}
colnames(shrub.model.summary)[2] <- "AICc Value"
colnames(shrub.model.summary)[3] <- "Pseudo-F"

temp <- shrub.model.summary[ 1:10, c(6, 3, 4, 2, 7)] 
temp <- temp[ order(temp$`Pseudo-F`, decreasing = T), ]

knitr::kable(temp,
       digits = c(1,3,5,2,2),
       caption = "PERMANOVA model summary comparing multivariate models of shrub community composition.",
             format = "pandoc", escape = TRUE,
       format.args = list(decimal.mark = '.', big.mark = ","))

```


Together, these results support my hypothesis that development and landscaping actions impact vegetation communities on office developments. 
They also agree with some residential researchers who found that homeowner attitudes and actions were more important than socio-economic descriptors [@shakeel2014].
Together with flexible beta clustering results, this suggests that for each development a suite of decisions is made that results in either retaining more trees and planting native shrubs or retaining fewer trees and planting ornamental trees and shrubs.

However, development and landscaping outcomes are the end point of economic decision making processes poorly studied in urban ecology.
Though the coarse socio-economic variables examined here were not significant, developer and landowner motivations and decision making were not considered explicitly, only their outcomes. 
To reach these end points, developers may consider ease of construction based on site conditions, relative cost of different construction approaches, preferences of the landowner and customer specifications, previous company experience or company aesthetic, and development regulations [which impose costs on developers; @conway2016; @grimes2015; @hakkinen2011; @nappi2006; @dorney1984].
The intended audience of prospective and existing tenants may influence both development and landscaping decisions [@laverne2003; @levy2013].
These considerations may influence financing available to developers, financial risk, and the appeal of and thus demand for the completed project [@laverne2003].
Further, when considering multiple competing options---such as different landscaping choices---developers and landowners may satisfice [@mohamed2009].
That is, they search through alternatives until one meets an acceptability threshold and choose that option.

## Implications for urban habitat quality and quantity

The woody vegetation communities I observed on office developments suggest that integrating habitat conservation during and following development is possible using currently existing developments as models. 
As local biological communities are largely determined by vegetation, sites with more trees preserved and a greater abundance of native conifers likely provide higher habitat quality and quantity to other organisms [@wittig2010; @mckinney2002; @faeth2011; @avolio2018; @faeth2005] as native vegetation is more likely to support native insects and native birds than ornamental plantings [@belaire2014; @chong2014; @pennington2011; @paker2014; @narango2018; @burghardt2009; @dysonPhD].
One estimate suggests native vegetation volume must be above 70% in order to maintain populations of native insectivorous bird species [@narango2018]; sites with high numbers of trees preserved may already hit this target.

We can point to actions and policies more likely to support high quality habitat and benefit other trophic levels, including tree preservation policies, promoting native tree and shrub planting, and removing policy barriers to native vegetation [@dysonPhD; @threlfall2016; @leroux2014]. 
However, the motivations driving exemplary adoption of these actions are currently opaque. 
Anecdotes shared during fieldwork suggest owner-occupied office space, cost, and personal values and connections to nature may be important factors in determining development and landscaping actions, as with homeowners [@goddard2013; @nassauer1993; @kiesling2010; @kendal2012; @beumer2018; @helfand2006; @peterson2012].

## Implications for future urban ecology research

Observed within land use heterogeneity and between land use differences in socio-economic variable importance both have implications for urban ecology research. 
Within land use heterogeneity results in vegetation distributions that are non-normal, with likely kurtosis and heteroscedasticity (e.g. Figure \@ref(fig:FIG-extrapolated-dist)). 
Therefore, choice of sampling design and statistical method can result in inaccurate conclusions, particularly in conjunction with small sample size [@mcintyre2000urban]. 
Sampling across important gradients, as with vegetation composition in this study, is particularly important. 
Additionally, potential solutions already extant on the landscape may be overlooked.
This provides support for stratified sampling designs, larger sample sizes, and choosing analysis methods robust to broken assumptions of normality of the sampled population [e.g. @dewinter2013]. 

Researchers should choose their sampling strategy carefully based on research questions and the underlying distribution of key variables in urban contexts with long environmental gradients [@mcintyre2000urban; @ellis1997; @telford2011].
If the phenomenon of interest is related to the vegetation community, researchers should attempt to better understand and sample the vegetation gradient (e.g. via stratified sampling) instead of sampling only along a measure of the built gradient [e.g. housing density; @lerman2011]. 
Here, collecting vegetation information for the population of office developments prior to sampling improved sampling design.


```{r FIG-extrapolated-dist, fig.cap = '(ref:extrapolate)', dpi=300}

vc.mean <- mutate(sample.covariates,
                       tot.abund = rowSums(matrify.tree.sp.only),
                       tot.density = rowSums(dens.matrify.tree.sp.only)) %>%
    group_by(vegetation_class) %>%
    summarise_at(.vars = c("tot.abund"),
                 .funs = list(treeab.mean = mean, 
                              treeab.sd = sd))

pop.tree.extrap <- merge(population.covariates, vc.mean, 
                         by = "vegetation_class", sort = F)
# need to set the order of factors

pop.tree.extrap$vegetation_class <- factor(pop.tree.extrap$vegetation_class,
                                           levels = c("HH", "MC", "MD", "MM", "LL"))



library(ggplot2)

pop.tree.extrap$jittered <- jitter(pop.tree.extrap$treeab.mean,
                                   pop.tree.extrap$treeab.sd)
pop.tree.extrap$jittered[pop.tree.extrap$jittered < 1] <- 0

ggplot(pop.tree.extrap, aes(jittered)) +
  geom_histogram(aes(fill = vegetation_class), binwidth = 10) +
  scale_fill_manual(values = c("#008837", "#a6dba0", "#ffffbf", "#c2a5cf", "#7b3294"), name = "") +
  labs(x = "Tree Abundance per Site",
       y = "Number of Sites",
       title = "Tree Abundance on Office Developments") + scale_x_continuous(limits = c(-1, max(pop.tree.extrap$jittered) + 5))


```

(ref:extrapolate) Extrapolated distribution of the number of trees on office developments based on observed mean and standard deviations for each vegetation class used in sampling (HH, MC, MD, MM, LL). Note heavy right tail from HH, MC, and MD sites (kurtosis); each vegetation class also has a different variance (heteroscedasticity).

Between land use differences in socio-economic variable importance suggests that creating vegetation models of land use within a city is likely inaccurate if all land uses are assumed to respond equivalently to a given variable. 
Researchers cannot assume that vegetation gradients and socio-economic gradients are parallel; these gradients may also interact resulting in heteroscedasticity.
Decision pathways to support carbon sequestration and habitat models need to be constructed based on research for each land use separately.

# Conclusion

Humans control woody vegetation communities in urban ecosystems, influencing ecosystem service provision and habitat quality and quantity. Commercial land uses, including office developments, have largely been overlooked in studies of urban woody vegetation composition and studies examining how these communities are assembled.
I filled this gap by examining tree and shrub communities on office developments in Redmond and Bellevue, Washington, USA. 

I found that the vegetation communities on these developments are heterogenous, with distinct groups of sites characterized by Native and Ornamental Tree and Shrub vegetation communities. 
I also found that socioeconomic measures aggregated to the neighborhood scale and at the parcel scale were less important in explaining variation in community composition than variables describing the outcomes of development and landscaping choices.

This research contributes to our understanding of vegetation communities outside of municipal parks and residential land uses. 
It is also one of few studies that uses site surveys where the unit of measurement is based on how management decisions are made, instead of methods derived from wildlands vegetation research [e.g. transects and plots; @clarke2013; @bourne2013] or remote-sensing [e.g. @luck2009]. 

The observed heterogeneity in vegetation communities suggests that different ecosystem functions and habitat quantity and quality are provided on office developments.
Greater provision of these functions is possible using currently existing developments as models.
Further, within land use heterogeneity suggests that urban ecology research must more carefully consider sampling design, particularly sampling along key environmental gradients. The observed differences in variable importance between office developments and residential land uses suggests that future research and models of the urban ecosystem must account for land uses’ different decision pathways.

Going forward, researchers should examine other commercial land uses, commercial land use in additional ecotypes, and particularly the decision pathways followed by actors on commercial and other land uses.
This research agrees with other studies suggesting that specific actions are more important than aggregated socio-economic variables [@shakeel2014].
Additional research is needed to link decision makers' personal values and aesthetic preferences, economic motivations, and social norms with tree and shrub community composition on commercial land following work on residential property [ @cook2012; @shakeel2014].
Needed studies include interviews to better understand tree preservation and planting motivations [@hakkinen2011; @conway2016]; aesthetic preference studies as on residential developments [@larson2009; @harris2012]; and tracing decision making pathways based on previous land use [@yang2017].
A better understanding of these processes may improve habitat quality and quantity on commercial property [@uren2015]. 
Finally, research is also needed to determine if vegetation inequity observed on residential properties [@heynen2006] is perpetuated on commercial properties.
The No Vegetation type excluded from analysis here was often adjacent to retail use, where worker compensation is generally less than in medical/dental, software, and other white collar jobs in office developments.




# Supplemental Table 1: Trees

```{r TAB-all-trees}

# add common names, italicize??, asterisk for PNW natives


alltree <- merge(descriptive.stats.tree.sp, tree.native,
              by.x = "species.name", by.y = "tree.species") %>%
        arrange(., tree.scientific)

alltree$species.range <- paste0(alltree$species.min, "--", alltree$species.max)
alltree$tree.origin.update <- ifelse(alltree$tree.origin.update %in% c("non.native", "non-native"),
                           "Non-native", 
                           ifelse(alltree$tree.origin.update == "native",
                                  "Native",
                                  ifelse(alltree$tree.origin.update == "ambiguous",
                                         "Ambiguous",
                                         alltree$tree.origin.update)))



PA.tree.sp.only <- PA.tree.sp.only[ , match(alltree$species.name, colnames(PA.tree.sp.only))]
    #colnames(PA.tree.sp.only) == alltree$species.name

alltree$tree.incidence <- colSums(PA.tree.sp.only)


knitr::kable(alltree[ , c("tree.scientific",
                          "tree.common",
                          "tree.origin.update",
                          "tree.incidence",
                           "species.range",
                           "species.total.abundance",
                           "species.mean.abundance",
                           "species.abundance.sd")], 
             col.names = c("Taxa", 
                           "Common Name",
                           "Origin",
                           "No. Sites Found",
                           "Abundance Range",
                           "Total Abundance",
                           "Mean Abundance",
                           "Abundance SD"),
             row.names = FALSE,
             caption = "All trees observed in site surveys. Abundance is count of individuals belonging to each taxonomic group. Ambiguous indicate both native, non-native, and hybrids used in horticulture.",
             format = "pandoc", escape = TRUE, 
             format.args = list(decimal.mark = '.', big.mark = ",", digits = 2))

```

# Supplemental Table 2: Shrubs

```{r TAB-all-shrubs}

allshrub <- merge(descriptive.stats.shrub.sp, shrub.native,
              by.x = "species.name", by.y = "shrub.scientific.update") %>%
        arrange(., species.name) %>%
        distinct(species.name, .keep_all = TRUE)

allshrub$species.range <- paste0(allshrub$species.min, "--", allshrub$species.max)
allshrub$shrub.origin.update <- ifelse(allshrub$shrub.origin.update == "non-native",
                           "Non-native",
                           ifelse(allshrub$shrub.origin.update == "native",
                                  "Native",
                                  ifelse(allshrub$shrub.origin.update == "ambiguous",
                                         "Ambiguous",
                                         allshrub$shrub.origin.update)))

frequency.shrub <- tibble(species.name = names(colSums(PA.shrub.site)),
                          frequency.shrub = colSums(PA.shrub.site) )

#frequency.shrub$species.name == allshrub$species.name


allshrub$incidence <- frequency.shrub$frequency.shrub


knitr::kable(allshrub[ c(1:72,80:85,73:79),
                       c("species.name",
                         "shrub.common.pretty",
                         "shrub.origin.update",
                         "incidence",
                         "species.range",
                      "species.total.abundance",
                      "species.mean.abundance",
                      "species.abundance.sd")],
             col.names = c("Taxa",
                           "Common Name",
                           "Origin",
                           "No. Sites Found",
                           "Abundance Range",
                           "Total Abundance",
                           "Mean Abundance",
                           "Abundance SD"),
             row.names = FALSE, 
             caption = "All shrubs observed in site surveys. Abundance is count of individuals belonging to each taxonomic group. Ambiguous indicate both native, non-native, and hybrids used in horticulture.",
             format = "pandoc", escape = TRUE, 
             format.args = list(decimal.mark = '.', big.mark = ","))
```

# Supplemental Table 3: All PERMANOVA Results

```{r TAB-all-PERMANOVA}

# need to include Pseudo-F, p-value, AICc

AICc.all.shrub$Model <- c("Area (acres)", "Percent Foreign-Born", "Town", "Median Income", "Impervious w/in 500 m (%)", "Building Age (in 2017)", "Building Quality", "Appraised Land Value per Acre", "Small and Medium Height Vegetation Cover", "Tall Vegetation Cover", "Tree cluster", "Dead Wood (count)", "Median Height of Dominant Conifer", "Stands Predate Development", "Density of Native COnifers") 
AICc.all.shrub$`Adjusted p-val` <- p.adjust(AICc.all.shrub$`p-value`, "holm")

AICc.all.tree$Model <- c("Area (acres)", "Percent Foreign-Born", "Town", "Median Income", "Impervious w/in 500 m (%)", "Building Age (in 2017)", "Building Quality", "Appraised Land Value per Acre", "Small and Medium Height Vegetation Cover", "Tall Vegetation Cover", "Dead Wood (count)")
AICc.all.tree$`Adjusted p-val` <- p.adjust(AICc.all.tree$`p-value`, "holm")


AICc.all.tree <- rbind(AICc.all.tree, rep(NA,9), rep(NA,9), rep(NA,9), rep(NA,9))


socio.economic.results <- cbind(AICc.all.shrub[ , c(6,3,4,9,2,7)], AICc.all.tree[ , c(3,4,9,2,7)])
socio.economic.results[ , 2:10] <- round(socio.economic.results[ , 2:10], 3)

colnames(socio.economic.results) <- c("Model", "Shrub Pseudo-_F_","Shrub p-value",
                                      "Shrub Adjusted p-val", "Shrub AICc","Shrub Delta AICc",
                                      "Tree Pseudo-_F_","Tree p-value",
                                      "Tree Adjusted p-val", "Tree AICc","Tree Delta AICc")
options(knitr.kable.NA = '')
knitr::kable(socio.economic.results,
             caption = "All simple multivariate PERMANOVA results for tree and shrub communities.",
             format = "pandoc", escape = TRUE, 
             digits = 2, 
             format.args = list(decimal.mark = '.', big.mark = ","))

```



# Works Cited