I am trying to fit a “non-linear” model. Note that it’s not really non-linear, but the goal is to fit a change-score model (a consulting request) with missing values.
Conceptually, what I want is:
outcome ~ alpha + (mi(time2) - mi(time1))*beta. That is - There are missings in time2 and time 1 that I want to model from other predictors; after filling those in, I want to compute a difference score and linearly predict the outcome from it.
The following does not work (because
mi() is not a function in Stan). (Variables changed to be generic).
cs.bf <- brmsformula(outcome ~ alpha + (mi(time2) - mi(time1))*beta,nl=TRUE,beta ~ 1, alpha ~ 1) +
brmsformula(time1|mi() ~ z1 + z2) +
brmsformula(time2|mi() ~ z1 + z2)
Is there a way to specify this model?
mi only works in linear formulas not in non-linear ones. So you need to replace the parts of the non-linear formula with some non-linear parameter on which you then specify a linear formula including the mi term(s).
I assumed as such.
Would it be as simple as:
outcome ~ alpha + (t2 - t1)*beta, t2 ~ mi(time2), time2|mi() ~ z1 + z2, ...
That seems incorrect, since t2 isn’t /modeled/ from mi(time2), but defined as mi(time2). e.g., it should conceptually be
t2 ~ 0 + 1*mi(time2) (no intercept estimated, no coefficient estimated).
outcome ~ alpha + (mi(time2) - mi(time1))*beta
would translate to
outcome ~ alpha + (t2 - t1)*beta
t2 ~ 0 + mi(time2)
t1 ~ 0 + mi(time1)
Ok, I see your point. Please disregard my former suggestions.
And that wouldn’t try to estimate a coefficient for mi(time2) on t2?
Edit: Nvm. You’re too quick :p.
Ways to solve this could be potentially be
(1) to allow to fix parameters to certain values. This may be a feature implemented in brms 3 but I haven’t thought about it in too much detail yet (though should be possible although painful for me).
(2) allow deterministic formulas that is formulas that don’t include any coefficients at all (i.e., set them to 1; or 0 for the intercept) but otherwise behave like linear formulas. Not sure if this is a something reasonable in the grant scheme of things though.
In the mean time, I assume I can generate the stan code using what you suggested, then just manually fix the coefficient(s) to be 1. That’d be the easiest method, right? (That is, easiest without breaking other brms tools).
I got it to work; thanks Paul. I know the solution isn’t ideal, but you at least confirmed that it can’t be (currently) be done within the syntax.
For anyone (including myself) in the future:
- Specify the formulas as Paul mentioned (e.g.,
alpha + beta*(t2 - t1), t1 ~ 0 + mi(time1), t2 ~ 0 + mi(time2), time1|mi() ~ preds, time2|mi() ~ preds)
- Generate the stan code and stan data via brms.
- Modify the stan code, move
vector[Ksp_cHCCscale_t1] bsp_cHCCscale_t1 = ';
// special effects coefficients
vector[Ksp_cHCCscale_t2] bsp_cHCCscale_t2 = ';
to transformed parameters (and set them to 1, obviously, as above).
4) remove the priors for those.
5) Initialize the brms model using the brm formula, with iter = 1 or something. We just need it to build the structure.
6) Sample using stan on the modified code and generated stan data.
7) Replace brmFit$fit with stanFit.
8) brmFit <- brms:::rename_pars(brmFit)
In short: Modify the code to implement the constraint; initialize a brms object; replace its fit with the stan fit on modified model code; rename the brm pars internally.