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.

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

get_response(<brms>) fails with aterms #779

Closed
Chinca-lang opened this issue Jun 20, 2023 · 6 comments · Fixed by #781 or #785
Closed

get_response(<brms>) fails with aterms #779

Chinca-lang opened this issue Jun 20, 2023 · 6 comments · Fixed by #781 or #785
Assignees

Comments

@Chinca-lang
Copy link

Hi,
I'm trying to use describe_posterior with a model fitted with brm but I get this error:

describe_posterior(modelSex2, ci=0.89)
Error in [.data.frame(model_data, , rn, drop = FALSE) :
undefined columns selected

My model is:

modelSex2 <- brm(formula= ALI | trunc(lb = 0, ub = 10)~Sex*as.factor(date) + (1|ID),
chains = 10, iter = 5000, warmup = 1000,
refresh = 0,
data = data_as_cal, family = poisson)

I suspect that trunc function is causing some trouble as I've fitted the same model without trunc and describe_posterior is working well.
Do you have any suggestions?

@Chinca-lang
Copy link
Author

Chinca-lang commented Jun 21, 2023

just a quick update
If I set the rope_range it works:

describe_posterior(modelSex2, ci=0.89, rope_range = rope_range)

where rope_range is:

rope_value <- 0.1*sd(data_as_cal$ALI)
rope_range <- c(-rope_value, rope_value)

I don't know why it was difficult to calculate automatically the default rope range. But this fixed the issue.

@mattansb
Copy link
Member

This is actually an issue in insight::get_response().

library(brms)
library(bayestestR)

fit4 <- brm(count | trunc(ub = 104) ~ zBase * Trt,
            data = epilepsy, family = poisson(),
            backend = "cmdstanr", cores = 4)

insight::get_response(fit4) # source = "environment"
#> Error in `[.data.frame`(model_data, , rn, drop = FALSE) : 
#>   undefined columns selected

insight::get_response(fit4, source = "mf")
#> Error in `[.data.frame`(model_data, , rn, drop = FALSE) : 
#>   undefined columns selected

@mattansb mattansb reopened this Jun 21, 2023
@mattansb mattansb transferred this issue from easystats/bayestestR Jun 21, 2023
@mattansb mattansb changed the title describe_posterior in brms get_response(<brms>) fails with a truncated response Jun 21, 2023
@strengejacke strengejacke self-assigned this Jun 21, 2023
@mattansb
Copy link
Member

The issue is a general issue with brms's aterms.
Another example: easystats/performance#596

@mattansb mattansb changed the title get_response(<brms>) fails with a truncated response get_response(<brms>) fails with aterms Jun 21, 2023
strengejacke added a commit that referenced this issue Jun 21, 2023
* `get_response(<brms>)` fails with a truncated response
Fixes #779

* add test

* news, version

* fix tests

* fix inaccuracy in test
@mattansb
Copy link
Member

When the aterm has a constant (not a data column), you get the error:

library(brms)

fit4 <- brm(am | trials(1) ~ hp,
            data = mtcars, family = binomial(),
            backend = "cmdstanr", cores = 4)

insight::get_response(fit4) # source = "environment"
#> Error in `[.data.frame`(model_data, , rn, drop = FALSE) : 
#>   undefined columns selected

insight::get_response(fit4, source = "mf")
#> Error in `[.data.frame`(model_data, , rn, drop = FALSE) : 
#>   undefined columns selected

This also means that your fix at #781 is wrong for this bug:

  1. If the aterm uses a variable from the data, that should be returned (regardless of the aterm type).
  2. If is it a constant, don't return it.

Originally posted by @mattansb in easystats/performance#596 (comment)

@mattansb mattansb reopened this Jun 21, 2023
strengejacke added a commit that referenced this issue Jun 21, 2023
strengejacke added a commit that referenced this issue Jun 21, 2023
strengejacke added a commit that referenced this issue Jun 21, 2023
* Fix aterms in brms-responses
Fixes #779

* fix, add tests

* spelling, lintr

* styler

* minor
@SethMusker
Copy link

I'm still getting this error (along with other related[?] ones) with a multivariate brms model with missing data imputation.

Example based on the marginaleffects brms vignette (https://marginaleffects.com/vignettes/brms.html)

library(marginaleffects)
library(brms)
dat <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/carData/TitanicSurvival.csv")
dat$survived <- ifelse(dat$survived == "yes", 1, 0)
dat$woman <- ifelse(dat$sex == "female", 1, 0) |> as.factor()
sum(is.na(dat$age))  # <-- 263 missing values

mod_miss <- brm(bf(survived ~ woman*mi(age) + passengerClass, family=bernoulli(link = "logit")) +
              bf(age | mi() ~ passengerClass + woman) +
              set_rescor(FALSE),
           data = dat,
           backend = "cmdstanr",
           cores=4)

#--------------------------------------
## Example 1
plot_slopes(mod_miss, 
   variables = "woman", 
   condition = "age")
# Error in `[.data.frame`(model_data, , rn, drop = FALSE) : 
#   undefined columns selected

#--------------------------------------
## Example 2
slopes(mod_miss,
    newdata = datagrid(
    woman = 1,
    passengerClass = "1st"))
# Error: The following variables can neither be found in 'data' nor in 'data2':
#   'age'
# In addition: Warning message:
#   These variables were not found: survived, age, .  Try specifying the `newdata` argument explicitly and make sure
# the missing variable is included.

#--------------------------------------
## Example 3
slopes(mod_miss,
    variables = "age",
    newdata = datagrid(woman = 0:1)) |>
    posterior_draws()
# Error: There is no valid predictor variable. Please change the `variables` argument or supply a new data frame to the
# `newdata` argument.
# In addition: Warning message:
#   These variables were not found: miage.  Try specifying the `newdata` argument explicitly and make sure the
# missing variable is included. 

And the tracebacks.

#--------------------------------------
## Example 1
8: stop("undefined columns selected")
7: `[.data.frame`(model_data, , rn, drop = FALSE)
6: model_data[, rn, drop = FALSE]
5: get_response.default(model)
4: insight::get_response(model)
3: sanitize_condition(model, condition, variables, modeldata = modeldata)
2: plot_comparisons(model, variables = variables, condition = condition, 
       by = by, newdata = newdata, type = type, vcov = vcov, conf_level = conf_level, 
       wts = wts, draw = draw, rug = rug, gray = gray, comparison = slope, 
       ...)
1: plot_slopes(mod_miss, variables = "woman", condition = "age")

#--------------------------------------
## Example 2
7: stop(format_message(string = string, ..., line_length = line_length, 
       indent = indent), call. = call)
6: format_alert(..., type = "error")
5: insight::format_error(pred_both[["error"]][["message"]])
4: get_contrasts(structure(list(formula = structure(list(forms = list(
       survived = structure(list(formula = survived ~ woman * mi(age) + 
           passengerClass, pforms = list(), pfix = list(), family = structure(list(
           family = "bernoulli", link = "logit", linkfun = function (mu) 
           link(mu, link = slink), linkinv = function (eta) 
           inv_link(eta, link = slink), dpars = "mu", type = "int", 
           ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c("weights", 
           "subset", "index"), specials = c("binary", "sbi_logit"
           )), class = c("brmsfamily", "family")), resp = "survived", 
           mecor = TRUE), class = c("brmsformula", "bform")), age = structure(list(
           formula = age | mi() ~ passengerClass + woman, pforms = list(), 
           pfix = list(), resp = "age", family = structure(list(
               family = "gaussian", link = "identity", linkfun = function (mu) 
               link(mu, link = slink), linkinv = function (eta) 
               inv_link(eta, link = slink), dpars = c("mu", "sigma"
               ), type = "real", ybounds = c(-Inf, Inf), closed = c(NA, 
               NA), ad = c("weights", "subset", "se", "cens", "trunc", 
               "mi", "index"), normalized = c("_time_hom", "_time_het", 
               "_lagsar", "_errorsar", "_fcor"), specials = c("residuals", 
               "rescor")), class = c("brmsfamily", "family")), mecor = TRUE), class = c("brmsformula", 
    ...
3: do.call("get_contrasts", args)
2: comparisons(model, newdata = newdata, variables = variables, 
       vcov = vcov, conf_level = conf_level, type = type, wts = wts, 
       hypothesis = hypothesis, equivalence = equivalence, df = df, 
       p_adjust = p_adjust, by = by, eps = eps, numderiv = numderiv, 
       comparison = slope, cross = FALSE, internal_call = TRUE, 
       ...)
1: slopes(mod_miss, newdata = datagrid(woman = 1, passengerClass = "1st"))

#--------------------------------------
## Example 3
7: stop(format_message(string = string, ..., line_length = line_length, 
       indent = indent), call. = call)
6: format_alert(..., type = "error")
5: insight::format_error(msg)
4: sanitize_variables(model = model, newdata = newdata, modeldata = modeldata, 
       variables = variables, cross = cross, by = by, comparison = comparison, 
       eps = eps)
3: comparisons(model, newdata = newdata, variables = variables, 
       vcov = vcov, conf_level = conf_level, type = type, wts = wts, 
       hypothesis = hypothesis, equivalence = equivalence, df = df, 
       p_adjust = p_adjust, by = by, eps = eps, numderiv = numderiv, 
       comparison = slope, cross = FALSE, internal_call = TRUE, 
       ...)
2: slopes(mod_miss, variables = "age", newdata = datagrid(woman = 0:1))
1: posterior_draws(slopes(mod_miss, variables = "age", newdata = datagrid(woman = 0:1)))

sessionInfo()

R version 4.3.0 (2023-04-21 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=English_South Africa.utf8  LC_CTYPE=English_South Africa.utf8    LC_MONETARY=English_South Africa.utf8
[4] LC_NUMERIC=C                          LC_TIME=English_South Africa.utf8    

time zone: Africa/Johannesburg
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] brms_2.21.0            Rcpp_1.0.12            marginaleffects_0.21.0

loaded via a namespace (and not attached):
 [1] bridgesampling_1.1-2 tensorA_0.36.2.1     utf8_1.2.4           generics_0.1.3       stringi_1.8.4       
 [6] lattice_0.21-8       magrittr_2.0.3       grid_4.3.0           mvtnorm_1.2-5        jsonlite_1.8.8      
[11] Matrix_1.5-4         processx_3.8.4       pkgbuild_1.4.4       backports_1.5.0      ps_1.7.7            
[16] gridExtra_2.3        Brobdingnag_1.2-9    fansi_1.0.6          QuickJSR_1.2.2       scales_1.3.0        
[21] codetools_0.2-19     abind_1.4-5          cli_3.6.3            rlang_1.1.4          cmdstanr_0.8.1      
[26] munsell_0.5.1        withr_3.0.0          StanHeaders_2.32.9   tools_4.3.0          rstan_2.32.6        
[31] inline_0.3.19        parallel_4.3.0       rstantools_2.4.0     checkmate_2.3.1      coda_0.19-4.1       
[36] dplyr_1.1.4          colorspace_2.1-0     ggplot2_3.5.1        vctrs_0.6.5          posterior_1.6.0     
[41] R6_2.5.1             matrixStats_1.3.0    stats4_4.3.0         lifecycle_1.0.4      stringr_1.5.1       
[46] insight_0.20.1.13    pkgconfig_2.0.3      RcppParallel_5.1.7   pillar_1.9.0         gtable_0.3.5        
[51] loo_2.8.0            data.table_1.15.4    glue_1.7.0           collapse_2.0.15      xfun_0.45           
[56] tibble_3.2.1         tidyselect_1.2.1     knitr_1.47           rstudioapi_0.16.0    bayesplot_1.11.1    
[61] nlme_3.1-162         compiler_4.3.0       distributional_0.4.0

@strengejacke
Copy link
Member

Maybe related to #896? Not sure, will check later.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
4 participants