options ls=72 ps=72;

/*------------------------------------------------------*/
/*                                                      */
/* Proc Glimmix to fit subject-specific (random effect) */
/* model  for the indonesian children infection disease */
/* data                                                 */
/*                                                      */
/*------------------------------------------------------*/

data indon; 
  infile "indon.dat";
  input id infect intercep age xero cosv sinv sex height stunted
  visit baseage season visitsq;
  time = age - baseage;
run;

title "Random intercept model for infectious disease data using proc GLIMMIX";
title2 "assuming conditional Bernoulli model";
proc glimmix data=indon method=MSPL;
  class id;
  model infect = season xero age time sex height 
    / dist=bin link=logit s;
  random int / subject=id;
run;

title "Random intercept model for infectious disease data using proc NLMIXED";
title2 "assuming conditional Bernoulli model";
proc nlmixed data=indon; * qpoints=15;
  parms beta0=-2 beta1=-0.04 beta2=0.6 beta3=-0.03  beta4=0 
        beta5=-0.4 beta6=-0.05 theta=0.5;
  bound theta>0;

  eta = beta0 + beta1*season + beta2*xero + beta3*age + beta4*time + beta5*sex 
        + beta6*height + b;
  expeta = exp(eta);
  p = expeta/(1+expeta);

  model infect ~ binary(p);
  random b ~ normal(0, theta) subject=id;
run;