Unable to fit nonlinear robit model to new data

Interfaces brms
1

I fit a nonlinear robit model in brms that was set up as follows:

stan_inv_robit <- "
  real inv_robit(real y, real nu) {
    return(student_t_cdf(y, nu, 0, (nu - 2) / nu));
  }
"
stanvar_inv_robit <- stanvar(scode = stan_inv_robit, block = "functions")

bf_1 <- bf(Accuracy_rev ~ inv_robit(eta, nu),
           nu ~ 1,
           b + m ~ 1,
           bl + tau ~ 1 + (1|Subset + User + Token),
           nlf(eta ~ (tau-((bl + repetition)^.1)*(model_time_z^-(b + m * stability_z)))/.1),
           nl = T,
           family = bernoulli('identity'))

prior_1 <-
  prior(gamma(2, 0.1), nlpar = "nu", lb = 2) +
  prior(normal(0, .5), nlpar = 'b') +
  prior(normal(0, .5), nlpar = 'm') +
  prior(normal(0, .5), nlpar = 'bl') +
  prior(normal(0, .5), nlpar = 'tau')

fit_1 <- brm(bf_1,
             prior = prior_1,
             stanvars = stanvar_inv_robit,
             data = filter(dat_prepped, `Train or Test` == 1),
             chains = 4)

I can obtain model fits for the data set that produced the parameter estimates, but when I try to get fitted values for a new set of data (i.e., fitted(fit_1, newdata = test, allow_new_levels = T) I get the following error:

Error: Error in (function (y, nu, pstream__ = <pointer: 0x000000000211ab40>) : Exception: student_t_cdf: Random variable is nan, but must not be nan! (in 'unknown file name' at line 5)

At first I thought my test set might have been returning NA values of eta that are passed on to the inv_robit function, but I was able to calculate non-NA values for every observation using the estimated parameter values. So I’m not sure what the issue is. Thanks for any advice!

1 Like
2

Can you show all the post processing code please? It seems you only showed parts of the code that is relevant.

3

Sure – let me know if there’s anything else that might be helpful.

# Expose inv_robit function
expose_functions(fit_1, vectorize = T)

# Get training observations from main data set
train <- na.omit(filter(dat_prepped, `Train or Test` == 1))
train$fitted <- fitted(fit_1)[,1] # Runs fine

# Get test observations from main data set
test <- na.omit(filter(dat_prepped, `Train or Test` == 0))
test$fitted <- fitted(fit_1, newdata = test, allow_new_levels = T)[,1] # Generates an error

# Get values of eta from model parameters and see if any are NA
test <- mutate(test, b = fitted(fit_1, newdata = test, allow_new_levels = T, nlpar = 'b')[,1],
               m = fitted(fit_1, newdata = test, allow_new_levels = T, nlpar = 'm')[,1],
               bl = fitted(fit_1, newdata = test, allow_new_levels = T, nlpar = 'bl')[,1],
               tau = fitted(fit_1, newdata = test, allow_new_levels = T, nlpar = 'tau')[,1],
               eta = (tau-((bl + repetition)^.1)*(model_time_z^-(b + m * stability_z)))/.1)

length(which(is.na(test$eta))) # Returns 0
4

Hmm unfortunately nothing in that code to identify the problem right away. Could you make a minimal reproducible example so that I can check myself?

5

While in the process of making up an example, I realized that I was getting NaN values for eta if I used summary = F in the fitted function. And that seems to be because my priors weren’t reflecting the fact that I’m using the identity link. So this issue seems to have been completely my fault. I’ll re-run with better priors and let you know if anything pops up.