James-Thorson-NOAA / VAST

Spatio-temporal analysis of univariate or multivariate data, e.g., standardizing data for multiple species or stages
http://www.FishStats.org
GNU General Public License v3.0
124 stars 54 forks source link

Summarize() function to output results dataframe not working #211

Closed afredston closed 4 years ago

afredston commented 4 years ago

I've been trying to figure out how to extract a results dataframe from a fitted VAST single-species model, with columns for knot position (northings and eastings), year, and density. I can't figure out where this actually lives in the fitted model (D_gcy?) or how to translate it into a dataframe. It looks like there were past efforts to create a summary function for this purpose (#2), but Summarize() isn't working in VAST or FishStatsUtils anymore. Any ideas?

aallyn commented 4 years ago

Hi Alexa, by no means a seasoned VAST user and not sure if this is the best way. Working off of the previous "Summarize" function and given an object returned by the VAST::fit_model function, this seems to work for me...

Summarize_VAST<- function(fit.model){

  # For debugging
  if(FALSE) {
    fit.model<- fit.basicfore
  }

  # Get object
  obj<- fit.model$tmb_list[["Obj"]]

  # Extract
  report<- obj$report()
  parhat<- obj$env$parList()

  # Slots worth describing -- density at each knot? Can project this after to get east/northing, or could pull east/northing out of the fit.model$extrapolation_list$Data_Extrap object...
  dens.df<- data.frame("Lon" = fit.model$spatial_list$latlon_g[,"Lon"],
                       "Lat" = fit.model$spatial_list$latlon_g[,"Lat"],
                       "Density" = report$D_gcy)

  # Bundle and return
  Return<- list("parhat"=parhat, "PredDens" = dens.df)
  return(Return)
}

vastmod.summary<- Summarize_VAST(your.model)
afredston commented 4 years ago

This is great, thanks! I've edited your function down into a shorter version to just return the density data in a "long" format. Sharing here in case it's useful for future users. @aallyn what does the first part of your function ("for debugging") do?

get_density <- function(fit.model){

  # Get object
  obj<- fit.model$tmb_list[["Obj"]]

  # Extract
  report<- obj$report()

  # Slots worth describing -- density at each knot? Can project this after to get east/northing, or could pull east/northing out of the fit.model$extrapolation_list$Data_Extrap object...
  dens.df<- data.frame("lon" = fit.model$spatial_list$latlon_g[,"Lon"],
                       "lat" = fit.model$spatial_list$latlon_g[,"Lat"],
                       "density" = report$D_gcy)

  dens.df <- reshape2::melt(dens.df, id.vars=c('lon','lat'))
  colnames(dens.df) <- c('lat','lon','year',
                     'density')

  dens.df$year <- gsub("density.","", dens.df$year)

  return(dens.df)
}
aallyn commented 4 years ago

That's awesome! The "for dugging" section is just to help walk through the function. You can manually run the lines within the "if(FALSE){}" section, with those lines corresponding to the function arguments which supply everything the function needs to run. When you call the function, though, what's inside the "if(FALSE)" isn't executed. Probably a more elegant way of doing it