License
Description
Setup
Get and prepare monitoring data
- Check for outliers
Get Spatial data and merge
Perform MCMC modeling
- Functions
- mcmc_calc function
Run model
Summary

License

Copyright (c) 2020 Geosyntec Consultants, Inc. Mozilla Public License Version 2.0

This software is provided “as is”, without warranty of any kind, express or implied, including but not limited to the warranties of merchantability, fitness for a particular purpose and noninfringement. In no event shall the authors or copyright holders be liable for any claim, damages or other liability, whether in an action of contract, tort or otherwise, arising from, out of or in connection with the software or the use or other dealings in the software.

Description

This notebook is the second step in the watershed regression task. It uses the relationships from Part 1 (linear models) to run generalized linear mixed models with censored dependent variables. This approach uses Markov chain Monte Carlo techniques and relies heavily on the MCMCglmm package (Hadfield 2009).

Here are the steps in the code below:

Code set up
Get and prepare monitoring data
Merge spatial data with monitoring data
Remove mulitcolinear predictors
Peform model selection
Select the best model (these results are then used in the Bayesian model in step 2)

Setup

clear workspace and load libraries. (code not shown in html file )

# library('renv')
library(caret)
library(readr)
library(dplyr)

setwd("~/repos/stormwaterheatmap/R-scripts/WatershedRegression")

Get and prepare monitoring data

The stormwater outfall data is available from the Department of Ecology at https://data.wa.gov/Natural-Resources-Environment/Municipal-Stormwater-Permit-Outfall-Data/d958-q2ci.

A .csv file is saved in WatershedRegression/data/S8_data.csv

all.S8.data <- read_csv("data/S8_data.csv", col_types = cols(field_collection_end_date = col_date(format = "%m/%d/%Y"), 
    field_collection_start_date = col_date(format = "%m/%d/%Y")))

# filter out rejected data
all.S8.data <- (filter(all.S8.data, !result_data_qualifier %in% "REJ"))

# filter out replicates
all.S8.data <- (filter(all.S8.data, !sample_replicate_flag %in% "Y"))

# change nondetect warnings to detects
warnings <- all.S8.data$nondetect_flag == "WARNING"
all.S8.data$nondetect_flag[warnings] <- FALSE

# Change NA to detect
all.S8.data$nondetect_flag[is.na(all.S8.data$nondetect_flag)] <- FALSE

# Change season to factor
all.S8.data$season <- as.factor(all.S8.data$season)

The chunk below makes a list of parameters we are interested in.

#Select Parameters
params <- c('Zinc - Water - Total',
 'Copper - Water - Total',
 'Nitrite-Nitrate - Water - Dissolved',
 'Lead - Water - Total',
 'Total Phosphorus - Water - Total',
 'Total Suspended Solids - Water - Total',
 'Total Phthalate - Water - Total',
'Total PAH - Water - Total',
#'Chrysene - Water - Total',
'CPAH - Water - Total',
'HPAH - Water - Total' 
#'Total Kjeldahl Nitrogen - Water - Total',
#'Total PCB - Water - Total'
)

#save a list of all the parameters in case we want to use mor. 
params.all <- data.frame(unique(all.S8.data$parameter))
s8data <- all.S8.data 
#

Clean up extracted data and rename columns:

s8data <- all.S8.data %>% 
dplyr::select(study_name, location_id, parameter, type, season, new_result_value, 
    nondetect_flag, study_id, access_id, field_collection_end_date, field_collection_start_date, 
    type)


# rename some columns
colnames(s8data)[colnames(s8data) == "location_id"] <- "Location"
colnames(s8data)[colnames(s8data) == "new_result_value"] <- "concentration"
s8data$nondetect_flag <- as.logical(s8data$nondetect_flag)
s8data$concentration <- as.numeric(s8data$concentration)

Check for outliers

Set up a jitter plot of all the data to look for outliers:

scatter_cocs <- function(df.coc,title) {
 p <- ggplot(df.coc, aes(1, concentration)) + geom_jitter() + labs(
  title = title,
  subtitle = "Data collected 2009-2013",
  caption =
    " Data source: Ecology, 2015",
  x = "Observations"
  #y = y.title
)  
p + facet_wrap( ~ parameter, ncol=3, scales = 'free') 
}

scatter_cocs(s8data[which(s8data$parameter %in% params),],'All Observations')

outliers are apprent for TSS, TP, and no2-no3. Remove these.

# remove and replot

outlierParams <- c("Total Suspended Solids - Water - Total", "Total Phosphorus - Water - Total", 
    "Nitrite-Nitrate - Water - Dissolved")

# This removes the highest values
outlierVals <- top_n(group_by(s8data[which(s8data$parameter %in% outlierParams), 
    ], parameter), 1, concentration)$concentration

s8data <- s8data %>% group_by(parameter) %>% slice(which(!(parameter %in% outlierParams & 
    concentration %in% outlierVals)))

scatter_cocs(s8data[which(s8data$parameter %in% params), ], "All Observations - Outliers Removed")

Looks better. Move on to the next Chunk.

Get Spatial data and merge

# Spatial predcitors have been extracted and saved as a csv file.

spatial_data <- read_csv("C:/Users/cnilsen/Google Drive/repos/spatialPredictors.csv")
# RES and COM are compositional data. Change to a ratio
spatial_data$LU_ratio = spatial_data$COM/spatial_data$RES
spatial_data <- dplyr::select(spatial_data, -c(RES, COM, .geo))
# merge spatial predictors with monitoring data
s8data.wPredictors <<- merge(s8data, spatial_data) %>% dplyr::select(-c(depSplusN))

# spatial_data<-read_csv('data/spatialPredictors42.csv', col_types = cols(X1 =
# col_skip())) #merge spatial predictors with monitoring data s8data.wPredictors
# <- merge(s8data, spatial_data)

Perform MCMC modeling

Functions

Some helper functions to help

Function to add a survival object to the S8 dataframe

add_surv <- function(df) {
    df$cenMin <- ifelse(df$nondetect_flag, -Inf, (df$concentration))
    df$cenMax <- (df$concentration)
    df$cenMin_log <- ifelse(df$nondetect_flag, -Inf, log(df$concentration))
    df$cenMax_log <- log(df$concentration)
    
    return(df)
}

Function to return a chart of predictions from the model

scatter_predict <- function(model_df,predictions) {
# model_to_predict <- CuModel 
# coc = params[2]
# df <- (subset(s8data.wPredictors, parameter == coc)) %>%
#   add_surv()
# predictions <- predict(model_to_predict, newdata=df, 
#          type="response", interval="none", level=0.95, it=NULL, 
#          posterior="all", verbose=FALSE, approx="numerical")
# 
 obs <- log(model_df$concentration)
 
 ggstatsplot::ggscatterstats(
  data =tibble(p = predictions, obs = obs,L = model_df$Location),
  x = p,
  y = obs,
  type = "bf",
  point.width.jitter = 0.02,
  #point.height.jitter = 0.1,
  marginal = FALSE,
  xlab = "Predicted log(µg/L)",
  ylab = "Observed  log(µg/L)",
  title = coc,
  results.subtitle = FALSE,
  subtitle = "Predictions vs. Observations",
  smooth.line.args = list(size = 1, color = "blue"),
  messages = FALSE
)
}

add tidy and glance functions to the mcmcglmm objects since they don’t play nicely with tidyverse objects.

# add custom functions to extract estimates (tidy) and goodness-of-fit (glance)
# information
tidy.MCMCglmm <- function(object, ...) {
    s <- summary(object, ...)
    ret <- tibble::tibble(term = row.names(s$solutions), estimate = s$solutions[, 
        1], conf.low = s$solutions[, 2], conf.high = s$solutions[, 3])
    ret
}
glance.MCMCglmm <- function(object, ...) {
    ret <- tibble::tibble(dic = object$DIC, n = nrow(object$X))
    ret
}

# estimate a simple model model <- MCMCglmm(PO ~ 1 + plate, random = ~ FSfamily,
# data = PlodiaPO, verbose=FALSE, pr=TRUE)

mcmc_calc function

This is the main funciton to run the mcmc model. It does the following:
1. subsets a dataframe to include only the parameter we want to predict
2. adds a survival object to handle censored data
3. sets up a simiple prior structure
4. performs mcmc modeling on either the log-transformed or non-log transformed responses.
5. returns the results

mcmc_calc <- function(coc.local, fixed_list, lhs) {
    
    df <- (subset(s8data.wPredictors, parameter == coc.local))
    
    data <- df %>% add_surv()
    # make the prior_structures
    prior.1 <- list(R = list(V = 1, fix = 1), G = list(G1 = list(V = 1, nu = 0.002)))
    prior.2 <- list(R = list(V = 2, fix = 1), G = list(G1 = list(V = 1, nu = 0)))
    
    
    if (lhs == "log") {
        mcmc_formula <- as.formula(paste("cbind(cenMin_log, cenMax_log) ~ ", paste0(fixed_list, 
            collapse = "+")))
    } else {
        mcmc_formula <- as.formula(paste("cbind(cenMin, cenMax) ~ ", paste0(fixed_list, 
            collapse = "+")))
    }
    
    mcmc_results <- MCMCglmm(mcmc_formula, random = ~Location, data = data, family = "cengaussian", 
        verbose = FALSE, prior = prior.1, singular.ok = TRUE, nitt = 60000, thin = 13, 
        burnin = 10000)
    return((mcmc_results))
}

function that returns bayesian plots

# Do some predictive checks
library(bayesplot)
color_scheme_set("blue")
bayes_plots <- function(fit, coc, df) {
    # fit <- TSSModel
    
    # coc = 'Total Suspended Solids - Water - Total'
    
    # df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
    yrep_c <- predict(fit, newdata = df, type = "response", interval = "confidence", 
        level = 0.9, it = NULL, posterior = "all", verbose = FALSE, approx = "numerical")
    yrep_p <- predict(fit, newdata = df, type = "response", interval = "prediction", 
        level = 0.9, it = NULL, posterior = "all", verbose = FALSE, approx = "numerical")
    
    # show uncertainty intervals under esimated posterior density curves
    plot.1 <- mcmc_areas(fit$Sol, prob = 0.8, pprob_outer = 0.95, point_est = "mean") + 
        ggplot2::labs(title = coc, subtitle = "Posterior distributions with medians and 80% intervals")
    
    # generate scatter plot of predictions
    
    colnames(yrep_p) <- c("fit.p", "lwr.p", "upr.p")
    scatterdata <- cbind(df, yrep_c, yrep_p)
    
    # generate scatter plot of predictions
    plot.2 <- ggplot(scatterdata) + geom_ribbon(aes(ymin = lwr.p, ymax = upr.p, x = fit), 
        fill = "grey", alpha = 0.5) + geom_ribbon(aes(ymin = lwr, ymax = upr, x = fit), 
        fill = "grey", alpha = 0.8) + geom_line(aes(x = fit, y = fit), color = "blue", 
        linetype = 5) + geom_point(aes(x = fit, y = log(concentration)), alpha = 0.5) + 
        ggplot2::labs(x = "yrep", y = "fit", title = coc, subtitle = "Scatter plot of observed data vs simulated", 
            caption = "dark shade: confidence intervals \n light shade: prediction intervals")
    
    
    # simulate with 100 draws
    ysim <- (simulate(fit, nsim = 100))
    
    
    # overlay of predictions
    plot.3 <- ppc_dens_overlay(log(df$concentration), t(ysim)) + ggplot2::labs(x = "log concentration, µg/L ", 
        title = coc, subtitle = "Observed (y) vs. simulated draws (yrep)")
    
    
    return(list(plot.1, plot.2, plot.3))
}

# function for other diagnostic plots
diagnostic_plots <- function(chains, coc) {
    plotTrace(chains, axes = TRUE, same.limits = TRUE)
    plotDens(chains, main = paste("Posterior Distributions \n", coc), probs = c(0.05, 
        0.95), same.limits = FALSE)
}

Run model

For each chunk below, we run the model and output diagnostic and prediction plots.

Zinc

# 
coc <- params[1]
ZnModel <- mcmc_calc(coc, c("impervious"), "log")
mod <- ZnModel
summary(mod)


 Iterations = 10001:59999
 Thinning interval  = 13
 Sample size  = 3847 

 DIC: 928.197 

 G-structure:  ~Location

         post.mean l-95% CI u-95% CI eff.samp
Location    0.5703   0.1667    1.156     3563

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: cbind(cenMin_log, cenMax_log) ~ impervious 

            post.mean l-95% CI u-95% CI eff.samp  pMCMC    
(Intercept)    3.6401   3.2098   4.0554     3847 <3e-04 ***
impervious     0.7137   0.4177   1.0310     3847 <3e-04 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
bayes_plots(mod, coc, df)

[[1]]


[[2]]


[[3]]

Copper

coc = "Copper - Water - Total"
CuModel <- mcmc_calc(coc, c("rev_logTraffic", "impervious"), "log")
mod <- CuModel
summary(mod$Sol)


Iterations = 10001:59999
Thinning interval = 13 
Number of chains = 1 
Sample size per chain = 3847 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                 Mean     SD Naive SE Time-series SE
(Intercept)    5.4050 0.8871 0.014303       0.013878
rev_logTraffic 2.2066 0.5519 0.008899       0.008660
impervious     0.3038 0.1274 0.002054       0.002057

2. Quantiles for each variable:

                  2.5%    25%    50%    75%  97.5%
(Intercept)    3.64854 4.8462 5.4047 5.9571 7.2030
rev_logTraffic 1.11494 1.8633 2.1970 2.5454 3.3350
impervious     0.05399 0.2237 0.3023 0.3798 0.5694

summary(mod)


 Iterations = 10001:59999
 Thinning interval  = 13
 Sample size  = 3847 

 DIC: 946.9554 

 G-structure:  ~Location

         post.mean l-95% CI u-95% CI eff.samp
Location    0.3326  0.09296   0.7163     3847

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: cbind(cenMin_log, cenMax_log) ~ rev_logTraffic + impervious 

               post.mean l-95% CI u-95% CI eff.samp   pMCMC    
(Intercept)      5.40501  3.64710  7.19880     4086 < 3e-04 ***
rev_logTraffic   2.20657  1.02546  3.23138     4062 0.00156 ** 
impervious       0.30377  0.04928  0.55974     3837 0.02132 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
bayes_plots(mod, coc, df)

[[1]]


[[2]]


[[3]]

Nitrite-Nitrate

coc = "Nitrite-Nitrate - Water - Dissolved"
NN_model <- mcmc_calc(coc, c("LU_ratio"), "log")
mod <- NN_model
summary(mod$Sol)


Iterations = 10001:59999
Thinning interval = 13 
Number of chains = 1 
Sample size per chain = 3847 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

              Mean      SD Naive SE Time-series SE
(Intercept) 5.7152 0.19984 0.003222       0.003148
LU_ratio    0.1577 0.07815 0.001260       0.001260

2. Quantiles for each variable:

                2.5%    25%    50%    75%  97.5%
(Intercept) 5.310966 5.5877 5.7167 5.8458 6.1042
LU_ratio    0.002186 0.1084 0.1557 0.2082 0.3119

summary(mod)


 Iterations = 10001:59999
 Thinning interval  = 13
 Sample size  = 3847 

 DIC: 984.2904 

 G-structure:  ~Location

         post.mean l-95% CI u-95% CI eff.samp
Location    0.5093   0.1568     1.04     3847

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: cbind(cenMin_log, cenMax_log) ~ LU_ratio 

            post.mean  l-95% CI  u-95% CI eff.samp  pMCMC    
(Intercept)  5.715214  5.300888  6.087447     4029 <3e-04 ***
LU_ratio     0.157691 -0.005646  0.301965     3847 0.0489 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
bayes_plots(mod, coc, df)

[[1]]


[[2]]


[[3]]

Lead

coc <- params[4]
PbModel <- mcmc_calc(params[4], c("impervious"), "log")

mod <- PbModel
summary(mod$Sol)


Iterations = 10001:59999
Thinning interval = 13 
Number of chains = 1 
Sample size per chain = 3847 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

              Mean     SD Naive SE Time-series SE
(Intercept) 1.1980 0.3201 0.005161       0.005161
impervious  0.6475 0.2334 0.003764       0.003856

2. Quantiles for each variable:

              2.5%    25%    50%    75% 97.5%
(Intercept) 0.5589 0.9954 1.2052 1.3998 1.833
impervious  0.1711 0.4995 0.6485 0.7953 1.100

summary(mod)


 Iterations = 10001:59999
 Thinning interval  = 13
 Sample size  = 3847 

 DIC: 1039.448 

 G-structure:  ~Location

         post.mean l-95% CI u-95% CI eff.samp
Location     1.441   0.4831    2.843     4425

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: cbind(cenMin_log, cenMax_log) ~ impervious 

            post.mean l-95% CI u-95% CI eff.samp   pMCMC   
(Intercept)    1.1980   0.5906   1.8490     3847 0.00364 **
impervious     0.6475   0.1905   1.1112     3665 0.00884 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
bayes_plots(mod, coc, df)

[[1]]


[[2]]


[[3]]

## Cadium

coc <- "Cadmium - Water - Total"
CdModel <- mcmc_calc("Cadmium - Water - Total", c("impervious", "logPopulation"), 
    "nonlog")
mod <- CdModel
summary(mod$Sol)


Iterations = 10001:59999
Thinning interval = 13 
Number of chains = 1 
Sample size per chain = 3847 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                 Mean      SD  Naive SE Time-series SE
(Intercept)   -0.1434 0.08136 0.0013117      0.0013117
impervious     0.2948 0.05925 0.0009553      0.0010129
logPopulation -0.2527 0.05613 0.0009050      0.0008809

2. Quantiles for each variable:

                 2.5%     25%     50%      75%    97.5%
(Intercept)   -0.3076 -0.1949 -0.1442 -0.09283  0.01626
impervious     0.1733  0.2595  0.2952  0.33262  0.41112
logPopulation -0.3668 -0.2882 -0.2515 -0.21677 -0.13994

summary(mod)


 Iterations = 10001:59999
 Thinning interval  = 13
 Sample size  = 3847 

 DIC: 804.353 

 G-structure:  ~Location

         post.mean  l-95% CI u-95% CI eff.samp
Location   0.04611 0.0003358   0.1396     2494

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: cbind(cenMin, cenMax) ~ impervious + logPopulation 

              post.mean l-95% CI u-95% CI eff.samp   pMCMC    
(Intercept)    -0.14338 -0.30841  0.01522     3847 0.08058 .  
impervious      0.29481  0.17676  0.41300     3422 < 3e-04 ***
logPopulation  -0.25274 -0.37335 -0.14813     4061 0.00052 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
bayes_plots(mod, coc, df)

[[1]]


[[2]]


[[3]]

Total Phosphorus

coc = "Total Phosphorus - Water - Total"
TPModel <- mcmc_calc("Total Phosphorus - Water - Total", c("rev_logTraffic"), "log")
mod <- TPModel
summary(mod$Sol)


Iterations = 10001:59999
Thinning interval = 13 
Number of chains = 1 
Sample size per chain = 3847 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                Mean     SD Naive SE Time-series SE
(Intercept)    7.231 0.7644 0.012324       0.012324
rev_logTraffic 1.715 0.4762 0.007678       0.007678

2. Quantiles for each variable:

                 2.5%   25%   50%   75% 97.5%
(Intercept)    5.7023 6.758 7.243 7.706 8.722
rev_logTraffic 0.7735 1.416 1.715 2.024 2.654

summary(mod)


 Iterations = 10001:59999
 Thinning interval  = 13
 Sample size  = 3847 

 DIC: 952.3584 

 G-structure:  ~Location

         post.mean l-95% CI u-95% CI eff.samp
Location    0.2531  0.06908    0.526     3847

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: cbind(cenMin_log, cenMax_log) ~ rev_logTraffic 

               post.mean l-95% CI u-95% CI eff.samp   pMCMC    
(Intercept)       7.2309   5.7004   8.7220     3847 < 3e-04 ***
rev_logTraffic    1.7146   0.7724   2.6537     3847 0.00208 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
bayes_plots(mod, coc, df)

[[1]]


[[2]]


[[3]]

Total Kjeldahl Nitrogen

coc = "Total Kjeldahl Nitrogen - Water - Total"
TKNModel <- mcmc_calc("Total Kjeldahl Nitrogen - Water - Total", "rev_logTraffic", 
    "log")
mod <- TKNModel
df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
bayes_plots(mod, coc, df)

[[1]]


[[2]]


[[3]]

## TSS

coc = "Total Suspended Solids - Water - Total"
TSSModel <- mcmc_calc("Total Suspended Solids - Water - Total", "rev_logTraffic", 
    "log")
df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
bayes_plots(TSSModel, coc, df)

[[1]]


[[2]]


[[3]]

‘Total PAH - Water - Total’,

coc = "Total PAH - Water - Total"
PAHModel <- mcmc_calc("Total PAH - Water - Total", "LU_ratio", "log")
mod <- PAHModel
df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
bayes_plots(mod, coc, df)

[[1]]


[[2]]


[[3]]

‘CPAH - Water - Total’

coc = "CPAH - Water - Total"
CPAHModel <- mcmc_calc("CPAH - Water - Total", c("logPopulation", "rev_logTraffic"), 
    "log")
mod <- CPAHModel
df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
bayes_plots(mod, coc, df)

[[1]]


[[2]]


[[3]]

‘HPAH - Water - Total’

coc = "HPAH - Water - Total"
HPAHModel <- mcmc_calc("HPAH - Water - Total", c("LU_ratio"), "log")
mod <- HPAHModel
df <- (subset(s8data.wPredictors, parameter == coc)) %>% add_surv()
bayes_plots(mod, coc, df)

[[1]]


[[2]]


[[3]]

Summary

Summarize posterior results for use in heatmap.

# summarize
metals <- list()
metals[["Total Copper"]] <- (CuModel)
metals[["Total Zinc"]] <- (ZnModel)
metals[["Total Cadmium"]] <- (CdModel)
metals[["Total Lead"]] <- (PbModel)

others <- list()
others[["Total Phosphorus"]] <- TPModel
others[["Total Kjeldahl Nitrogen"]] <- TKNModel
others[["Total Suspended Sediment"]] <- TSSModel
# others['FC'] <- FCModel





msummary(metals, title = "Total Metals", statistic_vertical = TRUE, statistic = "conf.int", 
    conf_level = 0.95)

Total Metals

	Total Copper	Total Zinc	Total Cadmium	Total Lead
(Intercept)	5.405	3.640	-0.143	1.198
	[3.647, 7.199]	[3.210, 4.055]	[-0.308, 0.015]	[0.591, 1.849]
impervious	0.304	0.714	0.295	0.647
	[0.049, 0.560]	[0.418, 1.031]	[0.177, 0.413]	[0.191, 1.111]
rev_logTraffic	2.207
	[1.025, 3.231]
logPopulation			-0.253
			[-0.373, -0.148]
Num.Obs.	0	0	0	0
dic	946.955	928.197	804.353	1039.448
n	438.000	414.000	416.000	417.000

msummary(others, title = "Nutrients and TSS", statistic_vertical = TRUE, statistic = "conf.int", 
    conf_level = 0.95)

Nutrients and TSS

	Total Phosphorus	Total Kjeldahl Nitrogen	Total Suspended Sediment
(Intercept)	7.231	8.766	13.236
	[5.700, 8.722]	[7.932, 9.720]	[11.466, 14.909]
rev_logTraffic	1.715	1.396	1.970
	[0.772, 2.654]	[0.849, 1.977]	[0.902, 3.073]
Num.Obs.	0	0	0
dic	952.358	823.060	1183.863
n	424.000	397.000	415.000

Print out a messy list of results for the heatmap

sols <- list()
sols[["Copper"]] <- summary(CuModel)$solutions
sols[["Zinc"]] <- summary(ZnModel)$solutions
sols[["Cadmium"]] <- summary(CdModel)$solutions
sols[["Lead"]] <- summary(PbModel)$solutions
sols[["TP"]] <- summary(TPModel)$solutions
sols[["TKN"]] <- summary(TKNModel)$solutions
sols[["TSS"]] <- summary(TSSModel)$solutions

for (n in 1:length(sols)) {
    print(kable(sols[[n]], row.names = TRUE, caption = paste(names(sols)[n], " - Summary of MCMC model")))
}



Table: Copper  - Summary of MCMC model

                  post.mean    l-95% CI    u-95% CI   eff.samp       pMCMC
---------------  ----------  ----------  ----------  ---------  ----------
(Intercept)        5.405015   3.6471017   7.1988047   4085.710   0.0002599
rev_logTraffic     2.206573   1.0254637   3.2313848   4061.867   0.0015597
impervious         0.303770   0.0492773   0.5597377   3836.965   0.0213153


Table: Zinc  - Summary of MCMC model

               post.mean    l-95% CI   u-95% CI   eff.samp       pMCMC
------------  ----------  ----------  ---------  ---------  ----------
(Intercept)    3.6401423   3.2097896   4.055369       3847   0.0002599
impervious     0.7136934   0.4176574   1.030960       3847   0.0002599


Table: Cadmium  - Summary of MCMC model

                  post.mean     l-95% CI     u-95% CI   eff.samp       pMCMC
--------------  -----------  -----------  -----------  ---------  ----------
(Intercept)      -0.1433757   -0.3084146    0.0152214   3847.000   0.0805823
impervious        0.2948066    0.1767649    0.4130013   3421.656   0.0002599
logPopulation    -0.2527430   -0.3733480   -0.1481330   4060.526   0.0005199


Table: Lead  - Summary of MCMC model

               post.mean    l-95% CI   u-95% CI   eff.samp       pMCMC
------------  ----------  ----------  ---------  ---------  ----------
(Intercept)    1.1980017   0.5906497   1.848978   3847.000   0.0036392
impervious     0.6474842   0.1905426   1.111220   3665.281   0.0088381


Table: TP  - Summary of MCMC model

                  post.mean    l-95% CI   u-95% CI   eff.samp       pMCMC
---------------  ----------  ----------  ---------  ---------  ----------
(Intercept)        7.230941   5.7004380   8.721957       3847   0.0002599
rev_logTraffic     1.714623   0.7723525   2.653688       3847   0.0020795


Table: TKN  - Summary of MCMC model

                  post.mean    l-95% CI   u-95% CI   eff.samp       pMCMC
---------------  ----------  ----------  ---------  ---------  ----------
(Intercept)        8.765541   7.9324569   9.720447       3847   0.0002599
rev_logTraffic     1.395636   0.8487623   1.976888       3847   0.0002599


Table: TSS  - Summary of MCMC model

                  post.mean     l-95% CI    u-95% CI   eff.samp       pMCMC
---------------  ----------  -----------  ----------  ---------  ----------
(Intercept)       13.235739   11.4664881   14.908653       3847   0.0002599
rev_logTraffic     1.969567    0.9022194    3.073489       3847   0.0025994

Watershed Regression - Part 2: Censored MCMC Generalized Linear Mixed Models

Christian Nilsen