Effect.fit_model new error on latest version of VAST

[smormede]

I've uploaded VAST package 3.10.0. I now get the following error when I try to plot partial effects plots with code that used to work.

require(VAST) require(effects)

run a VAST model which works fine

{...}

plot partial effects

pred <- Effect.fit_model( mod=fit,focal.predictors = c("method"), which_formula = "Q1",xlevels = 100, pad_values=1)

please read ?Effect.fit_model for details Error in eval(substitute(expr), data, enclos = parent.frame()) : object 'prior.weights' not found

I cannot replicate the same error with the code in the wiki: https://github.com/James-Thorson-NOAA/VAST/wiki/Visualize-covariate-response (which still works although it now has a warning message).

Anyone has had this issue?

thanks

Sophie

[James-Thorson-NOAA]

Sophie,

I've found that the effects package is pretty unstable and I've been slowly working out how to switch to using marginaleffects. Do you wanna help me work through the switch-over?

I think I've got it worked out, but would appreciate the chance to do some back-and-forth on the dev branch (or just sourcing in the relevant functions)

[smormede]

Hi James,

Happy to help where I can, in particular testing etc. Maybe send me an email and we can look at it.

Sophie

[Charles-Adams-NOAA]

Has this been resolved? I'm getting a slightly different error:

pred = Effect.fit_model( fit,

•                      focal.predictors = c("bottemp"),
  

•                      which_formula = "X1",
  

•                      xlevels = 100,
  

•                      transformation = list(link=identity, inverse=identity) )
  

please read ?Effect.fit_model for details Error in Effect.fit_model(fit, focal.predictors = c("bottemp"), which_formula = "X1", : Please load covariate_data_full and catchability_data_full into global memory

[smormede]

You just need to add these two lines before your call (see wiki for more details)

covariate_data_full = fit$effects$covariate_data_full catchability_data_full = fit$effects$catchability_data_full

Sophie

[Charles-Adams-NOAA]

Already did that. And loaded effects package

[James-Thorson-NOAA]

effects package has been annoying, so I added a new option using `marginaleffects package. Do you want to explore that instead?

On Tue, Feb 21, 2023 at 12:34 PM Charles-Adams-NOAA < @.***> wrote:

Already did that. And loaded effects package

— Reply to this email directly, view it on GitHub https://github.com/James-Thorson-NOAA/VAST/issues/359#issuecomment-1439057225, or unsubscribe https://github.com/notifications/unsubscribe-auth/AL62VMSLRXWFPFOS7RDENM3WYUREJANCNFSM6AAAAAAT5LDKSY . You are receiving this because you commented.Message ID: @.***>

[Charles-Adams-NOAA]

Happy to give it a try Jim!

[James-Thorson-NOAA]

Could you try out this comparison of effects and marginaleffects?

I tried making it a gory example, with multivariate, different predictors, mapped off stuff, etc, so you might explore a bit and I'd love to hear if it works for you!

# Load packages
library(VAST)
library(splines)

# load data set
example = load_example( data_set="EBS_pollock" )
covariate_data = data.frame( "Lat"=0, "Lon"=0, "Year"=example$covariate_data[,'Year'],
  "CPE"=(example$covariate_data[,'AREA_SUM_KM2_LTE2']-100000)/100000)

# Add species
example$sampling_data = rbind(
  cbind( example$sampling_data, "Species"=1 ),
  cbind( example$sampling_data, "Species"=2 )
)
example$sampling_data$Catch_KG = ifelse( example$sampling_data$Species==1, example$sampling_data$Catch_KG^1.1, example$sampling_data$Catch_KG )

# Make settings (turning off bias.correct to save time for example)
settings = make_settings( n_x=100,
  Region=example$Region,
  purpose="index2",
  bias.correct = TRUE )

# Change Beta1 to AR1, to allow linear covariate effect
settings$RhoConfig['Beta1'] = 4

# Define formula.
X1_formula = ~ bs(CPE, degree=2, intercept=FALSE)
X2_formula = ~ 1

#
X1config_cp = array( c(0,1,1,0), dim=c(2,2))

# Run model
fit = fit_model( "settings" = settings,
  Lat_i = example$sampling_data[,'Lat'],
  Lon_i = example$sampling_data[,'Lon'],
  t_i = example$sampling_data[,'Year'],
  b_i = example$sampling_data[,'Catch_KG'],
  a_i = example$sampling_data[,'AreaSwept_km2'],
  c_i = example$sampling_data[,'Species']-1,
  X1_formula = X1_formula,
  X2_formula = X2_formula,
  X1config_cp = X1config_cp,
  covariate_data = covariate_data,
  test_fit = FALSE,
  getsd = TRUE,
  #Parameters = ParHat,
  newtonsteps = 0 )

#####################
# Effects package
#####################

library(effects)  # Used to visualize covariate effects

# Must add data-frames to global environment (hope to fix in future)
covariate_data_full = fit$effects$covariate_data_full
catchability_data_full = fit$effects$catchability_data_full

# Plot 1st linear predictor, but could use `transformation` to apply link function
pred = Effect.fit_model( fit,
  focal.predictors = c("CPE"),
  which_formula = "X1",
  xlevels = 100,
  transformation = list(link=identity, inverse=identity),
  category_number = 2 )
plot(pred)

#####################
# marginaleffects package
#####################

# Plot 1st linear predictor, but could use `transformation` to apply link function
quant = function(x) seq(min(x),max(x),length=21)
newdata = datagrid( newdata=fit$covariate_data[,'CPE',drop=FALSE], CPE=quant )
  pred = predictions( fit, newdata=newdata, covariate="X1" )

library(ggplot2)
library(gridExtra)
ggplot( pred, aes(CPE, predicted)) +
  geom_line( aes(y=predicted), color="blue", size=1 ) +
  geom_ribbon( aes( x=CPE, ymin=conf.low, ymax=conf.high), fill=rgb(0,0,1,0.2) ) +
  facet_wrap(vars(category), scales="free", ncol=2) +
  labs(y="Predicted response")

[James-Thorson-NOAA]

I should say too ... it requires VAST >= 3.10.0

[Charles-Adams-NOAA]

Got a warning that the model is likely not converged

[James-Thorson-NOAA]

But did it successfully make the marginaleffects plots?

On Wednesday, February 22, 2023, Charles-Adams-NOAA < @.***> wrote:

Got a warning that the model is likely not converged

— Reply to this email directly, view it on GitHub https://github.com/James-Thorson-NOAA/VAST/issues/359#issuecomment-1440144783, or unsubscribe https://github.com/notifications/unsubscribe-auth/AL62VMSNAJTWXF6HPLOSNRTWYYQMNANCNFSM6AAAAAAT5LDKSY . You are receiving this because you commented.Message ID: @.***>

-- James Thorson (he/him)

Program leader, Habitat and Ecological Processes Research (HEPR) Alaska Fisheries Science Center, NMFS Affiliate Faculty, University of Washington

[Charles-Adams-NOAA]

No, it did not:

library(marginaleffects)

Plot 1st linear predictor, but could use transformation to apply link function

quant = function(x) seq(min(x),max(x),length=21) newdata = datagrid( newdata=fit$covariate_data[,'CPE',drop=FALSE], CPE=quant ) pred = predictions( fit, newdata=newdata, covariate="X1" ) Error: Unable to compute predicted values with this model. You can try to supply a different dataset to the newdata argument. This error was also raised:

unused arguments (covariate = "X1", newdata = list(c(-0.93076, -0.7374245, -0.544089, -0.3507535, -0.157418, 0.0359175, 0.229253, 0.4225885, 0.615924, 0.8092595, 1.002595, 1.1959305, 1.389266, 1.5826015, 1.775937, 1.9692725, 2.162608, 2.3559435, 2.549279, 2.7426145, 2.93595), 1:21), type = "response")

Bug Tracker: https://github.com/vincentarelbundock/marginaleffects/issues In addition: Warning message: These arguments are not supported for models of class fit_model: covariate. Please file a request on Github if you believe that additional arguments should be supported: https://github.com/vincentarelbundock/marginaleffects/issues Loading required package: sp

[James-Thorson-NOAA]

Hmm. Annoyingly, it looks like marginaleffects has already been updated since last month when this code snippet was working, and I didn't record the packageVersion that I was using before. I'll add it to the list to investigate more.

[James-Thorson-NOAA]

I finally found time to track down why marginaleffects stopped working ... they changed the syntax for one of their functions in (I think) release 0.10.0.

@Charles-Adams-NOAA are you interested in trying that script again?

[Charles-Adams-NOAA]

Sure, I've got marginaleffects_0.10.0, so I assume I need to install VAST from the dev branch?

[Charles-Adams-NOAA]

Disregard last post, I see you've changed the code. Stand by

[Charles-Adams-NOAA]

Getting a different error message now:

quant = function(x) seq(min(x),max(x),length=21) newdata = datagrid( newdata=fit$covariate_data[,'CPE',drop=FALSE], CPE=quant ) pred = predictions( fit, newdata=newdata, covariate="X1" ) Error: Unable to compute predicted values with this model. You can try to supply a different dataset to the newdata argument. This error was also raised:

unused arguments (covariate = "X1", newdata = list(c(-0.93076, -0.7374245, -0.544089, -0.3507535, -0.157418, 0.0359175, 0.229253, 0.4225885, 0.615924, 0.8092595, 1.002595, 1.1959305, 1.389266, 1.5826015, 1.775937, 1.9692725, 2.162608, 2.3559435, 2.549279, 2.7426145, 2.93595), 1:21), type = "response")

Bug Tracker: https://github.com/vincentarelbundock/marginaleffects/issues In addition: Warning message: These arguments are not supported for models of class fit_model: covariate. Please file a request on Github if you believe that additional arguments should be supported: https://github.com/vincentarelbundock/marginaleffects/issues

[James-Thorson-NOAA]

what packageVersion("VAST") and packageVersion("FishStatsUtils") do you have?

On Thu, Mar 9, 2023 at 1:06 PM Charles-Adams-NOAA @.***> wrote:

Getting a different error message now:

quant = function(x) seq(min(x),max(x),length=21) newdata = datagrid( newdata=fit$covariate_data[,'CPE',drop=FALSE], CPE=quant ) pred = predictions( fit, newdata=newdata, covariate="X1" ) Error: Unable to compute predicted values with this model. You can try to supply a different dataset to the newdata argument. This error was also raised:

unused arguments (covariate = "X1", newdata = list(c(-0.93076, -0.7374245, -0.544089, -0.3507535, -0.157418, 0.0359175, 0.229253, 0.4225885, 0.615924, 0.8092595, 1.002595, 1.1959305, 1.389266, 1.5826015, 1.775937, 1.9692725, 2.162608, 2.3559435, 2.549279, 2.7426145, 2.93595), 1:21), type = "response")

Bug Tracker: https://github.com/vincentarelbundock/marginaleffects/issues In addition: Warning message: These arguments are not supported for models of class fit_model: covariate. Please file a request on Github if you believe that additional arguments should be supported: https://github.com/vincentarelbundock/marginaleffects/issues

— Reply to this email directly, view it on GitHub https://github.com/James-Thorson-NOAA/VAST/issues/359#issuecomment-1462821360, or unsubscribe https://github.com/notifications/unsubscribe-auth/AL62VMU6IG3ZSGCTAXQ2X5DW3JA5JANCNFSM6AAAAAAT5LDKSY . You are receiving this because you commented.Message ID: @.***>

[Charles-Adams-NOAA]

VAST_3.10.0 FishStatsUtils_2.12.0

[agruss2]

Hi everyone,

I have been digging into the functions of the "effects" package and Google to find something out, but I am still unclear about what I was looking for.

Let's say that I am using Obs_Model = c(2,1) so that the link function is the log function in both linear predictors.

Then, I if I use the following line of code:

pred_X2 <- Effect.fit_model( Fit, focal.predictors = "Bottom_depth", which_formula = "X2",
	transformation = list( link = identity, inverse = identity ), pad_values = 1 )

I am getting insights into the marginal effect of bottom depth on density or insights into the marginal effect of bottom depth on log-density?

Many thanks in advance!

[James-Thorson-NOAA]

I don't actually remember if Effect.fit_model uses the link space or the natural (inverse-link) space by default. However, you could do some sanity checks by fitting a linear covariate response in a log-linked linear model (like a Poisson-link GLMM). If the plot looks linear, then it must be doing a log-link. If it is curved, then it must be using natural-space.

[agruss2]

Thanks a lot Jim, this is a great idea!