/******************************************************************

  Fit a subject specific model to the respiratory status data
  
******************************************************************/

options ls=80 ps=59 nodate; run;

data resp; infile 'respstatus.dat';
  input id drug gender age week status;
  w=0; if week>0 then w=1;
  run;

/*****************************************************************

  Fit the logistic model using proc  glimmix, proc nlmixed, and the
  nlinmix macro
 
******************************************************************/

*  proc glimmix; 

proc glimmix data=resp method=mspl;
  class id;
  model status = w drug*w / dist=binomial link=logit solution;
  random intercept / subject=id;
run;

proc glimmix data=resp method=quad;
  class id;
  model status = w drug*w / dist=binomial link=logit solution;
  random intercept / subject=id;
run;



*  proc nlmixed;

proc nlmixed data=resp;
*  parms b1=-0.4 b2=-0.1 b3=1.9 su2=6.5;
*  parms b1=-0.4212 b2=-0.08342 b3=1.9452 su2=6.4913;
  parms b1=-0.3 b2=-0.05 b3=1.9 su2=6;
  num = exp(b1+b2*w+b3*w*drug+u);
  prob = num/(1+num);
  model status ~ binary(prob);
  random u ~ normal(0,su2) subject=id;
run;


*  glimmix macro;

%inc 'glmm800.sas'  / nosource;

%glimmix(data=resp,
    procopt=method=ml,
    stmts=%str(
     class id;
     model status = w drug*w  / solution;
     random intercept / subject=id type=un;
  ), 
  error=binomial,
  link=logit
)