`conformal_cv+` Fails Inside a Function When Using `formula=as.formula()`

[mreece13]

---

name: "conformal_cv+" Fails Inside a Function When Using formula=as.formula()
about: Use of "conformal_cv+" inside inferences()

Thanks a lot for the great package, and for beginning to implement some conformal prediction functions. I think I've found a bug related to the use of the "conformal_cv+" option in inferences(). In short, there seems to be a scoping issue when I use as.formula() in the "formula" option of a given model when nested inside a function. This is useful when needing to dynamically change covariates depending on some features of the data or the data inputted to a function. A minimal example is included below. I'd be happy to offer a PR, but I know very little about scoping in R packages.

library(marginaleffects)

# not the most rigorous split, just for illustration
iris = dplyr::slice_sample(iris, n=150) # shuffle to ensure overlap
train = iris[1:50,]
calib = iris[51:100,]
test = iris[101:150,]

a <- function(){
  lm(
    formula = Sepal.Length ~ Sepal.Width * Species, 
    data = train
  )
}
b <- function(){
  
  form <- "Sepal.Length ~ Sepal.Width * Species"
  
  lm(
    formula = as.formula(form), 
    data = train
  )
  
}

# works with defined function
a() |> 
  predictions() |> 
  inferences(
    method = "conformal_cv+",
    conformal_test = test,
    R = 10
  )
#> 
#>  Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 % Pred. 2.5 % Pred. 97.5 %
#>      6.00      0.126 47.6   <0.001   Inf  5.75   6.25        5.04         6.96
#>      5.35      0.293 18.3   <0.001 246.0  4.78   5.93        4.39         6.31
#>      5.85      0.158 37.0   <0.001 992.3  5.54   6.16        4.89         6.81
#>      6.63      0.108 61.3   <0.001   Inf  6.41   6.84        5.66         7.59
#>      6.22      0.153 40.7   <0.001   Inf  5.92   6.52        5.26         7.18
#> --- 40 rows omitted. See ?print.marginaleffects ---
#>      6.15      0.134 45.7   <0.001   Inf  5.89   6.41        5.19         7.11
#>      7.05      0.199 35.4   <0.001 909.0  6.66   7.44        6.09         8.01
#>      6.84      0.143 47.8   <0.001   Inf  6.56   7.12        5.88         7.80
#>      6.84      0.143 47.8   <0.001   Inf  6.56   7.12        5.88         7.80
#>      6.63      0.108 61.3   <0.001   Inf  6.41   6.84        5.66         7.59
#> Type: response

# fails on `as.formula()` use
b() |> 
  predictions() |> 
  inferences(
    method = "conformal_cv+",
    conformal_test = test,
    R = 10
  )
#> Error in eval(mf, parent.frame()): object 'form' not found

# works with "conformal_split"
b() |> 
  predictions() |> 
  inferences(
    method = "conformal_split",
    conformal_calibration = calib,
    conformal_test = test
  )
#> 
#>  Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 % Pred. 2.5 % Pred. 97.5 %
#>      6.00      0.126 47.6   <0.001   Inf  5.75   6.25        5.12         6.88
#>      5.35      0.293 18.3   <0.001 246.0  4.78   5.93        4.47         6.23
#>      5.85      0.158 37.0   <0.001 992.3  5.54   6.16        4.97         6.73
#>      6.63      0.108 61.3   <0.001   Inf  6.41   6.84        5.75         7.51
#>      6.22      0.153 40.7   <0.001   Inf  5.92   6.52        5.34         7.10
#> --- 40 rows omitted. See ?print.marginaleffects ---
#>      6.15      0.134 45.7   <0.001   Inf  5.89   6.41        5.27         7.03
#>      7.05      0.199 35.4   <0.001 909.0  6.66   7.44        6.17         7.93
#>      6.84      0.143 47.8   <0.001   Inf  6.56   7.12        5.96         7.72
#>      6.84      0.143 47.8   <0.001   Inf  6.56   7.12        5.96         7.72
#>      6.63      0.108 61.3   <0.001   Inf  6.41   6.84        5.75         7.51
#> Type: response

sessionInfo()
#> R version 4.5.0 (2025-04-11)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS 26.0
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] marginaleffects_0.30.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] vctrs_0.6.5       cli_3.6.5         knitr_1.50        rlang_1.1.6      
#>  [5] xfun_0.52         generics_0.1.4    data.table_1.17.8 glue_1.8.0       
#>  [9] backports_1.5.0   htmltools_0.5.8.1 rmarkdown_2.29    evaluate_1.0.4   
#> [13] tibble_3.3.0      fastmap_1.2.0     yaml_2.3.10       lifecycle_1.0.4  
#> [17] insight_1.4.0     compiler_4.5.0    dplyr_1.1.4       fs_1.6.6         
#> [21] pkgconfig_2.0.3   rstudioapi_0.17.1 digest_0.6.37     R6_2.6.1         
#> [25] reprex_2.1.1      tidyselect_1.2.1  pillar_1.11.0     magrittr_2.0.3   
#> [29] checkmate_2.3.2   tools_4.5.0       withr_3.0.2

^{Created on 2025-09-25 with reprex v2.1.1}

[vincentarelbundock]

what happens if you assign the formula to a temporary variable instead of trying to create it inline in the model spec?

[mreece13]

Same problem unfortunately. I assume you mean form <- as.formula("Sepal.Length ~ Sepal.Width * Species") inside b().

[vincentarelbundock]

I meant something like this. The <<- is an ugly hack that might break parallelism, but I fear it might be hard to get around that because we rely on insight::find_formula() to extract the formula from the call, and scoping issues are very hard to resolve robustly across model types.

library(marginaleffects)
iris = dplyr::slice_sample(iris, n = 150) # shuffle to ensure overlap
train = iris[1:50, ]
calib = iris[51:100, ]
test = iris[101:150, ]

b <- function(dv = "Sepal.Length") {
    form <- paste(dv, "~ Sepal.Width * Species")
    form <<- as.formula(form)
    model <- lm(form, data = train)
    predictions(model) |>
        inferences(
            method = "conformal_cv+",
            conformal_test = test,
            R = 10
        )
}
b()

[mreece13]

Hmm, I see, this does work yes. I can probably work around the issues with parallelism but you're right that it's a bit unfortunate. Perhaps for cases when I need parallelism I will fully redefine the model depending on the input instead of just the formula.