Open smormede opened 1 year ago
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)
Hi James,
Happy to help where I can, in particular testing etc. Maybe send me an email and we can look at it.
Sophie
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 loadcovariate_data_full
andcatchability_data_full
into global memory
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
Already did that. And loaded effects package
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: @.***>
Happy to give it a try Jim!
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")
I should say too ... it requires VAST >= 3.10.0
Got a warning that the model is likely not converged
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
No, it did not:
library(marginaleffects)
Plot 1st linear predictor, but could use
transformation
to apply link functionquant = 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
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.
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?
Sure, I've got marginaleffects_0.10.0, so I assume I need to install VAST from the dev branch?
Disregard last post, I see you've changed the code. Stand by
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
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: @.***>
VAST_3.10.0 FishStatsUtils_2.12.0
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!
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.
Thanks a lot Jim, this is a great idea!
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 foundI 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