Adding S3 objects

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: FishStatsUtils
 Type: Package
 Title: Utilities (shared code and data) for FishStats spatio-temporal modeling toolbox
-Version: 2.1.0
-Date: 2019-06-10
+Version: 2.2.0
+Date: 2019-07-25
 Authors@R: person("James","Thorson", email="[email protected]", role=c("aut","cre"))
 Maintainer: James Thorson <[email protected]>
 Description: FishStatsUtils contains utilities (shared code and data) used by multiple 
@@ -30,7 +30,6 @@ Imports:
     ThorsonUtilities,
     abind,
     corpcor,
-    TMBhelper,
     pander,
     formatR
 Depends:
@@ -39,8 +38,8 @@ Depends:
     R (>= 3.1.0)
 Suggests:
     testthat
-Remotes: james-thorson/utilities,
-    kaskr/TMB_contrib_R/TMBhelper
+Remotes: 
+    james-thorson/utilities
 License: GPL-3
 LazyData: yes
 BuildVignettes: yes

NAMESPACE

@@ -1,5 +1,12 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(plot,fit_model)
+S3method(plot,make_extrapolation_info)
+S3method(plot,make_spatial_info)
+S3method(print,fit_model)
+S3method(print,make_extrapolation_info)
+S3method(print,make_spatial_info)
+S3method(summary,fit_model)
 export(Calc_Anisotropic_Mesh)
 export(Calc_Kmeans)
 export(Calc_Polygon_Areas_and_Polygons_Fn)
@@ -30,6 +37,7 @@ export(Rotate_Fn)
 export(calc_cov)
 export(calculate_proportion)
 export(combine_extrapolation_info)
+export(combine_lists)
 export(convert_version_name)
 export(fit_model)
 export(format_covariates)
@@ -41,12 +49,14 @@ export(make_settings)
 export(make_spatial_info)
 export(map_hypervariance)
 export(match_strata_fn)
+export(plot)
 export(plot_anisotropy)
 export(plot_biomass_index)
 export(plot_cov)
 export(plot_data)
 export(plot_encounter_diagnostic)
 export(plot_factors)
+export(plot_index)
 export(plot_lines)
 export(plot_loadings)
 export(plot_maps)
@@ -55,4 +65,6 @@ export(plot_quantile_diagnostic)
 export(plot_range_index)
 export(plot_residuals)
 export(plot_results)
+export(plot_timeseries)
+export(rotate_factors)
 export(summarize_covariance)

R/PlotMap_Fn.R

@@ -20,7 +20,8 @@ PlotMap_Fn <-
 function(MappingDetails, Mat, PlotDF, MapSizeRatio=c('Width(in)'=4,'Height(in)'=4), Xlim, Ylim, FileName=paste0(getwd(),"/"), Year_Set,
          Rescale=FALSE, Rotate=0, Format="png", Res=200, zone=NA, Cex=0.01, textmargin="", add=FALSE, pch=15,
          outermargintext=c("Eastings","Northings"), zlim=NULL, Col=NULL,
-         Legend=list("use"=FALSE, "x"=c(10,30), "y"=c(10,30)), mfrow=c(1,1), plot_legend_fig=TRUE, land_color="grey", ignore.na=FALSE, ...){
+         Legend=list("use"=FALSE, "x"=c(10,30), "y"=c(10,30)), mfrow=c(1,1), plot_legend_fig=TRUE, land_color="grey", ignore.na=FALSE,
+         map_style="rescale", ...){
 
   # Check for problems
   if( length(Year_Set) != ncol(Mat) ){
@@ -67,12 +68,19 @@ function(MappingDetails, Mat, PlotDF, MapSizeRatio=c('Width(in)'=4,'Height(in)'=
       }else{
         # If not rotating:  Use simple plot
         if( Rotate==0 ){
-          Map = maps::map(MappingDetails[[1]], MappingDetails[[2]], plot=FALSE)
-          plot( 1, type="n", ylim=mean(Ylim)+c(-0.5,0.5)*diff(Ylim), xlim=mean(Xlim)+c(-0.5,0.5)*diff(Xlim), xaxt="n", yaxt="n", xlab="", ylab="" )
-          map( Map, add=TRUE )  # , col=land_color, fill=TRUE ->  Using land-fill color produces weird plotting artefacts
-          #map(MappingDetails[[1]], MappingDetails[[2]], ylim=mean(Ylim)+c(-0.5,0.5)*diff(Ylim), xlim=mean(Xlim)+c(-0.5,0.5)*diff(Xlim), plot=TRUE)
           Col_Bin = ceiling( f(Mat[Which,,drop=FALSE],zlim=zlim)[,tI]*(length(Col)-1) ) + 1
-          points(x=PlotDF[Which,'Lon'], y=PlotDF[Which,'Lat'], col=Col[Col_Bin], cex=Cex, pch=pch)
+          if( map_style=="rescale" ){
+            # Make plot size using plot(.)
+            plot( 1, type="n", ylim=mean(Ylim)+c(-0.5,0.5)*diff(Ylim), xlim=mean(Xlim)+c(-0.5,0.5)*diff(Xlim), xaxt="n", yaxt="n", xlab="", ylab="" )
+            points(x=PlotDF[Which,'Lon'], y=PlotDF[Which,'Lat'], col=Col[Col_Bin], cex=Cex, pch=pch)
+            Map = maps::map(MappingDetails[[1]], MappingDetails[[2]], plot=FALSE)
+            map( Map, add=TRUE ) #, col=land_color, fill=TRUE )  # ->  Using land-fill color produces weird plotting artefacts
+          }else{
+            # Make plot size using map(.)
+            map(MappingDetails[[1]], MappingDetails[[2]], ylim=mean(Ylim)+c(-0.5,0.5)*diff(Ylim), xlim=mean(Xlim)+c(-0.5,0.5)*diff(Xlim), plot=TRUE, fill=FALSE, mar=c(0.1,0.1,par("mar")[3],0), myborder=0 )
+            points(x=PlotDF[Which,'Lon'], y=PlotDF[Which,'Lat'], col=Col[Col_Bin], cex=Cex, pch=pch)
+            map(MappingDetails[[1]], MappingDetails[[2]], plot=TRUE, col=land_color, fill=TRUE, add=TRUE )
+          }
         }
         # If rotating:  Record all polygons; Rotate them and all points;  Plot rotated polygons;  Plot points
         if( Rotate!=0 ){
@@ -94,7 +102,7 @@ function(MappingDetails, Mat, PlotDF, MapSizeRatio=c('Width(in)'=4,'Height(in)'=
           tmpUTM_rotated <- maptools::elide( tmpUTM, rotate=Rotate)
           plot( 1, type="n", xlim=range(tmpUTM_rotated@coords[-c(1:nrow(Tmp1)),'x']), ylim=range(tmpUTM_rotated@coords[-c(1:nrow(Tmp1)),'y']), xaxt="n", yaxt="n" )
           Col_Bin = ceiling( f(tmpUTM_rotated@data[-c(1:nrow(Tmp1)),-c(1:2),drop=FALSE],zlim=zlim)[,tI]*(length(Col)-1) ) + 1
-          if( ignore.na==FALSE && any(Col_Bin<1 | Col_Bin>50) ) stop("zlim doesn't span the range of the variable")
+          if( ignore.na==FALSE && any(Col_Bin<1 | Col_Bin>length(Col)) ) stop("zlim doesn't span the range of the variable")
           points(x=tmpUTM_rotated@coords[-c(1:nrow(Tmp1)),'x'], y=tmpUTM_rotated@coords[-c(1:nrow(Tmp1)),'y'], col=Col[Col_Bin], cex=Cex, pch=pch)
           # Plot map features
           lev = levels(as.factor(tmpUTM_rotated@data$PID))

R/calculate_proportion.R

@@ -6,7 +6,7 @@
 #'
 #' @param Index output from \code{FishStatsUtils::plot_biomass_index}
 #' @inheritParams plot_biomass_index
-#' @param ... list of settings to pass to \code{sdreport}
+#' @param ... list of arguments to pass to \code{plot_index}
 #'
 #' @return Tagged list of output
 #' \describe{
@@ -15,8 +15,9 @@
 #' }
 #'
 #' @export
-calculate_proportion = function( TmbData, Index, Year_Set=NULL, Years2Include=NULL, strata_names=NULL, category_names=NULL, plot_legend=TRUE,
-  DirName=paste0(getwd(),"/"), PlotName="Proportion.png", interval_width=1, width=6, height=6, ... ){
+calculate_proportion = function( TmbData, Index, Year_Set=NULL, Years2Include=NULL, strata_names=NULL, category_names=NULL,
+  plot_legend=ifelse(TmbData$n_l>1,TRUE,FALSE), DirName=paste0(getwd(),"/"), PlotName="Proportion.png", PlotName2="Average.png",
+  interval_width=1, width=6, height=6, xlab="Category", ylab="Proportion", ... ){
 
   # Warnings and errors
   if( !all(TmbData[['FieldConfig']] %in% c(-2,-1)) ){
@@ -26,7 +27,7 @@ calculate_proportion = function( TmbData, Index, Year_Set=NULL, Years2Include=NU
   SE_Index_ctl = array(Index$Index_ctl[,,,'Std. Error'],dim=dim(Index$Index_ctl)[1:3])
 
   # Calculate proportions, and total biomass
-  Prop_ctl = Index_ctl / outer(rep(1,TmbData$n_c),apply(Index_ctl,MARGIN=2:3,FUN=sum))
+  Prop_ctl = Index_ctl / outer(rep(1,dim(Index_ctl)[1]),apply(Index_ctl,MARGIN=2:3,FUN=sum))
   Index_tl = apply(Index_ctl,MARGIN=2:3,FUN=sum)
   SE_Index_tl = sqrt(apply(SE_Index_ctl^2,MARGIN=2:3,FUN=sum,na.rm=TRUE))
 
@@ -46,30 +47,31 @@ calculate_proportion = function( TmbData, Index, Year_Set=NULL, Years2Include=NU
   # Median effective sample size across categories
   Neff_tl = apply(Neff_ctl, MARGIN=2:3, FUN=median, na.rm=TRUE)
 
-  # Fill in missing
-  if( is.null(Year_Set) ) Year_Set = 1:TmbData$n_t
-  if( is.null(Years2Include) ) Years2Include = 1:TmbData$n_t
-  if( is.null(strata_names) ) strata_names = 1:TmbData$n_l
-  if( is.null(category_names) ) category_names = 1:TmbData$n_c
+  # Plot
+  if( !is.na(PlotName) ){
+    plot_index( Index_ctl=Prop_ctl, sd_Index_ctl=sqrt(var_Prop_ctl), Year_Set=Year_Set, Years2Include=Years2Include,
+      strata_names=strata_names, category_names=category_names, plot_legend=plot_legend,
+      DirName=DirName, PlotName=PlotName, interval_width=interval_width, width=width, height=height,
+      xlab=xlab, ylab=ylab, scale="uniform", ... )
+  }
+
+  # Calculate weighted mean
+  sd_Mean_tl = Mean_tl = apply( Prop_ctl, MARGIN=2:3, FUN=function(vec){sum(vec*(1:length(vec)))} )
+  for( tI in 1:nrow(sd_Mean_tl)){
+  for( lI in 1:ncol(sd_Mean_tl)){
+    sd_Mean_tl[tI,lI] = sqrt(sum( var_Prop_ctl[,tI,lI] * (1:dim(var_Prop_ctl)[1] - Mean_tl[tI,lI])^2 ))
+  }}
 
   # Plot
-  Par = list( mar=c(2,2,1,0), mgp=c(2,0.5,0), tck=-0.02, yaxs="i", oma=c(2,2,0,0), mfrow=c(ceiling(sqrt(TmbData$n_t)),ceiling(TmbData$n_t/ceiling(sqrt(TmbData$n_t)))), ... )
-  png( file=paste0(DirName,"/",PlotName), width=width, height=height, res=200, units="in")
-    par( Par )
-    for( tI in 1:TmbData$n_t ){
-      # Calculate y-axis limits
-      Ylim = c(0, max(Prop_ctl[,tI,]%o%c(1,1) + sqrt(var_Prop_ctl[,tI,])%o%c(-interval_width,interval_width),na.rm=TRUE))
-      # Plot stuff
-      plot(1, type="n", xlim=range(category_names), ylim=1.05*Ylim, xlab="", ylab="", main=ifelse(TmbData$n_t>1,paste0("Year ",Year_Set[tI]),"") )
-      for(l in 1:TmbData$n_l){
-        FishStatsUtils::plot_lines( y=Prop_ctl[,tI,l], x=1:TmbData$n_c+seq(-0.1,0.1,length=TmbData$n_l)[l], ybounds=Prop_ctl[,tI,]%o%c(1,1) + sqrt(var_Prop_ctl[,tI,])%o%c(-interval_width,interval_width), type="b", col=rainbow(TmbData[['n_l']])[l], col_bounds=rainbow(TmbData[['n_l']])[l], ylim=Ylim)
-      }
-      if(plot_legend==TRUE) legend( "top", bty="n", fill=rainbow(TmbData[['n_l']]), legend=as.character(strata_names), ncol=2 )
-    }
-    mtext( side=1:2, text=c("Age","Proportion of biomass"), outer=TRUE, line=c(0,0) )
-  dev.off()
+  if( !is.na(PlotName2) ){
+    plot_index( Index_ctl=1%o%Mean_tl, sd_Index_ctl=1%o%sd_Mean_tl, Year_Set=Year_Set, Years2Include=Years2Include,
+      strata_names=strata_names, category_names=category_names, plot_legend=plot_legend,
+      DirName=DirName, PlotName=PlotName2, interval_width=interval_width, width=width, height=height,
+      xlab=xlab, ylab="Category", scale="uniform", Yrange=c(NA,NA), ... )     # , Yrange=c(1,dim(var_Prop_ctl)[1])
+  }
 
   # Return stuff
-  Return = list("Prop_ctl"=Prop_ctl, "Neff_tl"=Neff_tl, "var_Prop_ctl"=var_Prop_ctl, "Index_tl"=Index_tl, "Neff_ctl"=Neff_ctl)
+  Return = list("Prop_ctl"=Prop_ctl, "Neff_tl"=Neff_tl, "var_Prop_ctl"=var_Prop_ctl, "Index_tl"=Index_tl, "Neff_ctl"=Neff_ctl,
+    "Mean_tl"=Mean_tl, "sd_Mean_tl"=sd_Mean_tl )
   return( invisible(Return) )
 }

R/combine_lists.R

@@ -0,0 +1,23 @@
+
+#' Combine lists
+#'
+#' \code{combine_lists} combines two lists with precident given to one over the other
+#'
+#' @param default default values for combined list
+#' @param input preferred values for combined list
+#'
+#' @return combined list.
+#'
+#' @export
+combine_lists = function( default, input ){
+  output = default
+  for( i in seq_along(input) ){
+    if( names(input)[i] %in% names(default) ){
+      output[[names(input)[i]]] = input[[i]]
+    }else{
+      output = c( output, input[i] )
+    }
+  }
+  return( output )
+}
+

R/data.R

@@ -16,3 +16,21 @@
 #' @usage data(stream_network_eel_example)
 #' @keywords data
 NULL
+
+#' Data to demonstrate and test covariate effects
+#'
+#' Data sufficient to demonstrate spatially varying coefficient models
+#'
+#' \itemize{
+#'   \item sampling_data data-frame of biological sampling data and associated covariate measurements
+#'   \item Region region for model demo
+#'   \item strata.limits user-specified stratification of results
+#'   \item X_xtp 3D array of covariates, specified at knots given the use of \code{fine_scale=FALSE}
+#' }
+#'
+#' @name GOA_pcod_covariate_example
+#' @docType data
+#' @author Dave McGowan
+#' @usage data(GOA_pcod_covariate_example)
+#' @keywords data
+NULL