stan code for a logistic regression spline model

[Pawel]

Hi, 

The data block (as generated by stancode) for a simple logistic regression spline model 

brm(y ~ s(x), data=list(y=y,x=x), family = "bernoulli")

is

data { 

  int<lower=1> N;  // total number of observations 

  int Y[N];  // response variable 

  int<lower=1> K;  // number of population-level effects 

  matrix[N, K] X;  // population-level design matrix 

  // data of smooth s(x)

  int nb_1;  // number of bases 

  int knots_1[nb_1]; 

  matrix[N, knots_1[1]] Zs_1_1; 

  int prior_only;  // should the likelihood be ignored? 

} 

How do I find what value of K was used?

How X is generated?

What is the value of nb_1?

and how the matrix Zs_1_1 is generated?

Thanks !

:)

pawel

[Travis Riddle]

Hey Pawel,

You can check the values for the data passed to a stan model like this with:

make_standata(y ~ s(x), data=list(y=y,x=x), family = "bernoulli", data=df)

X is the design matrix, and so is of dimension [n_observations, n_non_varying_parameters]. I'm not positive how it works for splines, but I think Z_1_1 are going to represent deflections for the splines. Again, you could check both of these with the make_standata function

[Pawel]

Thanks Travis !

This is exactly what I needed

:)
pawel

--
You received this message because you are subscribed to a topic in the Google Groups "brms-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/brms-users/BRgnn-7wSnM/unsubscribe.
To unsubscribe from this group and all its topics, send an email to brms-users+unsubscribe@googlegroups.com.
To post to this group, send email to brms-...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/5f4129c4-68a5-47d9-b8ff-295eaa8c4777%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Pawel]

Hey Travis, 

I would like to generate in stan's generated quantities block distribution of predicted values of s(x), where x is a vector of 100 evenly distributed values in the range of the original x variable. 

In the model block of the stan code generated by stancode for the model y ~ s(x), family = "bernoulli" the probability of y=1 in logit scale is 

mu = Xc * b + Zs_1_1 * s_1_1 + temp_Intercept;

so for the variable x2 of the uniformly distributed values in the range of the original x variable one would have to replace Xc and Zs_1_1 by the corresponding design matrices of x2. How to generate those matrices? 

Thanks !

:)

pawel

[Paul Buerkner]

Splines are complicated and everything related to spline data preparation in brms in handled by the mgcv package.

I wouldn't recommend doing this in Stan generated quantities block, yourself.

Rather, I recommend applying the marginal_smooths method to a fitted model object.

--
You received this message because you are subscribed to the Google Groups "brms-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+unsubscribe@googlegroups.com.

To post to this group, send email to brms-...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/97602082-f342-4ec3-ae3a-f64acd5fd5ce%40googlegroups.com.

[Pawel]

Hi Paul,

The reason I am trying to do this in stan is that I am working on a multi-level spline model that I don't know how to set up in brms. 

The simple version of the model is 

y ~ s(x), family = "bernoulli"

which can be trivially implemented in brms, where x is the log ratio of relative abundances of some bacteria in human body. Those log ratios are measured with some measurement error. In practice, one samples each bacteria type in a given biological sample and the number of bacteria of all types in the i-th sample is n[i], the number of bacteria 1 is x1[i] and bacteria 2 in sample i is x2[i]. Thus, lp1[i] and lp2[i] are the ML estimate of relative abundances of bacteria 1 and 2 in the i-th sample. These ML estimates of log relative abundandes are assumed to come from a normal distribution centered at true log relative abundances tlp1[i] and tlp2[i] for bug 1and 2, respectively. Standard deviations, denoted above as sigma[i], depend on ML estimates and are passed to the model as data. 

A spline model accounting for measurement errors has the following structure

 y ~ bernoulli_logit( mu );

for( i in 1:N )

{

    x1[i] ~ binomial(n[i], exp(lp1[i]));

    lp1[i] ~ normal(tlp1[i], sigma[i]);

    x2[i] ~ binomial(n[i], exp(lp2[i]));

    lp2[i] ~ normal(tlp2[i], sigma[i]);

}

x = tlp1 - tlp2;

mu = "spline expression for x"

Is there a way to set it up in brms?

Thanks !

:)

pawel

To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+...@googlegroups.com.

[Paul Buerkner]

I understand, thanks for the explanation. Unfortunately, this cannot be specified in brms, which requires fixed data for the spline covariates.

Maybe the following case study explaining to fit b-splines in Stan can be of help:

https://github.com/milkha/Splines_in_Stan/blob/master/splines_in_stan.pdf

To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+unsubscribe@googlegroups.com.

To post to this group, send email to brms-...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/9d0522ef-6120-4fa3-ae2b-0da84df57c6b%40googlegroups.com.