Error message when sampling from a multinomial distribution (hierarchical model)

Modeling
1

Hi everyone,

I am having difficulty using a multinomial distribution as part of a hierarchical meta-analysis model.

The line the issue pertains to is: x1[i,1:4] ~ multinomial(disease[i], p1[i,1:4]);. The error I get is: No matches for: Available argument signatures for multinomial: real return type required for probability function.

Can anyone help with this please? Thank you in advance.

data {
  int<lower = 0> Ns;
  int<lower = 0> x1[Ns, 4];
  int<lower = 0> x0[Ns, 4];
  int<lower = 0> disease[Ns];
  int<lower = 0> nodisease[Ns];
}

parameters {
  real rr;
  vector[6] b;
  vector<lower = 0, upper = 5>[6] tau;
  vector[6] z[Ns];
}

transformed parameters {
  matrix[6,6] Tau;
  matrix[6,6] L;
  vector<lower = 0, upper = 1>[4] p1[Ns];
  vector<lower = 0, upper = 1>[4] p0[Ns];
  real<lower = 0, upper = 1> p1_dot1[Ns];
  real<lower = 0, upper = 1> p1_1dot[Ns];
  real<lower = 0, upper = 1> p0_dot0[Ns];
  real<lower = 0, upper = 1> p0_0dot[Ns];
  vector[6] mu[Ns];
  real<lower = -1, upper = 1> rho;
  rho = rr*2 - 1;
  for (j in 1:6) {
    Tau[j,j] = tau[j]^2;
    for (k in (j+1):6) {
      Tau[j,k] = tau[j] * tau[k] * rho;
      Tau[k,j] = Tau[j,k];
    }
  }
  L = cholesky_decompose(Tau);
  for (i in 1:Ns) {
    //between-studies model
    mu[i] = b + L*z[i];
    //diseased group
    p1[i,4] = inv_logit(mu[i,3]);
    p1_dot1[i] = inv_logit(mu[i,2]);
    p1_1dot[i] = inv_logit(mu[i,1]);
    p1[i,2] = p1_dot1[i] - p1[i,4];
    p1[i,3] = p1_1dot[i] - p1[i,4];
    p1[i,1] = 1 - p1[i,2] - p1[i,3] - p1[i,4];
    //non-diseased group
    p0[i,1] = inv_logit(mu[i,6]);
    p0_dot0[i] = inv_logit(mu[i,5]);
    p0_0dot[i] = inv_logit(mu[i,4]);
    p0[i,2] = p0_0dot[i] - p0[i,1];
    p0[i,3] = p0_dot0[i] - p0[i,1];
    p0[i,4] = 1 - p0[i,1] - p0[i,2] - p0[i,3]; }
}

model {
  //priors
  b ~ normal(0,10);
  tau ~ uniform(0,2);
  rr ~ beta(1.3,1.7);
  for (i in 1:Ns) {
    z[i] ~ std_normal();
  //likelihoods
  x1[i,1:4] ~ multinomial(disease[i], p1[i,1:4]);
  x0[i,1:4] ~ multinomial(nodisease[i], p0[i,1:4]); }
}
2

The beta distribution is constraint. We have to constrain the parameter too:
real<lower=0, upper=1> rr;
The constraints of tau are in the range of [0,5]. (No need to explicitly same from uniform, can be removed from the model)
tau ~ uniform(0,5);

However the recommended way to build L in Stan is is to sample from lkj_corr_cholesky with parameter 1.0 in your case and diag_pre_multiply by Tau. See Stan User manual for details.

I haven’t check the p0, p1 part, verify that they form a simplex.

3

Thank you for the quick response. I have implemented the changes you have suggested and although this has improved my code, the error message persists.

I note in the Stan Functions Reference that the multinomial distribution sampling statement is specified as: y ~ multinomial(theta), rather than y ~ multinomial(N, theta), where N is an integer, as I have written. Indeed, when I change the problematic line to x1[i,1:4] ~ multinomial(p1[i,1:4]); my model now runs.

I am concerned however that this integer, in my example the number of patients in the diseased group, is no longer incorporated in the sampling statement. It is my understanding that this value is a required parameter to sample from the multinomial distribution, similar to y ~ binomial(N, theta) for the binomial case?

4

From User Manual Stan multinomial

\sum_{k=1}^K{y_k} = N

In your Stan model N is in sum(x1[i]), sum(x0[i]). So it’s ok.
I don’t understand why rho is real in your model. You may have a good reason that the correlation between your y_k are same though.