Annual model - No Regions
Tanya Rogers
4/22/2022
library(dplyr)
library(lavaan)
library(DiagrammeR)
library(tidyr)
library(ggplot2)
source("functions/table_funcs.R")
# For saving SEM diagrams:
library(purrr)
library(DiagrammeRsvg)
library(rsvg)
library(png)
library(grid)
library(ggpubr)Import data
combined=read.csv("data/annual_averages/annual_data_compiled_noregions.csv")
cnames=read.csv("analysis/column_names.csv", stringsAsFactors = F)
dsub=filter(combined, Year>=1975)
focaldata=dsub[,cnames$Datacolumn]
fvars=cnames$Shortname
colnames(focaldata)=fvars
focaldata = focaldata %>%
mutate(tzoop=hcope+clad+mysid+pcope+rotif_m,
tzoop_e=hcope_e+clad_e+mysid_e+pcope_e+rotif_e,
hzoop=hcope+clad+rotif_m,
hzoop_e=hcope_e+clad_e+rotif_e,
pzoop=mysid+pcope,
pzoop_e=mysid_e+pcope_e,
turbid=-secchi)
fvars=c(fvars,"tzoop","tzoop_e",
"hzoop","hzoop_e",
"pzoop","pzoop_e","turbid")
cnames=rbind(cnames,data.frame(Longname = c("Total zooplankton biomass",
"Total zooplankton energy",
"Herbivorous zooplankton biomass",
"Herbivorous zooplankton energy",
"Predatory zooplankton biomass",
"Predatory zooplankton energy",
"Turbidity"),
Shortname=c("tzoop","tzoop_e",
"hzoop","hzoop_e",
"pzoop","pzoop_e","turbid"),
Diagramname=c("total zooplankton",
"total zooplankton\nenergy",
"herbivorous\nzooplankton",
"herbivorous\nzooplankton\nenergy",
"predatory\nzooplankton",
"predatory\nzooplankton\nenergy",
"turbidity"),
Datacolumn=NA,Log=c(rep("yes",6),"no"),
Color=c("black","black","#ED7D31","#ED7D31","#7030A0",
"#7030A0","#4472C4"),
Definition = c("summed zooplankton biomass",
"summed zooplankton energy",
"summed herbivorous zooplankton biomass",
"summed herbivorous zooplankton energy",
"summed predatory zooplankton biomass",
"summed predatory zooplankton energy",
"negative secchi depth")))
#focal variables
varnames=c("temp","flow","turbid","chla","hzoop","pzoop","potam","corbic","estfish","estfish_bsmt","estfish_stn")
source("analysis/semDiagramFunctions.r")Data prep
Log transform, scale
#log transform
logvars=fvars[cnames$Log=="yes"]
logtrans=function(x) {
x2=x[which(!is.na(x))]
if(any(x2==0)) {log(x+min(x2[which(x2>0)],na.rm=T))}
else {log(x)}
}
focaldatalog = focaldata %>%
mutate_at(logvars,logtrans)
#scale data
fd0=focaldatalog
tvars=fvars[-1]
fd=fd0 %>%
#lag
mutate_at(tvars,list("1"=lag)) %>%
#scale
mutate_at(-1,scale) %>%
as.data.frame()
#detrended
fd_dtr=fd0 %>%
mutate_at(tvars,function(x) { #detrend
x2=x
x2[x2==0]=NA
res=residuals(lm(x2~fd$year))
out=x
out[which(!is.na(x2))]=res
return(out)
}) %>%
#lag
mutate_at(tvars,list("1"=lag)) %>%
#scale
mutate_at(-1,scale) %>%
as.data.frame()Time series plots
Other useful plots
Breakdown of total zooplankton biomass.
## # A tibble: 5 × 3
## Var Val prop
## <chr> <dbl> <dbl>
## 1 clad 67998. 0.0342
## 2 hcope 723390. 0.364
## 3 mysid 747951. 0.376
## 4 pcope 72802. 0.0366
## 5 rotif_m 376025. 0.189## Warning: Removed 5 rows containing non-finite values (`stat_align()`).
Similarity of fish indices.
## Warning: attributes are not identical across measure variables; they will be
## dropped## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Removed 7 rows containing missing values (`geom_line()`).
## Warning: attributes are not identical across measure variables; they will be dropped
## Removed 7 rows containing missing values (`geom_line()`).
Cross-correlation matrices
(only sig correlations shown… no correction for multiple comparisons)
## Warning: `expand_scale()` was deprecated in ggplot2 3.3.0.
## ℹ Please use `expansion()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Make a definitions table
Make a table similar to table 1 in Mac Nally et al. 2010 in Eco. Apps.
## `summarise()` has grouped output by 'Shortname', 'Longname'. You can override
## using the `.groups` argument.| Variable | Years (missing) | Definition |
|---|---|---|
| Ammonia | 1980‒2020 (0) | from the Discrete Environmental Monitoring Program (EMP) at DWR - year-round |
| Amphipods catch | 1975‒2020 (1) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Amphipods mass | 1975‒2020 (1) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Phytoplankton | 1980‒2020 (0) | from the Discrete Environmental Monitoring Program (EMP) at DWR - year-round |
| Cladocera | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Cladocera catch | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Cladocera energy | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Corbicula | 1975‒2020 (1) | from the Environmental Monitoring Program (EMP) Benthic Survey at DWR - year-round |
| Dissolved Orthophos | 1980‒2020 (0) | from the Discrete Environmental Monitoring Program (EMP) at DWR - year-round |
| Estuarine fishes FMWT | 1975‒2020 (1) | fall (September - December) - midwater trawl - biomass of estuarine pelagic forage fishes |
| Estuarine fishes BSMT | 1980‒2020 (1) | year-round - midwater trawl - biomass of estuarine pelagic forage fishes |
| Estuarine fishes BSOT | 1980‒2020 (0) | year-round - otter trawl - biomass of estuarine pelagic forage fishes |
| Estuarine fishes STN | 1975‒2021 (0) | summer (June - August) - townet - biomass of estuarine pelagic forage fishes |
| Flow | 1975‒2020 (0) | year-round - mean Delta outflow (water leaving the Delta to the Bay) |
| Herbivorous copepods | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Herbivorous copepods catch | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Herbivorous copepods energy | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Herbivorous zooplankton biomass | 1975‒2020 (0) | summed herbivorous zooplankton biomass |
| Herbivorous zooplankton energy | 1975‒2020 (0) | summed herbivorous zooplankton energy |
| Longfin smelt FMWT | 1975‒2020 (1) | fall (September - December) - midwater trawl - biomass |
| Longfin smelt BSMT | 1980‒2020 (1) | year-round - midwater trawl - biomass |
| Longfin smelt BSOT | 1980‒2020 (0) | year-round - otter trawl - biomass |
| Longfin smelt STN | 1975‒2021 (0) | summer (June - August) - townet - biomass |
| Marine fishes FMWT | 1975‒2020 (1) | fall (September - December) - midwater trawl - biomass |
| Marine fishes BSMT | 1980‒2020 (1) | year-round - midwater trawl - biomass |
| Marine fishes BSOT | 1980‒2020 (0) | year-round - otter trawl - biomass |
| Marine fishes STN | 1975‒2021 (0) | summer (June - August) - townet - biomass |
| Mysids | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Mysids catch | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Mysids energy | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Nitrate and Nitrite | 1980‒2020 (0) | from the Discrete Environmental Monitoring Program (EMP) at DWR - year-round |
| Predatory copepods | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Predatory copepods catch | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Predatory copepods energy | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Potamocorbula | 1975‒2020 (1) | from the Environmental Monitoring Program (EMP) Benthic Survey at DWR - year-round |
| Predatory zooplankton biomass | 1975‒2020 (0) | summed predatory zooplankton biomass |
| Predatory zooplankton energy | 1975‒2020 (0) | summed predatory zooplankton energy |
| Rotifers catch | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Rotifers energy | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Rotifers mass | 1975‒2020 (0) | from 5 different sources - year-round - see Bashevkin et al. 2022 |
| Striped bass FMWT | 1975‒2020 (1) | fall (September - December) - midwater trawl - biomass of age 0 individuals |
| Striped bass BSMT | 1980‒2020 (1) | year-round - midwater trawl - biomass of age 0 individuals |
| Striped bass BSOT | 1980‒2020 (0) | year-round - otter trawl - biomass of age 0 individuals |
| Striped bass STN | 1975‒2021 (0) | summer (June - August) - townet - biomass of age 0 individuals |
| Secchi | 1980‒2020 (0) | from the Discrete Environmental Monitoring Program (EMP) at DWR - year-round |
| Delta smelt FMWT | 1975‒2020 (1) | fall (September - December) - midwater trawl - biomass |
| Delta smelt BSMT | 1980‒2020 (1) | year-round - midwater trawl - biomass |
| Delta smelt BSOT | 1980‒2020 (0) | year-round - otter trawl - biomass |
| Delta smelt STN | 1975‒2021 (0) | summer (June - August) - townet - biomass |
| Mississippi silverside DJFMP | 1976‒2020 (0) | year-round - beach seines - biomass |
| Temperature | 1980‒2020 (0) | from the Discrete Environmental Monitoring Program (EMP) at DWR - year-round |
| Turbidity | 1980‒2020 (0) | negative secchi depth |
| Total zooplankton biomass | 1975‒2020 (0) | summed zooplankton biomass |
| Total zooplankton energy | 1975‒2020 (0) | summed zooplankton energy |
SEM model
With and without detrending.
#1
# model1='zoop=~hcope+clad+mysid
# fish=~estfish_bsmt+estfish_bsot
# zoop~chla+potam+flow
# chla~potam+flow
# fish~zoop+flow
# '
#2
model1='chla~potam+flow+temp
hzoop~chla+potam+flow+temp
pzoop~chla+potam+flow+hzoop+temp
fish~hzoop+pzoop+flow+turbid+temp+potam
fish=~estfish+estfish_stn+estfish_bsmt
'
modfit1=sem(model1, data=fd)
modfit1_dtr=sem(model1, data=fd_dtr)
summary(modfit1, standardized=T, rsq=T)## lavaan 0.6.13 ended normally after 34 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 27
##
## Used Total
## Number of observations 40 47
##
## Model Test User Model:
##
## Test statistic 27.481
## Degrees of freedom 18
## P-value (Chi-square) 0.070
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## fish =~
## estfish 1.000 0.875 0.894
## estfish_stn 0.858 0.104 8.230 0.000 0.751 0.877
## estfish_bsmt 0.850 0.137 6.223 0.000 0.744 0.761
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## chla ~
## potam -0.342 0.181 -1.885 0.059 -0.342 -0.341
## flow 0.055 0.175 0.314 0.754 0.055 0.055
## temp -0.036 0.166 -0.213 0.831 -0.036 -0.035
## hzoop ~
## chla 0.562 0.089 6.296 0.000 0.562 0.667
## potam -0.259 0.107 -2.425 0.015 -0.259 -0.306
## flow -0.080 0.099 -0.810 0.418 -0.080 -0.095
## temp 0.045 0.094 0.473 0.636 0.045 0.053
## pzoop ~
## chla 0.574 0.086 6.643 0.000 0.574 0.610
## potam -0.029 0.078 -0.374 0.709 -0.029 -0.031
## flow -0.189 0.068 -2.761 0.006 -0.189 -0.200
## hzoop 0.422 0.108 3.899 0.000 0.422 0.379
## temp -0.067 0.065 -1.032 0.302 -0.067 -0.071
## fish ~
## hzoop 0.311 0.137 2.265 0.023 0.355 0.299
## pzoop 0.143 0.117 1.224 0.221 0.163 0.153
## flow 0.118 0.075 1.564 0.118 0.135 0.134
## turbid 0.368 0.086 4.278 0.000 0.420 0.398
## temp 0.055 0.070 0.788 0.431 0.063 0.063
## potam -0.282 0.098 -2.893 0.004 -0.323 -0.321
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .estfish 0.191 0.056 3.415 0.001 0.191 0.200
## .estfish_stn 0.170 0.047 3.630 0.000 0.170 0.231
## .estfish_bsmt 0.401 0.096 4.159 0.000 0.401 0.420
## .chla 0.846 0.189 4.472 0.000 0.846 0.848
## .hzoop 0.270 0.060 4.472 0.000 0.270 0.380
## .pzoop 0.127 0.028 4.472 0.000 0.127 0.144
## .fish 0.057 0.037 1.542 0.123 0.074 0.074
##
## R-Square:
## Estimate
## estfish 0.800
## estfish_stn 0.769
## estfish_bsmt 0.580
## chla 0.152
## hzoop 0.620
## pzoop 0.856
## fish 0.926summary(modfit1_dtr, standardized=T, rsq=T)## lavaan 0.6.13 ended normally after 34 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 27
##
## Used Total
## Number of observations 40 47
##
## Model Test User Model:
##
## Test statistic 23.871
## Degrees of freedom 18
## P-value (Chi-square) 0.159
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## fish =~
## estfish 1.000 0.577 0.611
## estfish_stn 1.084 0.319 3.396 0.001 0.625 0.624
## estfish_bsmt 1.296 0.332 3.900 0.000 0.748 0.750
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## chla ~
## potam 0.251 0.180 1.398 0.162 0.251 0.258
## flow 0.344 0.186 1.848 0.065 0.344 0.342
## temp -0.062 0.166 -0.370 0.711 -0.062 -0.061
## hzoop ~
## chla 0.628 0.128 4.889 0.000 0.628 0.631
## potam -0.068 0.149 -0.453 0.651 -0.068 -0.070
## flow 0.007 0.158 0.047 0.963 0.007 0.007
## temp 0.064 0.135 0.474 0.636 0.064 0.064
## pzoop ~
## chla 0.645 0.099 6.533 0.000 0.645 0.655
## potam 0.052 0.091 0.566 0.572 0.052 0.054
## flow -0.192 0.096 -2.003 0.045 -0.192 -0.194
## hzoop 0.314 0.096 3.262 0.001 0.314 0.317
## temp -0.088 0.082 -1.071 0.284 -0.088 -0.089
## fish ~
## hzoop 0.068 0.092 0.742 0.458 0.118 0.117
## pzoop 0.234 0.102 2.290 0.022 0.406 0.398
## flow 0.393 0.115 3.413 0.001 0.681 0.674
## turbid 0.079 0.077 1.032 0.302 0.138 0.134
## temp 0.134 0.076 1.776 0.076 0.233 0.230
## potam -0.127 0.080 -1.586 0.113 -0.221 -0.226
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .estfish 0.558 0.134 4.178 0.000 0.558 0.626
## .estfish_stn 0.613 0.148 4.149 0.000 0.613 0.611
## .estfish_bsmt 0.435 0.121 3.584 0.000 0.435 0.437
## .chla 0.900 0.201 4.472 0.000 0.900 0.907
## .hzoop 0.593 0.133 4.472 0.000 0.593 0.605
## .pzoop 0.219 0.049 4.472 0.000 0.219 0.228
## .fish 0.023 0.043 0.549 0.583 0.070 0.070
##
## R-Square:
## Estimate
## estfish 0.374
## estfish_stn 0.389
## estfish_bsmt 0.563
## chla 0.093
## hzoop 0.395
## pzoop 0.772
## fish 0.930# residuals(modfit1)
# modificationindices(modfit1, sort=T, maximum.number=20)Nice plots
Without covariances
Original units
Detrended
With covariances
Original units
Detrended
Save updated SEM diagrams
modfit1_grobs <- map(list(plot_modfit1), ~convert_html_to_grob(.x, 2000))
modfit1_figure <- ggarrange(plotlist=modfit1_grobs) %>%
annotate_figure(top = text_grob("Annual SEM (whole estuary)",
color = "black",
face = "bold",
size = 8))
ggsave('./fig_output/sem_annual_noregions.png',modfit1_figure, width=3, height=3, dpi=300, bg = "white")Save figure object for combining with annual regional figure
# Give the object an informative name for when it is loaded into the annual regional file:
sem_annual_noregions <- plot_modfit1
sem_annual_noregions_dtr <- plot_modfit2
save(sem_annual_noregions, sem_annual_noregions_dtr,
file='./fig_output/sem_annual_noregions.RData')
focaldata_we <- focaldata
save(focaldata_we, file = "./fig_output/focaldata_we.Rdata")
# save model data for cross-correlation plots
save(fd, fd_dtr, file = "./fig_output/data_annual_noregions.Rdata")Save model coefficients for manuscript table
coeftable_we<-coef_tabler(modfit1, name="")
save(coeftable_we, file = "./fig_output/coeftable_we.Rdata")