====== Manuscript ======
:lab:people:bayesianmixturemissing_submission3.pdf|

====== Simulation Study ======

==== Design ====

Generating Data with parameters: 

==== CRC settings ====

{{:lab:people:sim-crc.zip}}

which includes 2 folders: (1) codes (runsim.sas, script.txt, model.txt, inits.txt, data.txt, qsub.sh. Most of them need edit)
(2) OpenBUGS307 (lib:OpenBUGS.sh, OpenBUGSCli.sh; bin)

  * FTP file transfer: connect: newcell.crc.nd.edu; go to: /dscratch/zlu; mkdir: Private; and then upload 2 folders downloaded; Disconnect.

  * FTP ssh client: 
  - ssh newcell (psw); 
  - cd /dscratch/zlu/Private/OPenBUGS307/bin ;  ls ;  
  - chmod 755 *.* ;  
  - ./OpenBUGS.sh ; 
  - modelQuit()  ; 
  - cd ../..  (or cd /dscratch/zlu/Private/codes) ;  ls ;  
  - nano runsimu.sas (and edit other files. I can also edit these files before uploading.);  
  - qsub qsub.sh (to submit files; checking use: qstat -u zlu; delete using: qdel  ); 
  - exit (or logout).


emacs:
  * create a file: ^x^f 
  * save: ^x^s
  * quit: ^x^c

1. R:
<code>
#!/bin/csh
#$ -M zlu@nd.edu
#$ -m abe
#$ -r y

module load R
R CMD BATCH  input_R_file.r  output_R_file.out
</code>


2. SAS:
<code>
#!/bin/csh
#$ -M zlu@nd.edu
#$ -m abe

module load sas
sas job.sas
</code>

==== Results ====

====== Real Data Analysis ======

==== Data Normality Transformation ====

  * Using power transformation to make it relatively normal: Y^p 

  * After transformation, if variances are too large to converge, there are three ways to make it converge:
  - change initial values
  - (Y/10)^p   (here we choose this way and p=1.5)
  - (Y^p)-E(Y^p)+E(Y) 



==== OpenBUGS code ====

<code>
########################
##         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]) 
              # p[i,t] <- phi(gamma0[t]*(2-z[i])+gamma1[t]*(z[i]-1))
          }
      } 
      ######################## 
      #      priors parameters 
      ######################## 
      # --------- etamu=beta --------- 
      # Intercept (put constrain:etamu[1,1] > etamu[2,1]):etamu[1,1]~dnorm(0,0.001)I(etamu[2,1],)
      etamu[1,1]~dnorm(0, 0.001)
      etamu[2,1]~dnorm(0, 0.001)  
      # Slope (put constrain: etamu[1,2] > etamu[2,2]): etamu[1,2]~dnorm(0, 0.001)I(etamu[2,2],)
      etamu[1,2]~dnorm(0, 0.001)
      etamu[2,2]~dnorm(0, 0.001)  
      # --------- 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]
                  rho[1]   <- Cov[1,1,2]/sqrt(Cov[1,1,1]*Cov[1,2,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] 
                rho[2]    <- Cov[2,1,2]/sqrt(Cov[2,1,1]*Cov[2,2,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, 0.001) 
      zeta1 ~dnorm(0, 0.001) 
      # --------- 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] 
      }  
} 


list(etamu=structure(.Data=c(1,1,1,1), .Dim=c(2,2)), pre.etaphi=structure(.Data=c(1, 0, 0, 1, 1, 0, 0, 1), .Dim=c(2,2,2)), pre.phi=c(.5,.5), zeta0=0, zeta1=1, gamma0=c(0,0,0,0,0), gamma1=c(0,0,0,0,0), gamma2=c(0,0,0,0,0) )

list( N = 1510 , Time = 5 , y = structure(.Data = c( 19.72351, 27.45125, 
26.10502, 30.21048, 30.67885, NA, 25.66131, 26.10502, 27.45125, 27.90498, NA, 
...
</code>

==== Results ====