Fix another bug in make_covariates(.) ...

... which failed with some columns orders for `covariate_data`

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: FishStatsUtils
 Type: Package
 Title: Utilities (shared code and data) for FishStats spatio-temporal modeling toolbox
-Version: 2.3.0
-Date: 2019-09-18
+Version: 2.3.2
+Date: 2019-10-10
 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 

R/load_example.R

@@ -122,7 +122,7 @@ load_example = function( data_set="EBS_pollock" ){
     strata.limits = data.frame('STRATA'="All_areas")
   }
   if( tolower(data_set) %in% tolower(c("covariate_example","GOA_pcod_covariate_example")) ){
-    data( GOA_Pcod_covariate_example, package="FishStatsUtils" )
+    data( GOA_pcod_covariate_example, package="FishStatsUtils" )
     sampling_data = GOA_Pcod_covariate_example$sampling_data
     strata.limits = data.frame('STRATA'="All_areas")
     X_xtp = GOA_Pcod_covariate_example$X_xtp

R/make_covariates.R

@@ -19,13 +19,14 @@
 make_covariates = function( formula=~0, covariate_data, Year_i, spatial_list, extrapolation_list ){
 
   # Errors
+  if( !is.data.frame(covariate_data) ) stop("Please ensure that `covariate_data` is a data frame")
   if( !all(c("Lat","Lon","Year") %in% names(covariate_data)) ){
     stop( "`data` in `make_covariates(.)` must include columns `Lat`, `Lon`, and `Year`" )
   }
-  if( !is.data.frame(covariate_data) ) stop("Please ensure that `covariate_data` is a data frame")
 
-  #
+  # transform data inputs
   sample_data = data.frame( "Year"=Year_i, "Lat"=spatial_list$latlon_i[,'Lat'], "Lon"=spatial_list$latlon_i[,'Lon'] )
+  covariate_names = setdiff( names(covariate_data), names(sample_data) )
 
   # set of years needed
   Year_Set = min(Year_i):max(Year_i)
@@ -35,15 +36,14 @@ make_covariates = function( formula=~0, covariate_data, Year_i, spatial_list, ex
 
   # Create data frame of necessary size
   DF_zp = NULL
-  DF_ip = cbind( sample_data, covariate_data[rep(1,nrow(sample_data)),-match(names(sample_data),colnames(covariate_data))] )
-  DF_ip[,-match(c("Year","Lat","Lon"),colnames(covariate_data))] = NA
+  DF_ip = cbind( sample_data, covariate_data[rep(1,nrow(sample_data)),covariate_names] )
+  DF_ip[,covariate_names] = NA
 
   # Loop through data and extrapolation-grid
   for( tI in seq_along(Year_Set) ){
-  #for( tI in 1:14 ){
 
     # Subset to same year
-    tmp_covariate_data = covariate_data[ which(Year_Set[tI]==covariate_data[,'Year'] | is.na(covariate_data[,'Year'])), ]
+    tmp_covariate_data = covariate_data[ which(Year_Set[tI]==covariate_data[,'Year'] | is.na(covariate_data[,'Year'])), , drop=FALSE]
     if( nrow(tmp_covariate_data)==0 ){
       stop("Year ", Year_Set[tI], " not found in `covariate_data` please specify covariate values for all years" )
     }
@@ -53,17 +53,18 @@ make_covariates = function( formula=~0, covariate_data, Year_i, spatial_list, ex
     if( length(Which) > 0 ){
       NN = RANN::nn2( data=tmp_covariate_data[,c("Lat","Lon")], query=sample_data[Which,c("Lat","Lon")], k=1 )
       # Add to data-frame
-      newcolumns = tmp_covariate_data[ NN$nn.idx[,1], -match(c("Lat","Lon","Year"),colnames(tmp_covariate_data)), drop=FALSE ]
-      DF_ip[Which, -match(c("Lat","Lon","Year"),colnames(covariate_data))] = newcolumns
+      nearest_covariates = tmp_covariate_data[ NN$nn.idx[,1], covariate_names, drop=FALSE ]
+      DF_ip[Which, covariate_names] = nearest_covariates
     }
 
     # Do nearest neighbors to define covariates for extrapolation grid, including years without observations
     NN = RANN::nn2( data=tmp_covariate_data[,c("Lat","Lon")], query=latlon_g[,c("Lat","Lon")], k=1 )
     # Add rows
-    newcolumns = tmp_covariate_data[ NN$nn.idx[,1], -match(names(sample_data),colnames(tmp_covariate_data)), drop=FALSE ]
-    newrows = cbind("Year"=Year_Set[tI], latlon_g, newcolumns )
+    nearest_covariates = tmp_covariate_data[ NN$nn.idx[,1], covariate_names, drop=FALSE ]
+    newrows = cbind("Year"=Year_Set[tI], latlon_g, nearest_covariates )
     DF_zp = rbind( DF_zp, newrows )
   }
+  if( any(is.na(DF_ip)) ) stop("Problem with `DF_ip` in `make_covariates(.)")
 
   # Convert to dimensions requested
   DF = rbind( DF_ip, DF_zp )
@@ -89,7 +90,7 @@ make_covariates = function( formula=~0, covariate_data, Year_i, spatial_list, ex
   }
 
   # return stuff
-  Return = list( "X_gtp"=X_gtp, "X_itp"=X_itp )
+  Return = list( "X_gtp"=X_gtp, "X_itp"=X_itp, "covariate_names"=covariate_names )
   return( Return )
 }
 

R/plot_maps.r

@@ -107,7 +107,8 @@ function(plot_set=3, Report, PlotDF, Sdreport=NULL, TmbData=NULL, projargs='+pro
   for(plot_num in plot_set){
 
     # Extract elements
-    plot_code <- c("encounter_prob", "pos_catch", "density", "", "", "epsilon_1", "epsilon_2", "linear_predictor_1", "linear_predictor_2", "density_CV", "covariates", "total_density", "covariate_effects_1", "covariate_effects_2")[plot_num]
+    Array_xct = NULL
+    plot_code <- c("encounter_prob", "pos_catch", "density", "", "", "epsilon_1", "epsilon_2", "linear_predictor_1", "linear_predictor_2", "density_CV", "covariates", "total_density", "covariate_effects_1", "covariate_effects_2", "omega_1", "omega_2")[plot_num]
     #if( missing(textmargin) ){
     #  textmargin <- c("Probability of encounter", "Density, ln(kg. per square km.)", "Density, ln(kg. per square km.)", "", "", "", "", "", "", "CV of density (dimensionless)", "Covariate value", "Density, ln(kg. per square km.)", "", "")[plot_num]
     #}
@@ -155,7 +156,7 @@ function(plot_set=3, Report, PlotDF, Sdreport=NULL, TmbData=NULL, projargs='+pro
     }
     if(plot_num==6){
       # Epsilon for presence/absence ("Eps_Pres")
-      if( quiet==FALSE ) message(" # Plotting spatio-temporal variation in 1st linear predictor")
+      if( quiet==FALSE ) message(" # Plotting spatio-temporal effects (Epsilon) in 1st linear predictor")
       if("D_xt"%in%names(Report)) Array_xct = Report$Epsilon1_st
       if("D_xct"%in%names(Report)) Array_xct = Report$Epsilon1_sct
       if("D_xcy"%in%names(Report)) Array_xct = Report$Epsilon1_sct
@@ -164,7 +165,7 @@ function(plot_set=3, Report, PlotDF, Sdreport=NULL, TmbData=NULL, projargs='+pro
     }
     if(plot_num==7){
       # Epsilon for positive values ("Eps_Pos")
-      if( quiet==FALSE ) message(" # Plotting spatio-temporal variation in 2nd linear predictor")
+      if( quiet==FALSE ) message(" # Plotting spatio-temporal effects (Epsilon) in 2nd linear predictor")
       if("D_xt"%in%names(Report)) Array_xct = Report$Epsilon2_st
       if("D_xct"%in%names(Report)) Array_xct = Report$Epsilon2_sct
       if("D_xcy"%in%names(Report)) Array_xct = Report$Epsilon2_sct
@@ -220,7 +221,7 @@ function(plot_set=3, Report, PlotDF, Sdreport=NULL, TmbData=NULL, projargs='+pro
     }
     if(plot_num==12){
       # Total density ("Dens")
-      if( quiet==FALSE ) message(" # Plotting covariate effects for 1st linear predictor")
+      if( quiet==FALSE ) message(" # Plotting total density")
       if("D_xt"%in%names(Report)) Array_xct = log(Report$D_xt)
       if("D_xct"%in%names(Report)) Array_xct = log(apply(Report$D_xct, FUN=sum, MARGIN=c(1,3)))
       if("D_xcy"%in%names(Report)) Array_xct = log(apply(Report$D_xcy, FUN=sum, MARGIN=c(1,3)))
@@ -249,14 +250,35 @@ function(plot_set=3, Report, PlotDF, Sdreport=NULL, TmbData=NULL, projargs='+pro
       if("dhat_ktp" %in% names(Report)) stop()
       if("dpred_ktp" %in% names(Report)) stop()
     }
+    #if(plot_num==15){
+    #  # Spatial effects for probability of encounter
+    #  if( quiet==FALSE ) message(" # Plotting spatial effects (Omega) for 1st linear predictor")
+    #  if("D_xt"%in%names(Report)) stop()
+    #  if("D_xct"%in%names(Report)) stop()
+    #  if("D_xcy"%in%names(Report)) Array_xct = Report$Omega1_sc %o% 1
+    #  if("D_gcy"%in%names(Report)) Array_xct = Report$Omega1_gc %o% 1
+    #  if("dhat_ktp" %in% names(Report)) stop()
+    #  if("dpred_ktp" %in% names(Report)) stop()
+    #}
+    #if(plot_num==16){
+    #  # Spatial effects for positive catch rates
+    #  if( quiet==FALSE ) message(" # Plotting spatial effects (Omega) for 2nd linear predictor")
+    #  if("D_xt"%in%names(Report)) stop()
+    #  if("D_xct"%in%names(Report)) stop()
+    #  if("D_xcy"%in%names(Report)) Array_xct = Report$Omega2_sc %o% 1
+    #  if("D_gcy"%in%names(Report)) Array_xct = Report$Omega2_gc %o% 1
+    #  if("dhat_ktp" %in% names(Report)) stop()
+    #  if("dpred_ktp" %in% names(Report)) stop()
+    #}
+    if( is.null(Array_xct)) stop("Problem with `plot_num` in `plot_maps(.)")
 
     # Plot for each category
     if( tolower(Panel)=="category" ){
       if(length(dim(Array_xct))==2) Nplot = 1
       if(length(dim(Array_xct))==3) Nplot = dim(Array_xct)[2]
       for( cI in 1:Nplot){
         if(length(dim(Array_xct))==2) Return = Mat_xt = Array_xct
-        if(length(dim(Array_xct))==3) Return = Mat_xt = Array_xct[,cI,]
+        if(length(dim(Array_xct))==3) Return = Mat_xt = array(as.vector(Array_xct[,cI,]),dim=dim(Array_xct)[c(1,3)])
 
         # Do plot
         #if( is.null(mfrow)) mfrow = c(ceiling(sqrt(length(Years2Include))), ceiling(length(Years2Include)/ceiling(sqrt(length(Years2Include)))))

R/plot_results.R

@@ -85,7 +85,7 @@ plot_results = function( fit, settings=fit$settings, plot_set=3, working_dir=pas
   # Plot densities
   message("\n### Making plots of spatial predictions")
   plot_maps_args = list(...)
-  plot_maps_args = combine_lists( input=plot_maps_args, default=list(plot_set=plot_set, category_names=category_names,
+  plot_maps_args = combine_lists( input=plot_maps_args, default=list(plot_set=plot_set, category_names=category_names, TmbData=fit$data_list,
     Report=fit$Report, Sdreport=fit$parameter_estimates$SD, PlotDF=map_list[["PlotDF"]], MapSizeRatio=map_list[["MapSizeRatio"]],
     working_dir=working_dir, Year_Set=year_labels, Years2Include=years_to_plot, legend_x=map_list[["Legend"]]$x/100, legend_y=map_list[["Legend"]]$y/100) )
   Dens_xt = do.call( what=plot_maps, args=plot_maps_args )

R/plot_variable.R

@@ -140,7 +140,7 @@ function( Y_gt, map_list, panel_labels, projargs='+proj=longlat', map_resolution
     yt = (1-legend_y[2])*par('usr')[3] + (legend_y[2])*par('usr')[4]
     if( diff(legend_y) > diff(legend_x) ){
       align = c("lt","rb")[2]
-      gradient = c("x","y")[1]
+      gradient = c("x","y")[2]
     }else{
       align = c("lt","rb")[1]
       gradient = c("x","y")[1]