/* Setup */
options nocenter nolabel nodate formchar="|----|+|---+=|-/\<>";
filename tmpsim url 'http://statistics.ats.ucla.edu/stat/data/simu.csv';
DATA Simu;
  INFILE tmpsim DLM=',' firstobs=2;
  input Y Ystar Ycensl Ycensr Xc1 Xc2 Xc3 Xc4 Xi1 Xi2 cat1 cat2 ID;
RUN;

/* Model set 1 */
*true model;
proc mixed data = Simu noclprint;
  class ID;
  model Ystar = Xc1 Xc2 Xi1 / solution;
  random int / subject=ID;
run;
*naive model;
proc mixed data = Simu noclprint;
  class ID;
  model Ycensr = Xc1 Xc2 Xi1 / solution;
  random int / subject=ID;
run;
*tobit model (right censored);
proc nlmixed data=Simu XTOL=1E-12 method=GAUSS qpoints=100;
parms sigma2_u=.7 sigma2=75 beta0=7.3 beta1=-.21 beta2=.96 beta3=3;
bounds sigma2_u sigma2 >= 0;
pi = constant('pi');
mu = beta0 + b_0j + beta1*Xc1 + beta2*Xc2 + beta3*Xi1;
if Ycensr < 16 then ll = (1 / (sqrt(2*pi*sigma2))) * exp( -(Ycensr-mu)**2 / (2*sigma2) );
if Ycensr >= 16 then ll = 1 - probnorm( (Ycensr - mu) / sqrt(sigma2) );
L=log(ll);
model Ycensr ~ general(L);
random b_0j ~ normal(0, sigma2_u) subject=ID;
run;
quit;


/* Model set 2 */
*true model;
proc mixed data = Simu noclprint;
  class ID;
  model Ystar = Xc1 Xc2 Xi1 / solution;
  random int Xc2 / subject=ID type = un;
run;
*naive model;
proc mixed data = Simu noclprint;
  class ID;
  model Ycensr = Xc1 Xc2 Xi1 / solution;
  random int Xc2 / subject=ID type = un;
run;
*tobit model;
proc nlmixed data=Simu XTOL=1E-8 method=GAUSS qpoints=20 TECH=TRUREG MAXITER=300;
parms sigma2_u1=.75 sigma2_u2=.38 sigma2_u3=.20 sigma2=74 beta0=3.79 beta1=.20 beta2=.98 beta3=3.36;
bounds sigma2_u1 sigma2_u3 sigma2 >= 0;
pi = constant('pi');
mu = beta0 + b_0j + beta1*Xc1 + beta2*Xc2 + beta3*Xi1 + b_1j*Xc2;
if Ycensr < 16 then ll = (1 / (sqrt(2*pi*sigma2))) * exp( -(Ycensr-mu)**2 / (2*sigma2) );
if Ycensr >= 16 then ll = 1 - probnorm( (Ycensr - mu) / sqrt(sigma2) );
L=log(ll);
model Ycensr ~ general(L);
random b_0j b_1j ~ normal([0, 0], [sigma2_u1, sigma2_u2, sigma2_u3]) subject=ID;
run;
quit;



/* Model set 3 */
*true model;
proc mixed data = Simu noclprint;
  class ID;
  model Ystar = Xc1 Xc2 Xi1 / solution;
  random int / subject=ID;
run;
*naive model;
proc mixed data = Simu noclprint;
  class ID;
  model Ycensl = Xc1 Xc2 Xi1 / solution;
  random int / subject=ID;
run;
*tobit model (right censored);
proc nlmixed data=Simu XTOL=1E-12 method=GAUSS qpoints=100;
parms sigma2_u=.7 sigma2=75 beta0=7.3 beta1=-.21 beta2=.96 beta3=3;
bounds sigma2_u sigma2 >= 0;
pi = constant('pi');
mu = beta0 + b_0j + beta1*Xc1 + beta2*Xc2 + beta3*Xi1;
if Ycensl > -5 then ll = (1 / (sqrt(2*pi*sigma2))) * exp( -(Ycensl-mu)**2 / (2*sigma2) );
if Ycensl <= -5 then ll = probnorm( (Ycensl - mu) / sqrt(sigma2) );
L=log(ll);
model Ycensl ~ general(L);
random b_0j ~ normal(0, sigma2_u) subject=ID;
run;
quit;