brms icon indicating copy to clipboard operation
brms copied to clipboard

Ability to use stats::filter or other distributed lag models in brms formulas

Open DoktorMike opened this issue 5 years ago • 5 comments

As always, your package is awesome. There's something I've been wondering about for quite some time. I often find myself using the stats::filter function in R to do parameterized lag structures.

Below is an example where I fit a simple linear regression model with and without the transformation. What I do in a Bayesian setting is that I put a prior distribution on lambda and fit it along with the other parameters.

tmpdf <- structure(list(x = c(
    130076, 136354, 115825, 104377, 117773,
    125095, 128427, 118909, 122596, 135917, 178711, 186966, 231921,
    252732, 238187, 221672, 205662, 234370, 233703, 234875
), y = c(
    14176,
    20609, 18791, 20271, 22575, 21892, 21898, 25324, 23767, 25016,
    31272, 34117, 36663, 45698, 47814, 48442, 44758, 44492, 43979,
    47550
)), row.names = c(NA, -20L), class = c(
    "tbl_df", "tbl",
    "data.frame"
))
summary(m1 <- lm(y ~ x, data = tmpdf))
lambda <- 0.5
summary(m2 <- lm(y ~ stats::filter(x, lambda, method = "recursive"), data = tmpdf))

However, naively trying to do this with brms leads to the following

summary(m3 <- brm(bf(y ~ b * stats::filter(x, la, method = "recursive"), la + b ~ 1, nl = TRUE),
    prior = prior(normal(0, 1), nlpar = "b") +
        prior(beta(1, 1), nlpar = "la", ub=1, lb=0),
    data = tmpdf
))

which does not compile

Compiling Stan program...
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Variable "stats" does not exist.
 error in 'model1069758d87d2_b3d538277eab8ac3823d67aca8534dc6' at line 35, column 30
  -------------------------------------------------
    33:     for (n in 1:N) {
    34:       // compute non-linear predictor values
    35:       mu[n] = nlp_b[n] * stats::filter(C_1[n] , nlp_la[n] , method="recursive");
                                     ^
    36:     }
  -------------------------------------------------

Error in stanc(file = file, model_code = model_code, model_name = model_name,  :
  failed to parse Stan model 'b3d538277eab8ac3823d67aca8534dc6' due to the above error.

So is there a way for me to do this? Or is there another way to do convolution using brms?

DoktorMike avatar Mar 20 '21 17:03 DoktorMike

Some more info on this. I'm using stats::filter to implement a Distributed Lag Model here where the Model looks like this

image

where

image

where y and x are the same as in the code I wrote.

DoktorMike avatar Mar 20 '21 19:03 DoktorMike

This guy is doing something similar here where he's implementing a convolution. https://discourse.mc-stan.org/t/discrete-convolution-by-direct-summation-can-this-be-made-more-efficient/969

DoktorMike avatar Mar 20 '21 20:03 DoktorMike

This looks very similar to an auto-regressive structure. Have you checked out ar() terms in brms?

stats::filter will probably work in brms if it works in lm. but only in linear formulas. You tried with non-linear formulas.

In the future I recommend asking questions on https://discourse.mc-stan.org/ if you are unsure whether or not you actually have a feature request.

paul-buerkner avatar Mar 24 '21 09:03 paul-buerkner

This looks very similar to an auto-regressive structure. Have you checked out ar() terms in brms?

stats::filter will probably work in brms if it works in lm. but only in linear formulas. You tried with non-linear formulas.

In the future I recommend asking questions on https://discourse.mc-stan.org/ if you are unsure whether or not you actually have a feature request.

As you can see in the references I posted. This is indeed not possible today (as a nonlinear formula would be needed which doesn't work as I posted.) if you want to sample the lag parameter. It's not an AR model because it's not defined on the response but rather the covariates. I'll post this issue in discourse and you can reopen later if you deem it interesting enough. :)

Another reference of a discussion of a similar topic is here https://discourse.mc-stan.org/t/autoregressive-distributive-lag-adl-model/11650/2

DoktorMike avatar Mar 24 '21 12:03 DoktorMike

I see. I didn't mean to close this issue actually. Just hit the wrong button.

paul-buerkner avatar Mar 24 '21 12:03 paul-buerkner

I will close this low priority issue to reduce the load of the brms issue tracker. If somebody is motivated to work on it, please write here and I am happy to reopen this issue.

paul-buerkner avatar Jan 27 '24 23:01 paul-buerkner