Watershed Regression - Part 1: Linear Mixed Models (adapted to GLS Survival Analysis)

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 uses spatial predictors and outfall monitoring data to develop predictive linear regression relationships. This represents the first step in the watershed regression for the stormwater heatmap. The second step is a censored bayesian model, which adds some capabilities that simple linear models do not have. This first step is necessary to select the most predictive parameters in an efficient way.

Here are the steps in the code below:

1. Code set up
   
2. Get and prepare monitoring data
   
3. Merge spatial data with monitoring data
   
4. Remove mulitcolinear predictors
   
5. Peform model selection
   
6. 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 )

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",
                        stringsAsFactors = FALSE )

#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 
kable(params,col.names = "Constituents of concern")

Constituents of concern
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
CPAH - Water - Total
HPAH - Water - Total

#

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:

#make a function for scatter plots 
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
)+
  theme_ipsum_rc()  
p + facet_wrap( ~ parameter, scales = 'free')+theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) 
}

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("data/spatialPredictors_4_13.csv", col_types = cols(X1 = col_skip()))
spatial_data <- read_csv("data/spatialPredictors_5_14.csv")

Parsed with column specification:
cols(
  `system:index` = col_character(),
  COM = col_double(),
  Location = col_character(),
  RES = col_double(),
  change_year_index = col_double(),
  depSplusN = col_double(),
  impervious = col_double(),
  logPopulation = col_double(),
  logTraffic = col_double(),
  nighttime_lights = col_double(),
  pm25 = col_double(),
  roadDensity = col_double(),
  soil = col_double(),
  .geo = col_character()
)

#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))

kable(head(s8data.wPredictors),caption = "S8 Monitoring Data")

Location	study_name	parameter	type	season	concentration	nondetect_flag	study_id	access_id	field_collection_end_date	field_collection_start_date	system:index	change_year_index	impervious	logPopulation	logTraffic	nighttime_lights	pm25	roadDensity	soil	LU_ratio
KICCOMS8D_OUT	King County Phase I Municipal Stormwater Permit	Copper - Water - Total	COM	4	22.800	FALSE	WAR044501_S8D	36149	11/2/2011	11/2/2011	0000000000000000000e	-0.084	-0.039	-0.564	-0.002	-1.68	-0.75	0.033	-1.37	-0.335
KICCOMS8D_OUT	King County Phase I Municipal Stormwater Permit	LPAH - Water - Total	COM	1	0.060	FALSE	WAR044501_S8D	102225	NA	2/12/2011	0000000000000000000e	-0.084	-0.039	-0.564	-0.002	-1.68	-0.75	0.033	-1.37	-0.335
KICCOMS8D_OUT	King County Phase I Municipal Stormwater Permit	Dichlobenil - Water - Total	COM	4	0.047	TRUE	WAR044501_S8D	35675	11/1/2010	11/1/2010	0000000000000000000e	-0.084	-0.039	-0.564	-0.002	-1.68	-0.75	0.033	-1.37	-0.335
KICCOMS8D_OUT	King County Phase I Municipal Stormwater Permit	Cadmium - Water - Total	COM	4	0.065	FALSE	WAR044501_S8D	33689	10/19/2012	10/18/2012	0000000000000000000e	-0.084	-0.039	-0.564	-0.002	-1.68	-0.75	0.033	-1.37	-0.335
KICCOMS8D_OUT	King County Phase I Municipal Stormwater Permit	Fecal coliform - Water - Total	COM	4	530.000	FALSE	WAR044501_S8D	39755	12/8/2010	12/8/2010	0000000000000000000e	-0.084	-0.039	-0.564	-0.002	-1.68	-0.75	0.033	-1.37	-0.335
KICCOMS8D_OUT	King County Phase I Municipal Stormwater Permit	Bis(2-ethylhexyl) phthalate - Water - Total	COM	4	1.270	TRUE	WAR044501_S8D	38768	11/23/2012	11/23/2012	0000000000000000000e	-0.084	-0.039	-0.564	-0.002	-1.68	-0.75	0.033	-1.37	-0.335

getBaseFormula <- function(df.coc) {
  #Function to make a formula for prediction
  predictors <- df.coc %>%
    select_if(is.numeric) %>%
    dplyr::select(-c(concentration, access_id)) %>%
    colnames()
  return(as.formula(paste(
    "(concentration) ~",  (paste((predictors),  collapse = " + ")), " + (1|Location)"
  )))
  
}

Perform Linear Regression

Below, we perform linear regression one parameter at a time.

Functions for regression

We set up a series of functions to make the analysis easier.

Data dictionary for these functions

coc: constituent of concern
df.coc: filtered dataframe
model: basic linear mixed model
log_model: linear-mixed model with log-transformed concentrations
model_with_seasonality: linear-mixed model with seasonality as an added factor
model_with_seasonality_log: linear-mixed model with seasonality as an added factor with log-transformed concentrations

Plotting function

Basic function to plot results by location

plot_s8 <- function(coc) {
    df.coc <- (base::subset(s8data.wPredictors,
                parameter == coc))
  #plot the data to inspect 
  plot <- ggplot(data=(df.coc))+geom_point(aes(x=Location,y=concentration),alpha = 0.5, color = "#0088a3")+scale_y_log10()+labs(
      y = "Concentration (µg/L)",
      x = "Location",
      title =  "Measured Concentrations",
      subtitle = coc ) + theme_ipsum_rc()
  return(plot)
}

Function to remove mulitcolinear predictors in base model

To address multicolinearity, we calculate the variance inflation factor (VIF) and iteratively remove paramters with the highest VIF. We keep removing parameters one at a time until all VIF values are below 5.0.

# #Calculate variance-inflation factors. 
# 
# calc_vif_base <- function(coc) {
#   df.coc <- (base::subset(s8data.wPredictors,
#                 parameter == coc))
#   base_formula <- getBaseFormula(df.coc)#returns a formula with all predictors
#   model.1 <<- lmer(base_formula, data = df.coc, na.action = na.omit)
#   v <- sort(vif(model.1),decreasing=TRUE)
#   return(v)
#   #kable(v,caption="Variance Inflation Factors - not filtered") #display the vif of the dataset
# }

Function to iteratively remove mulitcolinear predictors

Same as above, expect multicolinear predictors are removed.

## *Check VIF
 
check_vif <- function(coc) {
  df.coc <- (base::subset(s8data.wPredictors,
                parameter == coc))
  base_formula <- getBaseFormula(df.coc)#returns a formula with all predictors
  model.1 <- lmer(base_formula, data = df.coc, na.action = na.omit) #make into a lmer object 
  v <- sort(vif(model.1),decreasing=TRUE)
  
  #if the VIF of the highest ranked predictor is >10 then iteratively remove
  model_object <- model.1 #start with model object as the base model (all predictors included)
  
  
  for (i in 1:20) {
    interim_v <- sort(vif(model_object), decreasing = TRUE)
    if (max(interim_v) < 10) {
      break
    }
    predictor_to_drop = as.name(names(interim_v)[which(interim_v == max(interim_v))])
    model_object <-
      stats::update(model_object, paste(".~ . -", predictor_to_drop))
    }
  
  m1Terms <- (labels(terms(model.1)))
  m2Terms <- labels(terms(model_object))
  
  #compare the terms to get a list of the dropped terms
  droppedTerms <- setdiff(m1Terms, m2Terms)
  
  #make a list of selected predictors
  predictors <- m2Terms#colnames(model.frame(model_object)) 

  
  #filter df.coc to remove dropped terms.
  df.coc = dplyr::select(df.coc, -(droppedTerms))
  return(list("vif" = interim_v,"dropped" = droppedTerms,"predictors" = predictors))
  #kable(droppedTerms,caption = "These terms were dropped")
  
  #kable(interim_v,caption="Variance Inflation Factors - multicolinear factors dropped")
}

Function to perform stepwise selection and return a series of best models

## Stepwise Selection

#forward_selection(TRUE,coc,model_info$predictors)

# Extract the model that step found:

#perform forward selection on model parameters. First for non-transformed data, then for log-transformed data
forward_selection <- function(seasonal.bin, coc,predictors) {
  #seasonal.bin = binary (T/F) if seasonal model should be used
  library(lmerTest)
  #make this a lmer object 
  df.coc <- (base::subset(s8data.wPredictors,
                parameter == coc)) 
  model_object_formula <- as.formula(paste(
    "concentration ~",  (paste((predictors),  collapse = " + ")), " + (1|Location)"))
  
  model_object <- lmer(model_object_formula,data=df.coc)
  
  step.2 <-  lmerTest::step(model_object,reduce.random=FALSE,data=df.coc)
  step.2.log <- lmerTest::step(stats::update(model_object, log(concentration)~.))
  
  #extract the models 
  model.3 <- get_model(step.2)
  model.3.log <- get_model(step.2.log)
  
  #perform forward selection on model parameters, this time add seasonality . First for non-transformed data, then for log-transformed data
  step.4 <- lmerTest::step(stats::update(model_object,.~.+season))
  step.4.log <- lmerTest::step(stats::update(model_object,log(concentration) ~.+season))
  model.4 <- get_model(step.4)
  model.4.log <- get_model(step.4.log)
  
  #get formulas 
  model.3.formula <- as.formula(model.3@call$formula)
  model.3.log.formula <- as.formula(model.3.log@call$formula)
  model.4.formula <- as.formula(model.4@call$formula)
  model.4.log.formula <- as.formula(model.4.log@call$formula)
  
  #detach lmer test and remove the models. Keep the formulas. 
  detach("package:lmerTest", unload=TRUE)
  rm(model.3,model.3.log,model.4,model.4.log)
  
  #use lmer for performing modeling
  #calculate base model 
  df.coc <- (base::subset(s8data.wPredictors,
                parameter == coc)) 
  model.base <- lmer(model_object_formula,data=df.coc)
  model <- lmer(model.3.formula,data=df.coc)
  log_model<-  lmer(model.3.log.formula,data=df.coc)
  
  if(seasonal.bin) {
      #if seasonal model switch is on, calc seasonal models
    model_with_seasonality <- lmer(model.4.formula,data=df.coc)
    model_with_seasonality_log <- lmer(model.4.log.formula,data=df.coc)
      #add to list 
    modelList <- c(model,log_model,model_with_seasonality,model_with_seasonality_log)
    
    }
    else {
    modelList <-c(model,log_model)
    modelLables <-c('linear','log-linear') 
    }
   

 
 return(modelList)#,show.aic = TRUE)# = FALSE, title=coc, dv.labels = modelLables) 
  # #make a table of coefficients 
  # tab_model(modelList,
  #   model_log,
  #   model_with_seasonality,
  #   model.with_seasonality_log,
  #   show.aic = TRUE,auto.label = FALSE, title=coc, dv.labels = c('linear','log-linear','linear seasonal','log-linear seasonal'))#file=paste0("results/",coc,".html"))
}

Run through model selection functions

Now we use the functions from above to return model tables. We wrap a helper function to call the others.

results = c()
plots = c()
vifs = c()
tabs = c()
for (i in 1:length(params)){
coc = params[i]
df.coc <- (base::subset(s8data.wPredictors,
                parameter == coc)) 
plots[[coc]] <- plot_s8(coc)
model_info <- check_vif(coc)
vifs[[coc]] <- model_info$vif

results[[coc]] = forward_selection(TRUE,coc,model_info$predictors)}
modelLabels <-c('linear','log-linear','linear seasonal','log-linear seasonal') 

Wrapper function for displaying results for individual cocs.

show_results <- function(j){
#get parameter label 
lab = names(results)[j]

#plot observed 
print(plots[[j]])

#show variance inflation factors 
print(kable(vifs[[j]],caption = paste("Variance inflation factors",lab),col.names = c("vif")))

#show raw summary of results 
models <- results[[j]]
for (k in 2:2){
  print((summary(models[[k]])))
  #plot(models[[k]],,main=paste(lab,"\n","Resididuals"))
  print(qqmath(models[[k]],main=paste(lab,"\n",modelLabels[k],"\n","QQ plot of resididuals")))
}
  

#formatted table of results 
tabs = (tab_model(models,title=lab,show.aic = TRUE,dv.labels =modelLabels))

kable(tabs$knitr) #displays the table inline

}

Now we can specify one parameters at a time and display the results. Generally, the model with the lowest AIC is used in Step 2 (Bayesian modeling).

Zinc

(show_results(1))

	vif
impervious	9.21
change_year_index	8.94
nighttime_lights	4.00
roadDensity	3.60
logPopulation	2.70
LU_ratio	2.63

Linear mixed model fit by REML ['lmerMod']
Formula: log(concentration) ~ change_year_index + (1 | Location)
   Data: df.coc

REML criterion at convergence: 912

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.886 -0.524 -0.038  0.466  3.302 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.468    0.684   
 Residual             0.472    0.687   
Number of obs: 414, groups:  Location, 14

Fixed effects:
                  Estimate Std. Error t value
(Intercept)          4.042      0.198   20.37
change_year_index    1.013      0.189    5.36

Correlation of Fixed Effects:
            (Intr)
chng_yr_ndx 0.341 

x

|

Zinc - Water - Total

	linear			log-linear			linear seasonal			log-linear seasonal
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	92.44	67.34 – 117.55	<0.001	4.04	3.65 – 4.43	<0.001	77.16	48.87 – 105.45	<0.001	4.01	3.62 – 4.41	<0.001
impervious	64.81	40.18 – 89.43	<0.001				63.37	39.35 – 87.38	<0.001
change_year_index				1.01	0.64 – 1.38	<0.001				1.00	0.64 – 1.37	<0.001
season [2]							4.22	-25.91 – 34.34	0.784	-0.06	-0.26 – 0.13	0.521
season [3]							114.53	78.88 – 150.17	<0.001	0.39	0.16 – 0.62	0.001
season [4]							7.93	-16.74 – 32.60	0.529	0.00	-0.16 – 0.16	0.981
Random Effects
σ2	12496.95			0.47			11366.26			0.46
τ00	1809.55 Location			0.47 Location			1739.85 Location			0.46 Location
ICC	0.13			0.50			0.13			0.50
N	14 Location			14 Location			14 Location			14 Location
Observations	414			414			414			414
Marginal R2 / Conditional R2	0.229 / 0.326			0.484 / 0.741			0.290 / 0.384			0.491 / 0.745
AIC	5095.101			920.051			5038.048			920.701

show_results(2)

	vif
impervious	9.29
change_year_index	9.26
nighttime_lights	3.99
roadDensity	3.79
LU_ratio	2.64
logPopulation	2.61

Linear mixed model fit by REML ['lmerMod']
Formula: concentration ~ roadDensity + (1 | Location)
   Data: df.coc

REML criterion at convergence: 3141

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.124 -0.474 -0.107  0.175  5.633 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 47.5     6.89    
 Residual             70.5     8.40    
Number of obs: 438, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)    20.30       2.31    8.80
roadDensity    15.11       2.74    5.52

Correlation of Fixed Effects:
            (Intr)
roadDensity 0.574 
Linear mixed model fit by REML ['lmerMod']
Formula: log(concentration) ~ roadDensity + (1 | Location)
   Data: df.coc

REML criterion at convergence: 932

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-5.397 -0.532 -0.001  0.552  3.051 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.313    0.559   
 Residual             0.444    0.667   
Number of obs: 438, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)    2.559      0.187   13.68
roadDensity    1.263      0.222    5.69

Correlation of Fixed Effects:
            (Intr)
roadDensity 0.574 

Nitrite-Nitrate

	vif
impervious	9.27
change_year_index	9.24
nighttime_lights	3.98
roadDensity	3.79
LU_ratio	2.63
logPopulation	2.60

show_results(3)

Linear mixed model fit by REML ['lmerMod']
Formula: concentration ~ (1 | Location)
   Data: df.coc

REML criterion at convergence: 6108

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.944 -0.408 -0.154  0.239  6.027 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 143431   379     
 Residual             185115   430     
Number of obs: 406, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)      480        104    4.63
Linear mixed model fit by REML ['lmerMod']
Formula: log(concentration) ~ (1 | Location)
   Data: df.coc

REML criterion at convergence: 963

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.181 -0.615 -0.017  0.591  3.238 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.571    0.755   
 Residual             0.558    0.747   
Number of obs: 406, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)    5.689      0.206    27.6

Error in eval(parse(text = x, keep.source = FALSE)[[1L]]) : 
  object 'logseason' not found

Lead

show_results(4)

	vif
impervious	9.22
change_year_index	9.00
nighttime_lights	3.99
roadDensity	3.63
logPopulation	2.67
LU_ratio	2.63

Linear mixed model fit by REML ['lmerMod']
Formula: concentration ~ logPopulation + roadDensity + (1 | Location)
   Data: df.coc

REML criterion at convergence: 3233

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.342 -0.453 -0.105  0.026  9.890 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept)  19.6     4.43   
 Residual             132.5    11.51   
Number of obs: 417, groups:  Location, 14

Fixed effects:
              Estimate Std. Error t value
(Intercept)      13.48       1.65    8.19
logPopulation    -4.93       1.55   -3.17
roadDensity      11.77       2.66    4.43

Correlation of Fixed Effects:
            (Intr) lgPplt
logPopulatn -0.247       
roadDensity  0.561 -0.684
Linear mixed model fit by REML ['lmerMod']
Formula: log(concentration) ~ logPopulation + roadDensity + (1 | Location)
   Data: df.coc

REML criterion at convergence: 1077

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.732 -0.627 -0.016  0.577  3.301 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.858    0.926   
 Residual             0.689    0.830   
Number of obs: 417, groups:  Location, 14

Fixed effects:
              Estimate Std. Error t value
(Intercept)      2.039      0.317    6.44
logPopulation   -0.764      0.298   -2.56
roadDensity      2.129      0.500    4.26

Correlation of Fixed Effects:
            (Intr) lgPplt
logPopulatn -0.247       
roadDensity  0.574 -0.683

Total Phosphorus

show_results(5)

	vif
impervious	9.27
change_year_index	9.25
nighttime_lights	3.98
roadDensity	3.78
LU_ratio	2.64
logPopulation	2.60

Linear mixed model fit by REML ['lmerMod']
Formula: concentration ~ roadDensity + (1 | Location)
   Data: df.coc

REML criterion at convergence: 5118

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.922 -0.524 -0.213  0.278  7.274 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 3852     62.1    
 Residual             9811     99.1    
Number of obs: 424, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)    183.0       21.1    8.66
roadDensity     88.3       25.1    3.52

Correlation of Fixed Effects:
            (Intr)
roadDensity 0.572 
Linear mixed model fit by REML ['lmerMod']
Formula: log(concentration) ~ roadDensity + (1 | Location)
   Data: df.coc

REML criterion at convergence: 905

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.637 -0.669  0.003  0.673  3.338 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.257    0.507   
 Residual             0.449    0.670   
Number of obs: 424, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)    4.875      0.171   28.57
roadDensity    0.708      0.203    3.49

Correlation of Fixed Effects:
            (Intr)
roadDensity 0.574 

Total Suspended Solids

show_results(6)

	vif
impervious	9.27
change_year_index	9.20
nighttime_lights	3.98
roadDensity	3.77
LU_ratio	2.64
logPopulation	2.64

Linear mixed model fit by REML ['lmerMod']
Formula: concentration ~ roadDensity + (1 | Location)
   Data: df.coc

REML criterion at convergence: 10253

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.071 -0.457 -0.266  0.048  6.526 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 6.76e+08 26004   
 Residual             3.26e+09 57092   
Number of obs: 415, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)    64942       9143    7.10
roadDensity    28652      10909    2.63

Correlation of Fixed Effects:
            (Intr)
roadDensity 0.569 
Linear mixed model fit by REML ['lmerMod']
Formula: log(concentration) ~ roadDensity + (1 | Location)
   Data: df.coc

REML criterion at convergence: 1234

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.196 -0.584  0.016  0.557  3.276 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.263    0.513   
 Residual             1.064    1.031   
Number of obs: 415, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)   10.553      0.178   59.13
roadDensity    0.883      0.213    4.15

Correlation of Fixed Effects:
            (Intr)
roadDensity 0.570 

Total Phthalate

show_results(7)

	vif
impervious	9.21
change_year_index	9.04
nighttime_lights	3.99
roadDensity	3.65
logPopulation	2.66
LU_ratio	2.63

Linear mixed model fit by REML ['lmerMod']
Formula: concentration ~ change_year_index + impervious + (1 | Location)
   Data: df.coc

REML criterion at convergence: 2333

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.501 -0.305 -0.150 -0.003  9.027 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept)  2.06    1.44    
 Residual             16.21    4.03    
Number of obs: 412, groups:  Location, 14

Fixed effects:
                  Estimate Std. Error t value
(Intercept)          2.239      0.524    4.27
change_year_index   -2.551      1.099   -2.32
impervious           3.167      1.044    3.03

Correlation of Fixed Effects:
            (Intr) chng__
chng_yr_ndx  0.546       
impervious  -0.478 -0.911
Linear mixed model fit by REML ['lmerMod']
Formula: log(concentration) ~ (1 | Location)
   Data: df.coc

REML criterion at convergence: 1008

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.095 -0.717 -0.114  0.485  4.008 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.580    0.762   
 Residual             0.604    0.777   
Number of obs: 412, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)    0.418      0.208    2.01

Total PAH

show_results(8)

	vif
impervious	9.22
change_year_index	9.07
nighttime_lights	3.99
roadDensity	3.68
logPopulation	2.65
LU_ratio	2.63

Linear mixed model fit by REML ['lmerMod']
Formula: concentration ~ (1 | Location)
   Data: df.coc

REML criterion at convergence: 1267

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.648 -0.311 -0.159  0.016  9.223 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.271    0.521   
 Residual             1.182    1.087   
Number of obs: 412, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)    0.592      0.151    3.93
Linear mixed model fit by REML ['lmerMod']
Formula: log(concentration) ~ (1 | Location)
   Data: df.coc

REML criterion at convergence: 1193

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.676 -0.448 -0.183  0.542  4.240 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.566    0.752   
 Residual             0.960    0.980   
Number of obs: 412, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)   -1.417      0.208   -6.82

Error in eval(parse(text = x, keep.source = FALSE)[[1L]]) : 
  object 'logseason' not found

CPAH

show_results(9)

	vif
impervious	9.19
change_year_index	8.91
nighttime_lights	4.00
roadDensity	3.56
logPopulation	2.70
LU_ratio	2.64

Linear mixed model fit by REML ['lmerMod']
Formula: concentration ~ (1 | Location)
   Data: df.coc

REML criterion at convergence: -528

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.064 -0.266 -0.122 -0.046 10.387 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.00242  0.0492  
 Residual             0.01512  0.1230  
Number of obs: 412, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)   0.1164     0.0147    7.94
Linear mixed model fit by REML ['lmerMod']
Formula: log(concentration) ~ (1 | Location)
   Data: df.coc

REML criterion at convergence: 856

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.461 -0.274 -0.056  0.129  4.487 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.589    0.768   
 Residual             0.412    0.642   
Number of obs: 412, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)   -2.518      0.208   -12.1

HPAH

show_results(10)

	vif
impervious	9.19
change_year_index	8.90
nighttime_lights	4.00
roadDensity	3.55
logPopulation	2.71
LU_ratio	2.65

Linear mixed model fit by REML ['lmerMod']
Formula: concentration ~ logPopulation + roadDensity + (1 | Location)
   Data: df.coc

REML criterion at convergence: 757

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.493 -0.248 -0.119  0.064 11.596 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.0346   0.186   
 Residual             0.3454   0.588   
Number of obs: 412, groups:  Location, 14

Fixed effects:
              Estimate Std. Error t value
(Intercept)     0.1854     0.0721    2.57
logPopulation   0.2016     0.0679    2.97
roadDensity    -0.2976     0.1172   -2.54

Correlation of Fixed Effects:
            (Intr) lgPplt
logPopulatn -0.248       
roadDensity  0.557 -0.683
Linear mixed model fit by REML ['lmerMod']
Formula: log(concentration) ~ (1 | Location)
   Data: df.coc

REML criterion at convergence: 1162

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.548 -0.362 -0.128  0.249  4.158 

Random effects:
 Groups   Name        Variance Std.Dev.
 Location (Intercept) 0.497    0.705   
 Residual             0.892    0.944   
Number of obs: 412, groups:  Location, 14

Fixed effects:
            Estimate Std. Error t value
(Intercept)   -2.140      0.195     -11

Total Kjeldahl Nitrogen - Water - Total

	vif
impervious	9.26
change_year_index	8.69
nighttime_lights	4.00
roadDensity	3.35
logPopulation	2.80
LU_ratio	2.73

show_results(12)

Error in results[[j]] : subscript out of bounds