# Part 4: Deconstruction Wood Sensitivity Analysis

### 2018 Deconstruction Data Net Impact Analysis using OSU average softwood density factor

#### Andey Nunes, MS

#### Research Analyst 2

#### Oregon DEQ

#### October 22, 2018

# Introduction, Project Goals, and Notebook Setup

The 2018 Deconstruction Data Analysis uses prepared and cleaned project data from residential single family homes removed under a City of Portland Deconstruction permit. The goal of the data analysis is to quantify the net environmental benefits resulting from avoided disposal of materials due to salvage/reuse (measured as Global Warming Potential and Primary Energy Demand impacts). This document is part four of the project where previously constructed data sets are used to determine and graphically summarize impacts under the deconstruction and demolition scenario and the net benefits of deconstructing residential homes over the hypothetical demolition of those structures.

The analysis is intended to be reproducible using R coding procedures and employs the following R Packages:

package name
corrplot
ggpubr
ggrepel
ggthemes
grDevices
knitr
tidyverse

Prior to executing code and producing a report, the packages must be installed and accessed, and basic features of the document set up by specifying code chunk defaults.

    ## Loading required package: corrplot

    ## corrplot 0.84 loaded

    ## Loading required package: ggpubr

    ## Loading required package: ggplot2

    ## Loading required package: magrittr

    ## Loading required package: ggrepel

    ## Loading required package: ggthemes

    ## Loading required package: knitr

    ## Loading required package: tidyverse

    ## -- Attaching packages ----------------------------------------------------------------------------------------------------------------------- tidyverse 1.2.1 --

    ## v tibble  1.4.2     v purrr   0.2.5
    ## v tidyr   0.8.1     v dplyr   0.7.6
    ## v readr   1.1.1     v stringr 1.3.1
    ## v tibble  1.4.2     v forcats 0.3.0

    ## -- Conflicts -------------------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
    ## x tidyr::extract()   masks magrittr::extract()
    ## x dplyr::filter()    masks stats::filter()
    ## x dplyr::lag()       masks stats::lag()
    ## x purrr::set_names() masks magrittr::set_names()

    ## [[1]]
    ## [1] TRUE
    ##
    ## [[2]]
    ## [1] TRUE
    ##
    ## [[3]]
    ## [1] TRUE
    ##
    ## [[4]]
    ## [1] TRUE
    ##
    ## [[5]]
    ## [1] TRUE
    ##
    ## [[6]]
    ## [1] TRUE
    ##
    ## [[7]]
    ## [1] TRUE

Specify parameters for saving graphics.

    # designate a folder name for saving graphics and the default size specification
    graph_path <- file.path("graphs/")
    width <- 8
    height <- 5

## Data Import & consolidation

Previously developed data sets are imported and cleaned/manipulated as needed, the index values that R assigns when saving csv files is dropped, units of the values are formatted to be consistent across all data tables, and contractor names are replaced with a random number for anonymity in the final report.

For worker and equipment impacts that are associated with projects but not a specific material category (which is what all observations in our `decon` and `demo` data sets are up to this point), we have to assign compatible names within the existing variable categories in order to create observations of each per-project based transport impact (i.e. workers and equipment). For this the `worker_equipment_impacts` data values will be assigned the `site_transport_impact` variable name, since job-site related transportion based impacts is the best fit description for the inclusion of these values. Additional names are given in other variables in order to distinguish these values from the material impacts and material related transport impacts (including site and EOL based transport). The `worker_equipment_impacts` table will have the following variables and assigned variable states (values) added prior to being combined with the `decon` and `demo` scenario data frames:

variable name | value (assigned variable state)
--- | ---
`decon_equipment` | 0
`SimpleEOLname` | per_project_transport
`disposition` | none
`quantity` | 1
`quantity_units` | per_project
`material_impacts` | 0
`EOL_transport_impact` | 0

    # cannot use read_csv for reading in these files because of parsing error on some of the EOL transport dropbox values
    decon <- read.csv("intermediary/final_deconstruction_scenario_with_transport.csv", stringsAsFactors = F)
    demo <- read.csv("intermediary/final_demolition_scenario_with_transport.csv", stringsAsFactors = F)
    worker_equipment_impacts <- read.csv("intermediary/per_project_worker_equipment_transport_impacts.csv", stringsAsFactors = F)
    house_weight <- read.csv("intermediary/house_weight.csv", stringsAsFactors = F)
    house_weight$project <- as.character(house_weight$project)
    house_weight$contractor <- as.character(house_weight$contractor)

    # standardize impact_units
    decon$impact_units <- str_replace(decon$impact_units, "kgCO2e", "kg CO2e")

    # clear index columns and add the activity scenario name
    decon <- decon[,-1] %>%
       add_column(scenario = "decon", .before = 1)

    demo <- demo[,-1] %>%
       add_column(scenario = "demo", .before = 1)

    # adjust worker_equipment_impacts to fit the structure of the decon and demo activity scenario data frames
    worker_equipment_impacts <- worker_equipment_impacts %>%
       add_column(decon_equipment = 0, .before = 7) %>%
       add_column(SimpleEOLname = "per_project_transport") %>%
       add_column(disposition = "none") %>%
       add_column(quantity = 1) %>%
       add_column(quantity_units = "per_project") %>%
       add_column(material_impacts = 0) %>%
       add_column(EOL_transport_impact = 0) %>%
       gather(`decon_workers`, `decon_equipment`, `demo_workers`, `demo_equipment`, key = "scenario_deconMaterialName", value = "site_transport_impact") %>%
       separate(scenario_deconMaterialName, into = c("scenario", "deconMaterialName"), sep = "_") %>%
       select(c("scenario", "house_age", "house_size", "project", "contractor", "deconMaterialName", "SimpleEOLname", "disposition", "quantity", "quantity_units", "material_impacts","site_transport_impact", "EOL_transport_impact", "impact_units"))

    # join all 3 data frames to create a master data set that contains all material and transport impacts from both scenarios and save it
    final_complete_data <- bind_rows(decon, demo, worker_equipment_impacts)
    final_complete_data$project <- as.character(final_complete_data$project)
    final_complete_data$contractor <- as.character(final_complete_data$contractor)

    write.csv(final_complete_data, "intermediary/final_complete_data_all_scenario_impacts.csv")

    # in some instances, it may be desirable to see material impacts of either just the Softwood lumber products, or everything else minus the Softwood lumber products
    softwood_lumber <- final_complete_data %>%
       filter(deconMaterialName == "softwood lumber")

    nonSoftwood_lumber <- anti_join(final_complete_data, softwood_lumber, by = c("scenario", "house_age", "house_size", "project", "contractor", "deconMaterialName", "SimpleEOLname", "disposition", "quantity", "quantity_units", "material_impacts", "site_transport_impact", "EOL_transport_impact", "impact_units"))

    # remove comment to write/save these two tables separately as needed
    # write.csv(softwood_lumber, "intermediary/softwood lumber data.csv")
    # write.csv(nonSoftwood_lumber, "intermediary/nonSoftwood lumber data.csv")

The resulting `worker_equipment_impacts` data frame has 288 rows and the same 14 varible field columns as the `decon` and `demo` data tables. Now these three data tables (with exactly same column order and column names) can be combined into a single `final_complete_data` table and save as `final_complete_data_all_scenario_impacts.csv`. This final table contains 4308 rows of observations, and retains information specific to project, material weight quantity and material type category, disposition, impact category (designated by the `impact_units` variable), and impact values associated with the materials, site-transport (which includes worker and equipment per-project-transport impact values), and EOL-transport. This structure allows for the summary and break out of material based impacts, and all forms of transport impacts.

## Data Set Description

    ave_house_size <- summarise(house_weight, avg_sq_ft = mean(house_size))
    total_house_size <- summarise(house_weight, total_sq_ft = sum(house_size)) # total sq_Ft from all projects

FIGURE X: Relation of disposal material weight to house size

There was request from the project team to explore the weight ratio of salvaged materials to the dropbox disposal. This is depicted in FIGURE 3 where the dropbox total weight is on the x-axis, the salvaged material weight is on the y-axis, each project is a point along the x and y where the size of the point is relative to the size of the house, contractors are assigned to the point color, and each point is labeled with its respective percent salvaged material numerical value. In this situation, given the wide range and variation in the material weight values from each project, the median value is a more appropriate measure of central tendency to express an average value because it is not influenced by extreme values. For this reason, median values are included in the plot in order to provide context for the data points in terms of the percent salvaged materials and the total salvage and disposal weights. The median crosshairs have been added to visually assist in differentiating the better than average performers in the upper left quadrant (with higher than average salvage weight and lower than average dropbox weight), from the lower than average performers in the lower right quadrant (lower than average salvage weight and higher than average dropbox weight)

    ggplot(house_weight, aes(total_dropbox_quantity, total_salvage_quantity, color = contractor, size = house_size)) +
       geom_jitter() +
       #scale_size("house_size (sqft)") +
       geom_text_repel(aes(label = as.integer(percent_salvaged)), seed = 19, show.legend = F) +
       #geom_rug(aes(size = 0.25), sides = "bl", position = "jitter", show.legend = F) +
       geom_vline(xintercept = median(house_weight$total_dropbox_quantity), color = "grey30" ) +
       geom_hline(yintercept = median(house_weight$total_salvage_quantity), color = "grey30" ) +
       scale_radius("house size (sqft)", range = c(2,6)) +
       scale_color_calc() +
       labs(x = "kg dropbox materials", y = "kg salvaged materials") #+

       #ggtitle("Ratio of salvaged materials to dropbox materials by house size")

    #ggsave(filename = "ratio_salvage_dropbox.png", device = "png", path = graph_path, units = "in")
    # no need to save it here since the next document produces and saves the same graphic

FIGURE _ : Salvaged weight to dropbox weight ratio

# Methods

The work process for the project is illustrated by the following figure. Data frames developed from the cleaned raw data summary (`01_decon_material_weight_conversion.Rmd`) and the data preparation notebook (`02_decon_data_prep2018.Rmd`) are used with the impact factors from LCA data and transport data collected from Oregon DEQ and Oregon Metro (`03_Transport.Rmd`) in the development of this data analysis notebook. This notebook provides a detailed account of the analysis and is used to produce custom R Markdown reports for different stakeholders based on the stakeholder information needs.

Project Workflow

## Data Variable Distribution tests

Where we are interested in statistics of the variables in the data, the distributions of the data must be tested for assumptions before we can apply linear modeling methods. The Shapiro-Wilk’s method for normality test is used on the `house_size`, `house_age`, and salvage and disposal weight variables.

    shapiro.test(house_weight$house_age)

    ##
    ##  Shapiro-Wilk normality test
    ##
    ## data:  house_weight$house_age
    ## W = 1, p-value = 0.8

    shapiro.test(house_weight$house_size)

    ##
    ##  Shapiro-Wilk normality test
    ##
    ## data:  house_weight$house_size
    ## W = 0.9, p-value = 0.01

    shapiro.test(house_weight$total_salvage_quantity)

    ##
    ##  Shapiro-Wilk normality test
    ##
    ## data:  house_weight$total_salvage_quantity
    ## W = 1, p-value = 0.2

    shapiro.test(house_weight$total_dropbox_quantity)

    ##
    ##  Shapiro-Wilk normality test
    ##
    ## data:  house_weight$total_dropbox_quantity
    ## W = 1, p-value = 0.2

    ggqqplot(house_weight$house_age) + ggtitle("house_age qqplot for normality")

    ggqqplot(house_weight$house_size) + ggtitle("house_size qqplot for normality")

    ggqqplot(house_weight$total_salvage_quantity) + ggtitle("total_salvage_quantity qqplot for normality")

    ggqqplot(house_weight$total_dropbox_quantity) + ggtitle("total_dropbox_quantity qqplot for normality")

# Results

## Data Variable Relations and Linear Correlations

    # add lm of dropbox total kg ~ house size w/ r^2
    # house_weight$contractor <- as_factor(house_weight$contractor)
    house_in_dropbox <- lm(total_dropbox_quantity ~ house_size, data = house_weight)
    summary(house_in_dropbox)

    ##
    ## Call:
    ## lm(formula = total_dropbox_quantity ~ house_size, data = house_weight)
    ##
    ## Residuals:
    ##    Min     1Q Median     3Q    Max
    ##  -8625  -5423    557   3695  11647
    ##
    ## Coefficients:
    ##              Estimate Std. Error t value Pr(>|t|)
    ## (Intercept) 14068.973   2968.785    4.74 0.000037 ***
    ## house_size     -0.864      2.400   -0.36     0.72
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Residual standard error: 5490 on 34 degrees of freedom
    ## Multiple R-squared:  0.0038, Adjusted R-squared:  -0.0255
    ## F-statistic: 0.13 on 1 and 34 DF,  p-value: 0.721

    #plot(house_in_dropbox, 1:6)

    house_salvaged <- lm(total_salvage_quantity ~ house_size + house_age, data = house_weight)
    summary(house_salvaged)

    ##
    ## Call:
    ## lm(formula = total_salvage_quantity ~ house_size + house_age,
    ##     data = house_weight)
    ##
    ## Residuals:
    ##    Min     1Q Median     3Q    Max
    ##  -3506  -1437    -34   1424   4361
    ##
    ## Coefficients:
    ##             Estimate Std. Error t value Pr(>|t|)
    ## (Intercept) 2255.227   3762.597    0.60    0.553
    ## house_size     2.381      0.878    2.71    0.011 *
    ## house_age     -2.270     33.724   -0.07    0.947
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ## Residual standard error: 1980 on 33 degrees of freedom
    ## Multiple R-squared:  0.186,  Adjusted R-squared:  0.136
    ## F-statistic: 3.76 on 2 and 33 DF,  p-value: 0.0338

    #plot(house_salvaged, 1:6)

    num_house_weight <- house_weight %>%
       select("total_salvage_quantity", "total_dropbox_quantity", "house_age", "house_size")
    names(num_house_weight) <- c("salvage kg", "dropbox kg", "house age", "house size")

    house_weight_corr <- cor(num_house_weight, method = "pearson")
    print(house_weight_corr)

    ##            salvage kg dropbox kg house age house size
    ## salvage kg       1.00      -0.19      0.07       0.43
    ## dropbox kg      -0.19       1.00     -0.06      -0.06
    ## house age        0.07      -0.06      1.00       0.18
    ## house size       0.43      -0.06      0.18       1.00

    cor.test(house_weight$total_salvage_quantity, house_weight$house_size, method = "pearson")

    ##
    ##  Pearson's product-moment correlation
    ##
    ## data:  house_weight$total_salvage_quantity and house_weight$house_size
    ## t = 3, df = 30, p-value = 0.009
    ## alternative hypothesis: true correlation is not equal to 0
    ## 95 percent confidence interval:
    ##  0.1 0.7
    ## sample estimates:
    ## cor
    ## 0.4

    cor.test(house_weight$total_dropbox_quantity, house_weight$house_size, method = "pearson")

    ##
    ##  Pearson's product-moment correlation
    ##
    ## data:  house_weight$total_dropbox_quantity and house_weight$house_size
    ## t = -0.4, df = 30, p-value = 0.7
    ## alternative hypothesis: true correlation is not equal to 0
    ## 95 percent confidence interval:
    ##  -0.4  0.3
    ## sample estimates:
    ##   cor
    ## -0.06

    corrplot::corrplot(house_weight_corr, method = "color", type = "lower", addCoef.col = "black", tl.col = "black", tl.srt = 0, sig.level = 0.05, insig = "blank", diag = FALSE, mar = c(0,0,1,0))

House size is not a valid predictor of the total dropbox quantity since only a fraction of a percent of the variance in the house sizes explains variance in the dropbox quantity and the correlation is not statistically significant. A closer look at a correlation matrix between the salvage and dropbox quantities and the house characteristics reveals a modest positive correlation of 0.41 between the square foot size of the house and the salvage quantity, which when tested, turns out to be a statistically significant correlation at the 0.01 alpha level (a 0.05 p-value is conventionally accepted as an indicator of statistical significance).

For the following net impacts calculations, each impact is summarised by impact category and activity. The sets of `net` impacts are obtained by subtracting the `decon` activity from the `demo` activity: \[demo - decon = net\].

The `net` impacts are then expressed as a `net_per_sq_ft` value, which is `net` divided by the sum total of square feet of the projects (42365) and a `net_ave_project` value, which is `net` divided by the total number of projects (36): \[\frac{demo - decon} {total ~~ sq ~~ ft} = \frac{net}{sq~ft}\] \[\frac{demo - decon} {total ~~ number ~~ of ~~ projects} = average ~~ project ~~ net ~~ benefit\]

These calculations are completed for each impact category and saved as separate .csv files.

    # for the final data summaries, Jordan will want to see this broken out from the site transport leg
    # # create a table showing demo and decon worker and equipment total impacts and later add per sq_ft and per project

    worker_equipment_summary <- worker_equipment_impacts %>%
       group_by(scenario, deconMaterialName, impact_units) %>%
       summarise(impact = sum(site_transport_impact)) %>%
       rename(impact_origin = deconMaterialName)

Reshape the `final_complete_data` table for use with graphics and table construction & final results development.

    gathered_impacts <- final_complete_data %>%
       gather(`material_impacts`, `site_transport_impact`, `EOL_transport_impact`, key = "impact_origin", value = "impact")

Calculate net impacts by impact category and save tables.

    net_impacts <- filter(gathered_impacts, SimpleEOLname != "per_project_transport") %>%
       group_by(scenario, impact_units, impact_origin) %>%
       summarise(impact = sum(impact)) %>%
       full_join(worker_equipment_summary, by = c("scenario", "impact_units", "impact_origin", "impact"))

    # add the per house column now and export data table for the report draft
    # per comments from Jordan, the next several code chunks showing net impacts will not be used,
    # instead, make new tables showing impacts by scenario per average home with net value
    final_impact_summary <- net_impacts %>%
       mutate(per_home_impact = impact / length(house_weight$project))

    write.csv(final_impact_summary, "intermediary/final_impact_summary.csv")

Carbon net impact summary:

    carbon_net_impacts <- filter(net_impacts, impact_units == "kg CO2e") %>%
       ungroup() %>%
       spread(key = scenario, value = impact)

    carbon_net_impacts <- carbon_net_impacts %>%
       add_row(impact_units = "kg CO2e",
               impact_origin = "total",
               decon = sum(carbon_net_impacts$decon),
               demo = sum(carbon_net_impacts$demo)) %>%
       mutate(net = demo - decon) %>%
       mutate(net_per_sqft = net/as.numeric(total_house_size)) %>%
       mutate(net_per_project = net/length(house_weight$project)) %>%
       add_column(order = c(3,5,1,2,4,6)) %>%
       arrange(order) %>%
       select("impact_origin", "decon",  "demo", "net", "net_per_sqft", "net_per_project", "impact_units")

    write.csv(carbon_net_impacts, "output/carbon_impact_breakdown.csv")

    kable(carbon_net_impacts, caption = "Carbon Impact Table")

Carbon Impact Table
impact_origin | decon | demo | net | net_per_sqft | net_per_project | impact_units
--- | --- | --- | --- | --- | --- | ---
material_impacts | -543758 | -258315 | 285443 | 6.7 | 7929 | kg CO2e
site_transport_impact | 6363 | 4244 | -2120 | -0.1 | -59 | kg CO2e
EOL_transport_impact | 16621 | 21045 | 4424 | 0.1 | 123 | kg CO2e
workers | 22216 | 2222 | -19995 | -0.5 | -555 | kg CO2e
equipment | 0 | 8683 | 8683 | 0.2 | 241 | kg CO2e
total | -498558 | -222122 | 276436 | 6.5 | 7679 | kg CO2e

Energy net impact summary:

    energy_net_impacts <- filter(net_impacts, impact_units == "MJ") %>%
       ungroup() %>%
       spread(key = scenario, value = impact)

    energy_net_impacts <- energy_net_impacts %>%
       add_row(impact_units = "MJ",
               impact_origin = "total",
               decon = sum(energy_net_impacts$decon),
               demo = sum(energy_net_impacts$demo)) %>%
       mutate(net = demo - decon) %>%
       mutate(net_per_sqft = net/as.numeric(total_house_size)) %>%
       mutate(net_per_project = net/length(house_weight$project)) %>%
       add_column(order = c(3,5,1,2,4,6)) %>%
       arrange(order) %>%
       select("impact_origin", "decon",  "demo", "net", "net_per_sqft", "net_per_project", "impact_units")

    write.csv(energy_net_impacts, "output/energy_impact_breakdown.csv")

    kable(energy_net_impacts, caption = "Energy Impact Table")

Energy Impact Table
impact_origin | decon | demo | net | net_per_sqft | net_per_project | impact_units
--- | --- | --- | --- | --- | --- | ---
material_impacts | -3822482 | -4674200 | -851718 | -20.1 | -23659 | MJ
site_transport_impact | 93026 | 62050 | -30976 | -0.7 | -860 | MJ
EOL_transport_impact | 243043 | 307729 | 64686 | 1.5 | 1797 | MJ
workers | 282635 | 28264 | -254371 | -6.0 | -7066 | MJ
equipment | 0 | 126576 | 126576 | 3.0 | 3516 | MJ
total | -3203779 | -4149582 | -945803 | -22.3 | -26272 | MJ

## Salvaged Materials Overview

## Transport Factors and Equipment Use

## Contribution Analysis

In order to show material flows and their impact contributions, the gathered data set needs to be summarised.

Then for each material, material type, impact category, and origin of impacts from materials and various transport factors, the sum of the impact values can be displayed graphically to show the specific levels of impact contribution to the total impact value.

    ggplot(filter(ggimpact, impact_units == "MJ"), aes(x = deconMaterialName, y = impact_sum, fill = SimpleEOLname) ) +
       geom_bar(stat = "identity", position = "stack") +
       geom_text_repel(aes(label = as.integer(impact_sum) )) +
       scale_fill_manual(values = DEQ_pal[5:9]) +
       coord_flip() +
       ggtitle("Energy impacts (MJ) by origin material type") +
       facet_grid(impact_origin ~ scenario)

    ggplot(filter(ggimpact, impact_units == "kg CO2e"), aes(x = deconMaterialName, y = impact_sum, fill = SimpleEOLname) ) +
       geom_bar(stat = "identity", position = "stack") +
       geom_text_repel(aes(label = as.integer(impact_sum) )) +
       scale_fill_manual(values = DEQ_pal[5:9]) +
       coord_flip() +
       ggtitle("Carbon impacts (kg CO2e) by origin material type") +
       facet_grid(impact_origin ~ scenario)

    # use the next two sets of ggplot calls to experiment with creating a graphical chart of all impact origins & material types
    ggplot(filter(ggimpact, impact_units == "MJ"), aes(x = deconMaterialName, y = impact_sum, fill = impact_origin) ) +
       geom_bar(stat = "identity", position = "stack") +
       geom_text_repel(aes(label = as.integer(impact_sum) )) +
       scale_fill_manual(values = DEQ_pal[12:14]) +
       coord_flip() +
       ggtitle("Energy impacts (MJ) by origin material type") +
       facet_grid(SimpleEOLname ~ scenario )

    ggplot(filter(ggimpact, impact_units == "kg CO2e"), aes(x = deconMaterialName, y = impact_sum, fill = impact_origin) ) +
       geom_bar(stat = "identity", position = "stack") +
       geom_text_repel(aes(label = as.integer(impact_sum) )) +
       scale_fill_manual(values = DEQ_pal[12:14]) +
       coord_flip() +
       ggtitle("Carbon impacts (kg CO2e) by origin material type") +
       facet_grid(SimpleEOLname ~ scenario)