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