Make default of `formula` independent of `sigma`

[athowes]

can we entirely avoid specifying sigma and thus make this portable across all families? I think it should build down to the same thing?

_Originally posted by @seabbs in https://github.com/epinowcast/epidist/pull/208#discussion_r1700385087_

The formula option of epidist currently has default option as bf(mu ~ 1, sigma ~ 1). But this is lognormal specific. Can we put it to just bf(mu ~ 1)? Doing this breaks some tests, but those tests are written by me rather than brms so probably just need altering not to catch something which isn't an issue.

[seabbs]

I'm fairly confident this will work based on prior brms usage and think it would be a nice improvement.

[athowes]

To solve this, as well as one of the bugs in https://github.com/epinowcast/epidist/issues/227, I think we should:

1. Create a new function

replace_brms_prior <- function(priors, new_prior) {
  # Check if new_prior is in priors (aside from a different distribution)
  # Replace priors entry with new_prior
  # If not, return priors
}

2. Change the workflow in epidist_prior to

epidist_prior <- function() {
  priors <- brms::default_prior()
  if(...) { replace_brms_prior(priors, new_prior) }
  if(...) { replace_brms_prior(priors, new_prior) }
  if(...) { replace_brms_prior(priors, new_prior) }
}

I know you had reservations about this @seabbs but this is

1. Using existing brms functionality (default_prior)
2. Otherwise I would have to start getting involved in extracting out bits of the formula and checking if we have an intercept for sigma and things like that -- easier if I can just let brms do that for me

Note to self: in these if statements need to be checking what link function user provides. Suggested change of prior depends on link functions.

[seabbs]

Have you just tired doing this: we put it to just bf(mu ~ 1)?

I don't think it makes sense it use if gates here regardless but I also don't understand how this proposal relates to this formula default question?

[athowes]

Yes I have tried doing that. The reason why that doesn't work is because then it tries to set a prior on something which doesn't exist (the intercept for sigma) => we need to make sure that the parameters we try to set priors on are actually in the model => using the default_prior approach above.

[seabbs]

sorry why doesn't it exist? This isn't setting to be a constant right so in the default brms mode where you don't explicitly write the intercept the intercept will still exist

[athowes]

That would be reasonable behaviour but currently what we have is not doing that:

> fit <- epidist(data = prep_obs, seed = 1)
ℹ The following priors have been set:
• normal(2, 0.5) on the intercept of distributional parameter mu
• normal(0, 0.5) on the intercept of distributional parameter sigma
To alter priors, or set priors on other parameters, see ?epidist_prior.
This message is displayed once every 8 hours.
Error: The following priors do not correspond to any model parameter: 
Intercept_sigma ~ normal(0, 0.5)
Function 'default_prior' might be helpful to you.

When I've looked at the default_prior for this model it is:

> default_prior <- brms::default_prior(epidist_formula, data = data)
> default_prior
                prior     class coef group resp dpar nlpar lb   ub  source
 student_t(3, 5, 2.5) Intercept                                    default
 student_t(3, 0, 2.5)     sigma                             0 <NA> default

=> there is the sigma parameter not the Intercept_sigma parameter => can't set a prior on Intercept_sigma

Further, when running this removing the epidist_prior from brms (i.e. using default priors) then we have the error:

> fit_no_prior <- brms::brm(
+   formula = epidist_formula, family = epidist_family,
+   stanvars = epidist_stancode, backend = "cmdstanr", data = data
+ )
Compiling Stan program...
Semantic error in '/var/folders/4r/hkp4v9fn3wx044_hk8qnsjxw0000gn/T/RtmpWq657o/model-6c960467d7d.stan', line 62, column 14 to column 76:
   -------------------------------------------------
    60:      vector[N] mu = rep_vector(0.0, N);
    61:      mu += Intercept;
    62:      target += latent_lognormal_lpdf(Y | mu, sigma, pwindow, swindow, vreal1);
                       ^
    63:    }
    64:    // priors including constants
   -------------------------------------------------

Ill-typed arguments supplied to function 'latent_lognormal_lpdf':
(vector, vector, real, vector, vector, array[] real)
Available signatures:
(vector, vector, vector, vector, vector, array[] real) => real
  The third argument must be vector but got real

make: *** [/var/folders/4r/hkp4v9fn3wx044_hk8qnsjxw0000gn/T/RtmpWq657o/model-6c960467d7d.hpp] Error 1

Error: An error occured during compilation! See the message above for more information.

I think when you don't set a model then it puts sigma directly into the model rather than an intercept for sigma (which is then linked to sigma).

real<lower=0> sigma;  // dispersion parameter
...
lprior += student_t_lpdf(sigma | 3, 0, 2.5)
    - 1 * student_t_lccdf(0 | 3, 0, 2.5);

=> creates both of the bugs we are observing (putting priors on an intercept for sigma which doesn't exist, putting a real into a vector slot).

[seabbs]

Right got you.

For the latter we can overload latent_lognormal with another method that is a real and then broadcast that to a vector and call the current method.

For default priors that is interesting behaviour and I see where they are coming from. Your solution for the priors makes sense as long as we are using a prior function that is family specific.

[athowes]

That'd be nice r.e. the overloading. So if in Stan I were to write another function latent_lognormal_real_lpdf then call to latent_lognormal_lpdf inside it, then it would successfully work? (Stan would know to use the version that accepts real input?

For the prior design, yes there needs to be something family specific in there, but it also needs to be link function specific (so not only family). If someone choses to use a weird link function then the priors we think of might not make sense (when we are designing for the "mainline" links). When you say "as long as we are using a prior function that is family specific" it sounds like you don't consider an if statement to be family specific? If so, what design would you suggest?

[athowes]

(Note that this overloading approach I think could also work to fix the issue in #227)

[seabbs]

For overloading it is by type so you would have the same name for the function (but with args with different types) and it would internally do the transform from a real to a vector and then call itself. That means you don't need to change what is passed to brms etc.

Isn't the link related to the family (i.e stored as part of the object). Anyway I think it is fine for now to assign the prior to either sigma or sigma intercept depending on what is there.

[athowes]

Yep. So there are some number of dpars. Right now they're passed like vector dpar1, vector dpar2, ....

I have some concern that we need all possible combinations of real and vector for potentially unlimited length(dpars).

Like (and this is assuming there are just 2 dpar [which we can't in general]):

1. vector dpar1, vector dpar2
2. vector dpar1, real dpar2
3. real dpar1, vector dpar2
4. real dpar1, real dpar2

Yes links are stored somewhere in family. I'm not saying it's not possible, just that we need to have two lots of logic: 1. for family 2. for link within family.

Anyway I think it is fine for now to assign...

I think that if it's sigma not sigma_Intercept then it would be a prior on a quantity that isn't log transformed. So it'd be a different prior, and I haven't thought about what it should be. Hence I suggest not touching the cases I haven't thought about and make a new issue to cover more possibilities with priors.

[seabbs]

I have some concernt that we need all possible combinations of real and vector for potentially unlimited

We pass this on the fly though don't we so we only ever need to create the ones we need?

[athowes]

We pass this on the fly though don't we so we only ever need to create the ones we need?

// Here the strings
// * family
// * dpars_A
// * dpars_B
// are/have been replaced using regex

real latent_family_lpdf(vector y, dpars_A, vector pwindow,
                           vector swindow, array[] real obs_t) {
  int n = num_elements(y);
  vector[n] d = y - pwindow + swindow;
  vector[n] obs_time = to_vector(obs_t) - pwindow;
  return family_lpdf(d | dpars_B) - family_lcdf(obs_time | dpars_B);
}

Yes can just create the ones we need. Just some complicated regex required.

[athowes]

Closed via https://github.com/epinowcast/epidist/pull/255