Intricate watershed characteristics and chemical properties govern methyl mercury concentrations within headwater streams of boreal forests in Ontario, Canada: a structural equation modelling approach
Tyler J. Harrow-Lyle, Wai Ying Lam, Erik J.S. Emilson, Robert W. Mackereth, Carl P.J. Mitchell, and Stephanie J. Melles
2022-10-23
Descriptive statistics
setwd("E:/Postdoc- Stephanie Melles/Paper 1 sem submission")
library(psych)
mt<-read.csv("summary statistics.csv")
describe(mt)## vars n mean sd median trimmed mad min max
## Site* 1 176 9.94 5.57 10.00 9.93 7.41 1.00 19.00
## Temp.C 2 165 12.83 4.94 13.10 12.84 5.49 1.80 25.40
## Mt.ppt 3 176 0.23 0.23 0.15 0.19 0.16 0.02 1.21
## DOC.mg.L 4 170 16.25 7.37 14.95 15.91 7.03 0.09 37.35
## Nitrate.ppm 5 103 0.16 0.18 0.11 0.13 0.05 0.05 1.64
## Sulphate.ppm 6 103 0.84 0.90 0.47 0.66 0.42 0.10 5.42
## Phosphate.ppm 7 103 0.37 0.22 0.17 0.37 0.01 0.16 0.62
## Iron 8 143 381.33 282.01 310.23 347.93 215.42 0.00 2134.70
## Precip_wk 9 176 3.43 3.17 2.64 2.91 2.37 0.00 15.06
## Yr_Dist. 10 176 7.12 9.09 2.00 6.06 2.97 0.00 24.00
## Geology 11 176 3.60 1.12 4.00 3.62 0.00 1.00 6.00
## Mixed 12 176 19.50 8.77 22.33 19.67 9.39 3.75 36.95
## Wet 13 176 35.75 13.58 35.17 34.93 10.28 12.61 65.00
## WV5 14 176 0.01 0.00 0.01 0.01 0.00 0.00 0.01
## Wet.Cover 15 176 3.62 3.47 2.36 3.12 3.30 0.01 11.66
## range skew kurtosis se
## Site* 18.00 0.02 -1.26 0.42
## Temp.C 23.60 -0.06 -0.74 0.38
## Mt.ppt 1.19 1.87 4.22 0.02
## DOC.mg.L 37.26 0.40 -0.22 0.57
## Nitrate.ppm 1.59 5.57 39.27 0.02
## Sulphate.ppm 5.31 2.36 7.03 0.09
## Phosphate.ppm 0.46 0.02 -2.02 0.02
## Iron 2134.70 2.23 9.55 23.58
## Precip_wk 15.06 1.86 3.86 0.24
## Yr_Dist. 24.00 0.84 -1.03 0.68
## Geology 5.00 -0.45 0.29 0.08
## Mixed 33.20 -0.19 -0.69 0.66
## Wet 52.39 0.55 -0.27 1.02
## WV5 0.01 -0.55 -0.14 0.00
## Wet.Cover 11.65 1.15 0.25 0.26#Scaled summary statistics for roughness and phosphorus
scaledstatrough<-read.csv("roughness scaled summary statistics.csv")
describe(scaledstatrough)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 176 8.95 2 9.3 9.06 1.63 4.43 12.08 7.65 -0.55 -0.14 0.15scaledstatphosphate<-read.csv("phosphate scaled summary statistics.csv")
describe(scaledstatphosphate)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 103 3.75 2.21 1.67 3.73 0.15 1.57 6.22 4.65 0.02 -2.02 0.22#Principal component analyses - Figure 2 and Figure S1.
setwd("E:/Postdoc- Stephanie Melles/Paper 1 sem submission")
#Figure 2.
library(FactoMineR)
library(factoextra)## Loading required package: ggplot2##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alpha## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBalibrary(ggplot2)
mercurypca1<-read.csv("pcafigure2.csv")
head(mercurypca1)## Geology Geology.1 Region
## 1 Basaltic and andesitic flows, iron formation A CF
## 2 Basaltic and andesitic flows, iron formation A CF
## 3 Basaltic and andesitic flows, iron formation A CF
## 4 Basaltic and andesitic flows, iron formation A CF
## 5 Basaltic and andesitic flows, iron formation A CF
## 6 Basaltic, dacitic, and andesitic flows, iron formations B CF
## mhix Temp.C Mt.ppt DOC.mg.L Nitrate.ppm Sulphate.ppm Phosphate.ppm
## 1 0.9418975 5.7 0.044 9.16 0.154 3.093 0.592
## 2 0.9294352 7.7 0.020 4.68 0.116 2.152 0.157
## 3 0.9449565 11.4 0.119 19.70 0.171 4.293 0.592
## 4 0.9489639 10.5 0.051 9.77 0.166 2.442 0.157
## 5 0.9433831 9.3 0.075 7.75 0.181 2.592 0.599
## 6 0.9442036 10.3 0.630 19.34 0.168 0.429 0.592
## Iron Precip_wk Yr_Dist. Wet WV5 Wet.Cover
## 1 247.4900 1.571429 0 55.71429 0.009663 0.01
## 2 232.9600 2.128571 0 55.71429 0.009663 0.01
## 3 742.3900 3.500000 0 55.71429 0.009663 0.01
## 4 240.8200 2.271429 0 55.71429 0.009663 0.01
## 5 352.9500 1.157143 0 55.71429 0.009663 0.01
## 6 457.2041 2.128571 0 18.18182 0.004431 0.01#Subsetting the data and then standardizing the data
presencelanduse.pca1 <- PCA(mercurypca1[,4:16], graph = FALSE,scale=TRUE)
fviz_contrib(presencelanduse.pca1, choice = "ind", axes = 1:2)
res.desc <- dimdesc(presencelanduse.pca1, axes = c(1,2), proba = 0.05)
res.desc$Dim.1## $quanti
## correlation p.value
## Wet.Cover 0.6518443 5.582980e-10
## Temp.C 0.4747458 2.515740e-05
## DOC.mg.L 0.4721889 2.819748e-05
## WV5 0.4315947 1.535403e-04
## Wet -0.5934005 3.946600e-08
## Sulphate.ppm -0.7770315 1.026552e-15
##
## attr(,"class")
## [1] "condes" "list"res.desc$Dim.2## $quanti
## correlation p.value
## WV5 0.7821655 4.993508e-16
## Wet.Cover 0.4734967 2.660232e-05
## Wet 0.3121768 7.593693e-03
## Sulphate.ppm 0.3062099 8.897551e-03
## Iron -0.4246162 2.012171e-04
## DOC.mg.L -0.4653176 3.814274e-05
## Mt.ppt -0.6226775 5.218960e-09
##
## attr(,"class")
## [1] "condes" "list"eig.val <- get_eigenvalue(presencelanduse.pca1)
eig.val## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 2.2005786 16.927528 16.92753
## Dim.2 1.8683581 14.371985 31.29951
## Dim.3 1.5486474 11.912672 43.21218
## Dim.4 1.4679059 11.291584 54.50377
## Dim.5 1.2394159 9.533969 64.03774
## Dim.6 1.0366538 7.974260 72.01200
## Dim.7 0.9121358 7.016429 79.02843
## Dim.8 0.7394327 5.687944 84.71637
## Dim.9 0.7008610 5.391239 90.10761
## Dim.10 0.5868503 4.514233 94.62184
## Dim.11 0.3441403 2.647233 97.26908
## Dim.12 0.2232037 1.716952 98.98603
## Dim.13 0.1318163 1.013972 100.00000#Turning this into a biplot
landuse.presence1<-fviz_pca_biplot(presencelanduse.pca1,
# Individuals
geom.ind = "point",
fill.ind = mercurypca1$Region, col.ind = "black",
pointshape = 21, pointsize = 2,
palette = "jco",
addEllipses = TRUE,
# Variables
col.var = "contrib",
gradient.cols = "Dark2",
legend.title = list(fill = "Groups", color = "Contrib",
alpha = "Contrib"),xlab = "PC1: 17.68%", ylab = "PC2:16.66%")+theme(legend.position="right")
landuse.presence1
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.5-7synoptic.dist <- vegdist(mercurypca1[,4:16],method='bray')
adonis(synoptic.dist~mercurypca1[,3])##
## Call:
## adonis(formula = synoptic.dist ~ mercurypca1[, 3])
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## mercurypca1[, 3] 4 0.6019 0.150486 1.9218 0.10292 0.058 .
## Residuals 67 5.2465 0.078307 0.89708
## Total 71 5.8485 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Figure S1.
presencelanduse.pca <- PCA(mercurypca1[,4:16], graph = FALSE,scale=TRUE)
fviz_contrib(presencelanduse.pca, choice = "ind", axes = 1:2)
res.desc <- dimdesc(presencelanduse.pca, axes = c(1,2), proba = 0.05)
res.desc$Dim.1## $quanti
## correlation p.value
## Wet.Cover 0.6518443 5.582980e-10
## Temp.C 0.4747458 2.515740e-05
## DOC.mg.L 0.4721889 2.819748e-05
## WV5 0.4315947 1.535403e-04
## Wet -0.5934005 3.946600e-08
## Sulphate.ppm -0.7770315 1.026552e-15
##
## attr(,"class")
## [1] "condes" "list"res.desc$Dim.2## $quanti
## correlation p.value
## WV5 0.7821655 4.993508e-16
## Wet.Cover 0.4734967 2.660232e-05
## Wet 0.3121768 7.593693e-03
## Sulphate.ppm 0.3062099 8.897551e-03
## Iron -0.4246162 2.012171e-04
## DOC.mg.L -0.4653176 3.814274e-05
## Mt.ppt -0.6226775 5.218960e-09
##
## attr(,"class")
## [1] "condes" "list"eig.val <- get_eigenvalue(presencelanduse.pca)
eig.val## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 2.2005786 16.927528 16.92753
## Dim.2 1.8683581 14.371985 31.29951
## Dim.3 1.5486474 11.912672 43.21218
## Dim.4 1.4679059 11.291584 54.50377
## Dim.5 1.2394159 9.533969 64.03774
## Dim.6 1.0366538 7.974260 72.01200
## Dim.7 0.9121358 7.016429 79.02843
## Dim.8 0.7394327 5.687944 84.71637
## Dim.9 0.7008610 5.391239 90.10761
## Dim.10 0.5868503 4.514233 94.62184
## Dim.11 0.3441403 2.647233 97.26908
## Dim.12 0.2232037 1.716952 98.98603
## Dim.13 0.1318163 1.013972 100.00000#Turning this into a biplot
landuse.presence<-fviz_pca_biplot(presencelanduse.pca,
# Individuals
geom.ind = "point",
fill.ind = mercurypca1$Geology.1, col.ind = "black",
pointshape = 21, pointsize = 2,
palette = "jco",
addEllipses = TRUE,
# Variables
col.var = "contrib",
gradient.cols = "Dark2",
legend.title = list(fill = "Groups", color = "Contrib",
alpha = "Contrib"),xlab = "PC1: 16.93%", ylab = "PC2:14.37%")+theme(legend.position="right")
landuse.presence
synoptic.dist <- vegdist(mercurypca1[,4:16],method='bray')
adonis(synoptic.dist~mercurypca1[,2])##
## Call:
## adonis(formula = synoptic.dist ~ mercurypca1[, 2])
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## mercurypca1[, 2] 5 1.0924 0.218481 3.0319 0.18678 0.004 **
## Residuals 66 4.7561 0.072062 0.81322
## Total 71 5.8485 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Structural equation model - Figure 3.
setwd("E:/Postdoc- Stephanie Melles/Paper 1 sem submission")
library(lavaan)## This is lavaan 0.6-10
## lavaan is FREE software! Please report any bugs.##
## Attaching package: 'lavaan'## The following object is masked from 'package:psych':
##
## cor2covlibrary(semPlot)
library(OpenMx)##
## Attaching package: 'OpenMx'## The following object is masked from 'package:psych':
##
## trlibrary(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.1.6 v dplyr 1.0.8
## v tidyr 1.2.0 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1
## v purrr 0.3.4## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x ggplot2::%+%() masks psych::%+%()
## x ggplot2::alpha() masks psych::alpha()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(knitr)
library(kableExtra)##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2sem<-read.csv("figure3 sem.csv")
#Specifying the model based on the apriori circulated previously for feedback
modeltestsem <-'
#Latent variables
Watershed=~ Wet+WV51000+Geology+Yr_Dist.+X.wetcoverlog+Mixed
#Correlations
Wet ~~ Yr_Dist.
WV51000 ~~ Mixed
WV51000 ~~ Temp.C
Yr_Dist. ~~ Sulphate.ppm
WV51000 ~~ Sulphate.ppm
X.wetcoverlog ~~ Mixed
Yr_Dist. ~~ DOC.mg.L
#Measurement
Wet ~ Sulphate.ppm
DOC.mg.L ~ Temp.C
DOC.mg.L~Sulphate.ppm
Precip_wk~Phosphate10
Precip_wk~Nitrate.ppm
Wet ~ X.wetcoverlog
Geology ~ X.wetcoverlog
WV51000 ~ Sulphate.ppm
Iron ~ WV51000
WV51000 ~ X.wetcoverlog
DOC.mg.L ~ X.wetcoverlog
Temp.C ~ X.wetcoverlog
Temp.C ~ Phosphate10
Temp.C ~ Sulphate.ppm
Sulphate.ppm ~ X.wetcoverlog
Temp.C ~ Nitrate.ppm
Mt.ppt10~Watershed+Precip_wk+Iron+Temp.C+Sulphate.ppm+DOC.mg.L+Wet+Phosphate10+Nitrate.ppm+X.wetcoverlog
'
#Fitting the model to the dataframe
fittestsem <- cfa(modeltestsem, data = sem, missing='fiml',fixed.x=FALSE,std.lv=TRUE,estimator="MLR")
summary(fittestsem,fit.measures = TRUE, standardized=T,rsquare=T)## lavaan 0.6-10 ended normally after 191 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 68
##
## Number of observations 176
## Number of missing patterns 12
##
## Model Test User Model:
## Standard Robust
## Test Statistic 44.031 45.372
## Degrees of freedom 51 51
## P-value (Chi-square) 0.745 0.696
## Scaling correction factor 0.970
## Yuan-Bentler correction (Mplus variant)
##
## Model Test Baseline Model:
##
## Test statistic 508.051 518.746
## Degrees of freedom 90 90
## P-value 0.000 0.000
## Scaling correction factor 0.979
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000 1.000
## Tucker-Lewis Index (TLI) 1.029 1.023
##
## Robust Comparative Fit Index (CFI) 1.000
## Robust Tucker-Lewis Index (TLI) 1.023
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -311.672 -311.672
## Scaling correction factor 1.293
## for the MLR correction
## Loglikelihood unrestricted model (H1) -289.656 -289.656
## Scaling correction factor 1.155
## for the MLR correction
##
## Akaike (AIC) 759.343 759.343
## Bayesian (BIC) 974.936 974.936
## Sample-size adjusted Bayesian (BIC) 759.597 759.597
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000 0.000
## 90 Percent confidence interval - lower 0.000 0.000
## 90 Percent confidence interval - upper 0.036 0.039
## P-value RMSEA <= 0.05 0.992 0.988
##
## Robust RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.038
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.051 0.051
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Watershed =~
## Wet 0.132 0.019 6.930 0.000 0.132 0.790
## WV51000 -0.267 0.056 -4.786 0.000 -0.267 -0.312
## Geology 0.112 0.013 8.676 0.000 0.112 0.913
## Yr_Dist. -0.189 0.047 -4.007 0.000 -0.189 -0.325
## X.wetcoverlog -0.132 0.071 -1.841 0.066 -0.132 -0.389
## Mixed -0.040 0.018 -2.154 0.031 -0.040 -0.161
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Wet ~
## Sulphate.ppm 0.211 0.082 2.570 0.010 0.211 0.208
## DOC.mg.L ~
## Temp.C 0.309 0.159 1.944 0.052 0.309 0.225
## Sulphate.ppm -0.711 0.171 -4.146 0.000 -0.711 -0.471
## Precip_wk ~
## Phosphate10 0.031 0.025 1.247 0.213 0.031 0.097
## Nitrate.ppm 0.792 0.375 2.111 0.035 0.792 0.153
## Wet ~
## X.wetcoverlog 0.205 0.110 1.868 0.062 0.205 0.414
## Geology ~
## X.wetcoverlog 0.109 0.073 1.484 0.138 0.109 0.299
## WV51000 ~
## Sulphate.ppm 4.198 1.920 2.187 0.029 4.198 0.806
## Iron ~
## WV51000 -0.129 0.038 -3.373 0.001 -0.129 -0.232
## WV51000 ~
## X.wetcoverlog 1.294 0.344 3.760 0.000 1.294 0.511
## DOC.mg.L ~
## X.wetcoverlog -0.217 0.048 -4.568 0.000 -0.217 -0.296
## Temp.C ~
## X.wetcoverlog 0.132 0.034 3.895 0.000 0.132 0.247
## Phosphate10 0.081 0.012 6.537 0.000 0.081 0.380
## Sulphate.ppm -0.167 0.107 -1.561 0.119 -0.167 -0.152
## Sulphate.ppm ~
## X.wetcoverlog -0.117 0.043 -2.724 0.006 -0.117 -0.240
## Temp.C ~
## Nitrate.ppm -0.693 0.354 -1.961 0.050 -0.693 -0.203
## Mt.ppt10 ~
## Watershed 0.828 0.189 4.379 0.000 0.828 1.165
## Precip_wk -0.382 0.147 -2.591 0.010 -0.382 -0.147
## Iron 0.127 0.086 1.483 0.138 0.127 0.085
## Temp.C 0.921 0.292 3.149 0.002 0.921 0.234
## Sulphate.ppm 0.322 0.418 0.771 0.441 0.322 0.075
## DOC.mg.L 0.331 0.184 1.802 0.072 0.331 0.116
## Wet -4.176 0.916 -4.557 0.000 -4.176 -0.984
## Phosphate10 0.017 0.049 0.345 0.730 0.017 0.020
## Nitrate.ppm 3.794 0.978 3.879 0.000 3.794 0.283
## X.wetcoverlog 0.688 0.668 1.030 0.303 0.688 0.328
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Wet ~~
## .Yr_Dist. 0.023 0.005 4.568 0.000 0.023 0.377
## .WV51000 ~~
## .Mixed 0.044 0.013 3.339 0.001 0.044 0.194
## .Temp.C -0.026 0.008 -3.046 0.002 -0.026 -0.183
## .Yr_Dist. ~~
## .Sulphate.ppm -0.022 0.007 -3.375 0.001 -0.022 -0.253
## .WV51000 ~~
## .Sulphate.ppm -0.094 0.053 -1.755 0.079 -0.094 -0.629
## .X.wetcoverlog ~~
## .Mixed 0.017 0.006 2.613 0.009 0.017 0.222
## .Yr_Dist. ~~
## .DOC.mg.L -0.018 0.009 -1.979 0.048 -0.018 -0.155
## Phosphate10 ~~
## Nitrate.ppm -0.003 0.005 -0.697 0.486 -0.003 -0.076
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Wet 1.372 0.069 19.769 0.000 1.372 8.189
## .WV51000 2.140 0.602 3.553 0.000 2.140 2.496
## .Geology 0.588 0.039 15.144 0.000 0.588 4.788
## .Yr_Dist. 0.556 0.044 12.667 0.000 0.556 0.955
## .X.wetcoverlog 0.542 0.026 21.227 0.000 0.542 1.600
## .Mixed 1.255 0.019 67.742 0.000 1.255 5.106
## .DOC.mg.L 1.128 0.190 5.927 0.000 1.128 4.535
## .Precip_wk 0.482 0.039 12.346 0.000 0.482 1.762
## .Iron 2.941 0.139 21.180 0.000 2.941 6.179
## .Temp.C 1.058 0.046 23.019 0.000 1.058 5.853
## .Sulphate.ppm 0.308 0.030 10.350 0.000 0.308 1.871
## .Mt.ppt10 5.032 1.272 3.957 0.000 5.032 7.082
## Phosphate10 0.776 0.064 12.202 0.000 0.776 0.920
## Nitrate.ppm 0.064 0.008 7.834 0.000 0.064 1.202
## Watershed 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Wet 0.012 0.003 4.348 0.000 0.012 0.426
## .WV51000 0.870 0.410 2.121 0.034 0.870 1.184
## .Geology 0.004 0.001 4.742 0.000 0.004 0.289
## .Yr_Dist. 0.303 0.024 12.426 0.000 0.303 0.894
## .X.wetcoverlog 0.097 0.020 4.975 0.000 0.097 0.849
## .Mixed 0.059 0.006 10.255 0.000 0.059 0.974
## .DOC.mg.L 0.043 0.009 4.886 0.000 0.043 0.700
## .Precip_wk 0.073 0.007 10.428 0.000 0.073 0.969
## .Iron 0.214 0.068 3.135 0.002 0.214 0.946
## .Temp.C 0.023 0.003 6.678 0.000 0.023 0.701
## .Sulphate.ppm 0.026 0.004 6.315 0.000 0.026 0.942
## .Mt.ppt10 0.124 0.057 2.196 0.028 0.124 0.246
## Watershed 1.000 1.000 1.000
## Phosphate10 0.712 0.045 15.683 0.000 0.712 1.000
## Nitrate.ppm 0.003 0.001 1.988 0.047 0.003 1.000
##
## R-Square:
## Estimate
## Wet 0.574
## WV51000 -0.184
## Geology 0.711
## Yr_Dist. 0.106
## X.wetcoverlog 0.151
## Mixed 0.026
## DOC.mg.L 0.300
## Precip_wk 0.031
## Iron 0.054
## Temp.C 0.299
## Sulphate.ppm 0.058
## Mt.ppt10 0.754