Human disturbance as an indicator of ecological integrity:
statewide assessment of Florida’s wetlands and coastal waters
Publication:
(Pending, contact corresponding author to cite)
Data Accessibilty
Data and code to execute below models can be downloaded from here: (Pending)
date()## [1] "Sun Jun 25 17:19:44 2017"Options
library(knitr)
opts_knit$set(verbose = FALSE)Load Packages:
suppressMessages(library(INLA))
suppressMessages(library(rgdal))
suppressMessages(library(reshape2))
suppressMessages(library(raster))
suppressMessages(library(rasterVis))
suppressMessages(library(rgeos))
suppressMessages(library(maptools))
suppressMessages(library(mapproj))
suppressMessages(library(xtable))
suppressMessages(library(viridis))
suppressMessages(library(corrplot))
suppressMessages(library(ggplot2))
suppressMessages(library(GISTools))
suppressMessages(library(lattice))
suppressMessages(library(Thermimage))
suppressMessages(library(gridExtra))
suppressMessages(library(dplyr))
suppressMessages(library(spatstat))
suppressMessages(library(classInt))Set Directory
setwd("D:/Dist031017")Get AEI Values
AEIs.r = raster("./nLDI.tif")Get State Boundary
nProj = "+proj=aea +lat_1=24 +lat_2=31.5 +lat_0=24 +lon_0=-84
+x_0=400000 +y_0=0 +ellps=GRS80 +units=km +no_defs"
LL84 = "+proj=longlat +datum=WGS84
+no_defs +ellps=WGS84 +towgs84=0,0,0"
Counties = readOGR(dsn = "./County",
layer = "FL_County",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "./County", layer: "FL_County"
## with 67 features
## It has 46 fields
## Integer64 fields read as strings: NO_FARMS07 AVG_SIZE07 CROP_ACR07Florida = gUnaryUnion(Counties)
Florida = spTransform(Florida, nProj)
#Florida.rp = rasterize(aggregate(Florida,
# fact = 4),
# AEIs.r,
# field = 0,
# background = NA)
#writeRaster(Florida.rp, "./Florida.rp.tif")
Florida.rp = raster("./Florida.rp.tif")
Florida.buf = readOGR(dsn = "./ArcMap/CountyB",
layer = "Florida15K",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "./ArcMap/CountyB", layer: "Florida15K"
## with 1 features
## It has 1 fieldsFlorida.buf = spTransform(Florida.buf, CRS(proj4string(Florida.rp)))
#Florida.bufr = aggregate(Florida.rp, fact = 4)
#Florida.bufr = extend(Florida.bufr, c(100,100))
#Florida.bufr = rasterize(Florida.buf,
# Florida.bufr,
# field = 0,
# background = NA)
#writeRaster(Florida.bufr, "./Florida.bufr.tif")
Florida.bufr = raster("./Florida.bufr.tif")Plot LDI
Domain.r = aggregate(AEIs.r, fact = 12, method = "ngb")
Domain2.r = projectRaster(Domain.r, crs = LL84, method = "ngb")
Domain3.r = disaggregate(Domain2.r, fact=4)
CountiesLL = spTransform(Counties, CRS(LL84))
FloridaLL = spTransform(Florida, CRS(LL84))
rng = seq(0, 41.1134, 0.01)
cr = colorRampPalette(c("darkblue", "blue", "lightblue", "cyan2",
"yellow", "orange", "red", "darkred"),
bias = 1.5, space = "rgb")
Cp = levelplot(Domain3.r,
margin = FALSE,
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 41.1134, 0.01),
labels=list(at=c(0, 10, 20, 30, 40),
labels=c("0", "10", "20", "30", "40"))),
panel = panel.levelplot, space = "right", par.strip.text = list(fontface='bold'),
col = cr, par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(col = "black", cex = 2))
Cp +
latticeExtra::layer(sp.polygons(CountiesLL, col = "black", lwd = 2)) +
latticeExtra::layer(sp.polygons(FloridaLL, col = "darkgray", lwd = 2))
Load 2011 Wetland Condition Assessment (WCA) Site Locations
Site information for EPA Wetland Condition Assessment. Data available here: https://www.epa.gov/national-aquatic-resource-surveys/nwca
Sites = read.csv("./nwca2011_siteinfo.csv")
Sites = Sites %>%
filter(STATE == "FL")
dd = cbind(Sites$AA_CENTER_LON, Sites$AA_CENTER_LAT)
names(dd) = c("Long", "Lat")
Site.pnts = SpatialPointsDataFrame(dd, Sites)
proj4string(Site.pnts) = proj4string(FloridaLL)
Site.pnts = spTransform(Site.pnts, nProj)
Site.pnts$Region = 1:nrow(Site.pnts)
length(Site.pnts$Region) #Number of Sites## [1] 82Load 2010 Coastal Assessment Site Locations
Site information for EPA Coastal Condition Assessment. Data available here: https://www.epa.gov/national-aquatic-resource-surveys/nwca
CSites = read.csv("D:/Dist031017/WCA_raw/Coastal2010/assessed_ncca2010_siteinfo.revised.06212016.csv")
CSites = CSites %>%
filter(STATE == "FL")
dd = cbind(CSites$ALON_DD, CSites$ALAT_DD)
names(dd) = c("Long", "Lat")
CSite.pnts = SpatialPointsDataFrame(dd, CSites)
proj4string(CSite.pnts) = proj4string(FloridaLL)
CSite.pnts = spTransform(CSite.pnts, nProj)
CSite.pnts$Region = 1:nrow(CSite.pnts)
length(CSite.pnts$Region) #Number of Sites## [1] 94Plot Field Validation Sites
Circles = Locations of WCA Field Sites
Triangles = Locations of Coastal Field Sites
Site.pntsLL = spTransform(Site.pnts, CRS(proj4string(CountiesLL)))
CSite.pntsLL = spTransform(CSite.pnts, CRS(proj4string(CountiesLL)))
#Florida.r = projectRaster(aggregate(Florida.rp, fact = 12), crs = LL84, method = "ngb")
#writeRaster(Florida.r, "./Florida.r.tif")
Florida.r = raster("./Florida.r.tif")
rng = seq(0, 10, 0.01)
cr = colorRampPalette(c("tan", "blue", "lightblue",
"yellow", "orange", "red", "darkred"),
bias = 1.5, space = "rgb")
Cp = levelplot(Florida.r,
margin = FALSE,
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = NULL,
panel = panel.levelplot, space = "right", par.strip.text = list(fontface='bold'),
col = cr, par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(col = "black", cex = 2))
Cp +
latticeExtra::layer(sp.polygons(CountiesLL, col = "darkgray", lwd = 1.5)) +
latticeExtra::layer(sp.polygons(Site.pntsLL, col = "red", cex = 1, pch=19)) +
latticeExtra::layer(sp.polygons(CSite.pntsLL, col = "red", cex = 1, pch=17))
Define WCA and Coastal Site Locations
Site.5k = rbind(Site.pnts@coords, CSite.pnts@coords)
Site.5k = SpatialPointsDataFrame(Site.5k, as.data.frame(Site.5k))
proj4string(Site.5k) = proj4string(CSite.pnts)
Site.5k$Region = 1:nrow(Site.5k)
grd.ppp = as.ppp(Site.5k)
Site.5k$NN = nndist(grd.ppp)
Site.5k = subset(Site.5k, NN > 4.25)
Site.5k = gEnvelope(gBuffer(Site.5k,
width=4.25,
byid=TRUE,
id = Site.5k$Region),
byid=TRUE,
id = Site.5k$Region)Point Grid
A dense point grid is created for each focal area.
for (i in 1:length(Site.5k)){
Tmp.poly = Site.5k[i,]
Tmp.grd = makegrid(Tmp.poly, cellsize = 0.15)
if( i == 1){Grid.df = Tmp.grd}
else{Grid.df = rbind(Grid.df, Tmp.grd)}
}
Grd.pnts = SpatialPointsDataFrame(Grid.df, Grid.df)
proj4string(Grd.pnts) = proj4string(Site.5k)
Grd.pnts@data = Grd.pnts@data %>%
mutate(ID = 1:nrow(Grd.pnts),
Long = Grd.pnts@coords[,1],
Lat = Grd.pnts@coords[,2])
Grd.pnts$M_Out = extract(Florida.rp, Grd.pnts, method = "simple")
Grd.pnts = subset(Grd.pnts, M_Out == 0 ) Construct Mesh
Constructing a “nested” mesh.
Grd.coord = as.data.frame(cbind(Grd.pnts$Long, Grd.pnts$Lat))
FloridaB = gConvexHull(Florida)
bdry = inla.sp2segment(FloridaB)
mesh0 = inla.mesh.2d(loc = Grd.coord,
boundary = bdry,
min.angle = 22,
max.edge = 5)
mesh0b = inla.mesh.2d(loc = mesh0$loc,
min.angle = 22,
max.edge = 100)
mesh0b$n## [1] 333712mesh.dom = mesh0bMesh Verticies to Spatial Points
dd2 = as.data.frame(cbind(mesh.dom$loc[,1],
mesh.dom$loc[,2]))
names(dd2) = c("Long", "Lat")
V.pnts = SpatialPointsDataFrame(dd2, dd2)
proj4string(V.pnts) = proj4string(Florida)
V.pnts = spTransform(V.pnts, CRS(proj4string(Florida)))Get LDI Values for Model Fitting
Fl.pnts = V.pnts
Fl.pnts$LDIs = extract(AEIs.r, Fl.pnts, method = "simple")
Fl.pnts$M_Out = extract(Florida.bufr, Fl.pnts, method = "simple")
Fl.pnts = subset(Fl.pnts, M_Out == 0 )
Fl.pnts$LDIs = ifelse(is.na(Fl.pnts$LDIs) == TRUE,
0,
Fl.pnts$LDIs)
range(Fl.pnts$LDIs)## [1] 0.00 41.74Nearest Road
Calculate distance to nearest road.
#FLRoads = readOGR(dsn = "./Roads_fl",
# layer = "Roads_fl",
# stringsAsFactors = FALSE)
#FLRoads = spTransform(FLRoads, CRS(proj4string(Fl.pnts)))
#Rd.psp = as.psp(FLRoads)
load("D:/Dist031017/Roads_fl/Rd_psp.RData") #Previously converted from SLDF to psp (time intensive)
P.pp = as.ppp(Fl.pnts)
Fl.pnts$nRd = nncross(P.pp, Rd.psp)[,"dist"]
Fl.pnts$nRd = scale(Fl.pnts$nRd, scale=TRUE, center=TRUE)
range(Fl.pnts$nRd)## [1] -0.7031401 28.0319948Population Density
Population Density
http://sedac.ciesin.columbia.edu/data/set/gpw-v3-population-density
For near-shore areas, assigning the nearest terrestrial population density estimate.
PopDns = raster("./PopDns/PopDns.tif")
dist = distance(PopDns)
direct = direction(PopDns, from=FALSE)
rna = is.na(PopDns)
na.x = init(rna, 'x')
na.y = init(rna, 'y')
co.x = na.x + dist * sin(direct)
co.y = na.y + dist * cos(direct)
co = cbind(co.x[], co.y[])
NAVals = raster::extract(PopDns, co, method='simple')
r.NAVals = rna
r.NAVals[] = NAVals
r.filled = cover(x=PopDns, y= r.NAVals)
Fl.pnts$Pd = extract(r.filled,
spTransform(Fl.pnts, CRS(proj4string(r.filled))),
method = "simple")
Fl.pnts$Pd = ifelse(is.na(Fl.pnts$Pd) == TRUE,
mean(Fl.pnts$Pd, na.rm=T), Fl.pnts$Pd)
Fl.pnts$Pd = scale(Fl.pnts$Pd, scale=TRUE, center=TRUE)
range(Fl.pnts$Pd)## [1] -0.4175205 9.6383535Well Data
Calculate distance to nearest water well.
Well data (Florida DEP): http://geodata.dep.state.fl.us/datasets/c3e9b9a203e24bb78c2567fa5b156600_4
Wells = read.csv("./Wells/Well_FL.csv")
levels(Wells$WELL_TYPE)## [1] ""
## [2] "AGRICULTURAL SUPPLY WELL"
## [3] "DRAINAGE WELL"
## [4] "GROUND WATER LEVEL OBSERVATION WELL"
## [5] "GROUND WATER OBSERVATION WELL"
## [6] "GROUND WATER QUALITY MONITORING WELL"
## [7] "GROUND WATER QUALITY OBSERVATION WELL"
## [8] "INDUSTRIAL SUPPLY WELL"
## [9] "IRRIGATION WELL"
## [10] "NOT YET DETERMINED"
## [11] "OTHER WELL"
## [12] "PRIVATE DRINKING WATER WELL"
## [13] "PUBLIC DRINKING WATER WELL"Wells = Wells %>%
filter(WELL_TYPE == "PUBLIC DRINKING WATER WELL" |
WELL_TYPE == "AGRICULTURAL SUPPLY WELL" |
WELL_TYPE == "INDUSTRIAL SUPPLY WELL" |
WELL_TYPE == "IRRIGATION WELL" |
WELL_TYPE == "PRIVATE DRINKING WATER WELL") %>%
select(X, Y, PK_STATION)
dd = as.data.frame(cbind(Wells$X,
Wells$Y))
Wells = SpatialPointsDataFrame(dd, Wells)
proj4string(Wells) = proj4string(Florida.r)
Wells = spTransform(Wells, CRS(proj4string(Florida)))
W.pp = as.ppp(Wells)
Fl.pnts$nW = nncross(P.pp, W.pp)[,"dist"]
Fl.pnts$nW = scale(Fl.pnts$nW, scale=TRUE, center=TRUE)
range(Fl.pnts$nW)## [1] -0.6370749 14.2005881Fl.pnts@data = as.data.frame(apply(Fl.pnts@data, 2, as.numeric))Plot Well Locations
WellsLL = spTransform(Wells, CRS(proj4string(Florida.r)))
rng = seq(0, 10, 0.01)
cr = colorRampPalette(c("tan", "blue", "lightblue",
"yellow", "orange", "red", "darkred"),
bias = 1.5, space = "rgb")
Cp = levelplot(Florida.r,
margin = FALSE,
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = NULL,
panel = panel.levelplot, space = "right", par.strip.text = list(fontface='bold'),
col = cr, par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(col = "black", cex = 2))
Cp +
latticeExtra::layer(sp.polygons(CountiesLL, col = "darkgray", lwd = 1.5)) +
latticeExtra::layer(sp.polygons(WellsLL, col = "red", cex = 0.1, pch=19))
Project Points to Mesh
locs = cbind(Fl.pnts$Long, Fl.pnts$Lat)
A = inla.spde.make.A(mesh.dom, loc=locs)
spde = inla.spde2.pcmatern(mesh.dom,
prior.range=c(5, 0.01),
prior.sigma=c(1, 0.01))
field0 = inla.spde.make.index("field0",
spde$n.spde)Data Stack
FE.df = Fl.pnts@data
FE.lst = list(c(field0,
list(intercept0 = 1)),
list(Pd = FE.df[,"Pd"],
nRd = FE.df[,"nRd"],
nW = FE.df[,"nW"]))
Stack0 = inla.stack(data = list(LDI = FE.df$LDIs),
A = list(A, 1),
effects = FE.lst,
tag = "Obs0")Multi-Colinearity
CorrCov = FE.df %>% select(Pd, nRd, nW)
CorCov.cv = cor(CorrCov)
suppressMessages(library(perturb))
CI = colldiag(CorrCov)
CI## Condition
## Index Variance Decomposition Proportions
## intercept Pd nRd nW
## 1 1.000 0.000 0.087 0.214 0.196
## 2 1.255 1.000 0.000 0.000 0.000
## 3 1.309 0.000 0.851 0.015 0.132
## 4 1.766 0.000 0.062 0.771 0.672detach("package:perturb", unload=TRUE)Model1
Frm1 = LDI ~ -1 + intercept0 +
f(field0, model=spde)
#theta1 = Model1$internal.summary.hyperpar$mean
theta1 = c(-3.3879059, 0.5474899, 2.5184320)
#darwini
Model1 = inla(Frm1,
data = inla.stack.data(Stack0, spde=spde),
family = "gaussian",
verbose = TRUE,
control.predictor = list(
A = inla.stack.A(Stack0),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta1),
control.inla = list(int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) Model1 Metrics
summary(Model1)##
## Call:
## c("inla(formula = Frm1, family = \"gaussian\", data = inla.stack.data(Stack0, ", " spde = spde), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack0), ", " compute = TRUE, link = 1), control.inla = list(int.strategy = \"eb\"), ", " control.results = list(return.marginals.random = TRUE, return.marginals.predictor = TRUE), ", " control.mode = list(restart = TRUE, theta = theta1))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 4.4444 6218.0231 101.7837 6324.2511
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept0 5.7131 0.0725 5.5707 5.7131 5.8553 5.7131 0
##
## Random effects:
## Name Model
## field0 SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.0338 0.0002 0.0335 0.0338
## Range for field0 1.7289 0.0093 1.7108 1.7289
## Stdev for field0 12.4092 0.0518 12.3073 12.4093
## 0.975quant mode
## Precision for the Gaussian observations 0.0341 0.0338
## Range for field0 1.7473 1.7289
## Stdev for field0 12.5111 12.4098
##
## Expected number of effective parameters(std dev): 106202.17(0.00)
## Number of equivalent replicates : 2.99
##
## Deviance Information Criterion (DIC) ...: 2082862.93
## Effective number of parameters .........: 106202.28
##
## Watanabe-Akaike information criterion (WAIC) ...: 2082027.21
## Effective number of parameters .................: 85043.30
##
## Marginal log-Likelihood: -1085677.06
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModels = "Model1"
Fixed = "W (Coarse spatial structure only)"
#Deviance information criteria
DICs = Model1$dic$dic
#Watanabe-Akaike information criteria
WAICs = Model1$waic$waic
Model1_mets = as.data.frame(cbind(Model = Models,
DIC = round(DICs, 2),
WAIC = round(WAICs, 2),
COVARIATES = Fixed))
rm("Model1")Model2
Frm2 = LDI ~ -1 + intercept0 +
f(field0, model=spde) +
Pd + nRd + nW
#theta2 = Model2$internal.summary.hyperpar$mean
theta2 = c(-3.4167965, 0.3898347, 2.3316264)
Model2 = inla(Frm2,
data = inla.stack.data(Stack0, spde=spde),
family = "gaussian",
verbose = TRUE,
control.predictor = list(
A = inla.stack.A(Stack0),
compute = TRUE,
link = 1),
control.mode = list(restart = TRUE, theta = theta2),
control.inla = list(int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = FALSE, waic = TRUE)) Model2 Metrics
summary(Model2)##
## Call:
## c("inla(formula = Frm2, family = \"gaussian\", data = inla.stack.data(Stack0, ", " spde = spde), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = FALSE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack0), ", " compute = TRUE, link = 1), control.inla = list(int.strategy = \"eb\"), ", " control.results = list(return.marginals.random = TRUE, return.marginals.predictor = TRUE), ", " control.mode = list(restart = TRUE, theta = theta2))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 4.3970 3636.7933 114.2837 3755.4741
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept0 6.6138 0.0683 6.4797 6.6138 6.7477 6.6138 0
## Pd 2.7998 0.0817 2.6393 2.7998 2.9602 2.7998 0
## nRd -0.9002 0.0436 -0.9859 -0.9002 -0.8146 -0.9002 0
## nW -0.3309 0.0508 -0.4307 -0.3309 -0.2313 -0.3309 0
##
## Random effects:
## Name Model
## field0 SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.0341 0.0002 0.0338 0.0341
## Range for field0 1.6000 0.0090 1.5827 1.5998
## Stdev for field0 11.7903 0.0531 11.6916 11.7881
## 0.975quant mode
## Precision for the Gaussian observations 0.0344 0.0341
## Range for field0 1.6183 1.5991
## Stdev for field0 11.9001 11.7801
##
## Expected number of effective parameters(std dev): 107454.86(0.00)
## Number of equivalent replicates : 2.955
##
## Deviance Information Criterion (DIC) ...: 2081390.36
## Effective number of parameters .........: 107454.97
##
## Watanabe-Akaike information criterion (WAIC) ...: 2081184.15
## Effective number of parameters .................: 86425.98
##
## Marginal log-Likelihood: -1084560.44
## Posterior marginals for linear predictor and fitted values computedModels = "Model2"
Fixed = "Pop + nRd + nW + W"
#Deviance information criteria
DICs = Model2$dic$dic
#Watanabe-Akaike information criteria
WAICs = Model2$waic$waic
Model2_mets = as.data.frame(cbind(Model = Models,
DIC = round(DICs, 2),
WAIC = round(WAICs, 2),
COVARIATES = Fixed))Model3 (Non-Spatial Model)
Frm3 = LDI ~ -1 + intercept0 +
Pd + nRd + nW
Model3 = inla(Frm3,
data = inla.stack.data(Stack0, spde=spde),
family = "gaussian",
verbose = TRUE,
control.predictor = list(
A = inla.stack.A(Stack0),
compute = TRUE,
link = 1),
control.inla = list(int.strategy = "eb"),
control.results = list(return.marginals.random = TRUE,
return.marginals.predictor = TRUE),
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) Model3 Metrics
summary(Model3)##
## Call:
## c("inla(formula = Frm3, family = \"gaussian\", data = inla.stack.data(Stack0, ", " spde = spde), verbose = TRUE, control.compute = list(dic = TRUE, ", " cpo = TRUE, waic = TRUE), control.predictor = list(A = inla.stack.A(Stack0), ", " compute = TRUE, link = 1), control.inla = list(int.strategy = \"eb\"), ", " control.results = list(return.marginals.random = TRUE, return.marginals.predictor = TRUE))" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.3751 561.9245 101.7705 666.0702
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## intercept0 8.4482 0.0203 8.4083 8.4482 8.4881 8.4482 0
## Pd 5.1039 0.0209 5.0628 5.1039 5.1449 5.1039 0
## nRd -3.0575 0.0211 -3.0989 -3.0575 -3.0162 -3.0575 0
## nW -0.1567 0.0207 -0.1972 -0.1567 -0.1162 -0.1567 0
##
## The model has no random effects
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.0076 0.00 0.0076 0.0076
## 0.975quant mode
## Precision for the Gaussian observations 0.0077 0.0076
##
## Expected number of effective parameters(std dev): 4.03(0.00)
## Number of equivalent replicates : 78776.19
##
## Deviance Information Criterion (DIC) ...: 2449456.22
## Effective number of parameters .........: 4.03
##
## Watanabe-Akaike information criterion (WAIC) ...: 2449466.71
## Effective number of parameters .................: 14.49
##
## Marginal log-Likelihood: -1224773.17
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModels = "Model3"
Fixed = "Pop + nRd + nW"
#Deviance information criteria
DICs = Model3$dic$dic
#Watanabe-Akaike information criteria
WAICs = Model3$waic$waic
Model3_mets = as.data.frame(cbind(Model = Models,
DIC = round(DICs, 2),
WAIC = round(WAICs, 2),
COVARIATES = Fixed))
rm("Model3")Model Results
Model Comparison
Model_mets = as.data.frame(rbind(Model1_mets, Model2_mets, Model3_mets))
print.data.frame(Model_mets, right = FALSE)## Model DIC WAIC COVARIATES
## 1 Model1 2082862.93 2082027.21 W (Coarse spatial structure only)
## 2 Model2 2081390.36 2081184.15 Pop + nRd + nW + W
## 3 Model3 2449456.22 2449466.71 Pop + nRd + nW#ModComp = xtable(Model_mets)
#print(ModComp, include.rownames = TRUE)Selected Model (Fixed Effects)
SelMod.tab = as.data.frame(Model2$summary.fixed)
SelMod.tab2 = cbind(SelMod.tab[,1:3], SelMod.tab[,5])
names(SelMod.tab2) = c("Mean", "sd", "2.5% Q", "97.5% Q")
print.data.frame(SelMod.tab2, right = FALSE, digits = 3)## Mean sd 2.5% Q 97.5% Q
## intercept0 6.614 0.0683 6.480 6.748
## Pd 2.800 0.0817 2.639 2.960
## nRd -0.900 0.0436 -0.986 -0.815
## nW -0.331 0.0508 -0.431 -0.231#SelMod.tab2 = xtable(SelMod.tab2, digits = 4)
#print(SelMod.tab2, digits = 4, include.rownames = TRUE)Spatial Random Effect
Point estimates and 95% credible intervals for the course random spatial effect.
mRF.df = as.data.frame(Model2$summary.random$field0)[,1:6]
names(mRF.df) = c("ID", "Mean", "sd", "Q025", "Q50", "Q975")
mRF.df = arrange(mRF.df, Mean)
mRF.df$Index = 1:nrow(mRF.df)
RedSample = function(x){x[seq(1, length(x), 10)]} #Thinning for plotting
mRF.df2 = as.data.frame(apply(mRF.df, 2, RedSample))
SpCor2 = ggplot(mRF.df2, aes(Index, Mean)) +
#coord_flip() +
geom_smooth(method = "loess",
col = "black",
span = 0.25,
se = FALSE) +
geom_line(data = mRF.df2,
aes(Index, Q025),
col = "grey") +
geom_line(data = mRF.df2,
aes(Index, Q975),
col = "grey") +
geom_hline(yintercept = 0,
linetype = "dotted",
col = "red",
size = 1) +
xlab("Spatial Index") +
ylab("Spatial Field") +
theme_classic() +
theme(axis.text=element_text(size=14),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
plot.title = element_text(hjust = 0.5))
SpCor2
Density of the Marginal Terms
(Fixed Effects)
Int.df = as.data.frame(Model2$marginals.fix[[1]])
Pop.df = as.data.frame(Model2$marginals.fix[[2]])
nRd.df = as.data.frame(Model2$marginals.fix[[3]])
nW.df = as.data.frame(Model2$marginals.fix[[4]])
CI.df = as.data.frame(Model2$summary.fixed)[,c(3,5)]
In.plt = ggplot(Int.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[1,1], col = "gray") +
geom_vline(xintercept = CI.df[1,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
ggtitle("A") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
Pop.plt = ggplot(Pop.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[2,1], col = "gray") +
geom_vline(xintercept = CI.df[2,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab(NULL) +
ggtitle("B") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
nrd.plt = ggplot(nRd.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[3,1], col = "gray") +
geom_vline(xintercept = CI.df[3,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab("Density") +
ggtitle("C") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
nW.plt = ggplot(nW.df , aes(x, y)) +
geom_line(size = 1, col = "black") +
geom_vline(xintercept = CI.df[4,1], col = "gray") +
geom_vline(xintercept = CI.df[4,2], col = "gray") +
geom_vline(xintercept = 0, col = "red") +
xlab(NULL) +
ylab(NULL) +
ggtitle("D") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2),
plot.title = element_text(size = 20, hjust = 0.5))
grid.arrange(In.plt, Pop.plt, nrd.plt, nW.plt, ncol = 2)
Prediction
Create Raster and Extract Point Estimates
Create point grid for prediction locations.
Domain.b = aggregate(Florida.bufr, fact = 3)
Grd.pnts = rasterToPoints(Domain.b, sp = TRUE)
Grd.pnts@data = Grd.pnts@data %>%
mutate(Long = Grd.pnts@coords[,1],
Lat = Grd.pnts@coords[,2])Sum of Linear Components
Calculate covariates and then find predicted value for each location of point grid.
#Pop Density
Grd.pnts$Pd = extract(r.filled, spTransform(Grd.pnts, CRS(proj4string(PopDns))), method = "simple")
Grd.pnts$Pd = ifelse(is.na(Grd.pnts$Pd) == TRUE, mean(Grd.pnts$Pd, na.rm=T), Grd.pnts$Pd)
Grd.pnts$Pd = scale(Grd.pnts$Pd, scale=TRUE, center=TRUE)
range(Grd.pnts$Pd)## [1] -0.3919112 13.3856169#Nearest Road
Grd.pp = as.ppp(Grd.pnts)
Grd.pnts$nRd = nncross(Grd.pp, Rd.psp)[,"dist"]
Grd.pnts$nRd = scale(Grd.pnts$nRd, scale=TRUE, center=TRUE)
range(Grd.pnts$nRd)## [1] -0.5370561 12.1547942#Nearest Well
Grd.pnts$nW = nncross(Grd.pp, W.pp)[,"dist"]
Grd.pnts$nW = scale(Grd.pnts$nW, scale=TRUE, center=TRUE)
range(Grd.pnts$nW)## [1] -0.5086363 7.2225919#GRF
pLoc = cbind(Grd.pnts$Long, Grd.pnts$Lat)
Ap = inla.spde.make.A(mesh.dom, loc=pLoc)
Grd.pnts$m2RE = drop(Ap %*% Model2$summary.random$field0$mean)
Fixed.m2 = as.data.frame(Model2$summary.fixed)[,"mean"]
m2Int = Fixed.m2[1]
m2Pdm = Fixed.m2[2]
m2nRdm = Fixed.m2[3]
m2nWm = Fixed.m2[4]
Grd.pnts@data = as.data.frame(apply(Grd.pnts@data, 2, as.numeric))
Grd.pnts@data = Grd.pnts@data %>%
mutate(M2Pred = m2RE + m2Int + m2Pdm*Pd + m2nRdm*nRd + m2nWm*nW)Visualize State Prediction
Including natural and anthropogenic land use/land cover.
Pred_st = rasterize(Grd.pnts,
Domain.b,
"M2Pred",
background = NA)
Pred_st2 = Pred_st
Pred_st2 = projectRaster(Pred_st,
crs = proj4string(Domain2.r),
method = "ngb")
Pred_st2[Pred_st2<0]=0
Florida.buf2 = spTransform(Florida.buf, proj4string(Pred_st2))
rng = seq(0, 40, 0.01)
cr = colorRampPalette(c("darkblue", "blue", "cyan2",
"yellow", "orange", "red", "darkred"),
bias = 2, space = "rgb")
M1RE = levelplot(Pred_st2,
margin = FALSE,
#sub = "NA",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 40, 0.01),
labels=list(at=c(0, 10, 20, 30, 40)),
labels=c("0", "10", "20", "30", "40")),
panel = panel.levelplot, space = "right", par.strip.text = list(fontface='bold'),
col = cr, par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(col = "black", cex = 2))
M1RE +
latticeExtra::layer(sp.polygons(CountiesLL, col = "black", lwd = 2)) +
latticeExtra::layer(sp.polygons(Florida.buf2, col = "darkgray", lwd = 2)) 
Visualize Predictions for Natural Areas Only
Tmp.pred = disaggregate(Pred_st, fact = 4,
method = "bilinear")
Tmp.ldi = AEIs.r
Tmp.ldi = extend(Tmp.ldi, Tmp.pred)
Tmp.ldi = resample(Tmp.ldi, Tmp.pred, method = "ngb")
NAreas.r = Tmp.ldi
NAreas.r[is.na(NAreas.r)] = 0
values(NAreas.r) = ifelse(values(NAreas.r) == 0,
values(Tmp.pred), NA)
NAreas.r[NAreas.r < 0] = 0
NAreas.r = projectRaster(NAreas.r,
crs = proj4string(Domain2.r),
method = "ngb")Plot Prediction for Natural Areas
rng = seq(0, 40, 0.01)
M1RE = levelplot(NAreas.r,
margin = FALSE,
#sub = "NA",
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(at=seq(0, 40, 0.01),
labels=list(at=c(0, 10, 20, 30, 40)),
labels=c("0", "10", "20", "30", "40")),
panel = panel.levelplot, space = "right", par.strip.text = list(fontface='bold'),
col = cr, par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(col = "black", cex = 2))
M1RE +
latticeExtra::layer(sp.polygons(CountiesLL, col = "black", lwd = 2)) +
latticeExtra::layer(sp.polygons(Florida.buf2, col = "darkgray", lwd = 2)) 
Model Validation
Compare to Field Collected Samples
Model predictions are compared to field data collected in conjunction with the 2011 National Wetland Condition Assessment (WCA) and the 2010 Coastal Assessment. https://www.epa.gov/national-aquatic-resource-surveys/nwca
WCA Assessment Sites
Calculate covariates and then find predicted value for each WCA field site location.
#Pop Density
Site.pnts$Pd = extract(r.filled, spTransform(Site.pnts, CRS(proj4string(r.filled))), method = "simple")
Site.pnts$Pd = ifelse(is.na(Site.pnts$Pd) == TRUE, mean(Site.pnts$Pd, na.rm=T), Site.pnts$Pd)
Site.pnts$Pd = as.numeric(scale(Site.pnts$Pd, scale=TRUE, center=TRUE))
range(Site.pnts$Pd)## [1] -0.5128392 3.6302929#Nearest Road
Grd.pp = as.ppp(Site.pnts)
Site.pnts$nRd = nncross(Grd.pp, Rd.psp)[,"dist"]
Site.pnts$nRd = as.numeric(scale(Site.pnts$nRd, scale=TRUE, center=TRUE))
range(Site.pnts$nRd)## [1] -0.8872327 3.5463884#Nearest Well
Site.pnts$nW = nncross(Grd.pp, W.pp)[,"dist"]
Site.pnts$nW = as.numeric(scale(Site.pnts$nW, scale=TRUE, center=TRUE))
range(Site.pnts$nW)## [1] -0.4596036 5.9377740#GRF
Ap = inla.spde.make.A(mesh.dom, loc=Site.pnts)
Site.pnts$m2RE = drop(Ap %*% Model2$summary.random$field0$mean)
Fixed.m2 = as.data.frame(Model2$summary.fixed)[,"mean"]
m2Int = Fixed.m2[1]
m2Pdm = Fixed.m2[2]
m2nRdm = Fixed.m2[3]
m2nWm = Fixed.m2[4]
Site.pnts@data = Site.pnts@data %>%
mutate(M2Pred = m2RE + m2Int + m2Pdm*Pd + m2nRdm*nRd + m2nWm*nW)Vegetation Assessment
Match field results by site identifier (SITE_ID).
VegStress = read.csv("./WCA_raw/nwca2011_vegmmi.csv")
VegStress = VegStress[match(Site.pnts$SITE_ID, VegStress$SITE_ID, nomatch=0),]
VegStress = merge(VegStress, Site.pnts@data, by = "SITE_ID", nomatch=0)Select Indicators (Vegetation)
Acronyms:
VMMI = Vegetation Multimetric Index (MMI) score
N_TOL = Species Richness
Exp.VegStress = VegStress %>%
select(SITE_ID, M2Pred, VMMI, N_TOL)
Exp.VegStress1 = Exp.VegStress[,2:4]
ConCor = cor(Exp.VegStress1, use = "pairwise.complete.obs")
corrplot.mixed(ConCor)
Vegetation Indicators
for(i in 2:ncol(Exp.VegStress1)){
MyDF = Exp.VegStress1
MyDF$IndVar = Exp.VegStress1[,i]
M.tmp = lm(IndVar ~ M2Pred, data = MyDF)
Cor.tmp = cor(MyDF$IndVar, MyDF$M2Pred, use = "pairwise.complete.obs")
RSq.tmp = summary(M.tmp)$r.squared
ARSq.tmp = summary(M.tmp)$adj.r.squared
Fst.tmp = as.numeric(summary(M.tmp)$fstatistic[1])
Est.tmp = as.numeric(M.tmp$coefficients[2])
SE.tmp = (coef(summary(M.tmp))[2, 2])
Pval.tmp = (coef(summary(M.tmp))[2, 4])
tmp.res = as.data.frame(cbind(Est.tmp, SE.tmp, Cor.tmp, RSq.tmp,
ARSq.tmp, Fst.tmp, Pval.tmp))
rownames(tmp.res) = colnames(Exp.VegStress1)[i]
if(i == 2){VegIndx.df = tmp.res}
else(VegIndx.df = rbind(VegIndx.df, tmp.res))
}
names(VegIndx.df) = c("Estimate", "SE", "r", "R Square", "Adj R Sq", "F", "pval")
VegIndx.df$Lab = rownames(VegIndx.df)
VegIndx.df2 = VegIndx.df %>%
filter(pval <= 0.05)
VegIndx.df2## Estimate SE r R Square Adj R Sq F
## 1 -0.7992547 0.2161247 -0.3741852 0.14001459 0.12977667 13.676076
## 2 0.3913856 0.1523305 0.2699297 0.07286204 0.06182468 6.601403
## pval Lab
## 1 0.0003865108 VMMI
## 2 0.0119560910 N_TOLSel.tmp = levels(factor(as.character(VegIndx.df2$Lab)))
df.subset = Exp.VegStress1[, Sel.tmp]
df.subset$M2Pred = Exp.VegStress$M2Pred
df.subset.m = melt(df.subset, "M2Pred")
ggplot(df.subset.m , aes(x=value, y=M2Pred)) +
geom_point(size=2, pch=19, col = "red") +
geom_smooth(method = "lm", se=FALSE) +
facet_wrap(~ variable, scales="free") +
theme_classic() +
xlab(NULL) +
ylab(NULL) +
theme(axis.title.y = element_text(face="bold",
size=14),
axis.text.x = element_text(face="bold",
size=11))
Vegetation Metrics
Additional vegetative indices. Match by site identifier (SITE_ID).
VegMets = read.csv("./WCA_raw/nwca2011_vegmetrics.csv")
VegMets = VegMets[match(Site.pnts$SITE_ID, VegMets$SITE_ID, nomatch=0),]
VegMets = merge(VegMets, Site.pnts@data, by = "SITE_ID", nomatch=0)Select Metrics (Vegetation)
Acronyms: D_AC = Simpson Diversity Index - alien + cryptogenic species.
D_ALIEN = Simpson Diversity Index - alien species only.
D_ALL = Simpson Diversity Index - All species.
D_NAT = Simpson Diversity Index - native species only.
D_SANDT = Simpson’s Diversity - Heterogeneity of S&T Types in AA.
D_VASC_STRATA = Simpson’s Diversity - Heterogeneity of Vertical Vascular Structure in AA based on occurrence and relative cover of all strata in all plots.
H_AC = Shannon-Wiener Diversity Index - alien + cryptogenic species.
H_ALIEN = Shannon-Wiener Diversity Index - alien species.
H_ALL = Shannon-Wiener Diversity Index - All species.
H_NAT = Shannon-Wiener Diversity Index - native species only.
H_SANDT = Shannon-Wiener - Heterogeneity of S&T Types in AA.
H_VASC_STRATA = Shannon-Wiener - Heterogeneity of Vertical Vascular Structure in AA based on occurrence and relative cover of all strata in all plots.
MEDN_AC = Median Number of Alien + Cryptogenic Species/100-m2Plot
MEDN_ADVSPP = Median Number of Adventive Species across 100-m2Plots
MEDN_ALIENSPP = Median Number of Alien Species across 100-m2 Plots
MEDN_FAM = Median Number of Families across 100-m2 Plots
MEDN_GEN = Median Number of Genera across 100-m2 Plots
MEDN_INTRSPP = Median Number of Introduced Species across 100-m2Plots
MEDN_NATSPP = Median Number of Native Species across 100-m2Plots
MEDN_SPP = Median Number of Species across 100-m2Plots
TOTN_AC = Total Number of alien and cryptogenic Species across AA (based 500 m2 sample area - 5 100-m2 plots)
TOTN_ADVSPP = Total Number of Adventive Species across AA (based 500 m2 sample area, i.e, 5 100-m2 plots)
TOTN_ALIENSPP = Total Number of Alien Species in AA (based 500 m2 sample area - 5 100-m2 plots)
Exp.VegMets = VegMets %>%
select(SITE_ID, M2Pred, D_AC, D_ALIEN, D_ALL, D_NAT, D_SANDT, D_VASC_STRATA,
H_AC, H_ALIEN, H_ALL, H_NAT, H_SANDT, H_VASC_STRATA,
MEDN_AC, MEDN_ADVSPP, MEDN_ALIENSPP, MEDN_FAM, MEDN_GEN,
MEDN_INTRSPP, MEDN_NATSPP, MEDN_SPP, TOTN_AC, TOTN_ADVSPP,
TOTN_ALIENSPP)
Exp.VegMets1 = Exp.VegMets[,2:25]
ConCor = cor(Exp.VegMets1, use = "pairwise.complete.obs")
corrplot.mixed(ConCor, tl.cex = 0.5, number.cex = 0.75)
Vegetation Indicators
for(i in 2:ncol(Exp.VegMets1)){
MyDF = Exp.VegMets1
MyDF$IndVar = Exp.VegMets1[,i]
M.tmp = lm(IndVar ~ M2Pred, data = MyDF)
Cor.tmp = cor(MyDF$IndVar, MyDF$M2Pred, use = "pairwise.complete.obs")
RSq.tmp = summary(M.tmp)$r.squared
ARSq.tmp = summary(M.tmp)$adj.r.squared
Fst.tmp = as.numeric(summary(M.tmp)$fstatistic[1])
Est.tmp = as.numeric(M.tmp$coefficients[2])
SE.tmp = (coef(summary(M.tmp))[2, 2])
Pval.tmp = (coef(summary(M.tmp))[2, 4])
tmp.res = as.data.frame(cbind(Est.tmp, SE.tmp, Cor.tmp, RSq.tmp, ARSq.tmp, Fst.tmp, Pval.tmp))
rownames(tmp.res) = colnames(Exp.VegMets1)[i]
if(i == 2){VegMet.df = tmp.res}
else(VegMet.df = rbind(VegMet.df, tmp.res))
}
names(VegMet.df) = c("Estimate", "SE", "r", "R Square", "Adj R Sq", "F", "pval")
VegMet.df$Lab = rownames(VegMet.df)
VegMet.df2 = VegMet.df %>%
filter(pval <= 0.05)
VegMet.df2## Estimate SE r R Square Adj R Sq F
## 1 0.009114523 0.002373440 0.3864499 0.14934353 0.13921667 14.747265
## 2 0.007932402 0.002329539 0.3482709 0.12129262 0.11083182 11.594964
## 3 0.018614746 0.004236283 0.4323192 0.18689986 0.17722010 19.308308
## 4 0.015914644 0.004121635 0.3882478 0.15073638 0.14062610 14.909218
## 5 0.067992438 0.011985038 0.5263174 0.27700998 0.26840295 32.184176
## 6 0.058558116 0.011160376 0.4968338 0.24684385 0.23787770 27.530656
## 7 0.059378815 0.011059520 0.5054636 0.25549344 0.24663027 28.826408
## 8 0.141238415 0.024734108 0.5288030 0.27963266 0.27105686 32.607174
## 9 0.006667381 0.002969671 0.2379319 0.05661158 0.04538076 5.040737
## 10 0.131804094 0.024218612 0.5105704 0.26068210 0.25188070 29.618243
## pval Lab
## 1 2.376500e-04 D_AC
## 2 1.016252e-03 D_ALIEN
## 3 3.228484e-05 H_AC
## 4 2.209395e-04 H_ALIEN
## 5 1.943092e-07 MEDN_AC
## 6 1.140798e-06 MEDN_ALIENSPP
## 7 6.912491e-07 MEDN_INTRSPP
## 8 1.660769e-07 TOTN_AC
## 9 2.738499e-02 TOTN_ADVSPP
## 10 5.105557e-07 TOTN_ALIENSPPSel.tmp = levels(factor(as.character(VegMet.df2$Lab)))
df.subset = Exp.VegMets1[, Sel.tmp]
df.subset$M2Pred = Exp.VegMets$M2Pred
df.subset.m = melt(df.subset, "M2Pred")
ggplot(df.subset.m , aes(x=value, y=M2Pred)) +
geom_point(size=2, pch=19, col = "red") +
geom_smooth(method = "lm", se=FALSE) +
facet_wrap(~ variable, scales="free", ncol = 5) +
theme_classic() +
xlab(NULL) +
ylab(NULL) +
theme(axis.title.y = element_text(face="bold",
size=14),
axis.text.x = element_text(face="bold",
size=11))
Soil Chemistry
SoilStress = read.csv("./WCA_raw/nwca2011_soilchem.csv")
SoilStress = SoilStress[match(Site.pnts$SITE_ID, SoilStress$SITE_ID, nomatch=0),]
SoilStress = merge(SoilStress, Site.pnts@data, by = "SITE_ID", nomatch=0)Select Indicators (Soil Stress)
Acronyms: AG mg/kg = Silver trace element analysis
AL Percent = Aluminum by ammonium oxalate extraction
AS mg/kg = Arsenic trace element analysis
BA mg/kg = Barium trace element analysis
BE mg/kg = Beryllium trace element analysis
CA cmol(+)/kg = Calcium
CD mg/kg = Cadmium trace element analysis
CEC cmol(+)/kg = Overall cation exchange capacity
CO mg/kg = Cobalt trace element analysis
CR mg/kg = Chromium trace element analysis
CU mg/kg = Copper trace element analysis
EC mmhos/cm = Electrical conductivity
FE Percent = Iron by ammonium oxalate extraction
HG ug/kg = Mercury trace element analysis
K cmol(+)/kg = Potassium CEC
MG cmol(+)/kg = Magnesium CEC
MN mg/kg = Manganese by ammonium oxalate extraction
MO mg/kg = Molybdenum trace element analysis
NI mg/kg = Nickel trace element analysis
OLSEN_P mg P/kg = Olsen Phosphorus
P mg/kg = Phosphorus by ammonium oxalate extraction
P_T mg/kg = Phosphorus trace analysis
PB mg/kg = Lead trace element analysis
PH_CACL2 pH unit = pH using 1:2.0.01 M CaCl2
PH_H2O pH unit = pH using 1:1 H20
SB mg/kg = Antimony trace element analysis
SE ug/kg = Selenium trace element analysis
SI Percent = Silicon by ammonium oxalate extraction
SN mg/kg = Tin trace element analysis
SR mg/kg = Strontium trace element analysis
TOT_CARBON Percent = Total Carbon (total C is assumed to be organic)
TOT_NITROGEN Percent = Total Nitrogen
TOT_SULFUR Percent = Total Sulfur
V mg/kg = Vanadium trace element analysis
W mg/kg = Tungsten trace element analysis
ZN mg/kg = Zinc trace element analysis
Exp.SoilStress = SoilStress %>%
select(M2Pred, TOT_CARBON, TOT_NITROGEN, TOT_SULFUR, PH_H2O, PH_CACL2,
CEC, CA, K, AL, FE, MN, P, SI, EC, OLSEN_P, AG, AS, BA, BE, CD, CO, CR, CU, HG,
MN_T, MO, NI, P_T, PB, SB, SE, SN, SR, V, W, ZN)
SC.CorTab = cor(Exp.SoilStress, use = "pairwise.complete.obs")
corrplot.mixed(SC.CorTab, tl.cex = 0.5, number.cex = 0.75)
Soil Stress
Exp.SoilStress1 = Exp.SoilStress
for(i in 2:ncol(Exp.SoilStress1)){
MyDF = Exp.SoilStress1
MyDF$IndVar = Exp.SoilStress1[,i]
M.tmp = lm(IndVar ~ M2Pred, data = MyDF)
Cor.tmp = cor(MyDF$IndVar, MyDF$M2Pred, use = "pairwise.complete.obs")
RSq.tmp = summary(M.tmp)$r.squared
ARSq.tmp = summary(M.tmp)$adj.r.squared
Fst.tmp = as.numeric(summary(M.tmp)$fstatistic[1])
Est.tmp = as.numeric(M.tmp$coefficients[2])
SE.tmp = (coef(summary(M.tmp))[2, 2])
Pval.tmp = (coef(summary(M.tmp))[2, 4])
tmp.res = as.data.frame(cbind(Est.tmp, SE.tmp, Cor.tmp,
RSq.tmp, ARSq.tmp, Fst.tmp, Pval.tmp))
rownames(tmp.res) = colnames(Exp.SoilStress1)[i]
if(i == 2){SoilChem.df = tmp.res}
else(SoilChem.df = rbind(SoilChem.df, tmp.res))
}
names(SoilChem.df) = c("Estimate", "SE", "r", "R Square", "Adj R Sq", "F", "pval")
SoilChem.df$Lab = rownames(SoilChem.df)
SoilChem.df2 = SoilChem.df %>%
filter(pval <= 0.05)
SoilChem.df2## Estimate SE r R Square Adj R Sq F
## 1 9.931230e-03 3.068784e-03 0.3347291 0.11204359 0.10134532 10.473057
## 2 3.617999e+02 5.607308e+01 0.5779618 0.33403985 0.32601623 41.632082
## 3 1.945844e-03 3.254778e-04 0.5486372 0.30100276 0.29258110 35.741527
## 4 1.287975e+00 3.796821e-01 0.3508040 0.12306343 0.11236908 11.507333
## 5 4.280231e-03 1.240387e-03 0.3542091 0.12546406 0.11492748 11.907477
## 6 1.204682e-01 1.754982e-02 0.6017675 0.36212408 0.35443883 47.119351
## 7 7.982220e-02 2.921044e-02 0.2873025 0.08254275 0.07148905 7.467431
## 8 9.530720e-01 4.575710e-01 0.2228766 0.04967398 0.03822427 4.338448
## 9 3.354335e+00 9.147448e-01 0.3733902 0.13942026 0.12905183 13.446612
## 10 1.661994e+00 7.753049e-01 0.2290427 0.05246055 0.04104442 4.595298
## 11 4.207942e-01 1.274785e-01 0.3406507 0.11604288 0.10539279 10.895957
## 12 4.961185e+02 7.481178e+01 0.5885075 0.34634107 0.33846566 43.977535
## 13 1.195337e-02 3.838511e-03 0.3234400 0.10461344 0.09382565 9.697393
## 14 1.106587e+00 3.964899e-01 0.2929110 0.08579685 0.07478235 7.789449
## 15 2.972645e+00 4.723143e-01 0.5683894 0.32306648 0.31491065 39.611744
## pval Lab
## 1 1.740892e-03 FE
## 2 6.960412e-09 P
## 3 5.449748e-08 SI
## 4 1.069302e-03 OLSEN_P
## 5 8.815985e-04 AG
## 6 1.121066e-09 CD
## 7 7.675109e-03 CO
## 8 4.033886e-02 CR
## 9 4.319956e-04 CU
## 10 3.498604e-02 MN_T
## 11 1.422091e-03 NI
## 12 3.156791e-09 P_T
## 13 2.532581e-03 SB
## 14 6.519036e-03 V
## 15 1.393088e-08 ZNSel.tmp = levels(factor(as.character(SoilChem.df2$Lab)))
df.subset = Exp.SoilStress1[, Sel.tmp]
df.subset$M2Pred = Exp.SoilStress$M2Pred
df.subset.m = melt(df.subset, "M2Pred")
ggplot(df.subset.m , aes(x=value, y=M2Pred)) +
geom_point(size=2, pch=19, col = "red") +
geom_smooth(method = "lm", se=FALSE) +
facet_wrap(~ variable, scales="free", ncol = 4) +
theme_classic() +
xlab(NULL) +
ylab(NULL) +
theme(axis.title.y = element_text(face="bold",
size=14),
axis.text.x = element_text(face="bold",
size=11))## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
Combined Table for WCA
Validation Results Tbale for WCA
Validation.df = rbind(VegIndx.df, VegMet.df, SoilChem.df)
colnames(Validation.df) = c("Estimate", "SE", "R", "RSquare", "AdjRSq", "F", "pval", "Lab")
Validation.df2 = Validation.df %>%
filter(pval <= 0.05)
Validation.df2a = Validation.df2[,1:7]
Validation.df2a$pval = as.numeric(format(Validation.df2a$pval, digits = 4, scientific = FALSE))
RndFunc = function(x){round(x, 4)}
Validation.df2a = as.data.frame(apply(Validation.df2a, 2, RndFunc))
rownames(Validation.df2a) = Validation.df2$Lab
Validation.df2a## Estimate SE R RSquare AdjRSq F pval
## VMMI -0.7993 0.2161 -0.3742 0.1400 0.1298 13.6761 0.0004
## N_TOL 0.3914 0.1523 0.2699 0.0729 0.0618 6.6014 0.0120
## D_AC 0.0091 0.0024 0.3864 0.1493 0.1392 14.7473 0.0002
## D_ALIEN 0.0079 0.0023 0.3483 0.1213 0.1108 11.5950 0.0010
## H_AC 0.0186 0.0042 0.4323 0.1869 0.1772 19.3083 0.0000
## H_ALIEN 0.0159 0.0041 0.3882 0.1507 0.1406 14.9092 0.0002
## MEDN_AC 0.0680 0.0120 0.5263 0.2770 0.2684 32.1842 0.0000
## MEDN_ALIENSPP 0.0586 0.0112 0.4968 0.2468 0.2379 27.5307 0.0000
## MEDN_INTRSPP 0.0594 0.0111 0.5055 0.2555 0.2466 28.8264 0.0000
## TOTN_AC 0.1412 0.0247 0.5288 0.2796 0.2711 32.6072 0.0000
## TOTN_ADVSPP 0.0067 0.0030 0.2379 0.0566 0.0454 5.0407 0.0274
## TOTN_ALIENSPP 0.1318 0.0242 0.5106 0.2607 0.2519 29.6182 0.0000
## FE 0.0099 0.0031 0.3347 0.1120 0.1013 10.4731 0.0017
## P 361.7999 56.0731 0.5780 0.3340 0.3260 41.6321 0.0000
## SI 0.0019 0.0003 0.5486 0.3010 0.2926 35.7415 0.0000
## OLSEN_P 1.2880 0.3797 0.3508 0.1231 0.1124 11.5073 0.0011
## AG 0.0043 0.0012 0.3542 0.1255 0.1149 11.9075 0.0009
## CD 0.1205 0.0175 0.6018 0.3621 0.3544 47.1194 0.0000
## CO 0.0798 0.0292 0.2873 0.0825 0.0715 7.4674 0.0077
## CR 0.9531 0.4576 0.2229 0.0497 0.0382 4.3384 0.0403
## CU 3.3543 0.9147 0.3734 0.1394 0.1291 13.4466 0.0004
## MN_T 1.6620 0.7753 0.2290 0.0525 0.0410 4.5953 0.0350
## NI 0.4208 0.1275 0.3407 0.1160 0.1054 10.8960 0.0014
## P_T 496.1185 74.8118 0.5885 0.3463 0.3385 43.9775 0.0000
## SB 0.0120 0.0038 0.3234 0.1046 0.0938 9.6974 0.0025
## V 1.1066 0.3965 0.2929 0.0858 0.0748 7.7894 0.0065
## ZN 2.9726 0.4723 0.5684 0.3231 0.3149 39.6117 0.0000#ModComp = xtable(Validation.df2a, digits = 3)
#print(ModComp, include.rownames = TRUE)Validation for Coastal Areas
CSites = read.csv("D:/Dist031017/WCA_raw/Coastal2010/assessed_ncca2010_siteinfo.revised.06212016.csv")
CSites = CSites %>%
filter(STATE == "FL")
dd = cbind(CSites$ALON_DD, CSites$ALAT_DD)
names(dd) = c("Long", "Lat")
CSite.pnts = SpatialPointsDataFrame(dd, CSites)
proj4string(CSite.pnts) = proj4string(FloridaLL)
CSite.pnts = spTransform(CSite.pnts, nProj)
CSite.pnts$Region = 1:nrow(CSite.pnts)Coastal Assessment Sites
Locate the predicted value at each coastal field site.
CSite.pnts$Pd = extract(r.filled, spTransform(CSite.pnts, CRS(proj4string(r.filled))), method = "simple")
CSite.pnts$Pd = ifelse(is.na(CSite.pnts$Pd) == TRUE, mean(CSite.pnts$Pd, na.rm=T), CSite.pnts$Pd)
CSite.pnts$Pd = as.numeric(scale(CSite.pnts$Pd, scale=TRUE, center=TRUE))
range(CSite.pnts$Pd)## [1] -0.5604156 6.1159679#Nearest Road
Grd.pp = as.ppp(CSite.pnts)
CSite.pnts$nRd = nncross(Grd.pp, Rd.psp)[,"dist"]
CSite.pnts$nRd = as.numeric(scale(CSite.pnts$nRd, scale=TRUE, center=TRUE))
range(CSite.pnts$nRd)## [1] -0.925175 3.176584#Nearest Well
CSite.pnts$nW = nncross(Grd.pp, W.pp)[,"dist"]
CSite.pnts$nW = as.numeric(scale(CSite.pnts$nW, scale=TRUE, center=TRUE))
range(CSite.pnts$nW)## [1] -0.6757062 3.9690630#GRF
Ap = inla.spde.make.A(mesh.dom, loc=CSite.pnts)
CSite.pnts$m2RE = drop(Ap %*% Model2$summary.random$field0$mean)
Fixed.m2 = as.data.frame(Model2$summary.fixed)[,"mean"]
m2Int = Fixed.m2[1]
m2Pdm = Fixed.m2[2]
m2nRdm = Fixed.m2[3]
m2nWm = Fixed.m2[4]
CSite.pnts@data = CSite.pnts@data %>%
mutate(M2Pred = m2RE + m2Int + m2Pdm*Pd + m2nRdm*nRd + m2nWm*nW)Fish Tissue Contaminants Assessment
FTiss = read.csv("D:/Dist031017/WCA_raw/Coastal2010/assessed_ncca2010_ecological_fish_tissue_contaminant_data.csv")
FTiss = FTiss %>%
filter(STATE == "FL")
Test.param = levels(factor(FTiss$PARAMETER))
for (i in 1:length(Test.param)) {
FTiss1 = FTiss %>%
filter(PARAMETER == Test.param[i])
if (range(FTiss1$RESULT)[2] == 0) {next}
else {
FTiss1 = FTiss1[match(CSite.pnts$SITE_ID, FTiss1$SITE_ID, nomatch=0),]
FTiss1 = merge(FTiss1, CSite.pnts@data, by = "SITE_ID", nomatch=0)
M.tmp = lm(RESULT ~ M2Pred, data = FTiss1)
Cor.tmp = cor(FTiss1$RESULT, FTiss1$M2Pred, use = "pairwise.complete.obs")
RSq.tmp = summary(M.tmp)$r.squared
ARSq.tmp = summary(M.tmp)$adj.r.squared
Fst.tmp = as.numeric(summary(M.tmp)$fstatistic[1])
Est.tmp = as.numeric(M.tmp$coefficients[2])
SE.tmp = (coef(summary(M.tmp))[2, 2])
Pval.tmp = (coef(summary(M.tmp))[2, 4])
tmp.res = as.data.frame(cbind(Est.tmp, SE.tmp, Cor.tmp, RSq.tmp,
ARSq.tmp, Fst.tmp, Pval.tmp))
rownames(tmp.res) = Test.param[i]
if(i == 1){FishTiss.df = tmp.res}
else(FishTiss.df = rbind(FishTiss.df, tmp.res))}
}
FishTiss.df$Labs = rownames(FishTiss.df)
FishTiss.df2 = FishTiss.df %>%
filter(Pval.tmp <= 0.05)
rownames(FishTiss.df2) = FishTiss.df2$Labs
FishTiss.df2 = FishTiss.df2[,1:7]
names(FishTiss.df2) = c("Estimate", "SE", "R", "RSquare", "AdjRSq", "F", "pval")
FishTiss.df2## Estimate SE R RSquare AdjRSq
## ALPHACHL 0.015129866 0.006400663 0.2417739 0.05845464 0.04799302
## CISNONA 0.054570849 0.007207261 0.6237996 0.38912591 0.38233842
## CU 0.054381233 0.009640361 0.5110872 0.26121009 0.25300131
## FE 3.869070541 1.121702731 0.3417016 0.11675997 0.10694619
## GAMMACHL 0.013048793 0.006537805 0.2058791 0.04238622 0.03174606
## PB 0.015880420 0.003510253 0.4304353 0.18527458 0.17622207
## PCB101 0.094474273 0.023034672 0.3968280 0.15747247 0.14811105
## PCB110 0.048942858 0.013882898 0.3483364 0.12133822 0.11157531
## PCB118 0.060979525 0.017764667 0.3402432 0.11576544 0.10594061
## PCB119 0.005943876 0.002514546 0.2417739 0.05845464 0.04799302
## PCB138 0.192434220 0.034510588 0.5067232 0.25676838 0.24851025
## PCB141 0.005943876 0.002514546 0.2417739 0.05845464 0.04799302
## PCB149 0.063019932 0.018827868 0.3327201 0.11070266 0.10082158
## PCB153 0.281084890 0.067683013 0.4010191 0.16081633 0.15149207
## PCB174 0.007564933 0.003200331 0.2417739 0.05845464 0.04799302
## PCB18 0.015129866 0.006400663 0.2417739 0.05845464 0.04799302
## PCB180 0.079248861 0.015987386 0.4631027 0.21446412 0.20573595
## PCB183 0.012452098 0.005958388 0.2151308 0.04628127 0.03568439
## PCB187 0.055561908 0.019089405 0.2933115 0.08603164 0.07587643
## PCB194 0.005403523 0.002285951 0.2417739 0.05845464 0.04799302
## PCB28 0.005943876 0.002514546 0.2417739 0.05845464 0.04799302
## PCB44 0.011347399 0.004800497 0.2417739 0.05845464 0.04799302
## PCB52 0.139943192 0.016675529 0.6625708 0.43900009 0.43276675
## PCB66 0.012968456 0.005486283 0.2417739 0.05845464 0.04799302
## PPDDE 0.255517537 0.118283897 0.2220225 0.04929398 0.03873058
## TNONCHL 0.141197681 0.014063661 0.7268429 0.52830062 0.52305951
## F pval
## ALPHACHL 5.587535 2.024010e-02
## CISNONA 57.329870 3.085099e-11
## CU 31.820830 1.932360e-07
## FE 11.897556 8.572556e-04
## GAMMACHL 3.983610 4.896802e-02
## PB 20.466665 1.846938e-05
## PCB101 16.821435 8.992849e-05
## PCB110 12.428491 6.673055e-04
## PCB118 11.782948 9.051180e-04
## PCB119 5.587535 2.024010e-02
## PCB138 31.092803 2.550843e-07
## PCB141 5.587535 2.024010e-02
## PCB149 11.203496 1.192969e-03
## PCB153 17.247082 7.449253e-05
## PCB174 5.587535 2.024010e-02
## PCB18 5.587535 2.024010e-02
## PCB180 24.571470 3.346123e-06
## PCB183 4.367445 3.945312e-02
## PCB187 8.471680 4.545792e-03
## PCB194 5.587535 2.024010e-02
## PCB28 5.587535 2.024010e-02
## PCB44 5.587535 2.024010e-02
## PCB52 70.427832 6.307618e-13
## PCB66 5.587535 2.024010e-02
## PPDDE 4.666488 3.341406e-02
## TNONCHL 100.799487 2.360431e-16FishTiss.df2$pval = as.numeric(format(FishTiss.df2$pval, digits = 4, scientific = FALSE))
FishTiss.df2a = as.data.frame(apply(FishTiss.df2, 2, RndFunc))
FishTiss.df2a## Estimate SE R RSquare AdjRSq F pval
## ALPHACHL 0.0151 0.0064 0.2418 0.0585 0.0480 5.5875 0.0202
## CISNONA 0.0546 0.0072 0.6238 0.3891 0.3823 57.3299 0.0000
## CU 0.0544 0.0096 0.5111 0.2612 0.2530 31.8208 0.0000
## FE 3.8691 1.1217 0.3417 0.1168 0.1069 11.8976 0.0009
## GAMMACHL 0.0130 0.0065 0.2059 0.0424 0.0317 3.9836 0.0490
## PB 0.0159 0.0035 0.4304 0.1853 0.1762 20.4667 0.0000
## PCB101 0.0945 0.0230 0.3968 0.1575 0.1481 16.8214 0.0001
## PCB110 0.0489 0.0139 0.3483 0.1213 0.1116 12.4285 0.0007
## PCB118 0.0610 0.0178 0.3402 0.1158 0.1059 11.7829 0.0009
## PCB119 0.0059 0.0025 0.2418 0.0585 0.0480 5.5875 0.0202
## PCB138 0.1924 0.0345 0.5067 0.2568 0.2485 31.0928 0.0000
## PCB141 0.0059 0.0025 0.2418 0.0585 0.0480 5.5875 0.0202
## PCB149 0.0630 0.0188 0.3327 0.1107 0.1008 11.2035 0.0012
## PCB153 0.2811 0.0677 0.4010 0.1608 0.1515 17.2471 0.0001
## PCB174 0.0076 0.0032 0.2418 0.0585 0.0480 5.5875 0.0202
## PCB18 0.0151 0.0064 0.2418 0.0585 0.0480 5.5875 0.0202
## PCB180 0.0792 0.0160 0.4631 0.2145 0.2057 24.5715 0.0000
## PCB183 0.0125 0.0060 0.2151 0.0463 0.0357 4.3674 0.0395
## PCB187 0.0556 0.0191 0.2933 0.0860 0.0759 8.4717 0.0045
## PCB194 0.0054 0.0023 0.2418 0.0585 0.0480 5.5875 0.0202
## PCB28 0.0059 0.0025 0.2418 0.0585 0.0480 5.5875 0.0202
## PCB44 0.0113 0.0048 0.2418 0.0585 0.0480 5.5875 0.0202
## PCB52 0.1399 0.0167 0.6626 0.4390 0.4328 70.4278 0.0000
## PCB66 0.0130 0.0055 0.2418 0.0585 0.0480 5.5875 0.0202
## PPDDE 0.2555 0.1183 0.2220 0.0493 0.0387 4.6665 0.0334
## TNONCHL 0.1412 0.0141 0.7268 0.5283 0.5231 100.7995 0.0000#ModComp = xtable(FishTiss.df2a, digits = 3)
#print(ModComp, include.rownames = TRUE)Coastal Water Chemistry
WChem = read.csv("D:/Dist031017/WCA_raw/Coastal2010/assessed_ncca2010_waterchem.csv")
WChem = WChem %>%
filter(STATE == "FL")
Test.param = levels(factor(WChem$PARAMETER))
for (i in 1:length(Test.param)) {
WChem1 = WChem %>%
filter(PARAMETER == Test.param[i])
if (range(WChem1$RESULT)[2] == 0) {next}
else {
WChem1 = WChem1[match(CSite.pnts$SITE_ID, WChem1$SITE_ID, nomatch=0),]
WChem1 = merge(WChem1, CSite.pnts@data, by = "SITE_ID", nomatch=0)
M.tmp = lm(RESULT ~ M2Pred, data = WChem1)
Cor.tmp = cor(WChem1$RESULT, WChem1$M2Pred, use = "pairwise.complete.obs")
RSq.tmp = summary(M.tmp)$r.squared
ARSq.tmp = summary(M.tmp)$adj.r.squared
Fst.tmp = as.numeric(summary(M.tmp)$fstatistic[1])
Est.tmp = as.numeric(M.tmp$coefficients[2])
SE.tmp = (coef(summary(M.tmp))[2, 2])
Pval.tmp = (coef(summary(M.tmp))[2, 4])
tmp.res = as.data.frame(cbind(Est.tmp, SE.tmp, Cor.tmp,
RSq.tmp, ARSq.tmp, Fst.tmp, Pval.tmp))
rownames(tmp.res) = Test.param[i]
if(i == 1){WChem.df = tmp.res}
else(WChem.df = rbind(WChem.df, tmp.res))}
}
WChem.df$Labs = rownames(WChem.df)
WChem.df2 = WChem.df %>%
filter(Pval.tmp <= 0.05)
rownames(WChem.df2) = WChem.df2$Labs
WChem.df2 = WChem.df2[,1:7]
names(WChem.df2) = c("Estimate", "SE", "R", "RSquare", "AdjRSq", "F", "pval")
WChem.df2## Estimate SE R RSquare AdjRSq F
## DIN 0.006236587 0.0015126564 0.3967324 0.1573966 0.1481372 16.99861
## NH3 0.005408922 0.0008573623 0.5516213 0.3042861 0.2966408 39.80089
## pval
## DIN 8.248846e-05
## NH3 1.000187e-08WChem.df2$pval = as.numeric(format(WChem.df2$pval, digits = 4, scientific = FALSE))
WChem.df2a = as.data.frame(apply(WChem.df2, 2, RndFunc))
WChem.df2a## Estimate SE R RSquare AdjRSq F pval
## DIN 0.0062 0.0015 0.3967 0.1574 0.1481 16.9986 1e-04
## NH3 0.0054 0.0009 0.5516 0.3043 0.2966 39.8009 0e+00#ModComp = xtable(WChem.df2a, digits = 3)
#print(ModComp, include.rownames = TRUE)Coastal Sediment Chemistry
STox = read.csv("D:/Dist031017/WCA_raw/Coastal2010/assessed_ncca2010_sediment_chemistry.revised.06.21.2016.csv")
STox = STox %>%
filter(STATE == "FL")
Test.param = levels(factor(STox$PARAMETER))
for (i in 1:length(Test.param)) {
STox1 = STox %>%
filter(PARAMETER == Test.param[i])
if (range(STox1$RESULT)[2] == 0) {next}
else {
STox1 = STox1[match(CSite.pnts$SITE_ID, STox1$SITE_ID, nomatch=0),]
STox1 = merge(STox1, CSite.pnts@data, by = "SITE_ID", nomatch=0)
M.tmp = lm(RESULT ~ M2Pred, data = STox1)
Cor.tmp = cor(STox1$RESULT, STox1$M2Pred, use = "pairwise.complete.obs")
RSq.tmp = summary(M.tmp)$r.squared
ARSq.tmp = summary(M.tmp)$adj.r.squared
Fst.tmp = as.numeric(summary(M.tmp)$fstatistic[1])
Est.tmp = as.numeric(M.tmp$coefficients[2])
SE.tmp = (coef(summary(M.tmp))[2, 2])
Pval.tmp = (coef(summary(M.tmp))[2, 4])
tmp.res = as.data.frame(cbind(Est.tmp, SE.tmp, Cor.tmp,
RSq.tmp, ARSq.tmp, Fst.tmp, Pval.tmp))
rownames(tmp.res) = Test.param[i]
if(i == 1){STox.df = tmp.res}
else(STox.df = rbind(STox.df, tmp.res))}
}
STox.df$Labs = rownames(STox.df)
STox.df2 = STox.df %>%
filter(Pval.tmp <= 0.05)
rownames(STox.df2) = STox.df2$Labs
STox.df2 = STox.df2[,1:7]
names(STox.df2) = c("Estimate", "SE", "R", "RSquare", "AdjRSq", "F", "pval")
STox.df2## Estimate SE R RSquare AdjRSq
## ACENTHE 9.948396e+00 1.273203218 0.6357552 0.40418468 0.39756451
## ACENTHY 5.730889e-01 0.078932565 0.6077591 0.36937108 0.36236409
## ANTHRA 5.131563e+01 6.563112842 0.6360033 0.40450022 0.39788355
## BENANTH 1.595694e+02 20.397724768 0.6362022 0.40475325 0.39813939
## BENAPY 7.368317e+01 9.400190950 0.6369552 0.40571189 0.39910868
## BENEPY 5.961468e+01 7.604292598 0.6370099 0.40578165 0.39917922
## BENZOBFL 1.482578e+02 18.946439324 0.6363085 0.40488851 0.39827616
## BENZOKFL 7.726452e+01 9.857428260 0.6369420 0.40569510 0.39909171
## BENZOP 2.967480e+01 3.820515215 0.6334952 0.40131620 0.39466416
## BIPHENYL 2.315222e-01 0.040073288 0.5201358 0.27054121 0.26243611
## CHRYSENE 1.936634e+02 24.765973038 0.6360489 0.40455820 0.39794218
## CU 3.236619e+01 4.158088853 0.6343090 0.40234788 0.39570730
## DIBENZ 1.234053e+01 1.578725013 0.6359053 0.40437558 0.39775753
## DIBENZO 1.424687e+01 1.824836889 0.6354411 0.40378544 0.39716084
## DIMETH 5.394908e-01 0.076557471 0.5962973 0.35557052 0.34841019
## FLUORANT 5.538428e+02 70.874417650 0.6357916 0.40423091 0.39761125
## FLUORENE 1.493555e+01 1.912361275 0.6355771 0.40395830 0.39733561
## HG 1.842333e-02 0.003360588 0.5003387 0.25033880 0.24200923
## INDENO 2.996694e+01 3.842023374 0.6350808 0.40332768 0.39669799
## MENAP1 3.234181e-01 0.054792621 0.5282804 0.27908016 0.27106994
## MENAP2 2.966061e-01 0.076689465 0.3775162 0.14251845 0.13299088
## MEPHEN1 1.281369e+01 1.642842047 0.6350773 0.40332322 0.39669348
## PB 2.684930e+00 0.418746450 0.5599656 0.31356148 0.30593438
## PCB141 8.533721e-03 0.003366318 0.2581579 0.06664548 0.05627488
## PCB149 2.469454e-02 0.010445588 0.2418043 0.05846934 0.04800789
## PCB151 1.199570e-02 0.004991884 0.2455477 0.06029368 0.04985250
## PCB187 8.421642e-03 0.003552260 0.2424467 0.05878041 0.04832241
## PERYLENE 1.866207e+01 3.380968317 0.5029025 0.25291088 0.24460989
## PHENANTH 2.844433e+02 36.406262774 0.6357241 0.40414511 0.39752450
## PYRENE 4.262386e+02 54.549745868 0.6357593 0.40418985 0.39756974
## SB 2.840485e-02 0.005677422 0.4664804 0.21760396 0.20891068
## SN 5.348508e-01 0.092460844 0.5206051 0.27102969 0.26293002
## TRIMETH 5.723433e-01 0.074752207 0.6280441 0.39443940 0.38771094
## ZN 1.349293e+01 1.838446941 0.6118949 0.37441539 0.36746445
## F pval
## ACENTHE 61.053519 9.854662e-12
## ACENTHY 52.714673 1.324287e-10
## ANTHRA 61.133557 9.618963e-12
## BENANTH 61.197801 9.433949e-12
## BENAPY 61.441696 8.764070e-12
## BENEPY 61.459475 8.717183e-12
## BENZOBFL 61.232168 9.336476e-12
## BENZOKFL 61.437419 8.775389e-12
## BENZOP 60.329774 1.227365e-11
## BIPHENYL 33.379143 1.073013e-07
## CHRYSENE 61.148274 9.576253e-12
## CU 60.589275 1.134325e-11
## DIBENZ 61.101932 9.711396e-12
## DIBENZO 60.952369 1.016101e-11
## DIMETH 49.658415 3.571609e-10
## FLUORANT 61.065239 9.819782e-12
## FLUORENE 60.996146 1.002724e-11
## HG 30.054232 3.802923e-07
## INDENO 60.836560 1.052373e-11
## MENAP1 34.840510 6.225989e-08
## MENAP2 14.958526 2.074550e-04
## MEPHEN1 60.835432 1.052733e-11
## PB 41.111523 6.497286e-09
## PCB141 6.426383 1.297053e-02
## PCB149 5.589027 2.022389e-02
## PCB151 5.774603 1.831189e-02
## PCB187 5.620619 1.988411e-02
## PERYLENE 30.467555 3.242626e-07
## PHENANTH 61.043487 9.884618e-12
## PYRENE 61.054828 9.850759e-12
## SB 25.031258 2.775766e-06
## SN 33.461819 1.040279e-07
## TRIMETH 58.622614 2.068909e-11
## ZN 53.865431 9.167338e-11STox.df2$pval = as.numeric(format(STox.df2$pval, digits = 4, scientific = FALSE))
STox.df2a = as.data.frame(apply(STox.df2, 2, RndFunc))
STox.df2a## Estimate SE R RSquare AdjRSq F pval
## ACENTHE 9.9484 1.2732 0.6358 0.4042 0.3976 61.0535 0.0000
## ACENTHY 0.5731 0.0789 0.6078 0.3694 0.3624 52.7147 0.0000
## ANTHRA 51.3156 6.5631 0.6360 0.4045 0.3979 61.1336 0.0000
## BENANTH 159.5694 20.3977 0.6362 0.4048 0.3981 61.1978 0.0000
## BENAPY 73.6832 9.4002 0.6370 0.4057 0.3991 61.4417 0.0000
## BENEPY 59.6147 7.6043 0.6370 0.4058 0.3992 61.4595 0.0000
## BENZOBFL 148.2578 18.9464 0.6363 0.4049 0.3983 61.2322 0.0000
## BENZOKFL 77.2645 9.8574 0.6369 0.4057 0.3991 61.4374 0.0000
## BENZOP 29.6748 3.8205 0.6335 0.4013 0.3947 60.3298 0.0000
## BIPHENYL 0.2315 0.0401 0.5201 0.2705 0.2624 33.3791 0.0000
## CHRYSENE 193.6634 24.7660 0.6360 0.4046 0.3979 61.1483 0.0000
## CU 32.3662 4.1581 0.6343 0.4023 0.3957 60.5893 0.0000
## DIBENZ 12.3405 1.5787 0.6359 0.4044 0.3978 61.1019 0.0000
## DIBENZO 14.2469 1.8248 0.6354 0.4038 0.3972 60.9524 0.0000
## DIMETH 0.5395 0.0766 0.5963 0.3556 0.3484 49.6584 0.0000
## FLUORANT 553.8428 70.8744 0.6358 0.4042 0.3976 61.0652 0.0000
## FLUORENE 14.9355 1.9124 0.6356 0.4040 0.3973 60.9961 0.0000
## HG 0.0184 0.0034 0.5003 0.2503 0.2420 30.0542 0.0000
## INDENO 29.9669 3.8420 0.6351 0.4033 0.3967 60.8366 0.0000
## MENAP1 0.3234 0.0548 0.5283 0.2791 0.2711 34.8405 0.0000
## MENAP2 0.2966 0.0767 0.3775 0.1425 0.1330 14.9585 0.0002
## MEPHEN1 12.8137 1.6428 0.6351 0.4033 0.3967 60.8354 0.0000
## PB 2.6849 0.4187 0.5600 0.3136 0.3059 41.1115 0.0000
## PCB141 0.0085 0.0034 0.2582 0.0666 0.0563 6.4264 0.0130
## PCB149 0.0247 0.0104 0.2418 0.0585 0.0480 5.5890 0.0202
## PCB151 0.0120 0.0050 0.2455 0.0603 0.0499 5.7746 0.0183
## PCB187 0.0084 0.0036 0.2424 0.0588 0.0483 5.6206 0.0199
## PERYLENE 18.6621 3.3810 0.5029 0.2529 0.2446 30.4676 0.0000
## PHENANTH 284.4433 36.4063 0.6357 0.4041 0.3975 61.0435 0.0000
## PYRENE 426.2386 54.5497 0.6358 0.4042 0.3976 61.0548 0.0000
## SB 0.0284 0.0057 0.4665 0.2176 0.2089 25.0313 0.0000
## SN 0.5349 0.0925 0.5206 0.2710 0.2629 33.4618 0.0000
## TRIMETH 0.5723 0.0748 0.6280 0.3944 0.3877 58.6226 0.0000
## ZN 13.4929 1.8384 0.6119 0.3744 0.3675 53.8654 0.0000#ModComp = xtable(STox.df2a, digits = 3)
#print(ModComp, include.rownames = TRUE)sessionInfo()## R version 3.3.3 (2017-03-06)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] classInt_0.1-24 spatstat_1.51-0 rpart_4.1-10
## [4] nlme_3.1-131 dplyr_0.5.0 gridExtra_2.2.1
## [7] Thermimage_3.0.0 GISTools_0.7-4 MASS_7.3-45
## [10] ggplot2_2.2.1 corrplot_0.77 viridis_0.4.0
## [13] viridisLite_0.2.0 xtable_1.8-2 mapproj_1.2-4
## [16] maps_3.1.1 maptools_0.9-2 rgeos_0.3-23
## [19] rasterVis_0.41 latticeExtra_0.6-28 RColorBrewer_1.1-2
## [22] lattice_0.20-34 raster_2.5-8 reshape2_1.4.2
## [25] rgdal_1.2-7 INLA_0.0-1488620070 Matrix_1.2-8
## [28] sp_1.2-4 knitr_1.16
##
## loaded via a namespace (and not attached):
## [1] zoo_1.8-0 splines_3.3.3 colorspace_1.3-2
## [4] spatstat.utils_1.6-0 htmltools_0.3.6 mgcv_1.8-17
## [7] rlang_0.1.1 e1071_1.6-8 hexbin_1.27.1
## [10] foreign_0.8-67 DBI_0.6-1 plyr_1.8.4
## [13] stringr_1.2.0 munsell_0.4.3 gtable_0.2.0
## [16] evaluate_0.10 labeling_0.3 class_7.3-14
## [19] parallel_3.3.3 Rcpp_0.12.11 tensor_1.5
## [22] scales_0.4.1 backports_1.1.0 abind_1.4-5
## [25] deldir_0.1-14 digest_0.6.12 stringi_1.1.5
## [28] polyclip_1.6-1 grid_3.3.3 rprojroot_1.2
## [31] tools_3.3.3 magrittr_1.5 goftest_1.1-1
## [34] lazyeval_0.2.0 tibble_1.3.3 assertthat_0.2.0
## [37] rmarkdown_1.5 R6_2.2.1