Skip to content

Commit

Permalink
add'l script for report version of figure (multi-panel)
Browse files Browse the repository at this point in the history
  • Loading branch information
@milesmedina
milesmedina committed Oct 1, 2024
1 parent 5f55098 commit 10c1f5e
Show file tree
Hide file tree
Showing 6 changed files with 347 additions and 38 deletions.
21 changes: 10 additions & 11 deletions R/logistic_model_otb_85.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@ library(dplyr)
library(plyr)
library(lubridate)
library(ggplot2)
library(here)

load(here('data-clean/epcwq_clean.RData'))
load(here('data/loads.RData'))
Expand Down Expand Up @@ -120,16 +119,16 @@ dat <- inner_join( chl[,c('date','chl','chl_year_t1')],
y0 = plot_data$lwr$lwr[which(plot_data$TN_load==TN_targets[i])],
lty = 2, col = seg_colors[i]
)
text( x = rep(text_pos[i],3),
y = c( plot_data$median[which(plot_data$TN_load==TN_targets[i])],
plot_data$upr$upr[which(plot_data$TN_load==TN_targets[i])],
plot_data$lwr$lwr[which(plot_data$TN_load==TN_targets[i])]
),
labels = round( c( 100*plot_data$median[which(plot_data$TN_load==TN_targets[i])],
100*plot_data$upr$upr[which(plot_data$TN_load==TN_targets[i])],
100*plot_data$lwr$lwr[which(plot_data$TN_load==TN_targets[i])]
),0 ), cex = 0.7, offset = 0.1, pos = 3, col = seg_colors[i]
)
# text( x = rep(text_pos[i],3),
# y = c( plot_data$median[which(plot_data$TN_load==TN_targets[i])],
# plot_data$upr$upr[which(plot_data$TN_load==TN_targets[i])],
# plot_data$lwr$lwr[which(plot_data$TN_load==TN_targets[i])]
# ),
# labels = round( c( 100*plot_data$median[which(plot_data$TN_load==TN_targets[i])],
# 100*plot_data$upr$upr[which(plot_data$TN_load==TN_targets[i])],
# 100*plot_data$lwr$lwr[which(plot_data$TN_load==TN_targets[i])]
# ),0 ), cex = 0.7, offset = 0.1, pos = 3, col = seg_colors[i]
# )
}
legend( 'topright', bty = 'n', lty = 2, col = seg_colors,
legend = paste0(TN_targets," tons/month")
Expand Down
53 changes: 26 additions & 27 deletions R/logistic_model_otb_85_norm.R
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,10 @@

rm(list=ls(all=TRUE))
rm(list=ls(all=TRUE))

library(dplyr)
library(plyr)
library(lubridate)
library(ggplot2)
library(here)

load(here('data-clean/epcwq_clean.RData')) # WQ data
load(here('data/loads.RData')) # monthly TN loads
Expand All @@ -20,12 +19,12 @@ chl <- epcwq3[ which( epcwq3$param=="Chla" &
epcwq3$subseg %in% subseg ),
c('date','value') ]
colnames( chl ) <- c('date','chl')
chl <- chl |> group_by(date) |> dplyr::summarise(chl=mean(chl)) |>
chl <- chl |> group_by(date) |> dplyr::summarise(chl=mean(chl)) |>
data.frame()
chl$year <- year( chl$date )

# Annual mean chlorophyll
chl_yr <- chl |> group_by(year) |> dplyr::summarise(chl=mean(chl)) |>
chl_yr <- chl |> group_by(year) |> dplyr::summarise(chl=mean(chl)) |>
as.data.frame()
chl$chl_year_t1 <- mapvalues( chl$year, chl_yr$year, chl_yr$chl )

Expand Down Expand Up @@ -55,7 +54,7 @@ colnames( hydro ) <- c('year','date','hydro_load')
hydro_yr <- totanndat[ which(totanndat$year>=2000 &
totanndat$bay_segment=="Old Tampa Bay"),
c('year','hy_load')] |> as.data.frame()
# hydro_yr <- hydro |> group_by(year) |> dplyr::summarise(hydro_load=sum(hydro_load)) |>
# hydro_yr <- hydro |> group_by(year) |> dplyr::summarise(hydro_load=sum(hydro_load)) |>
# as.data.frame()
hydro$hydro_load_year_t1 <- mapvalues( hydro$year, hydro_yr$year, hydro_yr$hy_load )

Expand All @@ -72,20 +71,20 @@ deliv <- deliv[ ,c('date','ratio','ratio_year_t1') ]

# Input data
dat <- inner_join( chl[,c('date','chl','chl_year_t1')],
deliv[,c('date','ratio','ratio_year_t1')],
deliv[,c('date','ratio','ratio_year_t1')],
by = 'date' )


# Logistic model
# Specify chlorophyll target
chl_target <- 8.5

# Create binary response variable
dat$y <- ifelse( dat$chl > chl_target, 0, 1 )
# Fit the model
mod <- glm( y ~ ratio, #+ chl_year_t1 + ratio_year_t1,
data = dat, family = 'binomial' )

# Generate predictions
toprd <- data.frame( ratio = seq( 0, max(dat$ratio), 0.01 ) )
# toprd <- expand.grid(
Expand All @@ -95,21 +94,21 @@ dat <- inner_join( chl[,c('date','chl','chl_year_t1')],
# )
prds <- predict( mod, type = 'response',
newdata = toprd, se.fit = TRUE )
toplo <- toprd |>
toplo <- toprd |>
mutate(
prd = prds$fit,
hival = prds$fit + 1.96 * prds$se.fit,
loval = prds$fit - 1.96 * prds$se.fit
)

# Summarise predictions
plot_data <- toplo |> group_by(ratio) |>
dplyr::summarise(median=median(prd)) |> as.data.frame()
plot_data$upr <- toplo |> group_by(ratio) |>
dplyr::summarise(upr=median(hival)) |> as.data.frame() |> select(upr)
dplyr::summarise(upr=median(hival)) |> as.data.frame() |> select(upr)
plot_data$lwr <- toplo |> group_by(ratio) |>
dplyr::summarise(lwr=median(loval)) |> as.data.frame() |> select(lwr)
dplyr::summarise(lwr=median(loval)) |> as.data.frame() |> select(lwr)

# Plot predictions
png( "../figs/logistic_model_otb_85_norm.png", width = 7, height = 7, units = 'in', res = 600 )
par(mar=c(4.5,4,2,1))
Expand All @@ -124,10 +123,10 @@ dat <- inner_join( chl[,c('date','chl','chl_year_t1')],
polygon( x = c( plot_data$ratio, rev(plot_data$ratio) ),
y = c( plot_data$upr$upr, rev(plot_data$lwr$lwr) ),
col = rgb(0.3,0.7,1,0.2), border = rgb(0,0,0,0) )
abline( v = seq(0,2,0.1), col = rgb(0,0,0,0.1) )
abline( h = seq(0,1,0.1), col = rgb(0,0,0,0.1) )
lines( upr$upr ~ ratio, data = plot_data, col = rgb(0.3,0.7,1,1) )
lines( lwr$lwr ~ ratio, data = plot_data, col = rgb(0.3,0.7,1,1) )
abline( v = seq(0,2,0.1), col = rgb(0,0,0,0.1) )
abline( h = seq(0,1,0.1), col = rgb(0,0,0,0.1) )
lines( upr$upr ~ ratio, data = plot_data, col = rgb(0.3,0.7,1,1) )
lines( lwr$lwr ~ ratio, data = plot_data, col = rgb(0.3,0.7,1,1) )
# abline( v = 0 )
# Plot targets
deliv_targets <- c(1.08,1.71)
Expand All @@ -149,16 +148,16 @@ dat <- inner_join( chl[,c('date','chl','chl_year_t1')],
y0 = plot_data$lwr$lwr[which(plot_data$ratio==deliv_targets[i])],
lty = 2, col = seg_colors[i]
)
text( x = rep(text_pos[i],3),
y = c( plot_data$median[which(plot_data$ratio==deliv_targets[i])],
plot_data$upr$upr[which(plot_data$ratio==deliv_targets[i])],
plot_data$lwr$lwr[which(plot_data$ratio==deliv_targets[i])]
),
labels = round( c( 100*plot_data$median[which(plot_data$ratio==deliv_targets[i])],
100*plot_data$upr$upr[which(plot_data$ratio==deliv_targets[i])],
100*plot_data$lwr$lwr[which(plot_data$ratio==deliv_targets[i])]
),0 ), cex = 0.7, offset = 0.1, pos = 3, col = seg_colors[i]
)
# text( x = rep(text_pos[i],3),
# y = c( plot_data$median[which(plot_data$ratio==deliv_targets[i])],
# plot_data$upr$upr[which(plot_data$ratio==deliv_targets[i])],
# plot_data$lwr$lwr[which(plot_data$ratio==deliv_targets[i])]
# ),
# labels = round( c( 100*plot_data$median[which(plot_data$ratio==deliv_targets[i])],
# 100*plot_data$upr$upr[which(plot_data$ratio==deliv_targets[i])],
# 100*plot_data$lwr$lwr[which(plot_data$ratio==deliv_targets[i])]
# ),0 ), cex = 0.7, offset = 0.1, pos = 3, col = seg_colors[i]
# )
}
legend( 'topright', bty = 'n', lty = 2, col = seg_colors,
legend = paste0(deliv_targets," tons/Mm3")
Expand Down
Loading

0 comments on commit 10c1f5e

Please sign in to comment.