getBioStressorResponses: Biological Stressor Responses
In leppott/CASTfxn: Functions for Causal Assessment Screening Tool (CAST)

Description Usage Arguments Details Value Examples

View source: R/getBioStressorResponses.R

Get Biological (Algae or BMI) stressor responses.

getBioStressorResponses(TargetSiteID, stressors, BioResp,
  list.MatchBioData, LogTransf, ref.sites, biocomm = "bmi",
  dir_results = file.path(getwd(), "Results"),
  dir_sub = "StressorResponse", boo_pred_warn = TRUE)

`TargetSiteID`	Site ID
`stressors`	stressors
`BioResp`	Biological response variables. For example, BMI metrics or Algae metrics.
`list.MatchBioData`	list of matched biological (BMI or algae) and stressor data.
`LogTransf`	Value for if stressor variables should be log10 transformed; 1=TRUE, 0=FALSE.
`ref.sites`	Reference sites.
`biocomm`	Biological community; algae or BMI. Default = "BMI".
`dir_results`	Directory to save plots. Default = working directory and Results.
`dir_sub`	Subdirectory for outputs from this function. Default = "StressorResponse"
`boo_pred_warn`	Should warnings for prediction be suppressed. Default = TRUE.

Biological (Algae or BMI) stressor regressions.

A jpg in SiteID subfoler of the "Results" folder of working directory. And two tab-delimited text files; stressor correlations and scores.

## Not run: 
# Example 1, BMI
TargetSiteID <- "SRCKN001.61"
dir_results <- file.path(getwd(), "Results")
biocomm <- "bmi"

# datasets getSiteInfo
# data, example included with package
data.Stations.Info <- data_Sites       # need for getSiteInfo and getChemDataSubsets
data.SampSummary   <- data_SampSummary
data.303d.ComID    <- data_303d
data.bmi.metrics   <- data_BMIMetrics
data.algae.metrics <- data_AlgMetrics
data.mod           <- data_ReachMod

# Cluster based on elevation category  # need for getSiteInfo and getChemDataSubsets
elev_cat <- toupper(data.Stations.Info[data.Stations.Info[,"StationID_Master"]==TargetSiteID
                    , "ElevCategory"])
if(elev_cat=="HI"){
   data.cluster <- data_Cluster_Hi
} else if(elev_cat=="LO") {
   data.cluster <- data_Cluster_Lo
}

# Map data
# San Diego
#flowline <- rgdal::readOGR(dsn = "data_gis/NHDv2_Flowline_Ecoreg85", layer = "NHDv2_eco85_Project")
#outline <- rgdal::readOGR(dsn = "data_gis/Eco85", layer = "Ecoregion85")
# AZ
map_flowline  <- data_GIS_Flow_HI
map_flowline2 <- data_GIS_Flow_LO
if(elev_cat=="HI"){
   map_flowline <- data_GIS_Flow_HI
} else if(elev_cat=="LO") {
   map_flowline <- data_GIS_Flow_LO
}
map_outline   <- data_GIS_AZ_Outline
# Project site data to USGS Albers Equal Area
usgs.aea <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23
              +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83
              +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
# projection for outline
my.aea <- "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 
           +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
map_proj <- my.aea
#
dir_sub <- "SiteInfo"

# Run getSiteInfo
list.SiteSummary <- getSiteInfo(TargetSiteID, dir_results, data.Stations.Info
                                , data.SampSummary, data.303d.ComID
                                , data.bmi.metrics, data.algae.metrics
                                , data.cluster, data.mod
                                , map_proj, map_outline, map_flowline
                                , dir_sub=dir_sub)

# Data getChemDataSubsets
# data import, example 
# data.chem.raw <- read.delim(paste(myDir.Data,"data.chem.raw.tab",sep=""),na.strings = c(""," "))
# data.chem.info <- read.delim(paste(myDir.Data,"data.chem.info.tab",sep=""))
site.COMID <- list.SiteSummary$COMID
site.Clusters <- list.SiteSummary$ClustIDs
# data, example included with package
data.chem.raw <- data_Chem
data.chem.info <- data_ChemInfo

# Run getChemDataSubsets
list.data <- getChemDataSubsets(TargetSiteID, comid=site.COMID, cluster=site.Clusters
                                , data.cluster=data.cluster, data.Stations.Info=data.Stations.Info
                                , data.chem.raw=data.chem.raw, data.chem.info=data.chem.info)

# Data getStressorList
chem.info <- list.data$chem.info
cluster.chem <- list.data$cluster.chem
cluster.samps <- list.data$cluster.samps
ref.sites <- list.data$ref.sites
site.chem <- list.data$site.chem
dir_sub <- "CandidateCauses"

# set cutoff for possible stressor identification
probsLow <- 0.10
probsHigh <- 0.90 

# Run getStressorList
list.stressors <- getStressorList(TargetSiteID, site.Clusters, chem.info, cluster.chem
                                 , cluster.samps, ref.sites, site.chem
                                 , probsHigh, probsLow, biocomm, dir_results
                                 , dir_sub)
                                 
# Data getBioMatches, BMI
## remove "none"
stressors <- list.stressors$stressors[list.stressors$stressors != "none"]
stressors_logtransf <- list.stressors$stressors_LogTransf[list.stressors$stressors != "none"]
LogTransf <- stressors_logtransf
data.bio.metrics <- data_BMIMetrics

# Run getBioMatches
list.MatchBioData <- getBioMatches(stressors, list.data, list.SiteSummary, data.SampSummary
                                   , data.chem.raw, data.bio.metrics, biocomm)
  
# Data getBioStressorResponses, BMI 
BioResp <- c("IBI", "TotalTaxSPL_Sc", "DipTaxSPL_Sc"
             , "IntolTaxSPL_Sc", "HBISPL_Sc", "PlecoPct_Sc", "ScrapPctSPL_Sc"
             , "TrichTax_Sc", "EphemTax_Sc", "EphemPct_Sc", "Dom01PctSPL_Sc") 
dir_sub <- "StressorResponse"
             
# Run getBioStressorResponses, BMI               
getBioStressorResponses(TargetSiteID, stressors, BioResp, list.MatchBioData
                        , LogTransf, ref.sites, biocomm, dir_results, dir_sub)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Example 2, Algae
TargetSiteID <- "LCBEN002.57"
dir_results <- file.path(getwd(), "Results")
biocomm <- "algae"

# datasets getSiteInfo
# data, example included with package
data.Stations.Info <- data_Sites       # need for getSiteInfo and getChemDataSubsets
data.SampSummary   <- data_SampSummary
data.303d.ComID    <- data_303d
data.bmi.metrics   <- data_BMIMetrics
data.algae.metrics <- data_AlgMetrics
data.mod           <- data_ReachMod

#' # Cluster based on elevation category  # need for getSiteInfo and getChemDataSubsets
elev_cat <- toupper(data.Stations.Info[data.Stations.Info[,"StationID_Master"]==TargetSiteID
                    , "ElevCategory"])
if(elev_cat=="HI"){
   data.cluster <- data_Cluster_Hi
} else if(elev_cat=="LO") {
   data.cluster <- data_Cluster_Lo
}

# Map data
# San Diego
#flowline <- rgdal::readOGR(dsn = "data_gis/NHDv2_Flowline_Ecoreg85", layer = "NHDv2_eco85_Project")
#outline <- rgdal::readOGR(dsn = "data_gis/Eco85", layer = "Ecoregion85")
# AZ
map_flowline  <- data_GIS_Flow_HI
map_flowline2 <- data_GIS_Flow_LO
if(elev_cat=="HI"){
   map_flowline <- data_GIS_Flow_HI
} else if(elev_cat=="LO") {
   map_flowline <- data_GIS_Flow_LO
}
map_outline   <- data_GIS_AZ_Outline
# Project site data to USGS Albers Equal Area
usgs.aea <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23
              +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83
              +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
# projection for outline
my.aea <- "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 
           +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
map_proj <- my.aea
#
dir_sub <- "SiteInfo"

# Run getSiteInfo
list.SiteSummary <- getSiteInfo(TargetSiteID, dir_results, data.Stations.Info
                                , data.SampSummary, data.303d.ComID
                                , data.bmi.metrics, data.algae.metrics
                                , data.cluster, data.mod
                                , map_proj, map_outline, map_flowline
                                , dir_sub=dir_sub)

# Data getChemDataSubsets
# data, example included with package
data.chem.raw <- data_Chem
data.chem.info <- data_ChemInfo
site.COMID <- list.SiteSummary$COMID
site.Clusters <- list.SiteSummary$ClustIDs
dir_sub <- "CandidateCauses"

# Run getChemDataSubsets
list.data <- getChemDataSubsets(TargetSiteID, comid=site.COMID, cluster=site.Clusters
                                , data.cluster=data.cluster, data.Stations.Info=data.Stations.Info
                                , data.chem.raw=data.chem.raw, data.chem.info=data.chem.info)

# Data getStressorList
chem.info <- list.data$chem.info
cluster.chem <- list.data$cluster.chem
cluster.samps <- list.data$cluster.samps
ref.sites <- list.data$ref.sites
site.chem <- list.data$site.chem

# set cutoff for possible stressor identification
probsLow <- 0.10
probsHigh <- 0.90 
biocomm <- "algae"

# Run getStressorList
list.stressors <- getStressorList(TargetSiteID, site.Clusters, chem.info, cluster.chem
                                 , cluster.samps, ref.sites, site.chem
                                 , probsHigh, probsLow, biocomm, dir_results)
                                 
# Data getBioMatches, Algae
## remove "none"
stressors <- list.stressors$stressors[list.stressors$stressors != "none"]
stressors_logtransf <- list.stressors$stressors_LogTransf[list.stressors$stressors != "none"]
LogTransf <- stressors_logtransf
data.bio.metrics <- data_AlgMetrics

# Run getBioMatches, Algae
list.MatchBioData <- getBioMatches(stressors, list.data, list.SiteSummary, data.SampSummary
                                  , data.chem.raw, data.bio.metrics, biocomm)

# Data getBioStressorResponses, Algae
BioResp <- colnames(data.bio.metrics[6:52])
dir_sub <- "StressorResponse"

# Run getBioStressorResponses, Algae
getBioStressorResponses(TargetSiteID, stressors, BioResp, list.MatchBioData
                       , LogTransf, ref.sites, biocomm, dir_results, dir_sub)

## End(Not run)