jeudi 14 novembre 2019

Conditional distribution in WinBUGS using the step function

For the following data set, I would like to write a bug model to run from R.

 mydata
       ID        Y  x z
    1   1 5.302956  1 1
    2   1 3.358249  2 1
    3   1 4.976734  3 1
    4   1 4.290459  4 1
    5   1 0.000000  5 0
    6   2 5.975351  1 1
    7   2 6.620773  2 1
    8   2 8.045909  3 1
    9   2 7.378384  4 1
    10  2 6.908755  5 1
    11  2 8.672657  6 1
    12  2 8.284252  7 1
    13  2 8.455531  8 1
    14  2 7.415175  9 1
    15  2 8.634265 10 1
    16  3 7.356993  1 1
    17  3 6.607598  2 1
    18  3 0.000000  3 0
    19  3 0.000000  4 0
    20  3 0.000000  5 0
    21  3 0.000000  6 0
    22  3 0.000000  7 0
    23  3 0.000000  8 0
    24  3 6.398595  9 1
    25  3 6.580639 10 1
    26  4 5.525104  1 1
    27  4 2.326914  2 1
    28  4 6.744059  3 1
    29  4 6.710523  4 1
    30  4 6.876265  5 1

The Y values given z equal to 1 i,e Y|z=1 ~ LogNormal(mu_i, sig2). I am trying this model, Unfortunately it is not working

model{
  for (i in 1:N) { #N is the total obs
    z[i]~dbern(prob[i])
    logit(prob[i])<-alpha0+alpha1*x[i]+alpha2*pow(x[i],2)
    } #z  is for logistic regression 
    step(z[i]-1){        
    Y[i]~ dlnorm(mu[i], tau)
    mu[i] <- step(z[i]-1)*(zeta[ID[i]]+beta0[cls[ID[i]]]+beta1[cls[ID[i]]]*x[i]+beta2[cls[ID[i]]]*pow(x[i],2)
    ) #this loop is for Y values given z=1 follow LogN dist
  }      
  for ( k in 1:M){
    zeta[k]~dnorm(0, prec.zeta)
  }
  for ( j in 1:M){
    cls[j]~ dcat(p[])
  } #M is the number of unique ID 
  for ( b in 1:2){ #for two cluster, each beta will have two values
    beta0[b]~dnorm(0.0,1.0E-2)
    beta1[b]~dnorm(0.0, 1.0E-2)
    beta2[b]~dnorm(0.0, 1.0E-2)

  }
  alpha0~dnorm(0.0,1.0E-2)
  alpha1~dnorm(0.0,1.0E-2)
  alpha2~dnorm(0.0,1.0E-2)
  prec.zeta~dgamma(0.1,0.1)
  p[1:2]~ ddirch(alpha[])
  alpha[1] <- 1
  alpha[2] <- 1
  tau~dgamma(0.1,0.1)
  sig2<-1/tau
  sig2.zeta<-1/prec.zeta
}

z[i] follows Bernoulli and then Logistic regression. taking only Y values corresponding to z=1, Y follows LogNormal. This is a mixture of two quadratic models.
Any help is appreciated.

Aucun commentaire:

Enregistrer un commentaire