fixing some plotting errors for asap3 files with multiple fleets and …

…age comp plots for fleets with no age comp
timjmiller · Nov 2, 2022 · 83f23ff · 83f23ff
1 parent 6a896f4
commit 83f23ff
Show file tree

Hide file tree

Showing 6 changed files with 122 additions and 118 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -25,7 +25,7 @@ wham 1.0.6 (2022-04-08)
 
 ### Major improvements
 
-* One-step-ahead prediction quantile residuals for age compositon observations for most log-likelihoods (Trijoulet et al. in review) [80e3bba](https://github.com/timjmiller/wham/commit/80e3bbae0244bd1199f018dca53f09b08f5fd203).
+* One-step-ahead prediction quantile residuals for age compositon observations for most log-likelihoods (Trijoulet et al. 2023) [80e3bba](https://github.com/timjmiller/wham/commit/80e3bbae0244bd1199f018dca53f09b08f5fd203).
 * Can specify a prior distribution on fully-selected catchability (logit-scale) which is then estimated as a random effect [fcbc604](https://github.com/timjmiller/wham/commit/fcbc604068005cd7cce6c3f329376c2b4ef7b540).
 * Auto-regressive random effects for fully-selected catchability [fcbc604](https://github.com/timjmiller/wham/commit/fcbc604068005cd7cce6c3f329376c2b4ef7b540).
 * Environmental covariate effects on fully-selected catchability [fcbc604](https://github.com/timjmiller/wham/commit/fcbc604068005cd7cce6c3f329376c2b4ef7b540).

diff --git a/R/compare_wham_models.R b/R/compare_wham_models.R
@@ -686,10 +686,7 @@ plot.selectivity.compare <- function(x, plot.opts, type="fleet"){
       cat("Fleet selectivity blocks not identical, cannot produce comparison plot.")
       return(NULL)
     } else {
-      selblocks <- as.numeric(unique(x[[1]]$selblock_pointer_fleets))
-      print(selblocks)
       selblocks <- unique(as.integer(x[[1]]$selblock_pointer_fleets))
-      print(selblocks)
       yrs <- lapply(selblocks, function(y){
                                   tmp <- which(x[[1]]$selblock_pointer_fleets == y);
                                   tmp <- tmp %% length(x[[1]]$years);
@@ -699,10 +696,7 @@ plot.selectivity.compare <- function(x, plot.opts, type="fleet"){
   }
   if(type == 'indices'){
     if(!allSame(lapply(x, function(y) unname(y$selblock_pointer_indices)))) stop("Index selectivity blocks not identical, cannot produce comparison plot")
-    selblocks <- as.numeric(unique(x[[1]]$selblock_pointer_indices))
-      print(selblocks)
     selblocks <- unique(as.integer(x[[1]]$selblock_pointer_indices))
-      print(selblocks)
     yrs <- lapply(selblocks, function(y){
                                 tmp <- which(x[[1]]$selblock_pointer_indices == y);
                                 tmp <- tmp %% length(x[[1]]$years);

diff --git a/R/make_osa_residuals.R b/R/make_osa_residuals.R
@@ -1,14 +1,14 @@
 #' Calculate one-step-ahead residuals
 #' 
 #' Standard residuals are not appropriate for models with random effects. Instead, one-step-ahead (OSA) residuals
-#' can be used for evaluating model goodness-of-fit (\href{https://link.springer.com/article/10.1007/s10651-017-0372-4}{Thygeson et al. (2017)},
+#' can be used for evaluating model goodness-of-fit (\href{https://doi.org/10.1007/s10651-017-0372-4}{Thygeson et al. (2017)},
 #' implemented in \code{\link[TMB:oneStepPredict]{TMB::oneStepPredict}}). OSA residual options
 #' are passed to \code{\link[TMB:oneStepPredict]{TMB::oneStepPredict}} in a list \code{osa.opts}. Current options are method: 
 #' oneStepGaussianOffMode (default), oneStepGaussian, or oneStepGeneric, and parallel: TRUE/FALSE. 
 #' See \code{\link[TMB:oneStepPredict]{TMB::oneStepPredict}} for further details.
 #' It is not recommended to run this function (or \code{\link[TMB:oneStepPredict]{TMB::oneStepPredict}}) with any random effects and
 #' mvtweedie age composition likelihoods due to extensive computational demand. An error will be thrown in such cases. 
-#' See Trijoulet et al. In review for OSA methods for age composition OSA residuals.
+#' See \href{https://doi.org/10.1016/j.fishres.2022.106487}{Trijoulet et al. (2023)} for OSA methods for age composition OSA residuals.
 #'
 #' @param model A fit WHAM model, output from \code{\link{fit_wham}}.
 #'
@@ -47,6 +47,7 @@ make_osa_residuals = function(model,osa.opts = list(method="oneStepGaussianOffMo
 
   cat("Doing OSA residuals...\n");
   input = model$input
+  if(class(input$data$obs) == "list") input$data$obs <- as.data.frame(input$data$obs) #simulated data$obs will be a list
   model$osa = input$data$obs
   model$osa$residual = NA
   #first do continuous obs, condition on obs without osa (probably none)

diff --git a/R/read_asap3_fit.R b/R/read_asap3_fit.R
@@ -87,7 +87,8 @@ read_asap3_fit <- function(wd, asap.name, pSPR=40)  {
   logFmult1  <- unlist(a1$asap.std[ which(a1$asap.std$name=="log_Fmult_year1")  ,3] )
   logFmult_devs <- a1$asap.std[which(a1$asap.std$name=="log_Fmult_devs")  ,3]
   log_F <- matrix(NA, nrow=nyears*n.fleet, ncol=2)
-  log_F[1,1] <- logFmult1
+  #year 1 for each fleet
+  log_F[1+nyears*(0:(n.fleet-1)),] <- logFmult1
 
   log.fmult.rows <- which(substr(a1$asap.cor.names, 1, 9)=="log_Fmult")  #this is nyears*n.fleet
   log.fmult.cor <- a1$asap.cor.mat[log.fmult.rows, log.fmult.rows]  #this has dimension nyears*n.fleet x nyears*n.fleet
@@ -96,13 +97,17 @@ read_asap3_fit <- function(wd, asap.name, pSPR=40)  {
   log.fmult.cov <- log.fmult.cor*(log.fmult.std %o% log.fmult.std) #creates var-cov matrix
   var.logFmult <- rep(NA, nyears*n.fleet)
   sigma.logFmult <- rep(NA, nyears*n.fleet)
-  var.logFmult[1] <- log.fmult.cov[1,1] #need to modify for multifleet
+  #year 1 for each fleet
+  var.logFmult[1+nyears*(0:(n.fleet-1))] <- log.fmult.cov[1+nyears*(0:(n.fleet-1)),1+nyears*(0:(n.fleet-1))] #need to modify for multifleet
+  #var.logFmult[1] <- log.fmult.cov[1,1] #need to modify for multifleet
 
+  for(f in 1:n.fleet) for (y in 2:nyears ) {
 
-  for (y in 2:nyears ) {
-
-    log_F[y,1] <- log_F[y-1] + logFmult_devs[y-1]
-    var.logFmult[y] <- var.logFmult[(y-1)] + a1$asap.std[log.fmult.rows[y],4]*a1$asap.std[log.fmult.rows[y],4]+2*(log.fmult.cov[y,(y-1)])
+    ind <- y + nyears*(0:(f-1))
+    log_F[ind,1] <- log_F[ind-1] + logFmult_devs[ind-1]
+    var.logFmult[ind] <- var.logFmult[ind-1] + a1$asap.std[log.fmult.rows[ind],4]*a1$asap.std[log.fmult.rows[ind],4]+2*(log.fmult.cov[ind,(ind-1)])
+    #log_F[y,1] <- log_F[y-1] + logFmult_devs[y-1]
+    #var.logFmult[y] <- var.logFmult[(y-1)] + a1$asap.std[log.fmult.rows[y],4]*a1$asap.std[log.fmult.rows[y],4]+2*(log.fmult.cov[y,(y-1)])
   }
 
   sigma.logFmult <- sqrt(var.logFmult)

diff --git a/R/wham_plots_tables.R b/R/wham_plots_tables.R
@@ -1395,58 +1395,60 @@ plot.catch.age.comp.resids <- function(mod, ages, ages.lab, scale.catch.bubble2
 	for (i in fleets)
 	{
     yind = which(dat$use_catch_paa[,i] ==1)
-    if(osa){
-      df = subset(mod$osa, type == "catchpaa")
-      my.title <- "Age Comp OSA Quantile Residuals for Fleet "
-      fname = paste0("Catch_age_comp_osa_resids_fleet_",i)
-      resids = matrix(NA, nrow = dat$n_years_model, ncol = dat$n_ages)
-      vals = resids
-      for(j in yind){
-        tmp = subset(df, year == j & fleet == paste0("fleet_",i))
-        resids[j,tmp$age] = tmp$residual
-        vals[j,tmp$age] = tmp$val
-        if(dat$age_comp_model_fleets[i] %in% c(1:2,10)) vals[j,tmp$age]/sum(vals[j,tmp$age]) #obs are numbers not proportions
-      }
+    if(length(yind)){
+      if(osa & "catchpaa" %in% mod$osa$type){
+        df = subset(mod$osa, type == "catchpaa")
+        my.title <- "Age Comp OSA Quantile Residuals for Fleet "
+        fname = paste0("Catch_age_comp_osa_resids_fleet_",i)
+        resids = matrix(NA, nrow = dat$n_years_model, ncol = dat$n_ages)
+        vals = resids
+        for(j in yind){
+          tmp = subset(df, year == j & fleet == paste0("fleet_",i))
+          resids[j,tmp$age] = tmp$residual
+          vals[j,tmp$age] = tmp$val
+          if(dat$age_comp_model_fleets[i] %in% c(1:2,10)) vals[j,tmp$age]/sum(vals[j,tmp$age]) #obs are numbers not proportions
+        }
 
-      scale.resid.bubble.catch <- 2
-    } else{
-      acomp.obs = dat$catch_paa[i,,]
-      acomp.pred = mod$rep$pred_catch_paa[1:length(years),i,]
-      #acomp.pred = aperm(mod$rep$pred_catch_paa[1:length(years),,,drop=FALSE], c(2,1,3))[i,,] #biomass is accounted for on the cpp side
-      #acomp.pred = acomp.pred/apply(acomp.pred,1,sum)
-      my.title <- "Age Comp Residuals (Observed-Predicted) for Fleet "
-      resids <- acomp.obs - acomp.pred  # NOTE obs-pred
-      resids[dat$use_catch_paa[,i]==0,] = NA # don't plot residuals for catch paa not fit in model
-      fname = paste0("Catch_age_comp_resids_fleet_",i)
-      scale.resid.bubble.catch <- 25
-    }
-    if(do.tex) cairo_pdf(file.path(od, paste0(fname,".pdf")), family = fontfam, height = 10, width = 10)
-    if(do.png) png(filename = file.path(od, paste0(fname,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
-    par(mar=c(4,4,2,2), oma=c(1,1,1,1), mfrow=c(1,1))
-    z1 <- resids
-    if(any(!is.na(resids))) range.resids<-range(abs((as.vector(z1))), na.rm=T)
-    else range.resids = c(0,0)
-
-    z3 <- z1 * scale.resid.bubble.catch * scale.catch.bubble2
-    resid.col <- matrix(NA, nrow=nyrs, ncol=n_ages)   # set color for residual bubbles
-    resid.col <- ifelse(z3 > 0.0, pos.resid.col, neg.resid.col)
-    plot(ages, rev(ages),  xlim = c(1, n_ages+1), ylim = c(years[nyrs],(years[1]-2)), xlab = "Age", ylab = tylab,
-      type = "n", axes=F)
-    axis(1, at= ages, lab=ages.lab)
-    axis(2, at = rev(years), lab = rev(years), cex.axis=0.75, las=1)
-    box()
-    abline(h=years, col="lightgray")
-    segments(x0=ages, y0=rep(years[1],n_ages), x1=ages, y1=rep(years[nyrs],n_ages), col = "lightgray", lty = 1)
+        scale.resid.bubble.catch <- 2
+      } else{
+        acomp.obs = dat$catch_paa[i,,]
+        acomp.pred = mod$rep$pred_catch_paa[1:length(years),i,]
+        #acomp.pred = aperm(mod$rep$pred_catch_paa[1:length(years),,,drop=FALSE], c(2,1,3))[i,,] #biomass is accounted for on the cpp side
+        #acomp.pred = acomp.pred/apply(acomp.pred,1,sum)
+        my.title <- "Age Comp Residuals (Observed-Predicted) for Fleet "
+        resids <- acomp.obs - acomp.pred  # NOTE obs-pred
+        resids[dat$use_catch_paa[,i]==0,] = NA # don't plot residuals for catch paa not fit in model
+        fname = paste0("Catch_age_comp_resids_fleet_",i)
+        scale.resid.bubble.catch <- 25
+      }
+      if(do.tex) cairo_pdf(file.path(od, paste0(fname,".pdf")), family = fontfam, height = 10, width = 10)
+      if(do.png) png(filename = file.path(od, paste0(fname,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
+      par(mar=c(4,4,2,2), oma=c(1,1,1,1), mfrow=c(1,1))
+      z1 <- resids
+      if(any(!is.na(resids))) range.resids<-range(abs((as.vector(z1))), na.rm=T)
+      else range.resids = c(0,0)
+
+      z3 <- z1 * scale.resid.bubble.catch * scale.catch.bubble2
+      resid.col <- matrix(NA, nrow=nyrs, ncol=n_ages)   # set color for residual bubbles
+      resid.col <- ifelse(z3 > 0.0, pos.resid.col, neg.resid.col)
+      plot(ages, rev(ages),  xlim = c(1, n_ages+1), ylim = c(years[nyrs],(years[1]-2)), xlab = "Age", ylab = tylab,
+        type = "n", axes=F)
+      axis(1, at= ages, lab=ages.lab)
+      axis(2, at = rev(years), lab = rev(years), cex.axis=0.75, las=1)
+      box()
+      abline(h=years, col="lightgray")
+      segments(x0=ages, y0=rep(years[1],n_ages), x1=ages, y1=rep(years[nyrs],n_ages), col = "lightgray", lty = 1)
 
-    for (j in 1:nyrs) points(ages, rep(years[j], n_ages), cex=abs(z3[j,]), col="black", bg = resid.col[j,],  pch = 21)
+      for (j in 1:nyrs) points(ages, rep(years[j], n_ages), cex=abs(z3[j,]), col="black", bg = resid.col[j,],  pch = 21)
 
-    bubble.legend1 <- round(quantile(abs(resids), probs = c(0.05,0.5,0.95), na.rm = TRUE),3)
-    bubble.legend2 <- bubble.legend1 * scale.resid.bubble.catch*scale.catch.bubble2
-    legend("topright", xpd=T, legend=bubble.legend1, pch=rep(1, 3), pt.cex=bubble.legend2, horiz=T , col='black')
-    legend("topleft", xpd=T, legend=c("Neg.", "Pos."), pch=rep(21, 2), pt.cex=3, horiz=T, pt.bg=c(neg.resid.col, pos.resid.col), col="black")
-    legend("top", xpd = TRUE, legend = paste("Max(resid)=",round(range.resids[2],2), sep=""), horiz = TRUE)
-    title (paste0(my.title,i), outer=T, line=-1)
-    if(do.tex | do.png) dev.off() else par(origpar)
+      bubble.legend1 <- round(quantile(abs(resids), probs = c(0.05,0.5,0.95), na.rm = TRUE),3)
+      bubble.legend2 <- bubble.legend1 * scale.resid.bubble.catch*scale.catch.bubble2
+      legend("topright", xpd=T, legend=bubble.legend1, pch=rep(1, 3), pt.cex=bubble.legend2, horiz=T , col='black')
+      legend("topleft", xpd=T, legend=c("Neg.", "Pos."), pch=rep(21, 2), pt.cex=3, horiz=T, pt.bg=c(neg.resid.col, pos.resid.col), col="black")
+      legend("top", xpd = TRUE, legend = paste("Max(resid)=",round(range.resids[2],2), sep=""), horiz = TRUE)
+      title (paste0(my.title,i), outer=T, line=-1)
+      if(do.tex | do.png) dev.off() else par(origpar)
+    } #end if any age comp for fleet
 	}   #end loop n_fleets
   # par(origpar)
 }
@@ -1467,60 +1469,62 @@ plot.index.age.comp.resids <- function(mod, ages, ages.lab, scale.catch.bubble2
 	for (i in indices)
 	{
     yind = which(dat$use_index_paa[,i] ==1)
-    if(osa){
-      df = subset(mod$osa, type == "indexpaa")
-      my.title <- "Age Comp OSA Quantile Residuals for Index "
-      fname = paste0("Catch_age_comp_osa_resids_index_",i)
-      resids = matrix(NA, nrow = dat$n_years_model, ncol = dat$n_ages)
-      vals = resids
-      for(j in yind){
-        tmp = subset(df, year == j & fleet == paste0("index_",i))
-        resids[j,tmp$age] = tmp$residual
-        vals[j,tmp$age] = tmp$val
-        if(dat$age_comp_model_indices[i] %in% c(1:2,10)) vals[j,tmp$age]/sum(vals[j,tmp$age]) #obs are numbers not proportions
+    if(length(yind)){
+      if(osa & "indexpaa" %in% mod$osa$type){
+        df = subset(mod$osa, type == "indexpaa")
+        my.title <- "Age Comp OSA Quantile Residuals for Index "
+        fname = paste0("Catch_age_comp_osa_resids_index_",i)
+        resids = matrix(NA, nrow = dat$n_years_model, ncol = dat$n_ages)
+        vals = resids
+        for(j in yind){
+          tmp = subset(df, year == j & fleet == paste0("index_",i))
+          resids[j,tmp$age] = tmp$residual
+          vals[j,tmp$age] = tmp$val
+          if(dat$age_comp_model_indices[i] %in% c(1:2,10)) vals[j,tmp$age]/sum(vals[j,tmp$age]) #obs are numbers not proportions
+        }
+        scale.resid.bubble.catch <- 2
+      } else {
+        acomp.obs = dat$index_paa[i,,]
+        acomp.pred = mod$rep$pred_index_paa[1:length(years),i,]
+        #acomp.pred = aperm(mod$rep$pred_IAA[1:length(years),,,drop=FALSE], c(2,1,3))[i,,] #biomass is accounted for on the cpp side
+        #acomp.pred = acomp.pred/apply(acomp.pred,1,sum)
+        my.title <- "Age Comp Residuals (Observed-Predicted) for Index "
+        resids <- acomp.obs - acomp.pred  # NOTE obs-pred
+        resids[dat$use_index_paa[,i]==0,] = NA # don't plot residuals for index paa not fit in model
+        fname = paste0("Catch_age_comp_resids_index_",i)
+        scale.resid.bubble.catch <- 25
       }
-      scale.resid.bubble.catch <- 2
-    } else{
-      acomp.obs = dat$index_paa[i,,]
-      acomp.pred = mod$rep$pred_index_paa[1:length(years),i,]
-      #acomp.pred = aperm(mod$rep$pred_IAA[1:length(years),,,drop=FALSE], c(2,1,3))[i,,] #biomass is accounted for on the cpp side
-      #acomp.pred = acomp.pred/apply(acomp.pred,1,sum)
-      my.title <- "Age Comp Residuals (Observed-Predicted) for Index "
-      resids <- acomp.obs - acomp.pred  # NOTE obs-pred
-      resids[dat$use_index_paa[,i]==0,] = NA # don't plot residuals for index paa not fit in model
-      fname = paste0("Catch_age_comp_resids_index_",i)
-      scale.resid.bubble.catch <- 25
-    }
-    tylab <- "Year"
-
-    if(do.tex) cairo_pdf(file.path(od, paste0(fname,".pdf")), family = fontfam, height = 10, width = 10)
-    if(do.png) png(filename = file.path(od, paste0(fname,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
-    par(mar=c(4,4,2,2), oma=c(1,1,1,1), mfrow=c(1,1))
-    z1 <- resids
-    if(any(!is.na(resids))) range.resids<-range(abs((as.vector(z1))), na.rm=T)
-    else range.resids <- c(0,0)
-
-    z3 <- z1 * scale.resid.bubble.catch * scale.catch.bubble2
-    resid.col <- matrix(NA, nrow=nyrs, ncol=n_ages)   # set color for residual bubbles
-    resid.col <- ifelse(z3 > 0.0, pos.resid.col, neg.resid.col)
-    plot(ages, rev(ages),  xlim = c(1, n_ages+1), ylim = c(years[nyrs],(years[1]-2)), xlab = "Age", ylab = tylab,
-      type = "n", axes=F)
-    axis(1, at= ages, lab=ages.lab)
-    axis(2, at = rev(years), lab = rev(years), cex.axis=0.75, las=1)
-    box()
-    abline(h=years, col="lightgray")
-    segments(x0=ages, y0=rep(years[1],n_ages), x1=ages, y1=rep(years[nyrs],n_ages), col = "lightgray", lty = 1)
-
-    for (j in 1:nyrs) if(dat$use_index_paa[j,i] == 1)
-      points(ages, rep(years[j], n_ages), cex=abs(z3[j,]), col="black", bg = resid.col[j,],  pch = 21)
-
-    bubble.legend1 <- round(quantile(abs(resids), probs = c(0.05,0.5,0.95), na.rm = TRUE),3)
-    bubble.legend2 <- bubble.legend1 * scale.resid.bubble.catch*scale.catch.bubble2
-    legend("topright", xpd=T, legend=bubble.legend1, pch=rep(1, 3), pt.cex=bubble.legend2, horiz=T , col='black')
-    legend("topleft", xpd=T, legend=c("Neg.", "Pos."), pch=rep(21, 2), pt.cex=3, horiz=T, pt.bg=c(neg.resid.col, pos.resid.col), col="black")
-    legend("top", xpd = TRUE, legend = paste("Max(resid)=",round(range.resids[2],2), sep=""), horiz = TRUE)
-    title (paste0(my.title,i), outer=T, line=-1)
-    if(do.tex | do.png) dev.off() else par(origpar)
+      tylab <- "Year"
+
+      if(do.tex) cairo_pdf(file.path(od, paste0(fname,".pdf")), family = fontfam, height = 10, width = 10)
+      if(do.png) png(filename = file.path(od, paste0(fname,'.png')), width = 10*144, height = 10*144, res = 144, pointsize = 12, family = fontfam)
+      par(mar=c(4,4,2,2), oma=c(1,1,1,1), mfrow=c(1,1))
+      z1 <- resids
+      if(any(!is.na(resids))) range.resids<-range(abs((as.vector(z1))), na.rm=T)
+      else range.resids <- c(0,0)
+
+      z3 <- z1 * scale.resid.bubble.catch * scale.catch.bubble2
+      resid.col <- matrix(NA, nrow=nyrs, ncol=n_ages)   # set color for residual bubbles
+      resid.col <- ifelse(z3 > 0.0, pos.resid.col, neg.resid.col)
+      plot(ages, rev(ages),  xlim = c(1, n_ages+1), ylim = c(years[nyrs],(years[1]-2)), xlab = "Age", ylab = tylab,
+        type = "n", axes=F)
+      axis(1, at= ages, lab=ages.lab)
+      axis(2, at = rev(years), lab = rev(years), cex.axis=0.75, las=1)
+      box()
+      abline(h=years, col="lightgray")
+      segments(x0=ages, y0=rep(years[1],n_ages), x1=ages, y1=rep(years[nyrs],n_ages), col = "lightgray", lty = 1)
+
+      for (j in 1:nyrs) if(dat$use_index_paa[j,i] == 1)
+        points(ages, rep(years[j], n_ages), cex=abs(z3[j,]), col="black", bg = resid.col[j,],  pch = 21)
+
+      bubble.legend1 <- round(quantile(abs(resids), probs = c(0.05,0.5,0.95), na.rm = TRUE),3)
+      bubble.legend2 <- bubble.legend1 * scale.resid.bubble.catch*scale.catch.bubble2
+      legend("topright", xpd=T, legend=bubble.legend1, pch=rep(1, 3), pt.cex=bubble.legend2, horiz=T , col='black')
+      legend("topleft", xpd=T, legend=c("Neg.", "Pos."), pch=rep(21, 2), pt.cex=3, horiz=T, pt.bg=c(neg.resid.col, pos.resid.col), col="black")
+      legend("top", xpd = TRUE, legend = paste("Max(resid)=",round(range.resids[2],2), sep=""), horiz = TRUE)
+      title (paste0(my.title,i), outer=T, line=-1)
+      if(do.tex | do.png) dev.off() else par(origpar)
+    } #end if any age comp for index
 	}   #end loop n_fleets
   # par(origpar)
 }

diff --git a/src/age_comp_osa.hpp b/src/age_comp_osa.hpp
@@ -1,4 +1,4 @@
-//from Trijoulet et al. (in review). OSA residuals for composition data
+//from Trijoulet et al. (2023). OSA residuals for composition data
 //osa residual versions of dmultinom, ddirmultinom, dlogisticnormal
 //see repo at