% Warnings: - The computations are time-consuming % - MCMC estimation sometimes crashes due % numerical problems. % the dataset contains the following variables: % % id: patient id % st: state the process jumps to at the end of the current interval % beg: starting point of the current interval % end: end point of the current interal % tas: indicator for transition awake => sleep at the end of the interval % tsa: indicator for transition sleep => awake at the end of the interval % tnr: indicator for transition Non-REM => REM at the end of the interval % trn: indicator for transition REM => Non-REM at the end of the interval % cortisol: cortisol level % read the data and store it in the dataset object d dataset d d.infile using c:\temp\sleep.dat % create the binary indicator for high cortisol level d.generate cortisoldummy=0 d.replace cortisoldummy=1 if cortisol>60 % create regression objects remlreg r bayesreg b % change the delimiter to allow for linebreaks in commands delimiter=; % basic model without time-varying effects logopen, replace using c:\temp\basic.log; r.outfile = c:\temp\basic; r.mregress tas = end(baseline, nrknots=40, gridchoice=all) + id(random) : tsa = end(baseline, nrknots=40, gridchoice=all) + id(random) : tnr = end(baseline, nrknots=40, gridchoice=all) + id(random) : trn = end(baseline, nrknots=40, gridchoice=all) + id(random), family=multistate state=st lefttrunc=beg using d; logclose; % Model with time-varying effect for high cortisol level logopen, replace using c:\temp\reml.log; r.outfile = c:\temp\reml; r.mregress tas = end(baseline, nrknots=40, gridchoice=all) + id(random) : tsa = end(baseline, nrknots=40, gridchoice=all) + id(random) : tnr = end(baseline, nrknots=40, gridchoice=all) + cortisoldummy*end(baseline) + id(random) : trn = end(baseline, nrknots=40, gridchoice=all) + cortisoldummy*end(baseline) + id(random), family=multistate state=st lefttrunc=beg using d; logclose; % the same, estimated with MCMC logopen, replace using c:\temp\mcmc.log; b.outfile = c:\temp\mcmc; b.mregress tas = end(baseline, nrknots=40) + id(random) : tsa = end(baseline, nrknots=40) + id(random) : tnr = end(baseline, nrknots=40) + cortisoldummy*end(baseline, nrknots=40) + id(random) : trn = end(baseline, nrknots=40) + cortisoldummy*end(baseline, nrknots=40) + id(random), family=multistate state=st begin=beg noposteriormode iterations=30000 burnin=10000 step=20 using d; logclose; % restore usual delimiter delimiter = return;