##############################################
##    GMM with non-ignorable missingdata    ##   
##############################################
     
model {
      ########################
      ##         model        
      ########################
      for( i in 1:N )  {      
          # --------- Model the latent mixture, parameter: lambda.---------
          z[i] ~ dcat(lambda[i,1:2])
              lambda[i,1] <- phi(zeta0+zeta1*x[i])
              lambda[i,2] <- 1-lambda[i,1]
          # --------- Model the growth curve, parameters: L, S, pre.etasigma, pre.phi ---------
          for(t in 1:Time) {
              y[i,t] ~ dnorm(mu[i,t], pre.phi[z[i]])
              mu[i,t]<-eta[i,1]+(t-1)*eta[i,2]
          }
          eta[i,1:2]  ~ dmnorm(etamu[z[i],1:2], pre.etaphi[z[i],1:2,1:2])
          # --------- Model the missingness, parameters:gamma0, gamma1 ---------     
          for (t in 1:Time) {
              R[i,t] ~ dbern(p[i,t])
              p[i,t] <- phi(gamma0[t]+gamma1[t]*(z[i]-1)+gamma2[t]*x[i]) 
          }
      } 
      ######################## 
      #      priors parameters 
      ######################## 
      # --------- etamu=beta --------- 
      # Intercept 
      etamu[1,1]~dnorm(0, 1.0E-6)
      etamu[2,1]~dnorm(0, 1.0E-6)  
      # Slope 
      etamu[1,2]~dnorm(0, 1.0E-6)
      etamu[2,2]~dnorm(0, 1.0E-6)  
      # --------- D (pre.etaphi) --------- 
      pre.etaphi[1,1:2,1:2] ~ dwish(V[1:2,1:2],2) 
      pre.etaphi[2,1:2,1:2] ~ dwish(V[1:2,1:2],2) 
       V[1,1]<-1 
       V[2,2]<-1 
       V[1,2]<-0 
       V[2,1]<-V[1,2] 
       Cov[1,1:2,1:2] <-inverse(pre.etaphi[1,1:2,1:2]) 
                  varI[1]  <- Cov[1,1,1] 
                  varS[1] <- Cov[1,2,2] 
                  cov[1]  <- Cov[1,1,2]
       Cov[2,1:2,1:2] <-inverse(pre.etaphi[2,1:2,1:2]) 
                varI[2]    <- Cov[2,1,1] 
                varS[2]  <- Cov[2,2,2] 
                cov[2]   <- Cov[2,1,2] 
      # --------- pre.phi --------- 
      pre.phi[1] ~ dgamma(0.001, 0.001) 
      pre.phi[2] ~ dgamma(0.001, 0.001) 
      phi[1] <- 1/pre.phi[1] 
      phi[2] <- 1/pre.phi[2] 
      # --------- zeta0 & zeta1 --------- 
      zeta0 ~dnorm(0, 1.0E-6) 
      zeta1 ~dnorm(0, 1.0E-6) 
      # --------- gamma --------- 
      for (t in 1:Time) {
          gamma0[t] ~ dnorm(0,1.0E-6) 
          gamma1[t] ~ dnorm(0,1.0E-6) 
          gamma2[t] ~ dnorm(0,1.0E-6) 
      }
      # --------- Re-parametric --------- 
      par[1] <- etamu[1,1] 
      par[2] <- etamu[2,1] 
      par[3] <- etamu[1,2] 
      par[4] <- etamu[2,2] 
      par[5] <- varI[1] 
      par[6] <- varI[2] 
      par[7] <- varS[1] 
      par[8] <- varS[2] 
      par[9] <- cov[1] 
      par[10]<- cov[2] 
      par[11]<- phi[1] 
      par[12]<- phi[2] 
      par[13]<- zeta0  
      par[14]<- zeta1  
      for (j in 1:Time) { 
          par[14+j]          <- gamma0[j]    #par[15],[16],[17],[18],[19] 
          par[14+Time+j]    <- gamma1[j]    #par[20],[21],[22],[23],[24] 
          par[14+2*Time+j] <- gamma2[j]    #par[25],[26],[27],[28],[29] 
      }  
}