Sampling from a shifted lognormal distribution (rshifted_lnorm) produces negative values

I am simulating data according to a logshifted distribution in order to test my model before running it on real data.

As expected, the shifted lognormal family in brms only accepts positive values, in other words, the shift can only be positive.
However, the brms::rshifted_lnorm distribution can also generate negative values.

e.g. if I take estimated parameter values from a previous paper on reaction times and input them into the function:

rshifted_lnorm(1, meanlog = 0.14, sdlog = 0.7, shift = -0.85)` 

I get a consistent number of negative values.

I can get the kind of values I want by manually building a shifted lognormal:

exp(rnorm(1, mean = 0.14, sd = 0.7) - 0.85)

But I am wondering as to whether the inner workings of the shifted lognormal distribution in brms differ from this hand-coded solution (and therefore assessing parameter recovery would be useless, given the mismatch in the formula).

I don’t know what the shifted log normal is supposed to be doing but I can say that

exp(rnorm(1, mean = 0.14, sd = 0.7) - 0.85))

is equivalent to

exp(rnorm(1, mean = 0.14 - 0.85, sd = 0.7)))

So, I think what you want is

rshifted_lnorm(1, meanlog = 0.14 - 0.85, sdlog = 0.7, shift = 0)` 

However, if that is what you want, the shift in the shifted lognormal is superfluous. So I think what the lognormal is supposed to be is equivalent to

exp(rnorm(1, mean = meanlog, sd = sdlog)) - shift

edit: this should be

exp(rnorm(1, mean = meanlog, sd = sdlog)) + shift

thanks @jsocolar

1 Like

As @stijn suggests, a shifted lognormal is a lognormal that has been shifted, not a distribution whose logarithm is distributed as a shifted normal (which is again a normal; therefore this distribution is just a lognormal). If you want to see the inner workings of brms::rshifted_lnorm, these are just:

> brms::rshifted_lnorm
function (n, meanlog = 0, sdlog = 1, shift = 0) 
{
    args <- nlist(dist = "lnorm", n, shift, meanlog, sdlog)
    do_call(rshifted, args)
}
<bytecode: 0x7fa6d10de9b0>
<environment: namespace:brms>
> brms:::rshifted
function (dist, n, shift = 0, ...) 
{
    do_call(paste0("r", dist), list(n, ...)) + shift
}
<bytecode: 0x7fa5e3664358>
<environment: namespace:brms>

i.e. exp(rnorm(1, mean = meanlog, sd = sdlog)) + shift, which is exactly what @stijn had except for a + shift rather than a - shift

1 Like

You are both absolutely right. This was a silly mistake on my part.

My misunderstanding was due to the shift (ndt) parameter in brms being on a log scale, as far as I can see.

So, I used (all positive) empirical data (reaction times) to fit a model, and extracted the parameter values, which were inputed in my simulation. Shift/ndt was -0.85, which would - I guess - be a reasonable 400 ms of non-decision-time, but directly feeding it into rshifted_lnorm() made it generate a consistent chunk of negative values.

So I’m a bit baffled as to how to properly input brms output of the shifted_lognormal() into rshifted_lnorm(). Is it by simply exponentiating the shift parameter?

1 Like

From brms/families.R at master · paul-buerkner/brms · GitHub near line 556:

#' @rdname brmsfamily
#' @export
shifted_lognormal <- function(link = "identity", link_sigma = "log",
                              link_ndt = "log") {
  slink <- substitute(link)
  .brmsfamily("shifted_lognormal", link = link, slink = slink,
              link_sigma = link_sigma, link_ndt = link_ndt)
}

So you’re absolutely right that ndt is estimated by default on a log scale, which is useful when interpreting as a non-decision-time because it will never be estimated to be negative. If your domain knowledge is consistent with the possibility of negative shifts, then this obviously doesn’t work, but you can get around that if desired using shifted_lognormal(link_ndt = "identity"). This would yield a parameter that is directly interpretable as the shift.

As you say, the way to generate values from the shifted lognormal when estimating the ndt using a log link is to exponentiate the linear predictor. This isn’t specific to the shifted lognormal: in GLMs we always need to feed the response distribution the inverse-link of the linear predictor.

1 Like