R/Module_plotModelFit.R
In MichaelFolkes/forecastR_package: Forecast stuff
Defines functions
plotModelFit
Documented in
plotModelFit
# This function does various diagnostics of the model fits
# (NOT the FC, just the fits used to generate the FC)
# - some are generic (e.g. obs vs fitted, residuals qqplot)
# - others are model specific (e.g. time series plot of Kalman Filter time-varying parameter)
#' Plot a model fit...
#'
#' @param fit.obj
#' @param options
#' @param fc.add
#'
#' @return
#' @export
#'
#' @examples
plotModelFit <- function(fit.obj, options= list(plot.which = "all",age.which="all",plot.add=FALSE),fc.add=NULL,tracing = FALSE){
# fit.obj is a list object generated by fitFC()
# currently, this object contains 1 list element for each age class,
# with the following required components:
# "model.type" "formula" "est.fn"
# "obs.values" "fitted.values" "data"
# run.yrs (run years for response variable)
# and additional components depending on the model
# "coefficients" "coefficients.table" "residuals"
# "sigma" "df" "fstatistic"
# "r.squared" "adj.r.squared"
# plot.which = "all" -> do all plots
# = one of "fitted_ts" "resid_ts" "resid_hist" "resid_qq" "fitvsobs" "residvsfitted" "modeldiagnostic" # TBI: "curvefit"
# = "precheck.report" -> do the first 4 listed in previous line
# NOTE: if plot.add = TRUE, then the plots are just added to the existing plotting device (for mix and match plots)
# if fc.add = NULL, then it has no effect
# = if it contains the ouput from calcFC(), then it includes a pt.fc or distr.fc
# (this only applies to the "fitted_ts" plot option)
# WARNING: plot panelling only works up to 9 age classes for now
# need to make sure that this is robust
if(options$age.which %in% c("all","total")){
if(any(grepl("Age",names(fit.obj)))){ ages.list <- names(fit.obj)[grepl("Age",names(fit.obj))]}
if(!any(grepl("Age",names(fit.obj)))){ ages.list <- "Total" }
}
if(!(options$age.which %in% c("all","total"))){ages.list <- options$age.which}
if(length(ages.list)==1){mfrow.use <- c(1,1)}
if(length(ages.list) > 1 & length(ages.list) <=4 ){par(mfrow.use <- c(2,2))}
if(length(ages.list) > 4 & length(ages.list) <=9 ){par(mfrow.use <- c(3,3))}
# -----------------------------------------------------
# Plot 1: Fitted vs. obs scatterplot
if(options$plot.which %in% c("all","fitvsobs")){
if(tracing){print("starting Plot 1: Fitted vs. obs scatterplot")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
par(pty="s")
# only do plot if have fitted values for that age class
if("fitted.values" %in% names(fit.obj[[age.plot]]) ){
# fading colors
x <- fit.obj[[age.plot]]$obs.values
bg.fade <- c(rep(0,max(0,length(x)-35)),seq(0.045,1, length.out=min(35,length(x))))
bg.col <- rgb(1, 0, 0,bg.fade)
plot.lims <- c(0, max(fit.obj[[age.plot]]$obs.values,fit.obj[[age.plot]]$fitted.values,na.rm=TRUE))
plot(fit.obj[[age.plot]]$obs.values,fit.obj[[age.plot]]$fitted.values, bty="n",xlab="Observed", ylab="Fitted",
xlim=plot.lims,ylim=plot.lims,col="darkblue",bg=bg.col,cex=1.4)
abline(0,1,col="black",lwd=1)
abline(lm(fit.obj[[age.plot]]$fitted.values ~ fit.obj[[age.plot]]$obs.values),col="red",lwd=2)
points(fit.obj[[age.plot]]$obs.values,fit.obj[[age.plot]]$fitted.values,pch=21,col="darkblue",bg=bg.col,cex=1.4)
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
perc.over <- round(100*sum(fit.obj[[age.plot]]$fitted.values > fit.obj[[age.plot]]$obs.values) / length(fit.obj[[age.plot]]$fitted.values))
text(0,par("usr")[4]*0.9,labels=paste(perc.over,"% Fitted > Obs"),col="red",cex=1,adj=0)
}
} # end looping through age classes
} # end if doing fit vs. obs
# -----------------------------------------------------
# Plot 1b: Residual vs Fitted scatterplot
if(options$plot.which %in% c("all","residvsfitted")){
if(tracing){print("starting Plot 1b: Residual vs Fitted scatterplot")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
# only do plot if have fitted values for that age class
if("fitted.values" %in% names(fit.obj[[age.plot]]) ){
# fading colors
x <- fit.obj[[age.plot]]$fitted.values
bg.fade <- c(rep(0,max(0,length(x)-35)),seq(0.045,1, length.out=min(35,length(x))))
bg.col <- rgb(1, 0, 0,bg.fade)
plot(fit.obj[[age.plot]]$fitted.values, fit.obj[[age.plot]]$residuals, bty="n",xlab="Fitted", ylab="Residuals",
col="darkblue",bg=bg.col,cex=1.4)
abline(h=0,col="darkblue")
# add fitted line only if resid vs. fit is NOT a straight vertical line (i.e. from ARIMA(0,0,0))
# see https://github.com/SOLV-Code/forecastR-ServerApp/issues/31
#print("-------")
min.fit <- min(fit.obj[[age.plot]]$fitted.values,na.rm = TRUE)
max.fit <- max(fit.obj[[age.plot]]$fitted.values,na.rm = TRUE)
#print(max.fit - min.fit)
#print(min.fit != max.fit)
if(round(min.fit) != round(max.fit)){
abline(lm(fit.obj[[age.plot]]$residuals ~ fit.obj[[age.plot]]$fitted.values),col="red",lwd=2)
}
points(fit.obj[[age.plot]]$fitted.values,fit.obj[[age.plot]]$residuals,pch=21,col="darkblue",bg=bg.col,cex=1.4)
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
}
} # end looping through age classes
} # end if doing residuals vs. fitted
# -----------------------------------------------------
# Plot 2: Fitted vs. obs time series
if(options$plot.which %in% c("all","fitted_ts","precheck.report")){
if(tracing){print("starting Plot 2: Fitted vs. obs time series")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
# only do plot if have fitted values for that age class
if("fitted.values" %in% names(fit.obj[[age.plot]]) ){
plot.lims <- c(0, max(fit.obj[[age.plot]]$obs.values,fit.obj[[age.plot]]$fitted.values,na.rm=TRUE))
xlim.use <- range(fit.obj[[age.plot]]$run.yrs)
if(!is.null(fc.add)){
fc.yr <- as.numeric(gsub("FC","",dimnames(fc.add$pt.fc)[[1]]))
fc.val <- fc.add$pt.fc[,age.plot]
fc.val.display <- prettyNum(round(fc.val),big.mark=",")
last.obs.idx <- length(fit.obj[[age.plot]]$obs.values)
plot.lims <- range(plot.lims, fc.val)
xlim.use <- range(xlim.use, fc.yr+2)
}
plot(fit.obj[[age.plot]]$run.yrs,fit.obj[[age.plot]]$obs.values , bty="n", type="o",cex=0.8,pch=19,col="darkblue",
xlab="Run Years", ylab="Abundance (Fix!)", ylim=plot.lims,xlim=xlim.use)
if(!is.null(fc.add)){
segments(fit.obj[[age.plot]]$run.yrs[last.obs.idx],fit.obj[[age.plot]]$fitted.values[last.obs.idx],
fc.yr,fc.val, col="red",lty=3)
}
lines(fit.obj[[age.plot]]$run.yrs,fit.obj[[age.plot]]$fitted.values,type="o",pch=22,col="red",bg="white",cex=0.8)
if(!is.null(fc.add)){
points(fc.yr,fc.val,pch = 22, col="red", bg="red",cex=1.3)
text(fc.yr+1.5,fc.val,labels=fc.val.display,adj=0,col="red",cex=0.8,xpd=NA)
}
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
if(age.plot==ages.list[1] & !options$plot.add ){
if(length(ages.list) == 1){leg.coords <- c( par("usr")[1] - ((par("usr")[2]-par("usr")[1])* 0.1),
par("usr")[3] - ((par("usr")[4]-par("usr")[3])* 0.06))
}
if(length(ages.list) > 1){leg.coords <- c( par("usr")[2] + ((par("usr")[2]-par("usr")[1])* 0.1),
par("usr")[3] - ((par("usr")[4]-par("usr")[3])* 0.06))
}
par(xpd=NA)
if(is.null(fc.add)){
legend( leg.coords[1],leg.coords[2] ,legend=c("Observed","Fitted"),pch=c(19,22),col=c("darkblue","red"),cex=0.8,bty="n")
}
if(!is.null(fc.add)){
legend( leg.coords[1],leg.coords[2] ,legend=c("Observed","Fitted",paste0("FC",fc.yr)),pch=c(19,22,22),col=c("darkblue","red","red"),
cex=0.8, pt.cex=c(0.8,0.8,1.3),pt.bg=c(NA,"white","red"),bty="n")
}
par(xpd=FALSE)
}
}
} # end looping through age classes
} # end if doing fitted time series
# -----------------------------------------------------
# Plot 3: curve fit (where applicable)
if(options$plot.which %in% c("all","curvefit")){
if(tracing){print("starting Plot 3: curve fit (where applicable)")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
plot(1:5,1:5, type="n", bty="n",xlab="",ylab="",axes=FALSE)
text(2.5,2.5, labels="placeholder for fitted curve plot")
} # end looping through age classes
} # end if doing curve plot
# -----------------------------------------------------
# Plot 4: residual time series
if(options$plot.which %in% c("all","resid_ts","precheck.report")){
if(tracing){print("starting Plot 4: residual time series")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
# only do plot if have residuals for that age class
if("residuals" %in% names(fit.obj[[age.plot]]) ){
ymax <- max(abs(fit.obj[[age.plot]]$residuals),na.rm=TRUE)
ylim.use <- c(-ymax,ymax)
plot(fit.obj[[age.plot]]$run.yrs,fit.obj[[age.plot]]$residuals , bty="n", type="o",cex=0.8,pch=22,col="red",bg="white",
xlab="Run Years", ylab="Raw Residuals (Model-Specific Units)",ylim=ylim.use)
abline(h=0,col="black",lty=2)
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
}
} # end looping through age classes
} # end if doing residual time series
# -----------------------------------------------------
# Plot 5: residual histogram
if(options$plot.which %in% c("all","resid_hist","precheck.report")){
if(tracing){print("starting Plot 5: residual histogram")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
# only do plot if have residuals for that age class
if("residuals" %in% names(fit.obj[[age.plot]]) ){
hist(fit.obj[[age.plot]]$residuals , xlab= "Raw Residuals (Model-Specific Units)",main="",
breaks=10,col="lightgrey",border="darkblue")
abline(v=0,col="red", lwd=2, lty=1)
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
}
} # end looping through age classes
} # end if doing residual time series
###########################################################
# Plot 6: residual qqplot
if(options$plot.which %in% c("all","resid_qq","precheck.report")){
if(tracing){print("starting Plot 6: residual qqplot")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
# only do plot if have residuals for that age class
if("residuals" %in% names(fit.obj[[age.plot]]) ){
qqnorm(fit.obj[[age.plot]]$residuals,main="",bty="n", pch=21,col="darkblue",bg="lightgrey")
qqline(fit.obj[[age.plot]]$residuals,col="red")
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
}
} # end looping through age classes
} # end if doing residual qqplot
# -----------------------------------------------------
# Plot 7: Model-Specific Diagnostic Plot
# THIS DIFFERS BY MODEL TYPE
# FOR NOW HAVE ONLY KALMAN FILTER SIBREG AND LOGPOWER SIBREG
if(options$plot.which %in% c("all","modeldiagnostic","precheck.report")){
if(tracing){print("starting Plot 7: Model-Specific Diagnostic Plot")}
if(!options$plot.add){par(mfrow=mfrow.use)}
# sib reg only apply if have more than 1 age class
ages.check <- length(ages.list) > 1
if(!ages.check){
# for now don't have any diagnostics yet, so just produce empy plot
plot(1:10,1:10,bty="n",type="n",xlab="",ylab="",axes = FALSE);text(5,5,"N/A")
} # end if have no age data
if(ages.check){
# blank plots for others (none yet)
if(fit.obj[[2]]$model.type!="SibRegKalman"){
for(age.plot in ages.list){
plot(1:10,1:10,bty="n",type="n",xlab="",ylab="",axes = FALSE);text(5,5,"N/A")
} # end looping through age classes
} # end if doing other than kalman
# Diagnostic plot for kalman filter: time varying parameter
# check if second youngest has kalman filter
if(fit.obj[[2]]$model.type=="SibRegKalman"){
for(age.plot in ages.list){
kf.have <- "smoothe.mean.a" %in% names(fit.obj[[age.plot]]$fit.obj[[2]])
if(!kf.have){ plot(1:10,1:10,bty="n",type="n",xlab="",ylab="",axes=FALSE);text(5,5,"N/A")}
if(kf.have ){
plot(fit.obj[[age.plot]]$run.yrs, fit.obj[[age.plot]]$fit.obj[[2]]$smoothe.mean.a, bty="n",xlab="Year", ylab="a parameter",
col="darkblue",type="l",cex=1.4)
abline(h=mean(fit.obj[[age.plot]]$fit.obj[[2]]$smoothe.mean.a),col="red")
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
} # end if have kf
} # end looping through age classes
title(main="Time Varying Intercept",outer=TRUE,line=-1.5)
} # end if doing SIbRegKalman
} # end if have age classes
} # end if doing modeldiagnostic
} # end plotModelFit()
MichaelFolkes/forecastR_package documentation built on April 4, 2021, 5:14 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plotModelFit
Documented in plotModelFit
# This function does various diagnostics of the model fits
# (NOT the FC, just the fits used to generate the FC)
# - some are generic (e.g. obs vs fitted, residuals qqplot)
# - others are model specific (e.g. time series plot of Kalman Filter time-varying parameter)
#' Plot a model fit...
#'
#' @param fit.obj
#' @param options
#' @param fc.add
#'
#' @return
#' @export
#'
#' @examples
plotModelFit <- function(fit.obj, options= list(plot.which = "all",age.which="all",plot.add=FALSE),fc.add=NULL,tracing = FALSE){
# fit.obj is a list object generated by fitFC()
# currently, this object contains 1 list element for each age class,
# with the following required components:
# "model.type" "formula" "est.fn"
# "obs.values" "fitted.values" "data"
# run.yrs (run years for response variable)
# and additional components depending on the model
# "coefficients" "coefficients.table" "residuals"
# "sigma" "df" "fstatistic"
# "r.squared" "adj.r.squared"
# plot.which = "all" -> do all plots
# = one of "fitted_ts" "resid_ts" "resid_hist" "resid_qq" "fitvsobs" "residvsfitted" "modeldiagnostic" # TBI: "curvefit"
# = "precheck.report" -> do the first 4 listed in previous line
# NOTE: if plot.add = TRUE, then the plots are just added to the existing plotting device (for mix and match plots)
# if fc.add = NULL, then it has no effect
# = if it contains the ouput from calcFC(), then it includes a pt.fc or distr.fc
# (this only applies to the "fitted_ts" plot option)
# WARNING: plot panelling only works up to 9 age classes for now
# need to make sure that this is robust
if(options$age.which %in% c("all","total")){
if(any(grepl("Age",names(fit.obj)))){ ages.list <- names(fit.obj)[grepl("Age",names(fit.obj))]}
if(!any(grepl("Age",names(fit.obj)))){ ages.list <- "Total" }
}
if(!(options$age.which %in% c("all","total"))){ages.list <- options$age.which}
if(length(ages.list)==1){mfrow.use <- c(1,1)}
if(length(ages.list) > 1 & length(ages.list) <=4 ){par(mfrow.use <- c(2,2))}
if(length(ages.list) > 4 & length(ages.list) <=9 ){par(mfrow.use <- c(3,3))}
# -----------------------------------------------------
# Plot 1: Fitted vs. obs scatterplot
if(options$plot.which %in% c("all","fitvsobs")){
if(tracing){print("starting Plot 1: Fitted vs. obs scatterplot")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
par(pty="s")
# only do plot if have fitted values for that age class
if("fitted.values" %in% names(fit.obj[[age.plot]]) ){
# fading colors
x <- fit.obj[[age.plot]]$obs.values
bg.fade <- c(rep(0,max(0,length(x)-35)),seq(0.045,1, length.out=min(35,length(x))))
bg.col <- rgb(1, 0, 0,bg.fade)
plot.lims <- c(0, max(fit.obj[[age.plot]]$obs.values,fit.obj[[age.plot]]$fitted.values,na.rm=TRUE))
plot(fit.obj[[age.plot]]$obs.values,fit.obj[[age.plot]]$fitted.values, bty="n",xlab="Observed", ylab="Fitted",
xlim=plot.lims,ylim=plot.lims,col="darkblue",bg=bg.col,cex=1.4)
abline(0,1,col="black",lwd=1)
abline(lm(fit.obj[[age.plot]]$fitted.values ~ fit.obj[[age.plot]]$obs.values),col="red",lwd=2)
points(fit.obj[[age.plot]]$obs.values,fit.obj[[age.plot]]$fitted.values,pch=21,col="darkblue",bg=bg.col,cex=1.4)
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
perc.over <- round(100*sum(fit.obj[[age.plot]]$fitted.values > fit.obj[[age.plot]]$obs.values) / length(fit.obj[[age.plot]]$fitted.values))
text(0,par("usr")[4]*0.9,labels=paste(perc.over,"% Fitted > Obs"),col="red",cex=1,adj=0)
}
} # end looping through age classes
} # end if doing fit vs. obs
# -----------------------------------------------------
# Plot 1b: Residual vs Fitted scatterplot
if(options$plot.which %in% c("all","residvsfitted")){
if(tracing){print("starting Plot 1b: Residual vs Fitted scatterplot")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
# only do plot if have fitted values for that age class
if("fitted.values" %in% names(fit.obj[[age.plot]]) ){
# fading colors
x <- fit.obj[[age.plot]]$fitted.values
bg.fade <- c(rep(0,max(0,length(x)-35)),seq(0.045,1, length.out=min(35,length(x))))
bg.col <- rgb(1, 0, 0,bg.fade)
plot(fit.obj[[age.plot]]$fitted.values, fit.obj[[age.plot]]$residuals, bty="n",xlab="Fitted", ylab="Residuals",
col="darkblue",bg=bg.col,cex=1.4)
abline(h=0,col="darkblue")
# add fitted line only if resid vs. fit is NOT a straight vertical line (i.e. from ARIMA(0,0,0))
# see https://github.com/SOLV-Code/forecastR-ServerApp/issues/31
#print("-------")
min.fit <- min(fit.obj[[age.plot]]$fitted.values,na.rm = TRUE)
max.fit <- max(fit.obj[[age.plot]]$fitted.values,na.rm = TRUE)
#print(max.fit - min.fit)
#print(min.fit != max.fit)
if(round(min.fit) != round(max.fit)){
abline(lm(fit.obj[[age.plot]]$residuals ~ fit.obj[[age.plot]]$fitted.values),col="red",lwd=2)
}
points(fit.obj[[age.plot]]$fitted.values,fit.obj[[age.plot]]$residuals,pch=21,col="darkblue",bg=bg.col,cex=1.4)
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
}
} # end looping through age classes
} # end if doing residuals vs. fitted
# -----------------------------------------------------
# Plot 2: Fitted vs. obs time series
if(options$plot.which %in% c("all","fitted_ts","precheck.report")){
if(tracing){print("starting Plot 2: Fitted vs. obs time series")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
# only do plot if have fitted values for that age class
if("fitted.values" %in% names(fit.obj[[age.plot]]) ){
plot.lims <- c(0, max(fit.obj[[age.plot]]$obs.values,fit.obj[[age.plot]]$fitted.values,na.rm=TRUE))
xlim.use <- range(fit.obj[[age.plot]]$run.yrs)
if(!is.null(fc.add)){
fc.yr <- as.numeric(gsub("FC","",dimnames(fc.add$pt.fc)[[1]]))
fc.val <- fc.add$pt.fc[,age.plot]
fc.val.display <- prettyNum(round(fc.val),big.mark=",")
last.obs.idx <- length(fit.obj[[age.plot]]$obs.values)
plot.lims <- range(plot.lims, fc.val)
xlim.use <- range(xlim.use, fc.yr+2)
}
plot(fit.obj[[age.plot]]$run.yrs,fit.obj[[age.plot]]$obs.values , bty="n", type="o",cex=0.8,pch=19,col="darkblue",
xlab="Run Years", ylab="Abundance (Fix!)", ylim=plot.lims,xlim=xlim.use)
if(!is.null(fc.add)){
segments(fit.obj[[age.plot]]$run.yrs[last.obs.idx],fit.obj[[age.plot]]$fitted.values[last.obs.idx],
fc.yr,fc.val, col="red",lty=3)
}
lines(fit.obj[[age.plot]]$run.yrs,fit.obj[[age.plot]]$fitted.values,type="o",pch=22,col="red",bg="white",cex=0.8)
if(!is.null(fc.add)){
points(fc.yr,fc.val,pch = 22, col="red", bg="red",cex=1.3)
text(fc.yr+1.5,fc.val,labels=fc.val.display,adj=0,col="red",cex=0.8,xpd=NA)
}
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
if(age.plot==ages.list[1] & !options$plot.add ){
if(length(ages.list) == 1){leg.coords <- c( par("usr")[1] - ((par("usr")[2]-par("usr")[1])* 0.1),
par("usr")[3] - ((par("usr")[4]-par("usr")[3])* 0.06))
}
if(length(ages.list) > 1){leg.coords <- c( par("usr")[2] + ((par("usr")[2]-par("usr")[1])* 0.1),
par("usr")[3] - ((par("usr")[4]-par("usr")[3])* 0.06))
}
par(xpd=NA)
if(is.null(fc.add)){
legend( leg.coords[1],leg.coords[2] ,legend=c("Observed","Fitted"),pch=c(19,22),col=c("darkblue","red"),cex=0.8,bty="n")
}
if(!is.null(fc.add)){
legend( leg.coords[1],leg.coords[2] ,legend=c("Observed","Fitted",paste0("FC",fc.yr)),pch=c(19,22,22),col=c("darkblue","red","red"),
cex=0.8, pt.cex=c(0.8,0.8,1.3),pt.bg=c(NA,"white","red"),bty="n")
}
par(xpd=FALSE)
}
}
} # end looping through age classes
} # end if doing fitted time series
# -----------------------------------------------------
# Plot 3: curve fit (where applicable)
if(options$plot.which %in% c("all","curvefit")){
if(tracing){print("starting Plot 3: curve fit (where applicable)")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
plot(1:5,1:5, type="n", bty="n",xlab="",ylab="",axes=FALSE)
text(2.5,2.5, labels="placeholder for fitted curve plot")
} # end looping through age classes
} # end if doing curve plot
# -----------------------------------------------------
# Plot 4: residual time series
if(options$plot.which %in% c("all","resid_ts","precheck.report")){
if(tracing){print("starting Plot 4: residual time series")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
# only do plot if have residuals for that age class
if("residuals" %in% names(fit.obj[[age.plot]]) ){
ymax <- max(abs(fit.obj[[age.plot]]$residuals),na.rm=TRUE)
ylim.use <- c(-ymax,ymax)
plot(fit.obj[[age.plot]]$run.yrs,fit.obj[[age.plot]]$residuals , bty="n", type="o",cex=0.8,pch=22,col="red",bg="white",
xlab="Run Years", ylab="Raw Residuals (Model-Specific Units)",ylim=ylim.use)
abline(h=0,col="black",lty=2)
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
}
} # end looping through age classes
} # end if doing residual time series
# -----------------------------------------------------
# Plot 5: residual histogram
if(options$plot.which %in% c("all","resid_hist","precheck.report")){
if(tracing){print("starting Plot 5: residual histogram")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
# only do plot if have residuals for that age class
if("residuals" %in% names(fit.obj[[age.plot]]) ){
hist(fit.obj[[age.plot]]$residuals , xlab= "Raw Residuals (Model-Specific Units)",main="",
breaks=10,col="lightgrey",border="darkblue")
abline(v=0,col="red", lwd=2, lty=1)
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
}
} # end looping through age classes
} # end if doing residual time series
###########################################################
# Plot 6: residual qqplot
if(options$plot.which %in% c("all","resid_qq","precheck.report")){
if(tracing){print("starting Plot 6: residual qqplot")}
if(!options$plot.add){par(mfrow=mfrow.use)}
for(age.plot in ages.list){
# only do plot if have residuals for that age class
if("residuals" %in% names(fit.obj[[age.plot]]) ){
qqnorm(fit.obj[[age.plot]]$residuals,main="",bty="n", pch=21,col="darkblue",bg="lightgrey")
qqline(fit.obj[[age.plot]]$residuals,col="red")
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
}
} # end looping through age classes
} # end if doing residual qqplot
# -----------------------------------------------------
# Plot 7: Model-Specific Diagnostic Plot
# THIS DIFFERS BY MODEL TYPE
# FOR NOW HAVE ONLY KALMAN FILTER SIBREG AND LOGPOWER SIBREG
if(options$plot.which %in% c("all","modeldiagnostic","precheck.report")){
if(tracing){print("starting Plot 7: Model-Specific Diagnostic Plot")}
if(!options$plot.add){par(mfrow=mfrow.use)}
# sib reg only apply if have more than 1 age class
ages.check <- length(ages.list) > 1
if(!ages.check){
# for now don't have any diagnostics yet, so just produce empy plot
plot(1:10,1:10,bty="n",type="n",xlab="",ylab="",axes = FALSE);text(5,5,"N/A")
} # end if have no age data
if(ages.check){
# blank plots for others (none yet)
if(fit.obj[[2]]$model.type!="SibRegKalman"){
for(age.plot in ages.list){
plot(1:10,1:10,bty="n",type="n",xlab="",ylab="",axes = FALSE);text(5,5,"N/A")
} # end looping through age classes
} # end if doing other than kalman
# Diagnostic plot for kalman filter: time varying parameter
# check if second youngest has kalman filter
if(fit.obj[[2]]$model.type=="SibRegKalman"){
for(age.plot in ages.list){
kf.have <- "smoothe.mean.a" %in% names(fit.obj[[age.plot]]$fit.obj[[2]])
if(!kf.have){ plot(1:10,1:10,bty="n",type="n",xlab="",ylab="",axes=FALSE);text(5,5,"N/A")}
if(kf.have ){
plot(fit.obj[[age.plot]]$run.yrs, fit.obj[[age.plot]]$fit.obj[[2]]$smoothe.mean.a, bty="n",xlab="Year", ylab="a parameter",
col="darkblue",type="l",cex=1.4)
abline(h=mean(fit.obj[[age.plot]]$fit.obj[[2]]$smoothe.mean.a),col="red")
title(main= paste(fit.obj[[age.plot]]$model.type,age.plot,sep=" - "))
} # end if have kf
} # end looping through age classes
title(main="Time Varying Intercept",outer=TRUE,line=-1.5)
} # end if doing SIbRegKalman
} # end if have age classes
} # end if doing modeldiagnostic
} # end plotModelFit()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.