Non-linear specification (variable is not a valid distributional or non-linear parameter)

Interfaces brms
1

I am trying to use non-linear syntax to estimate several parameters in a model with brms where I am treating the outcome as a Beta proportion and trying to infer several other parameters.

The general specification is:

\begin{align*} S_{r,t} & \sim Beta(\mu, \kappa) \\ \mu & = I \\ I_{p,d} & = 1 + B_{p,d} \\ B_{p,d} & = a * R^b + c \\ \end{align*}

S_{r,t}, C, and R are observed, so I want to estimate I, B, a, b, and c

I first tried

bform <- bf(
  Srt ~ I,
  I ~ 1 + (C * B )+ (1|r) + (1|t),
  B ~ 1 + (a * R^b + c),
  a ~ 1,
  b ~ 1,
  c ~ 1,
  nl = TRUE
  )

But, that gives me a power error:
Error in terms.formula(as.formula(x)) : invalid power in formula

I assumed this was do to how R is interpreting the ^

So I eliminated b and assumed a 2

bform <- bf(
  seroprev ~ I,
  I ~ 1 + (C * B + (1|r) + (1|t),
  B ~ 1 + (a * R^2 + c),
  a ~ 1,
  c ~ 1,
  nl = TRUE
  )

This gives a different error:
Error: The parameter 'B' is not a valid distributional or non-linear parameter. Did you forget to set 'nl = TRUE'?

I have set nl = TRUE.

Is it possible to estimate multiple parameters this way? Or should I be “injecting” these parameters into the transformed parameters block somehow?

2

When you eliminate the b and assume a 2, the model is no longer non-linear (R is observed, so R^2 is a constant).

You may well be right that the original problem has to do with how R interprets ^ in a model formula (i.e. as a degree of interaction). However,^ works in the example A real-world non-linear model in the relevant brms vignette here:
https://cran.r-project.org/web/packages/brms/vignettes/brms_nonlinear.html

If you could post a reproducible example of the original failure that yields Error in terms.formula(as.formula(x)) : invalid power in formula (running the call to bf from your post doesn’t reproduce the error for me) then I might be able to help. But I’ll also go ahead and pass you up the food chain by paging @paul.buerkner.

3

Thank you very much!

Here’s is an (unprincipled) reprex:

# load packaged
library(dplyr)
library(brms)

# Some fake data
Data <- expand.grid(r = factor(c(1:5)),
                   t = factor(c(1:8))
                   ) %>%
        mutate(C = round(rnorm(nrow(.),500, 150)),
               D = round(rnorm(nrow(.),10000, 1000)),
               R = C/D,
               Srt = rnorm(nrow(.),0.5, 0.15))

# Non-linear formula
bform <- bf(
        Srt ~ I,
        I ~ 1 + (C * B) + (1|r) + (1|t),
        B ~ 1 + (a * R(...^b) + c),
        a ~ 1,
        b ~ 1,
        c ~ 1,
        nl = TRUE
)

# Get available priors
get_prior(bform,
          data = Data,
          family = Beta())

This gives the following error:
Error: The parameter 'B' is not a valid distributional or non-linear parameter. Did you forget to set 'nl = TRUE'?

4

I see 2 problems:

(a) you can’t mix linear and non-linear formula syntax, that is I ~ B + (1 | r) for example is invalid if B is supposed to be a non-linear parameter.

(b) non-linear formula except for the main one need to be wrapped in nlf() to work. that is for example nlf(B ~ b^2)

5

Thank you, @paul.buerkner !

Where is the correct place to enter grouping variables for the non-linear parameters? In this case I want to group I, C, and B by r and t and leave a, b, and c independent of r and t.

I tried:

bform <- bf(
        Srt ~ I + (1|r) + (1|t),
        nlf(I ~ C * B),
        nlf(B ~ a * R^b + c),
        nlf(a ~ 1),
        nlf(b ~ 1),
        nlf(c ~ 1),
        nl = TRUE
)

get_prior(bform,
          data = Data,
          family = Beta())

But this gives me:
Error in vector("list", nrow(re)) : invalid 'length' argument

6

group level terms need to be entered in linear formulas only.

7

and you need to distinguish between a non-linear parameter (e.g. c) and a non-linear formula. for example you don’t need nlf() around c ~ 1 because the formula itself is linear.

1 Like
8

Thank you very much, @paul.buerkner !

This seems to be approaching the correct syntax (at least no immediate errors):

bform <- bf(
        Srt ~ I + (1|r) + (1|t),
        nlf(I ~ C * B),
        nlf(B ~ a * R^b + c),
        a ~ 1,
        b ~ 1,
        c ~ 1,
        nl = TRUE
)

get_prior(bform,
          data = Data,
          family = Beta())
9

Is it possible to set priors on group-level effects in non-linear models?

Make some fake data:

Data <- expand.grid(r = factor(c(1:5)),
                   t = factor(c(1:8))
                   ) %>%
        mutate(C = round(rnorm(nrow(.),500, 150)),
               D = round(rnorm(nrow(.),10000, 1000)),
               R = C/D,
               Srt = rnorm(nrow(.),0.5, 0.15))

Model formula:

bform <- bf(
        Srt ~ I,
        I ~ (1|r) + (1|t),
        nlf(I ~ C * B),
        nlf(B ~ a * R^b + c),
        a ~ 1,
        b ~ 1,
        c ~ 1,
        nl = TRUE
)

Get available priors:

prior_list <- get_prior(bform,
                        data = Data,
                        family = Beta())

Set priors on nlpars:

nl_priors <- c(
        set_prior("normal(0, 10)",
                  class = "b",
                  nlpar = "a"),
        set_prior("normal(0, 10)",
                  class = "b",
                  nlpar = "b"),
        set_prior("normal(0, 10)",
                  class = "b",
                  nlpar = "c")
)

Compile and run model:

NL_Mod <- brm(bform,
                Data,
                family = Beta(), 
                prior = nl_priors,
                inits = "random",
                iter = 2000,
                warmup = 1000,
                chains = 4,
                cores = ncores,
                backend = "cmdstan",
                normalize = FALSE,
                save_pars = save_pars(all = TRUE),
                control = list(adapt_delta = 0.9,
                               max_treedepth = 12)
)

Am I still mis-specifying in some way?

11

@paul.buerkner , I have the model running now, at least with the fake data.

However, I only get estimates for population level-parameters (a, b, and c) and I is a primary estimate of interest. How might I specify the model to return the estimates of I?

12

You can use posterior_epred(..., nlpar = "I")

1 Like