/* The Heights of Schoolgirls
   ==========================
An example taken from section 12.7 of
Verbeke, G. and Molenberghs, G. (2000) Linear Mixed Models for Longitudinal Data.
New York: Springer-Verlag, pp. 183-187.
*/

%include '.\hetmixed11.sas';

data schgirls;
 input height child age group @@;
 agecat=age;
 int=1;
 cards;
111.0  1  6  1  116.4  1  7  1  121.7  1  8  1  126.3  1  9  1  130.5  1 10  1
110.0  2  6  1  115.8  2  7  1  121.5  2  8  1  126.6  2  9  1  131.4  2 10  1
113.7  3  6  1  119.7  3  7  1  125.3  3  8  1  130.1  3  9  1  136.0  3 10  1
114.0  4  6  1  118.9  4  7  1  124.6  4  8  1  129.1  4  9  1  134.0  4 10  1
114.5  5  6  1  122.0  5  7  1  126.4  5  8  1  131.2  5  9  1  135.0  5 10  1
112.0  6  6  1  117.3  6  7  1  124.4  6  8  1  129.2  6  9  1  135.2  6 10  1
116.0  7  6  2  122.0  7  7  2  126.6  7  8  2  132.6  7  9  2  137.6  7 10  2
117.6  8  6  2  123.2  8  7  2  129.3  8  8  2  134.5  8  9  2  138.9  8 10  2
121.0  9  6  2  127.3  9  7  2  134.5  9  8  2  139.9  9  9  2  145.4  9 10  2
114.5 10  6  2  119.0 10  7  2  124.0 10  8  2  130.0 10  9  2  135.1 10 10  2
117.4 11  6  2  123.2 11  7  2  129.5 11  8  2  134.5 11  9  2  140.0 11 10  2
113.7 12  6  2  119.7 12  7  2  125.3 12  8  2  130.1 12  9  2  135.9 12 10  2
113.6 13  6  2  119.1 13  7  2  124.8 13  8  2  130.8 13  9  2  136.3 13 10  2
120.4 14  6  3  125.0 14  7  3  132.0 14  8  3  136.6 14  9  3  140.7 14 10  3
120.2 15  6  3  128.5 15  7  3  134.6 15  8  3  141.0 15  9  3  146.5 15 10  3
118.9 16  6  3  125.6 16  7  3  132.1 16  8  3  139.1 16  9  3  144.0 16 10  3 
120.7 17  6  3  126.7 17  7  3  133.8 17  8  3  140.7 17  9  3  146.0 17 10  3
121.0 18  6  3  128.1 18  7  3  134.3 18  8  3  140.3 18  9  3  144.0 18 10  3
115.9 19  6  3  121.3 19  7  3  127.4 19  8  3  135.1 19  9  3  141.1 19 10  3
125.1 20  6  3  131.8 20  7  3  141.3 20  8  3  146.8 20  9  3  152.3 20 10  3
;run;

/* ----- THE TWO-COMPONENT MODEL ----- */
/* ----------------------------------- */

/* Compute estimates with A=10 */
/* --------------------------- */
%HetMixed(DATA = schgirls,
    	  SUBJECT = child,
	      REPEATED = agecat,
		  RESPONSE = height,
		  FIXED = ,
		  RANDOM = int age,
		  TYPEREP = simple,
		  TYPERAND = un,
          G = 2,
          AMIN = 10,
		  AMAX = 10,
		  ABY = 10,
		  DECISWIT = 1,
		  DECISBET = 1,
          STOPWIT = 0.000001,             
          STOPBET = 0.01, 
          MAXITER = 1000,
		  PPIEPS = 0.1,
		  ENDPOST = dvapo10);

/* Print the posterior probabilities after the last iteration of the EM algorithm. */ 
title 'Posterior probabilities after A = 10';
proc print data=dvapo10; run;

/* Compute estimates with max A = 80 */
/* --------------------------------- */
%include '.\hetmixed11.sas';
%HetMixed(DATA = schgirls,
    	  SUBJECT = child,
	      REPEATED = agecat,
		  RESPONSE = height,
		  FIXED = ,
		  RANDOM = int age,
		  TYPEREP = simple,
		  TYPERAND = un,
          G = 2,
          AMIN = 30,
		  AMAX = 80,
		  ABY = 10,
		  DECISWIT = 1,
		  DECISBET = 1,
          STOPWIT = 0.000001,             
          STOPBET = 0.0001, 
          MAXITER = 1000,
          INITPOST = dvapo10, 
		  ENDPOST = dvapo80);

/* Print the posterior probabilities after the last iteration of the EM algorithm. */ 
title 'Posterior probabilities after A = 80';
proc print data=dvapo80; run;

/* Compute estimates with max A = 160 */
/* ---------------------------------- */
%HetMixed(DATA = schgirls,
    	  SUBJECT = child,
	      REPEATED = agecat,
		  RESPONSE = height,
		  FIXED = ,
		  RANDOM = int age,
		  TYPEREP = simple,
		  TYPERAND = un,
          G = 2,
          AMIN = 80,
		  AMAX = 160,
		  ABY = 10,
		  DECISWIT = 1,
		  DECISBET = 1,
          STOPWIT = 0.000001,             
          STOPBET = 0.0001, 
          MAXITER = 1000,
          INITPOST = dvapo80, 
		  ENDPOST = dvapo160);

/* Print the posterior probabilities after the last iteration of the EM algorithm. */ 
title 'Posterior probabilities after A = 160';
proc print data=dvapo160; run;

/* Compute estimates with max A = 200 */
/* ---------------------------------- */
%HetMixed(DATA = schgirls,
    	  SUBJECT = child,
	      REPEATED = agecat,
		  RESPONSE = height,
		  FIXED = ,
		  RANDOM = int age,
		  TYPEREP = simple,
		  TYPERAND = un,
          G = 2,
          AMIN = 200,
		  AMAX = 200,
		  ABY = 10,
		  DECISWIT = 1,
		  DECISBET = 1,
          STOPWIT = 0.0001,             
          STOPBET = 0.0001, 
          MAXITER = 1000,
          INITPOST = dvapo160, 
		  ENDPOST = dvapo200);

/* Print the posterior probabilities after the last iteration of the EM algorithm. */ 
title 'Posterior probabilities after A = 200';
proc print data=dvapo200; run;

/* Compute one iteration with A = 240 */
/* ---------------------------------- */
%HetMixed(DATA = schgirls,
    	  SUBJECT = child,
	      REPEATED = agecat,
		  RESPONSE = height,
		  FIXED = ,
		  RANDOM = int age,
		  TYPEREP = simple,
		  TYPERAND = un,
          G = 2,
          AMIN = 240,
		  AMAX = 240,
		  ABY = 10,
		  DECISWIT = 1,
		  DECISBET = 1,
          STOPWIT = 0.0001,             
          STOPBET = 0.0001, 
          MAXITER = 1,
          INITPOST = dvapo200, 
		  ENDPOST = dvapo240a);

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

/* ----- THE THREE-COMPONENT MODEL ----- */
/* ------------------------------------- */
/* Initial posterior probabilities for the three-component model.
   They are subjectively based on posterior probabilities obtained
   from the two-component model.      */
data initdrie;
 input child post1 post2 post3;
 cards;
1  0.8  0.1  0.1
2  0.75 0.1  0.15
3  0.77 0.12 0.11
4  0.76 0.09 0.15
5  0.78 0.06 0.16
6  0.3  0.05 0.65
7  0.79 0.11 0.1
8  0.74 0.1  0.16
9  0.05 0.75 0.2
10 0.78 0.03 0.19
11 0.65 0.05 0.3
12 0.8  0.05 0.15
13 0.75 0.04 0.21
14 0.84 0.01 0.15
15 0.05 0.75 0.2
16 0.02 0.74 0.24
17 0.08 0.70 0.22
18 0.1  0.1  0.8
19 0.06 0.5  0.44
20 0.05 0.60 0.35
; run;

/* Compute estimates with A=10 */
/* --------------------------- */
%HetMixed(DATA = schgirls,
    	  SUBJECT = child,
	      REPEATED = agecat,
		  RESPONSE = height,
		  FIXED = ,
		  RANDOM = int age,
		  TYPEREP = simple,
		  TYPERAND = un,
          G = 3,
          AMIN = 10,
		  AMAX = 10,
		  ABY = 10,
		  DECISWIT = 1,
		  DECISBET = 1,
          STOPWIT = 0.000001,             
          STOPBET = 0.01, 
          MAXITER = 1000,
		  INITPOST = initdrie,
		  ENDPOST = tripo10);

/* Print the posterior probabilities after the last iteration of the EM algorithm. */ 
title 'Posterior probabilities after A = 10';
proc print data=tripo10; run;

/* Compute estimates with max A=80 */
/* ------------------------------- */
%HetMixed(DATA = schgirls,
    	  SUBJECT = child,
	      REPEATED = agecat,
		  RESPONSE = height,
		  FIXED = ,
		  RANDOM = int age,
		  TYPEREP = simple,
		  TYPERAND = un,
          G = 3,
          AMIN = 30,
		  AMAX = 80,
		  ABY = 10,
		  DECISWIT = 1,
		  DECISBET = 1,
          STOPWIT = 0.000001,             
          STOPBET = 0.0001, 
          MAXITER = 1000,
		  INITPOST = tripo10,
		  ENDPOST = tripo80);

/* Print the posterior probabilities after the last iteration of the EM algorithm. */ 
title 'Posterior probabilities after A = 80';
proc print data=tripo80; run;

/* Compute estimates with max A=160 */
/* -------------------------------- */
%HetMixed(DATA = schgirls,
    	  SUBJECT = child,
	      REPEATED = agecat,
		  RESPONSE = height,
		  FIXED = ,
		  RANDOM = int age,
		  TYPEREP = simple,
		  TYPERAND = un,
          G = 3,
          AMIN = 80,
		  AMAX = 160,
		  ABY = 10,
		  DECISWIT = 1,
		  DECISBET = 1,
          STOPWIT = 0.000001,             
          STOPBET = 0.0001, 
          MAXITER = 1000,
		  INITPOST = tripo80,
		  ENDPOST = tripo160);

/* Print the posterior probabilities after the last iteration of the EM algorithm. */ 
title 'Posterior probabilities after A = 160';
proc print data=tripo160; run;

/* Compute estimates with max A=200 */
/* -------------------------------- */
%HetMixed(DATA = schgirls,
    	  SUBJECT = child,
	      REPEATED = agecat,
		  RESPONSE = height,
		  FIXED = ,
		  RANDOM = int age,
		  TYPEREP = simple,
		  TYPERAND = un,
          G = 3,
          AMIN = 200,
		  AMAX = 200,
		  ABY = 10,
		  DECISWIT = 1,
		  DECISBET = 1,
          STOPWIT = 0.0001,             
          STOPBET = 0.0001, 
          MAXITER = 1000,
		  INITPOST = tripo160,
		  ENDPOST = tripo200);

/* Print the posterior probabilities after the last iteration of the EM algorithm. */ 
title 'Posterior probabilities after A = 200';
proc print data=tripo200; run;

/* Compute one iteration with A = 240 */
/* ---------------------------------- */
%HetMixed(DATA = schgirls,
    	  SUBJECT = child,
	      REPEATED = agecat,
		  RESPONSE = height,
		  FIXED = ,
		  RANDOM = int age,
		  TYPEREP = simple,
		  TYPERAND = un,
          G = 3,
          AMIN = 240,
		  AMAX = 240,
		  ABY = 10,
		  DECISWIT = 1,
		  DECISBET = 1,
          STOPWIT = 0.0001,             
          STOPBET = 0.0001, 
          MAXITER = 1,
          INITPOST = tripo200, 
		  ENDPOST = tripo240a);
data girls.tripo240a;
 set tripo240a;
run;