Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Effect.fit_model new error on latest version of VAST #359

Open
smormede opened this issue Jan 17, 2023 · 22 comments
Open

Effect.fit_model new error on latest version of VAST #359

smormede opened this issue Jan 17, 2023 · 22 comments

Comments

@smormede
Copy link

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
Copy link
Owner

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
Copy link
Author

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
Copy link

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
Copy link
Author

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
Copy link

Already did that. And loaded effects package

@James-Thorson-NOAA
Copy link
Owner

James-Thorson-NOAA commented Feb 21, 2023 via email

@Charles-Adams-NOAA
Copy link

Happy to give it a try Jim!

@James-Thorson-NOAA
Copy link
Owner

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
Copy link
Owner

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

@Charles-Adams-NOAA
Copy link

Got a warning that the model is likely not converged

@James-Thorson-NOAA
Copy link
Owner

James-Thorson-NOAA commented Feb 22, 2023 via email

@Charles-Adams-NOAA
Copy link

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
Copy link
Owner

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
Copy link
Owner

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
Copy link

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

@Charles-Adams-NOAA
Copy link

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

@Charles-Adams-NOAA
Copy link

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
Copy link
Owner

James-Thorson-NOAA commented Mar 9, 2023 via email

@Charles-Adams-NOAA
Copy link

VAST_3.10.0
FishStatsUtils_2.12.0

@agruss2
Copy link
Collaborator

agruss2 commented Jul 7, 2024

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
Copy link
Owner

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
Copy link
Collaborator

agruss2 commented Jul 16, 2024

Thanks a lot Jim, this is a great idea!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@smormede @agruss2 @Charles-Adams-NOAA @James-Thorson-NOAA