Dear all,
I’m using a truncated quantile (median) regression (using asym_laplace) of the type
bf(Y | trunc(lb = 0) ~ 1 + x + (1 | id1) + (1 + x | id2), quantile = 0.5),
and would like to get conditional effects via marginal_effects. brms tells me
Error: Fitted values on the respone scale not yet implemented for truncated 'asym_laplace' models.
My question: How difficult and effortful would it be (for me?) to properly implement that functionality and how would I approach this? And if it’s too difficult and effortful, is there a workaround I can use to create the kind of data marginal_effects is producing?
Thanks!
- Operating System: macOS 10.14.6
- brms Version: 2.10.0
What you need is mainly a formula for the mean of a truncated asymmetric laplace distribution.
I wasn’t able to find a closed-form formula online, but I found a numerical solution based on this paper: https://www.jstatsoft.org/article/view/v016c02 (see p. 3)
dtrunc <- function(x,
spec,
a = -Inf,
b = Inf,
...)
{
tt <- rep(0, length(x))
g <- get(paste("d", spec, sep = ""), mode = "function")
G <- get(paste("p", spec, sep = ""), mode = "function")
tt[x >= a & x <= b] <-
g(x[x >= a & x <= b], ...) / (G(b, ...) - G(a, ...))
return(tt)
}
extrunc <- function(spec, a = -Inf, b = Inf, ...)
{
f <- function(x)
x * dtrunc(x, spec, a = a, b = b, ...)
return(integrate(f, lower = a, upper = b)$value)
}
This then seems to work with brms’s dasym_laplace and pasym_laplace
extrunc("asym_laplace", mu = 1, a = 0)
This approach seems generic enough to support any kind of truncated distributions that brms already supports in a non-truncated way.
Looking at vignette("brms_customfamilies") it seems to me that the quick fix would be to add a function fitted_trunc_asym_laplace in fitted.R.
fitted_trunc_asym_laplace <- function(draws, lb, ub) {
###
}
Or would it make more sense to replace dasym_laplace, pasym_laplace, qasym_laplace, and rasym_laplace in distributions.R with versions that can handle truncation? Or is this irrelevant for postprocessing?
Computing the truncated density is not a problem and brms does that already for arbitary models
The remaining problem is the truncated mean. Of course numerical integration, as done in extrunc is always a brute force solution but remember that this needs to be done for each single posterior draws of each single observation, which will be highly inefficient.
OK. Then simply focussing on the use case of me creating a workaround (where I accept the inefficiency of the numerical approach), that is, not a general solution for brms, would it be sufficient to implement it in
fitted_trunc_asym_laplace <- function(draws, lb, ub) {
###
}
?
Yes it should be sufficient.
1 Like
Dear Paul,
Last October above we discussed about how to implement a brute-force method to calculate the expectation of the median of a truncated asymmetric Laplace distribution.
Based on the functions dtrunc and extrunc above (taken from p. 3 of https://www.jstatsoft.org/article/view/v016c02), I’ve defined
fnc_exp_median_trunc_asym_laplace <- function(mu, sigma, lb, ub) {
extrunc(
mu = mu,
sigma = sigma,
spec = "asym_laplace",
a = lb,
b = ub,
quantile = .5
)
}
and
posterior_epred_trunc_asym_laplace <- function(prep, lb, ub) {
mapply(fnc_exp_median_trunc_asym_laplace,
prep$dpars$mu,
prep$dpars$sigma,
lb,
ub)
}
which computes when I apply it.
But I realize that I do not yet properly handle or understand prep$dpar$mu and prep$dpar$sigma, since they may differ in length (or rather dimensions) depending on the model fitted.
Is there an another example in posterior_epred.R you can point me to better understand this and possibly even adapt to my use case?
I’m aware that I haven’t yet provided a MRE, so if you think one would be useful here, please let me know.
Thanks!
Perhaps the code used in https://github.com/paul-buerkner/brms/blob/master/R/posterior_predict.R, especially the use of the get_dpar function could be helpful in your setting.
1 Like
Thanks a lot for the pointer!