Zhiyong Zhang's Psychometric Website
 
Home Control Panel Login Register Search Submit Logout About me Contact me
My Research
Google
Web
This Site
Home > My Research > SAS WinBUGS > SAS WinBUGS for growth curve models
SAS WinBUGS for growth curve models
2008-05-16          Read: 6221 times
Cite this page: (2008). SAS WinBUGS for growth curve models. Retrieved December 16, 2018, from http://www.psychstat.org/us/article.php/75.htm.
SAS WinBUGS for growth curve models

Please cite this using

Zhang, Z., McArdle, J. J., Wang, L., & Hamagami, F. (2008). A SAS interface for Bayesian analysis with WinBUGS. Structural Equation Modeling, 15(4), 705?C728.

TITLE 'Run WinBUGS from SAS: A Growth Model Example by Zhang & McArdle';
TITLE2 'Generate the Growth Model Data';

DATA Sim_LinGM;
* setting mathematical parameters;
  MuL=10; MuS=5;
  Sigma_e=1; Sigma_L=2; Sigma_S=1; rho=.5;
 
* setting statistical parameters;
  N = 1000; seed = 20030930;

* need to setup arrays so we can have more variables;
ARRAY y_score{5} y1-y5;

* generating raw data;
  DO _N_ = 1 TO N;
        e_L=RANNOR(seed);
  e_S=rho*e_L+SQRT(1-rho**2)*RANNOR(seed);
      L_score=MuL+Sigma_L*e_L;
  S_score=MuS+Sigma_S*e_S;
* now the indicator variables ;
       DO t = 1 TO 5;
         y_score{t} = L_score + t * S_score + Sigma_e*RANNOR(seed);
 END;
   KEEP y1-y5 L_score S_score;
   OUTPUT;
   END;
RUN;

PROC CORR DATA=Sim_LinGM;
RUN;

PROC MEANS DATA=Sim_LinGM;
VAR y1-y5 L_score S_score;
RUN;

DATA model;
INPUT model \$80.;
CARDS;/*Start of the model scripts*/
model
{#Model
   for (i in 1:N){
       LS[i,1:2]~dmnorm(Mu[1:2], Inv_cov[1:2,1:2])
           for (t in 1:T){
                y[i,t]~dnorm(MuY[i,t], Inv_sig2_e)
                MuY[i,t]<-LS[i,1]+LS[i,2]*t
            }
    }
   #Prior
   Mu[1:2]~dmnorm(Mu0[1:2], Inv_cov0[1:2,1:2])
   Mu0[1]<-0
   Mu0[2]<-0
   Inv_cov0[1,1]<-1.0E-6
   Inv_cov0[2,2]<-1.0E-6
   Inv_cov0[2,1]<-Inv_cov0[1,2]
   Inv_cov0[1,2]<-0

   Inv_cov[1:2,1:2]~dwish(R[1:2,1:2], 2)
   R[1,1]<-1
   R[2,2]<-1
   R[2,1]<-R[1,2]
   R[1,2]<-0

   Inv_sig2_e~dgamma(.001,.001)
  
   #Transform of the parameters
   MuL<-Mu[1]
   MuS<-Mu[2]
   Cov[1:2,1:2]<-inverse(Inv_cov[1:2,1:2])
   Sig2_L<-Cov[1,1]
   Sig2_S<-Cov[2,2]
   rho<-Cov[1,2]/sqrt(Cov[1,1]*Cov[2,2])
   Sig2_e<-1/Inv_sig2_e
}
;
/*end of the model scripts*/
RUN;

DATA _NULL_;
  SET model;
  FILE 'C:\SASWinBUGS\GrowthModel.txt';
  PUT model;
RUN;


/*Transform the SAS data file to the WinBUGS format
   The SAS MACRO is called here. The data are saved into the file
'c:\SASWinBUGS\SAS\test\struc.txt'
*/
%_sexport(data=Sim_LinGM, file='c:\SASWinBUGS\GrowthData.txt',
var=y1-y5);

/*Create the initial values*/
DATA _NULL_;
FILE "c:\SASWinBUGS\InitValues.txt";
PUT "list(Mu=c(0,0), Inv_cov= structure(.Data = c(1,0,0,1),.Dim=c(2,2)), Inv_sig2_e=1)";
RUN;

/*Create a script file that includes the command to be used in WinBUGS*/
DATA _NULL_;
FILENAME scrip "C:\programs\WinBUGS14\bugscript.txt";
FILE scrip;
PUT // @@

#1 "display('log')"

/* Check model
    1. using "/" instead of "\"
*/
#2 "check('C:/SASWinBUGS/GrowthModel.txt')"

/*Read in the data */
#3 "data('C:/SASWinBUGS/GrowthData.txt')"
/* Compile the model */
#4 "compile(1)"

/* Give the initial values */
#5 "inits(1, 'C:/SASWinBUGS/InitValues.txt')"

/* use "gen.inits()" to generate initial values for these not specified in inital value file*/
#6 "gen.inits()"
/*Discard the first 5000*/
#7 "update(2000)"

/* set variable names to be monitored, one command line for each variable*/
#8 "set(MuL)"
#9 "set(MuS)"
#10 "set(Sig2_L)"
#11 "set(Sig2_S)"
#12 "set(rho)"
#13 "set(Sig2_e)"
/*Set DIC*/
#14 "dic.set()"
/*Update to get the estimations*/
#15 "thin.updater(10)"
#16 "update(2000)"
/* call for stats to be output, one command for each variable  */
#17 "dic.stats()"

/* coda file to save the coda. The data will be read into SAS later.
    The index data are saved in file outputIndex.txt
    The sample data are saved in file output1.txt for chain 1*/
#18 "coda(*,'C:/SASWinBUGS/output')"

/* BUGS log will be saved as bugslog.txt */
 
#19 "save('C:/SASWinBUGS/bugslog.txt')"
/*#20 "quit()"*/
;
RUN;
/*Creat a bat file
   Can run winbugs directly. This way can save a file for further use.*/
DATA _NULL_;
FILE "c:\SASWinBUGS\run.bat";
PUT '"C:\programs\WinBUGS14\WinBUGS14.exe" /PAR bugscript.txt';
PUT 'exit';
RUN;
/*Run WinBUGS*/
DATA _NULL_;
X "c:\SASWinBUGS\run.bat";
RUN;
QUIT;

/*Read in the log file to view the DIC*/
DATA log;
INFILE "c:\SASWinBUGS\bugslog.log" TRUNCOVER ;
INPUT log \$80.;
log=translate(log," ","09"x);
RUN;

PROC PRINT DATA=log;
RUN;

/*Read in the generate samples and do the analysis
   The data are saved in work.post1
*/
TITLE 'Run WinBUGS from SAS: A Growth Model Example by Zhang & McArdle';
TITLE2 'Results from the simulation study';
%coda2sas(out=post1, infile='c:\SASWinBUGS\outputIndex.txt',
chain='c:\SASWinBUGS\output1.txt', stats=1);
QUIT;

Submitted by: johnny
Add a comment View comment Add to my favorite Email to a friend Print
If you want to rate, please login first, or click here to register. Or you can use USERNAME: guest and PASSWORD: guest to log in.
Average score 0, based on 0 comments
1 2 3 4 5 6 7 8 9 10
Latest comment (Total: 2 ) Time Author Reply
'"c:\SASWinBUGS\run.bat"' 2009-05-14 20:42 pm guest 2
should $ be added? 2009-05-14 20:27 pm guest 1
More comments...

Copyright © 2003-13 Zhiyong Zhang \'s Psychometric Website
Maintained by Zhiyong Zhang