In Blocktt/MassIBItools: Used to calculate Massachusetts DEP IBI Metrics and scores, and produce output files

knitr::opts_chunk$set(results='asis', echo=FALSE, warning=FALSE, message = FALSE)
# needed for trouble shooting
boo_DEBUG <- FALSE
if(boo_DEBUG==TRUE){
  # myConfig <- file.path(system.file(package="ContDataQC"), "extdata", "config.ORIG.R")
  # source(myConfig)
}## IF ~ boo_DEBUG ~ END

INTRODUCTION

The following document is dynamic - it changes depending on the data input into the Mass DEP R Shiny Application.

Please refer to the "Data File Information" section for general properties of the input dataset used to genereate this report.

Histograms in the "Plots" section will change depending on the input dataset.

DATA FILE INFORMATION

  library(dplyr)
  # Results File
  df_metval <- df_metval
  df_metsc <- df_metsc

  # df_metval <- read.csv("C:/Users/Ben.Block/OneDrive - Tetra Tech, Inc/GitHub/MassIBI_results_metval.csv")
  # df_metsc <- read.csv("C:/Users/Ben.Block/OneDrive - Tetra Tech, Inc/GitHub/MassIBI_results_metsc.csv")

  # create general index region names
  df_metval$Region_Name <- ifelse(df_metval$INDEX_REGION == "KickIBI_CH_300ct", "CH",
                     ifelse(df_metval$INDEX_REGION == "KickIBI_CH_100ct", "CH",
                            ifelse(df_metval$INDEX_REGION == "KickIBI_WH_300ct", "WH",
                                   ifelse(df_metval$INDEX_REGION == "KickIBI_WH_100ct", "WH",
                                          "LowGrad"))))
  df_metsc$Region_Name <- ifelse(df_metval$INDEX_REGION == "KickIBI_CH_300ct", "CH",
                     ifelse(df_metval$INDEX_REGION == "KickIBI_CH_100ct", "CH",
                            ifelse(df_metval$INDEX_REGION == "KickIBI_WH_300ct", "WH",
                                   ifelse(df_metval$INDEX_REGION == "KickIBI_WH_100ct", "WH",
                                          "LowGrad"))))

  # Report Info
  myReportDate <- format(Sys.time(), '%Y-%m-%d %H:%M:%S')
  cat(paste("**Report Date:** " ,myReportDate, "\n\n", sep=""))
  myUser <- Sys.getenv("USERNAME")
  cat(paste("**Generated By:** ", myUser, "\n\n", sep=""))

  # Number of Samples
  myNumRecords <- nrow(df_metval)
  cat(paste("**Number of Records:** ", myNumRecords, "\n\n", sep=""))

  # Period of Record
  df_metval$COLLDATE <- as.Date(df_metval$COLLDATE, format = "%m/%d/%Y")
  POR <- paste(min(df_metval$COLLDATE)," to ", max(df_metval$COLLDATE), sep="")
  cat(paste("**Period of Record:** ", POR, sep="", collapse="\n\n"))

SUMMARY TABLES

The tables below summarize the data that were input into the shiny app.

Table 1. Number of sites by Index Narrative and Index Region.

library(knitr)
library(tidyr)
library(dplyr)
library(flextable)

options(dplyr.summarise.inform = FALSE)

df_scores <- df_metsc
df_scores$counts <- rep(1)

# create a vector in the desired order
row_ord <- c("Exceptional", "Satisfactory", "Moderately Degraded", "Severely Degraded")

table1 <- df_scores %>% 
  group_by(Index_Nar, INDEX_REGION) %>% 
  summarize(counts = sum(counts)) %>% # renames Index_Region to nice names
  mutate(INDEX_REGION = replace(INDEX_REGION, INDEX_REGION == "KickIBI_CH_100ct", "Central Hills (100 ct)")) %>% 
  mutate(INDEX_REGION = replace(INDEX_REGION, INDEX_REGION == "KickIBI_CH_300ct", "Central Hills (300 ct)")) %>% 
  mutate(INDEX_REGION = replace(INDEX_REGION, INDEX_REGION == "KickIBI_WH_100ct", "Western Highlands (100ct)")) %>% 
  mutate(INDEX_REGION = replace(INDEX_REGION, INDEX_REGION == "KickIBI_WH_300ct", "Western Highlands (300ct)")) %>% 
  mutate(INDEX_REGION = replace(INDEX_REGION, INDEX_REGION == "LowGradientIBI", "Low Gradient")) %>% 
  pivot_wider(names_from = INDEX_REGION,
              values_from = counts) %>% 
  mutate(Index_Nar = factor(Index_Nar, levels = row_ord)) %>% 
  arrange(Index_Nar) %>% 
  rename("Index Narrative" = Index_Nar)

# create flextable for output 
# https://davidgohel.github.io/flextable/articles/layout.html
ftable1 <- flextable(table1)
ftable1 <- autofit(ftable1)
ftable1 <- set_table_properties(ftable1, width = .5, layout = "autofit")
theme_booktabs(ftable1)

Table 2. Index scores summary statistics by Index Region. Rows that contain "Stress Sites" and "Reference Sites" are static and always appear regardless of input data set. Essentially, you can use these rows to compare the results of your input dataset to the reference/stress sites used in calibration of each IBI. Note, do not compare across Index Regions as this is not a correct interpretation (e.g., do not compare CH to WH or Low Gradient). Rather, only compare inputs to their equivalent index region results (e.g., CH to CH Stress/Reference Sites).

library(knitr)
library(tidyr)
library(readxl)
library(flextable)

scores_trim <- df_scores %>% 
  select(INDEX_REGION, Index, counts)

table2_temp <- scores_trim %>% 
  group_by(INDEX_REGION) %>% 
  summarize(N = sum(counts),
            Min = round(min(Index), 1),
            q25 = round(quantile(Index, 0.25, na.rm = TRUE), 1),
            Med = round(median(Index, na.rm = TRUE), 1),
            q75 = round(quantile(Index, 0.75, na.rm = TRUE), 1),
            Max = round(max(Index),1))%>% 
  rename("Index Region" = INDEX_REGION)

# state directories
table.dir <- "tables"
table.file <- "Instruction_Tables.xlsx"
IBI_tab.dir <- "IBI_data"

IBI_data <- read_excel(file.path(table.dir, table.file), sheet = IBI_tab.dir
                     , na = c("NA", ""), trim_ws = TRUE, skip = 0
                     , col_names = TRUE)
table2 <- rbind(table2_temp, IBI_data)

table2 <- table2 %>% # renames Index_Region to nice names
  mutate(`Index Region` = replace(`Index Region`, `Index Region` == "KickIBI_CH_100ct", "Cent. Hills (100 ct) - Input")) %>% 
  mutate(`Index Region` = replace(`Index Region`, `Index Region` == "KickIBI_CH_300ct", "Cent. Hills (300 ct) - Input")) %>% 
  mutate(`Index Region` = replace(`Index Region`, `Index Region` == "KickIBI_WH_100ct", "West. High. (100ct) - Input")) %>% 
  mutate(`Index Region` = replace(`Index Region`, `Index Region` == "KickIBI_WH_300ct", "West. High. (300ct) - Input")) %>% 
  mutate(`Index Region` = replace(`Index Region`, `Index Region` == "LowGradientIBI", "Low Gradient - Input"))

# create flextable for output 
# https://davidgohel.github.io/flextable/articles/layout.html
ftable2 <- flextable(table2)
ftable2 <- autofit(ftable2)
theme_booktabs(ftable2)

Table 3. Metrics in the Western Highlands 100-count IBI, with scoring formulas. Trend is the direction of metric response with increasing stress.

library(readxl)
library(knitr)
library(flextable)

# state directories
table.dir <- "tables"
table.file <- "Instruction_Tables.xlsx"
tab3.dir <- "Back_table1"

table3 <- read_excel(file.path(table.dir, table.file), sheet = tab3.dir
                     , na = c("NA", ""), trim_ws = TRUE, skip = 0
                     , col_names = TRUE)

# create flextable for output 
# https://davidgohel.github.io/flextable/articles/layout.html
ftable3 <- flextable(table3)
ftable3 <- autofit(ftable3)
ftable3 <- set_table_properties(ftable3, width = .5, layout = "autofit")
theme_booktabs(ftable3)

Table 4. Metrics in the Western Highlands 300-count IBI, with scoring formulas. Trend is the direction of metric response with increasing stress.

library(readxl)
library(knitr)
library(flextable)

# state directories
table.dir <- "tables"
table.file <- "Instruction_Tables.xlsx"
tab4.dir <- "Back_table2"

table4 <- read_excel(file.path(table.dir, table.file), sheet = tab4.dir
                     , na = c("NA", ""), trim_ws = TRUE, skip = 0
                     , col_names = TRUE)

# create flextable for output 
# https://davidgohel.github.io/flextable/articles/layout.html
ftable4 <- flextable(table4)
ftable4 <- autofit(ftable4)
ftable4 <- set_table_properties(ftable4, width = .5, layout = "autofit")
theme_booktabs(ftable4)

Table 5. Metrics in the Central Hills 100-count IBI, with scoring formulas. Trend is the direction of metric response with increasing stress.

library(readxl)
library(knitr)
library(flextable)

# state directories
table.dir <- "tables"
table.file <- "Instruction_Tables.xlsx"
tab5.dir <- "Back_table3"

table5 <- read_excel(file.path(table.dir, table.file), sheet = tab5.dir
                     , na = c("NA", ""), trim_ws = TRUE, skip = 0
                     , col_names = TRUE)

ftable5 <- flextable(table5)
ftable5 <- autofit(ftable5)
ftable5 <- set_table_properties(ftable5, width = .5, layout = "autofit")
theme_booktabs(ftable5)

Table 6. Metrics in the Central Hills 300-count IBI, with scoring formulas. Trend is the direction of metric response with increasing stress.

library(readxl)
library(knitr)
library(flextable)

# state directories
table.dir <- "tables"
table.file <- "Instruction_Tables.xlsx"
tab6.dir <- "Back_table4"

table6 <- read_excel(file.path(table.dir, table.file), sheet = tab6.dir
                     , na = c("NA", ""), trim_ws = TRUE, skip = 0
                     , col_names = TRUE)

ftable6 <- flextable(table6)
ftable6 <- autofit(ftable6)
ftable6 <- set_table_properties(ftable6, width = .5, layout = "autofit")
theme_booktabs(ftable6)

Table 7. Metrics in the Low Gradient IBI, with scoring formulas. Trend is the direction of metric response with increasing stress.

library(readxl)
library(knitr)
library(flextable)

# state directories
table.dir <- "tables"
table.file <- "Instruction_Tables.xlsx"
tab7.dir <- "Back_table5"

table7 <- read_excel(file.path(table.dir, table.file), sheet = tab7.dir
                     , na = c("NA", ""), trim_ws = TRUE, skip = 0
                     , col_names = TRUE)

ftable7 <- flextable(table7)
ftable7 <- autofit(ftable7)
ftable7 <- set_table_properties(ftable7, width = .5, layout = "autofit")
theme_booktabs(ftable7)

PLOTS

Map of Sites in Massachusetts

These are the sites that have been input into the R shiny application.

library(ggplot2)
##### Massachusetts Map of Sites
## Color palette
region.colors <- c("KickIBI_CH_300ct" = "#e41a1c", "KickIBI_WH_300ct"= "#377eb8",
                   "KickIBI_CH_100ct" = "#4daf4a", "KickIBI_WH_100ct"= "#984ea3"
                   , "LowGradientIBI" = "#ff7f00")
region.labels <- c("KickIBI_CH_300ct" = "Cent. Hills (300 ct)", "KickIBI_WH_300ct"= "West. High. (300ct)",
                   "KickIBI_CH_100ct" = "Cent. Hills (100 ct)", "KickIBI_WH_100ct"= "West. High. (100ct)"
                   , "LowGradientIBI" = "Low Gradient")

## Map of MA
m1 <- ggplot(data = subset(map_data("state"), region %in% c("massachusetts", "rhode island", "connecticut"))) + 
            geom_polygon(aes(x=long, y=lat, group=group), fill="light gray", color="black") +
            coord_fixed(1.3)+
          labs(x= "Longitude",
               y= "Latitude")+
          theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "right")

## Add Points
m1 <- m1 + geom_point(data=df_metsc
                      , aes(LONG, LAT, color=INDEX_REGION), size = 2)

## Change colors
m1 <- m1 + scale_colour_manual(name="Index Region"
                               , values=region.colors
                               , labels = region.labels)
## Print Plot
m1

Overall Index Scores by Index Region

# Packages
library(ggplot2)
library(dplyr)

# Plot
region.colors <- c("KickIBI_CH_300ct" = "#e41a1c", "KickIBI_WH_300ct"= "#377eb8",
                   "KickIBI_CH_100ct" = "#4daf4a", "KickIBI_WH_100ct"= "#984ea3"
                   , "LowGradientIBI" = "#ff7f00") # plot colors
region.labels <- c("KickIBI_CH_300ct" = "Cent. Hills (300 ct)", "KickIBI_WH_300ct"= "West. High. (300ct)",
                   "KickIBI_CH_100ct" = "Cent. Hills (100 ct)", "KickIBI_WH_100ct"= "West. High. (100ct)"
                   , "LowGradientIBI" = "Low Gradient")

# Western Highlands Scores 300ct
WH_300 <- df_metsc %>% 
  filter(INDEX_REGION == "KickIBI_WH_300ct") %>% 
  select_if(~sum(!is.na(.)) > 0)

if(dim(WH_300)[1]==0){
  cat("")
} else {
    p1 <- ggplot(WH_300, aes(x = Index, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(x= "Scores",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")
} # end if else statement

# Western Highlands Scores 100ct
WH_100 <- df_metsc %>% 
  filter(INDEX_REGION == "KickIBI_WH_100ct") %>% 
  select_if(~sum(!is.na(.)) > 0)

if(dim(WH_100)[1]==0){
  cat("")
} else {
    p2 <- ggplot(WH_100, aes(x = Index, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(x= "Scores",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")
} # end if else statement

# Central Hills Scores 300ct
CH_300 <- df_metsc %>% 
  filter(INDEX_REGION == "KickIBI_CH_300ct") %>% 
  select_if(~sum(!is.na(.)) > 0)

if(dim(CH_300)[1]==0){
  cat("")
} else {
    p3 <- ggplot(CH_300, aes(x = Index, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(x= "Scores",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")
} # end if else statement

# Central Hills Scores 100ct
CH_100 <- df_metsc %>% 
  filter(INDEX_REGION == "KickIBI_CH_100ct") %>% 
  select_if(~sum(!is.na(.)) > 0)

if(dim(CH_100)[1]==0){
  cat("")
} else {
    p4 <- ggplot(CH_100, aes(x = Index, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(x= "Scores",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")
} # end if else statement

# Low Gradient Scores
LowGrad_data <- df_metsc %>% 
  filter(INDEX_REGION == "LowGradientIBI") %>% 
  select_if(~sum(!is.na(.)) > 0)

if(dim(LowGrad_data)[1]==0){
  cat("")
} else {
    p5 <- ggplot(LowGrad_data, aes(x = Index, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(x= "Scores",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")
} # end if else statement

# Output
if(dim(WH_300)[1]==0){
  cat("")
} else {
  p1
} # end if statement

if(dim(WH_100)[1]==0){
  cat("")
} else {
  p2
} # end if statement

if(dim(CH_300)[1]==0){
  cat("")
} else {
  p3
} # end if statement

if(dim(CH_100)[1]==0){
  cat("")
} else {
  p4
} # end if statement

if(dim(LowGrad_data)[1]==0){
  cat("")
} else {
  p5
} # end if statement

Metric Values Distribution by Index Region.

region.colors <- c("KickIBI_CH_300ct" = "#e41a1c", "KickIBI_WH_300ct"= "#377eb8",
                   "KickIBI_CH_100ct" = "#4daf4a", "KickIBI_WH_100ct"= "#984ea3"
                   , "LowGradientIBI" = "#ff7f00") # plot colors
region.labels <- c("KickIBI_CH_300ct" = "Cent. Hills (300 ct)", "KickIBI_WH_300ct"= "West. High. (300ct)",
                   "KickIBI_CH_100ct" = "Cent. Hills (100 ct)", "KickIBI_WH_100ct"= "West. High. (100ct)"
                   , "LowGradientIBI" = "Low Gradient")
##### METRIC VALUES LOOP

#create table of metric names

col_metabbrev <- c("nt_total"
                 ,"pt_EPT"
                 ,"pi_EphemNoCaeBae"
                 ,"pi_ffg_filt"
                 ,"pt_ffg_pred"
                 ,"pt_tv_intol"
                 ,"pi_Pleco"
                 ,"pi_ffg_shred"
                 ,"pi_tv_intol"
                 ,"x_Becks"
                 ,"pi_OET"
                 ,"pt_NonIns"
                 ,"pt_POET"
                 ,"pt_tv_toler"
                 ,"pt_volt_semi")

col_metscore <- c("SC_nt_total"
                 ,"SC_pt_EPT"
                 ,"SC_pi_EphemNoCaeBae"
                 ,"SC_pi_ffg_filt"
                 ,"SC_pt_ffg_pred"
                 ,"SC_pt_tv_intol"
                 ,"SC_pi_Pleco"
                 ,"SC_pi_ffg_shred"
                 ,"SC_pi_tv_intol"
                 ,"SC_x_Becks"
                 ,"SC_pi_OET"
                 ,"SC_pt_NonIns"
                 ,"SC_pt_POET"
                 ,"SC_pt_tv_toler"
                 ,"SC_pt_volt_semi")

col_metnames <- c("Total Taxa Richness"
                  , "% Ephem., Plec., Trich. taxa"
                  , "% Ephem indiv (no Caenids or Baetids)"
                  , "% Filterer individuals"
                  , "% Predator taxa"
                  , "% Intolerant taxa"
                  , "% Plecoptera individuals"
                  , "% Shredder individuals"
                  , "% Intolerant individuals"
                  , "Becks Biotic Index"
                  , "% Odon., Ephem., Trich. indivs"
                  , "% non-insect taxa"
                  , "% Plec., Odon., Ephem., Trich. taxa"
                  , "% tolerant taxa"
                  , "% semivoltine taxa")

df_metnames <- as.data.frame(cbind(col_metabbrev,col_metscore, col_metnames))
df_metnames <- df_metnames %>% 
  rename("Met_Names" = col_metnames,
         "Met_Scores" = col_metscore,
         "Met_Abbrev" = col_metabbrev)

# western highlands metrics - 300ct
# for columns "ni_total" to ncol(df_metval)
WH_300_metval <- df_metval %>% 
  select(Region_Name, everything()) %>% 
  filter(INDEX_REGION == "KickIBI_WH_300ct")

#start column
if(dim(WH_300_metval)[1]==0){
  start_col <- 0
} else {
  start_col <- which(colnames(WH_300_metval)=="nt_total")
} # end start_col statement

#end column
if(dim(WH_300_metval)[1]==0){
  end_col <- 0
} else {
  end_col <- ncol(WH_300_metval)
} # end end_col statement

# CREATE PLOTS - only if data are present
if(start_col == 0 & end_col == 0){
  cat("")
} else {

  metval_cols <- colnames(WH_300_metval[,start_col:end_col])

for(i in metval_cols){
  subset.data <- WH_300_metval[,c("INDEX_REGION", i)]

  subset.data <- subset.data %>% 
    rename(Value = eval(i))

  subset.metnames <- df_metnames %>% 
    filter(Met_Abbrev == i)


  title <- paste0(subset.metnames$Met_Names)

  subset.data$INDEX_REGION <- as.factor(subset.data$INDEX_REGION)

  p6 <- ggplot(subset.data, aes(x = Value, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(title = title,
         x= "Metric Values",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")

  print(p6)

} # END of Loop
} # end if else statement

# western highlands metrics - 100ct
# for columns "ni_total" to ncol(df_metval)
WH_100_metval <- df_metval %>% 
  select(Region_Name, everything()) %>% 
  filter(INDEX_REGION == "KickIBI_WH_100ct")

#start column
if(dim(WH_100_metval)[1]==0){
  start_col <- 0
} else {
  start_col <- which(colnames(WH_100_metval)=="nt_total")
} # end start_col statement

#end column
if(dim(WH_100_metval)[1]==0){
  end_col <- 0
} else {
  end_col <- ncol(WH_100_metval)
} # end end_col statement

# CREATE PLOTS - only if data are present
if(start_col == 0 & end_col == 0){
  cat("")
} else {

  metval_cols <- colnames(WH_100_metval[,start_col:end_col])

for(i in metval_cols){
  subset.data <- WH_100_metval[,c("INDEX_REGION", i)]

  subset.data <- subset.data %>% 
    rename(Value = eval(i))

  subset.metnames <- df_metnames %>% 
    filter(Met_Abbrev == i)


  title <- paste0(subset.metnames$Met_Names)

  subset.data$INDEX_REGION <- as.factor(subset.data$INDEX_REGION)

  p7 <- ggplot(subset.data, aes(x = Value, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(title = title,
         x= "Metric Values",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")

  print(p7)

} # END of Loop
} # end if else statement

# central hills metrics - 300ct
# for columns "ni_total" to ncol(df_metval)
CH_300_metval <- df_metval %>% 
  select(Region_Name, everything()) %>% 
  filter(INDEX_REGION == "KickIBI_CH_300ct")

#start column
if(dim(CH_300_metval)[1]==0){
  start_col <- 0
} else {
  start_col <- which(colnames(CH_300_metval)=="nt_total")
} # end start_col statement

#end column
if(dim(CH_300_metval)[1]==0){
  end_col <- 0
} else {
  end_col <- ncol(CH_300_metval)
} # end end_col statement

# CREATE PLOTS - only if data are present
if(start_col == 0 & end_col == 0){
  cat("")
} else {

  metval_cols <- colnames(CH_300_metval[,start_col:end_col])

for(i in metval_cols){
  subset.data <- CH_300_metval[,c("INDEX_REGION", i)]

  subset.data <- subset.data %>% 
    rename(Value = eval(i))

  subset.metnames <- df_metnames %>% 
    filter(Met_Abbrev == i)


  title <- paste0(subset.metnames$Met_Names)

  subset.data$INDEX_REGION <- as.factor(subset.data$INDEX_REGION)

  p8 <- ggplot(subset.data, aes(x = Value, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(title = title,
         x= "Metric Values",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")

  print(p8)

} # END of Loop
} # end if else statement

# central hills metrics - 100ct
# for columns "ni_total" to ncol(df_metval)
CH_100_metval <- df_metval %>% 
  select(Region_Name, everything()) %>% 
  filter(INDEX_REGION == "KickIBI_CH_100ct")

#start column
if(dim(CH_100_metval)[1]==0){
  start_col <- 0
} else {
  start_col <- which(colnames(CH_100_metval)=="nt_total")
} # end start_col statement

#end column
if(dim(CH_100_metval)[1]==0){
  end_col <- 0
} else {
  end_col <- ncol(CH_100_metval)
} # end end_col statement

# CREATE PLOTS - only if data are present
if(start_col == 0 & end_col == 0){
  cat("")
} else {

  metval_cols <- colnames(CH_100_metval[,start_col:end_col])

for(i in metval_cols){
  subset.data <- CH_100_metval[,c("INDEX_REGION", i)]

  subset.data <- subset.data %>% 
    rename(Value = eval(i))

  subset.metnames <- df_metnames %>% 
    filter(Met_Abbrev == i)


  title <- paste0(subset.metnames$Met_Names)

  subset.data$INDEX_REGION <- as.factor(subset.data$INDEX_REGION)

  p9 <- ggplot(subset.data, aes(x = Value, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(title = title,
         x= "Metric Values",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")

  print(p9)

} # END of Loop
} # end if else statement


# low gradient metrics
# for columns "ni_total" to ncol(df_metval)
LowGrad_metval <- df_metval %>% 
  select(Region_Name, everything()) %>% 
  filter(INDEX_REGION == "LowGradientIBI")

#start column
if(dim(LowGrad_metval)[1]==0){
  start_col <- 0
} else {
  start_col <- which(colnames(LowGrad_metval)=="nt_total")
} # end start_col statement

#end column
if(dim(LowGrad_metval)[1]==0){
  end_col <- 0
} else {
  end_col <- ncol(LowGrad_metval)
} # end end_col statement

# CREATE PLOTS - only if data are present
if(start_col == 0 & end_col == 0){
  cat("")
} else {

  metval_cols <- colnames(LowGrad_metval[,start_col:end_col])

for(i in metval_cols){
  subset.data <- LowGrad_metval[,c("INDEX_REGION", i)]

  subset.data <- subset.data %>% 
    rename(Value = eval(i))

  subset.metnames <- df_metnames %>% 
    filter(Met_Abbrev == i)


  title <- paste0(subset.metnames$Met_Names)

  subset.data$INDEX_REGION <- as.factor(subset.data$INDEX_REGION)

  p10 <- ggplot(subset.data, aes(x = Value, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(title = title,
         x= "Metric Values",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")

  print(p10)

} # END of Loop
} # end if else statement

Metric Scores Distribution by Index Region.

Only metrics applied within each region are plotted.

region.colors <- c("KickIBI_CH_300ct" = "#e41a1c", "KickIBI_WH_300ct"= "#377eb8",
                   "KickIBI_CH_100ct" = "#4daf4a", "KickIBI_WH_100ct"= "#984ea3"
                   , "LowGradientIBI" = "#ff7f00") # plot colors
region.labels <- c("KickIBI_CH_300ct" = "Cent. Hills (300 ct)", "KickIBI_WH_300ct"= "West. High. (300ct)",
                   "KickIBI_CH_100ct" = "Cent. Hills (100 ct)", "KickIBI_WH_100ct"= "West. High. (100ct)"
                   , "LowGradientIBI" = "Low Gradient")


##### METRIC SCORES LOOP
# WESTERN HIGHLANDS - 300ct
# start column
if(dim(WH_300)[1]==0){
  start_col <- 0
} else {
  start_col <- which(colnames(WH_300)=="SC_nt_total")
} # end start_col statement

# end column
if(dim(WH_300)[1]==0){
  end_col <- 0
} else {
  end_col <- which(colnames(WH_300)=="SC_x_Becks")
} # end end_col statement

# CREATE PLOTS - only if data are present
if(start_col == 0 & end_col == 0){
  cat("")
} else {
  metsc_cols <- colnames(WH_300[,start_col:end_col])


for(i in metsc_cols){
  subset.data <- WH_300[,c("INDEX_REGION", i)]

  subset.data <- subset.data %>% 
    rename(Value = eval(i))

  subset.metnames <- df_metnames %>% 
    filter(Met_Scores == i)


  title <- paste0("Met Score - ", subset.metnames$Met_Names)

  subset.data$INDEX_REGION <- as.factor(subset.data$INDEX_REGION)

  subset.data <- subset.data %>% filter(!is.na(Value))

  p11 <- ggplot(subset.data, aes(x = Value, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(title = title,
         x= "Metric Scores",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")

  print(p11)

} # END of Loop
} # end if else statement

# WESTERN HIGHLANDS - 100ct
# start column
if(dim(WH_100)[1]==0){
  start_col <- 0
} else {
  start_col <- which(colnames(WH_100)=="SC_nt_total")
} # end start_col statement

# end column
if(dim(WH_100)[1]==0){
  end_col <- 0
} else {
  end_col <- which(colnames(WH_100)=="SC_x_Becks")
} # end end_col statement

# CREATE PLOTS - only if data are present
if(start_col == 0 & end_col == 0){
  cat("")
} else {
  metsc_cols <- colnames(WH_100[,start_col:end_col])


for(i in metsc_cols){
  subset.data <- WH_100[,c("INDEX_REGION", i)]

  subset.data <- subset.data %>% 
    rename(Value = eval(i))

  subset.metnames <- df_metnames %>% 
    filter(Met_Scores == i)


  title <- paste0("Met Score - ", subset.metnames$Met_Names)

  subset.data$INDEX_REGION <- as.factor(subset.data$INDEX_REGION)

  subset.data <- subset.data %>% filter(!is.na(Value))

  p12 <- ggplot(subset.data, aes(x = Value, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(title = title,
         x= "Metric Scores",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")

  print(p12)

} # END of Loop
} # end if else statement

# CENTRAL HILLS - 300ct
# start column
if(dim(CH_300)[1]==0){
  start_col <- 0
} else {
  start_col <- which(colnames(CH_300)=="SC_nt_total")
} # end start_col statement

# end column
if(dim(CH_300)[1]==0){
  end_col <- 0
} else {
  end_col <- which(colnames(CH_300)=="SC_pt_tv_intol")
} # end end_col statement

# CREATE PLOTS - only if data are present
if(start_col == 0 & end_col == 0){
  cat("")
} else {
  metsc_cols <- colnames(CH_300[,start_col:end_col])


for(i in metsc_cols){
  subset.data <- CH_300[,c("INDEX_REGION", i)]

  subset.data <- subset.data %>% 
    rename(Value = eval(i))

  subset.metnames <- df_metnames %>% 
    filter(Met_Scores == i)


  title <- paste0("Met Score - ", subset.metnames$Met_Names)

  subset.data$INDEX_REGION <- as.factor(subset.data$INDEX_REGION)

  subset.data <- subset.data %>% filter(!is.na(Value))

  p13 <- ggplot(subset.data, aes(x = Value, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(title = title,
         x= "Metric Scores",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")

  print(p13)

} # END of Loop
} # end if else statement

# CENTRAL HILLS - 100ct
# start column
if(dim(CH_100)[1]==0){
  start_col <- 0
} else {
  start_col <- which(colnames(CH_100)=="SC_nt_total")
} # end start_col statement

# end column
if(dim(CH_100)[1]==0){
  end_col <- 0
} else {
  end_col <- which(colnames(CH_100)=="SC_pt_tv_intol")
} # end end_col statement

# CREATE PLOTS - only if data are present
if(start_col == 0 & end_col == 0){
  cat("")
} else {
  metsc_cols <- colnames(CH_100[,start_col:end_col])


for(i in metsc_cols){
  subset.data <- CH_100[,c("INDEX_REGION", i)]

  subset.data <- subset.data %>% 
    rename(Value = eval(i))

  subset.metnames <- df_metnames %>% 
    filter(Met_Scores == i)


  title <- paste0("Met Score - ", subset.metnames$Met_Names)

  subset.data$INDEX_REGION <- as.factor(subset.data$INDEX_REGION)

  subset.data <- subset.data %>% filter(!is.na(Value))

  p14 <- ggplot(subset.data, aes(x = Value, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(title = title,
         x= "Metric Scores",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")

  print(p14)

} # END of Loop
} # end if else statement

# LOW GRADIENT scores
# start column
if(dim(LowGrad_data)[1]==0){
  start_col <- 0
} else {
  start_col <- which(colnames(LowGrad_data)=="SC_pi_OET")
} # end start_col statement

# end column
if(dim(LowGrad_data)[1]==0){
  end_col <- 0
} else {
  end_col <- which(colnames(LowGrad_data)=="SC_pt_volt_semi")
} # end end_col statement

# CREATE PLOTS - only if data are present
if(start_col == 0 & end_col == 0){
  cat("")
} else {
  metsc_cols <- colnames(LowGrad_data[,start_col:end_col])


for(i in metsc_cols){
  subset.data <- LowGrad_data[,c("INDEX_REGION", i)]

  subset.data <- subset.data %>% 
    rename(Value = eval(i))

  subset.metnames <- df_metnames %>% 
    filter(Met_Scores == i)


  title <- paste0("Met Score - ", subset.metnames$Met_Names)

  subset.data$INDEX_REGION <- as.factor(subset.data$INDEX_REGION)

  subset.data <- subset.data %>% filter(!is.na(Value))

  p14 <- ggplot(subset.data, aes(x = Value, fill = INDEX_REGION))+ 
    geom_histogram(alpha = 1, position = "dodge", bins = 20, color = "black")+
    labs(title = title,
         x= "Metric Scores",
         y= "Frequency",
         fill= "Index Region")+
    scale_fill_manual(values=region.colors, labels = region.labels)+
    theme(text = element_text(size = 12),
          axis.text = element_text(color = "black", size = 12),
          axis.text.x = element_text(angle = 0, hjust = 0.5),
          panel.background = element_rect(fill = "white"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          axis.line = element_line(color = "black"),
          legend.position = "top")

  print(p14)

} # END of Loop
} # end if else statement

Blocktt/MassIBItools documentation built on May 6, 2021, 11:46 a.m.