Help to write my first stan model

General
1

I would need help in order to write a specific Stan model.

The biological question

The idea of the model is modeling the number of Bones (NbBones : discret continuous count) according to Orders (3 different discret orders) and Time (continuous). And then show that The Orders 2 has a number of bones greater than the other Orders.

The Baysian model

I’m trying to write a GLM negative binomial model with zero inflated data since the mean of NbBones is not equal to the variance because of a huge number of zeros, and since within these zeros, I have a proportion p that is coming from a specific category (Bones was destroyed and then is not observable, and bones was not find because it never existed (but observable).

So far I wrote this model such as :

model
{
  for(i in 1 : N) 
  {
    observable[i] ~ dbern(p)
    NbBones[i] ~ dpois(nu[i] * theta[i] * observable[i])  
    theta[i] <-  intercept + Coef_Orders[Orders[i]] + Coef_Time*Time[i] 
    nu[i] ~ dgamma(1/alpha, 1/alpha)
  }
  p ~ dunif(0, 1)
  alpha ~ dunif(0, 10)
  Coef_Orders[1] <- 0
  Coef_Orders[2] ~ dexp(1)
  Coef_Orders[3] ~ dexp(1) 
  Coef_Time~ dexp(1)
  intercept ~ dexp(1)
}

Where in the negative binomial model with excess zeros, it is considered that only a proportion p of individuals are likely to have bones (observables). For observable individuals (p %), the number of bones follows a negative binomial distribution (allowing 0’s, as an observable individual may not have any bones). Unobservable individuals (1 -p %) have a number of bones following a Poisson distribution with an expectation of 0 (so their number of bones is necessarily zero).

My question

But i cannot manage to write it on stan, can someone help me ? (I found thid code :

stan_glm1 <- stan_glm(NbBones ~ Orders * Time,
                      data = data, family = neg_binomial_2,
                      prior = normal(0, 2.5),
                      prior_intercept = normal(0, 5),
                      seed = 12345)

But it does not include the zero inflated model…

Thanks a lot for your help and time :-)

Here are the data :

  Orders    Time NbBones
1      1 0.21771       0
2      1 0.08723       0
3      1 0.03911       0
4      1 0.08236       2
5      1 0.08744       6
6      1 0.00417       1

In dput format if you want a toy example :

So far I tried ;


    //Declare the data 
    data{
      int<lower=0> N;//
      real Time[N]; //
      real NbBones[N]; //
    }
    //Declare the parameters 
    parameters{
      real  theta; //continuous unconstrained variable 
      real<lower=0>observable;//constrained to not be negative 
    }
    //Model block - declare the likelihood and priors 
    model{
      //likelihood
      for(i in 1 : N){
        NbBones[i] ~poisson(theta[i] * observable[i]);
        theta[i] <-  intercept + Coef_orders[Orders[i]] + Coef_Time*Time[i];
        observable[i] ~ bernoulli(p);
      }
      //Priors 
      p ~ uniform(0, 1);
      alpha ~ uniform(0, 10);
      Coef_orders[1] <- 0;
      Coef_orders[2] ~ exponential(1);
      Coef_orders[3] ~ exponential(1) ;
      Coef_Time  ~ exponential(1);
      intercept ~ exponential(1);
    }



> dput(data)
structure(list(Orders = c("1", "1", "1", "1", "1", "1", "1", 
"1", "1", "3", "1", "2", "1", "2", "2", "2", "2", "2", "2", "2", 
"2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", 
"2", "3", "2", "2", "3", "2", "3", "3", "1", "3", "1", "1", "3", 
"3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
"3", "3", "3", "3", "3", "3", "1", "2", "1", "1", "1", "1", "2", 
"2", "3", "1", "3", "3", "1", "1", "2", "2", "2", "2", "2", "2", 
"2", "2", "2", "1", "2", "2", "2", "2", "2", "2", "2", "2", "1", 
"1", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", 
"2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", 
"3", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "3", "3", 
"2", "2", "3", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", 
"1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", 
"1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", 
"1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", 
"1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", 
"1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", 
"1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", 
"1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", 
"1", "1", "1", "1", "1", "1"), Time = c(0.21771, 0.08723, 0.03911, 
0.08236, 0.08744, 0.00417, 0.0162, 0.11631, 0.01712, 0.18276, 
0.0169, 0.01049, 0.04138, 0.03169, 0.05169, 0.01718, 0.11169, 
0.04654, 0.09138, 0.23242, 0.08843, 0.0185, 0.27725, 0.01363, 
0.01198, 0.06394, 0.0218, 0.05122, 0.09829, 0.07324, 0.07186, 
0.05843, 0.06693, 0.18545, 0.18037, 0.00655, 0.0122, 0.00572, 
0.01076, 0.07204, 0.00957, 0.00579, 0.02072, 0.01819, 0.02319, 
0.00948, 0.00263, 0.00698, 0.00364, 0.00502, 0.00597, 0.07171, 
0.09189, 0.01643, 0.00803, 0.00688, 0.00864, 0.01067, 0.01123, 
0.0078, 0.01015, 0.00944, 0.00755, 0.0227, 0.02841, 0.1087, 0.07557, 
0.04564, 0.03514, 0.00946, 0.011, 0.01385, 0.01837, 0.09937, 
0.10749, 0.01513, 0.01839, 0.01219, 0.01021, 0.01, 0.00844, 0.18295, 
0.22844, 0.12524, 0.14649, 0.03092, 0.02957, 0.05148, 0.08305, 
0.02372, 0.04403, 0.06083, 0.05597, 0.04202, 0.01476, 0.03361, 
0.04169, 0.00828, 0.01348, 0.06318, 0.06703, 0.02774, 0.03293, 
0.27072, 0.11012, 0.05687, 0.01636, 0.03637, 0.04128, 0.01157, 
0.01237, 0.13546, 0.01796, 0.07682, 0.0674, 0.188, 0.13317, 0.01142, 
0.01522, 0, 0, 0.01998, 0.07889, 0.02845, 0.08173, 0.01455, 0.06392, 
0.03808, 0.08862, 0.05304, 0.04987, 0.02237, 0.02612, 0.06138, 
0.05949, 0.00959, 0.01014, 0.00687, 0.00701, 0.21064, 0.01838, 
0.14369, 0.0136, 0.01749, 0.02122, 0.06394, 0.07276, 0.02882, 
0.01453, 0.05646, 0.0174, 0.0714, 0.05069, 0.03303, 0.01158, 
0.04024, 0.03759, 0.06414, 0.00566, 0.06868, 0.01201, 0.03361, 
0.03775, 0.03854, 0.00506, 0.01818, 0.03194, 0.04363, 0.01267, 
0.01589, 0.00328, 0.00965, 0.00229, 0.01062, 0.00908, 0.02851, 
0.02654, 0.00732, 0.01082, 0.02281, 0.02332, 0.04561, 0.01981, 
0.0493, 0.04335, 0.07554, 0.0074, 0.00729, 0.00631, 0.0631, 0.10638, 
0.03161, 0.02191, 0.00658, 0.00831, 0.03528, 0.01206, 0.02015, 
0.01668, 0.00712, 0.00762, 0.01905, 0.0512, 0.05175, 0.00956, 
0.02592, 0.00313, 0.00415, 0.0024, 0.00201, 0.00077, 0.00415, 
6e-04, 0.00275, 0.00109, 0.0023, 0.00156, 0.00601, 0.00223, 0.03884, 
0.00295, 0.04161, 0.02727, 0.00582, 0.00725, 0.01396, 0.00538, 
0.03965, 0.01708, 0.01833, 0.02333, 0.03832, 0.02266, 0.02942, 
0.00625, 0.0192, 0.00577, 0.01567, 0.00193, 0.01579, 0.00174, 
0.01351, 0.00787, 0.00599, 0.00326, 0.00456, 0.00482), NbBones = c(0L, 
0L, 0L, 2L, 6L, 1L, 0L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 3L, 3L, 2L, 5L, 
4L, 4L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 2L, 1L, 0L, 1L, 0L, 
1L, 0L, 0L, 0L, 2L, 2L, 2L, 2L, 4L, 3L, 0L, 0L, 0L, 2L, 2L, 3L, 
0L, 1L, 0L, 0L, 1L, 0L, 2L, 3L, 1L, 0L, 1L, 3L, 2L, 3L, 2L, 0L, 
2L, 0L, 4L, 2L, 2L, 2L, 1L, 3L, 1L, 1L, 2L, 6L, 6L, 1L, 5L, 4L, 
2L, 1L, 0L, 0L, 0L, 1L, 4L, 0L, 2L, 0L, 3L, 1L, 0L, 2L, 2L, 0L, 
5L, 4L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 
3L, 5L, 3L, 5L, 3L, 1L, 0L, 1L, 4L, 3L, 3L, 0L, 0L, 0L, 0L, 2L, 
0L, 0L, 2L, 0L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 0L, 3L, 0L, 1L, 0L, 
2L, 2L, 1L, 0L, 0L, 0L, 1L, 0L, 3L, 0L, 1L, 0L, 2L, 0L, 0L, 0L, 
2L, 1L, 0L, 1L, 1L, 0L, 2L, 2L, 0L, 0L, 0L, 0L, 3L, 4L, 0L, 0L, 
1L, 0L, 1L, 0L, 1L, 0L, 2L, 0L, 4L, 4L, 0L, 0L, 0L, 0L, 1L, 1L, 
0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 2L, 1L, 
0L, 0L, 0L, 1L, 2L, 1L, 2L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L)), row.names = c("1", "2", "3", "4", "5", 
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", 
"17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", 
"28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", 
"39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", 
"50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60", 
"61", "62", "63", "64", "65", "66", "67", "68", "69", "70", "71", 
"72", "73", "74", "75", "76", "77", "78", "79", "80", "81", "82", 
"83", "84", "85", "86", "87", "88", "89", "90", "91", "92", "93", 
"94", "95", "96", "97", "98", "99", "100", "101", "102", "103", 
"104", "105", "106", "107", "108", "109", "110", "111", "112", 
"113", "114", "115", "116", "117", "118", "119", "120", "121", 
"122", "123", "124", "125", "126", "127", "128", "129", "130", 
"131", "132", "133", "134", "135", "136", "137", "138", "139", 
"140", "141", "142", "143", "144", "145", "146", "147", "148", 
"149", "150", "151", "152", "153", "154", "155", "156", "157", 
"158", "159", "160", "161", "162", "163", "164", "165", "166", 
"167", "168", "169", "170", "171", "172", "173", "174", "175", 
"176", "177", "178", "179", "180", "181", "182", "183", "184", 
"185", "186", "187", "188", "189", "190", "191", "192", "193", 
"194", "195", "196", "197", "198", "199", "200", "201", "202", 
"203", "204", "205", "206", "207", "208", "209", "210", "211", 
"212", "213", "214", "215", "216", "217", "218", "219", "220", 
"221", "222", "223", "224", "225", "226", "227", "228", "229", 
"230", "231", "232", "233", "234", "235", "236", "237", "238", 
"239", "240", "241", "242", "243", "244", "245", "246", "247"
), class = "data.frame")
2 Likes
2

I found a solution by using the package brms

and the following code :

posterior_glm1 = brm(NbBones~  Orders * Time  , data = data,
                     family = zero_inflated_negbinomial(),
                     iter  = 6000, chains = 2, 
                     seed  = 1234,
                     cores = 5,thin=5)

If I understand it well it allows taking into account two kind of zero, those not observable and those observable, and give you at the end a zi posteriori probability which is the expected % number of unobservable zeros.

1 Like
3

Hello,

I wanted to check in, did your brms solution work?