Skip to content

Commit

Permalink
couple small changes for multiple Ecov
Browse files Browse the repository at this point in the history
  • Loading branch information
brianstock-NOAA committed Dec 7, 2020
1 parent 60f1358 commit 9e46b48
Show file tree
Hide file tree
Showing 2 changed files with 11 additions and 9 deletions.
10 changes: 6 additions & 4 deletions R/prepare_wham_input.R
Expand Up @@ -677,7 +677,7 @@ without changing ASAP file, specify M$initial_means.")
data$Ecov_model <- sapply(ecov$process_model, match, c("rw", "ar1"))

# data$n_Ecov_pars <- c(1,3)[data$Ecov_model] # rw: 1 par (sig), ar1: 3 par (phi, sig)
if(!ecov$where %in% c('recruit','M')){
if(!any(ecov$where %in% c('recruit','M'))){
stop("Sorry, only ecov effects on recruitment and M currently implemented.
Set ecov$where = 'recruit' or 'M'.")
}
Expand Down Expand Up @@ -735,14 +735,14 @@ Ecov years: ", data$year1_Ecov, " to ", end_Ecov,"

if(data$Ecov_where[i] == 1){ # recruitment
cat(paste0("Ecov ",i,": ",ecov$label[i],"
",c('*NO*','Controlling','Limiting','Lethal','Masking','Directive')[data$Ecov_how+1]," (",ecov_str[[i]],") effect on: ", c('recruitment','M')[data$Ecov_where[i]],"
",c('*NO*','Controlling','Limiting','Lethal','Masking','Directive')[data$Ecov_how[i]+1]," (",ecov_str[[i]],") effect on: ", c('recruitment','M')[data$Ecov_where[i]],"
In model years:
"))
}
if(data$Ecov_where[i] == 2){ # M
cat(paste0("Ecov ",i,": ",ecov$label[i],"
",ecov_str[[i]]," effect on: ", c('recruitment','M')[data$Ecov_where[i]],"
",c('*NO*',ecov_str[[i]])[data$Ecov_how[i]+1]," effect on: ", c('recruitment','M')[data$Ecov_where[i]],"
In model years:
"))
Expand All @@ -753,8 +753,10 @@ lastyr <- tail(years,1)
cat(paste0("Lag: ",data$Ecov_lag[i],"
Ex: ",ecov$label[i]," in ",years[1]," affects ", c('recruitment','M')[data$Ecov_where[i]]," in ",years[1+data$Ecov_lag[i]],"
",ecov$label[i]," in ",lastyr," affects ", c('recruitment','M')[data$Ecov_where[i]]," in ",lastyr+data$Ecov_lag[i],"
"))
}
data$Ecov_label <- list(data$Ecov_label)
} # end load Ecov

# add vector of all observations for one step ahead residuals ==========================
Expand Down Expand Up @@ -1003,7 +1005,7 @@ Ex: ",ecov$label[i]," in ",years[1]," affects ", c('recruitment','M')[data$Ecov_
for(j in 1:data$n_Ecov){
tmp[,j] <- ifelse(data$Ecov_how[j]==0, NA, 0)
for(i in 1:max.poly){
if(i > data$Ecov_poly) tmp[i,j] = NA
if(i > data$Ecov_poly[j]) tmp[i,j] = NA
}
}
ind.notNA <- which(!is.na(tmp))
Expand Down
10 changes: 5 additions & 5 deletions R/wham_plots_tables.R
Expand Up @@ -202,7 +202,7 @@ plot.osa.residuals <- function(mod, do.tex=FALSE, do.png=FALSE, fontfam="", res=
abline(h=0, col=plot.colors[f], lwd=2)

# 2. trend vs. fitted val
plot(tmp$pred, tmp$residual, type='p', col=plot.colors[f], pch=19, xlab=paste0("Predicted ", mod$env$data$Ecov_label[f]), ylab="OSA Residuals",
plot(tmp$pred, tmp$residual, type='p', col=plot.colors[f], pch=19, xlab=paste0("Predicted ", mod$env$data$Ecov_label[[1]][f]), ylab="OSA Residuals",
ylim=ylims)
abline(h=0, col=plot.colors[f], lwd=2)

Expand Down Expand Up @@ -231,7 +231,7 @@ plot.osa.residuals <- function(mod, do.tex=FALSE, do.png=FALSE, fontfam="", res=
lines(z, upper, lty=2, col=plot.colors[f])
lines(z, lower, lty=2, col=plot.colors[f])

title (paste0("OSA residual diagnostics: Ecov ",f," (",mod$env$data$Ecov_label[f],")"), outer=T, line=-1)
title (paste0("OSA residual diagnostics: Ecov ",f," (",mod$env$data$Ecov_label[[1]][f],")"), outer=T, line=-1)
if(do.tex | do.png) dev.off() else par(origpar)
}
}
Expand Down Expand Up @@ -264,7 +264,7 @@ fit.summary.text.plot.fn <- function(mod){
text(5,nl <- nl-0.5, paste0("Index Age Comp Models: ", paste(acm[mod$env$data$age_comp_model_indices], collapse = ", ")))
text(5,nl <- nl-0.5,paste0("Recruitment model: ", recs[mod$env$data$recruit_model]))
if(!all(mod$env$data$Ecov_model == 0)){
for(ec in 1:mod$env$data$n_Ecov) text(5,nl <- nl-0.5, paste0("Environmental effect ", ec,": ", mod$env$data$Ecov_label[ec]," (",env.mod[mod$env$data$Ecov_model[ec]],") on ",env.where[mod$env$data$Ecov_where[ec]], " (", env.how[mod$env$data$Ecov_how[ec]],")"))
for(ec in 1:mod$env$data$n_Ecov) text(5,nl <- nl-0.5, paste0("Environmental effect ", ec,": ", mod$env$data$Ecov_label[[1]][ec]," (",env.mod[mod$env$data$Ecov_model[ec]],") on ",env.where[mod$env$data$Ecov_where[ec]], " (", env.how[mod$env$data$Ecov_how[ec]],")"))
} else {
text(5,nl <- nl-0.5, "Environmental effects: none")
}
Expand Down Expand Up @@ -612,7 +612,7 @@ plot.all.stdresids.fn = function(mod, do.tex = FALSE, do.png = FALSE, fontfam=""
}
xe$row = xe$Label
xe$Label = factor(xe$Label)
levels(xe$Label) = mod$env$data$Ecov_label
levels(xe$Label) = mod$env$data$Ecov_label[[1]]
xe$type = "Ecov"
}

Expand Down Expand Up @@ -723,7 +723,7 @@ plot.ecov.stdresids.fn = function(mod, years, do.tex = FALSE, do.png = FALSE, fo
x <- rbind(x, td)
}
x$Ecov = factor(x$Ecov)
levels(x$Ecov) = mod$env$data$Ecov_label
levels(x$Ecov) = mod$env$data$Ecov_label[[1]]
names(plot.colors) = levels(x$Ecov)
ggp = ggplot2::ggplot(x, ggplot2::aes(x=Year, y = Stdres, color=Ecov)) +
ggplot2::geom_line(size=1.1) +
Expand Down

0 comments on commit 9e46b48

Please sign in to comment.