/*************************************************************************
This program contains analyses in "Multivariate Tests for Phase Capacity"
by Taylor Lewis in Survey Research Methods Vol 13, Issue 1.

Last updated: 3.13.19.

Questions can be directed to author at thlewis@rti.org.
**************************************************************************/


*** code to reproduce Figure 1;
proc iml;
delta={0.21, 0.53, 0.44, 0.19, 0.20, 0.29, 0.14, 0.00, 0.17};
X={1 0 0 0 0 0,
   1 0 0 1 0 0,
   1 0 0 2 0 0,
   0 1 0 0 0 0,
   0 1 0 0 1 0,
   0 1 0 0 2 0,
   0 0 1 0 0 0,
   0 0 1 0 0 1,
   0 0 1 0 0 2};

* print outcome vector/matrix to verify created properly;
print delta;
print X;

* estimate model parameters;
beta_hat=INV(X`*X)*X`*delta;

print beta_hat;
/* output of beta_hat vector:
0.2783333
0.1766667
0.0966667
0.115
0.05
0.015
*/
quit;

data toplot;
 input Y X;
if _N_=1 then y_hat=0.2783333;
if _N_=2 then y_hat=0.2783333 + 1*0.1766667;
if _N_=3 then y_hat=0.2783333 + 2*0.1766667;
if _N_=4 then y_hat=0.0966667;
if _N_=5 then y_hat=0.0966667 + 1*0.115;
if _N_=6 then y_hat=0.0966667 + 2*0.115;
if _N_=7 then y_hat=0.05;
if _N_=8 then y_hat=0.05 + 1*0.015;
if _N_=9 then y_hat=0.05 + 2*0.015;
if 1 le _N_ le 3 then estimate='Item 1';
  else if 4 le _N_ le 6 then estimate='Item 2';
  else if 7 le _N_ le 9 then estimate='Item 3';
if 1 le _N_ le 3 then pattern='solid      ';
  else if 4 le _N_ le 6 then pattern='mediumdash';
  else if 7 le _N_ le 9 then pattern='shortdash';
if 1 le _N_ le 3 then pattern='solid      ';
  else if 4 le _N_ le 6 then pattern='mediumdash';
  else if 7 le _N_ le 9 then pattern='shortdash';
if 1 le _N_ le 3 then GROUPMS='1';
  else if 4 le _N_ le 6 then GROUPMS='2';
  else if 7 le _N_ le 9 then GROUPMS='3';
datalines;
0.21 0
0.53 1
0.44 2
0.19 0
0.20 1
0.29 2
0.14 0
0.00 1
0.17 2
;
run;

proc sgplot data=toplot;
  series  Y=y_hat X=X / group=estimate grouplp=pattern lineattrs=(color=black);
  scatter Y=Y X=X     / group=estimate markerchar=GROUPMS;
label estimate='Point Estimate Relative Percent Change Trend Line';
xaxis label='Wave' values=(0 1 2) valuesdisplay=("k - 2" "k - 1" "k");
yaxis label='Point Estimate Relative Percent Change' values=(0 to 0.7 by 0.1);
run;
quit;





****** BEGIN: simulation study ***************************************************************************************;
*** create condition-specific, supplemental data sets housing EMPID and simulated wave for each of the 1000 iterations;
* 1) create simulated waves for the condition where wave is unrelated to outcome;
data waves_cond1;
  set work (keep=EMPID);
  do iteration=1 to 1000;
  wave=rantbl(54432,0.250719466,0.17451364,0.15010935,0.11039484,0.07056521,0.05870841,0.05053528,0.04385863,0.04696673);
  output;
  end;
run;
proc sort data=waves_cond1; by EMPID; run;


* 2) create simulated waves for the condition where wave is related to outcome;
* first step is to assign a group variable within each class based on the ranked sum of positive responses for items in JS index;
* read in the raw responses;
data combined (keep=EMPID sum_pos);
  set responses (where=(agency in('AGY1' 'AGY2' 'AGY3')));
attrib _ALL_ label=' ';
format _ALL_;
sum_pos=sum(input(Q4,1.),input(Q5,1.),input(Q13,1.),input(Q63,1.),input(Q67,1.),input(Q69,1.),input(Q70,1.)) 
        + ranuni(12331); * <-- fudge factor to ensure no ties during the pending ranking step;
run;
proc sort data=combined; by EMPID; run;
proc sort data=work;     by EMPID; run;
data junk;
  merge work     (in=a)
        combined (in=b);
  by EMPID;
if a;
run;
proc sort data=junk; by class; run;
proc rank data=junk out=junk groups=2;
  by    class; *  control on weighting class prior to ranking;
  var   sum_pos;
  ranks group;
run;
proc sort data=junk; by class group; run;
/*
* inspect the distribution of group within class to ensure it is about a 50/50 split;
proc freq data=junk;
  by class;
  table group;
run;
* looks good;
*/
data waves_cond2;
  set junk (keep=EMPID group);
  do iteration=1 to 1000;
  * assign such that less positive sample units tend to respond sooner;
  if group=0 then wave=rantbl(8722,0.345343617,0.207206170,0.115114539,0.092091631,0.046045816,0.046045816,0.036836652,0.034534362,0.039138943);
		     else wave=rantbl(7668,0.156095315,0.141821112,0.185104179,0.128698055,0.095084609,0.071371014,0.064233913,0.053182917,0.054794521);
  output;
  end;
run;
proc sort data=waves_cond2; by EMPID; run;


/*
* verify wave distributions are marginally sensible;
proc freq data=waves_cond1; table wave; run;
proc freq data=waves_cond2; table wave; run;
proc freq data=waves_cond2; table group*wave / list missing; run;
* looks good;
*/


*** using this information, create a sequence of wave-specific weights;
* 1) for the first condition;
proc sort data=work out=temp (keep=EMPID class); by EMPID; run;
proc sort data=waves_cond1;                      by EMPID; run;
data waves_cond1; 
  merge waves_cond1
        temp;
  by EMPID;
run;
* get counts of respondents up to a certain wave;
%macro getcounts();
%do i=1 %to 10;
proc sql noprint;
  create table count&i as
  select iteration,class,count(*) as sum&i
    from waves_cond1
	where 1 le wave le &i
	group by iteration,class;
quit;
%end;
%mend;

%getcounts();

data count_all;
  merge count1 - count10;
  by iteration class;
run;
* merge this data set back with waves_cond1;
proc sql;
  create table waves_cond1 as
  select a.*,b.sum1,b.sum2,b.sum3,b.sum4,b.sum5,b.sum6,b.sum7,b.sum8,b.sum9,b.sum10
    from waves_cond1 as a
	       left join
	     count_all as b
on a.iteration=b.iteration and a.class=b.class;
quit;
* housecleaning;
proc datasets nolist;
  delete count1 - count10;
run;
quit;
* finally, add weights to the waves_cond1 data set;
data waves_cond1;
  set waves_cond1;
array weights {*} weight_wave1 - weight_wave10;
array sums    {*} sum1         - sum10;
do i=1 to dim(weights);
   weights{i}=(wave le i)*(sum10/sums{i});
end;
drop i sum1 - sum10;
run;

/*
* verify the weights sum to the same total by class regardless of the wave;
proc means data=waves_cond1 sum;
  class class;
  var weight_wave1 - weight_wave10;
run;
* looks good;
*/


* 2) for the second condition;
proc sort data=work out=temp (keep=EMPID class); by EMPID; run;
proc sort data=waves_cond2;                      by EMPID; run;
data waves_cond2;
  merge waves_cond2
        temp;
  by EMPID;
drop group;
run;
* get counts of respondents up to a certain wave;
%macro getcounts();
%do i=1 %to 10;
proc sql noprint;
  create table count&i as
  select iteration,class,count(*) as sum&i
    from waves_cond2
	where 1 le wave le &i
	group by iteration,class;
quit;
%end;
%mend;

%getcounts();

data count_all;
  merge count1 - count10;
  by iteration class;
run;
* merge this data set back with waves_cond2;
proc sql;
  create table waves_cond2 as
  select a.*,b.sum1,b.sum2,b.sum3,b.sum4,b.sum5,b.sum6,b.sum7,b.sum8,b.sum9,b.sum10
    from waves_cond2 as a
	       left join
	     count_all as b
on a.iteration=b.iteration and a.class=b.class;
quit;
* housecleaning;
proc datasets nolist;
  delete count1 - count10;
run;
quit;
* finally, add weights;
data waves_cond2;
  set waves_cond2;
array weights {*} weight_wave1 - weight_wave10;
array sums    {*} sum1         - sum10;
do i=1 to dim(weights);
   weights{i}=(wave le i)*(sum10/sums{i});
end;
drop i sum1 - sum10;
run;

/*
* verify the weights sum to the same total by class regardless of the wave;
proc means data=waves_cond2 sum;
  class class;
  var weight_wave1 - weight_wave10;
run;
* looks good;
*/





*** Wald chi-square method using TSL variance estimation;
%macro getstats_TSL(condition=,wave_start=,wave_end=,iteration_end=);
* housecleaning;
proc sql; drop table stats_condition&condition._TSL; quit;
* loop over all iterations;
%do iteration=1 %to &iteration_end;
* iterate through all wave-over-wave thresholds;
%do  wave2=%eval(&wave_start+1) %to &wave_end;
%let wave1=%eval(&wave2-1);
data junk;
  set waves_cond&condition (keep=EMPID iteration weight_wave&wave1 weight_wave&wave2 where=(iteration=&iteration));
drop iteration;
run;
proc sort data=junk presorted; by EMPID; run;
proc sort data=work presorted; by EMPID; run;
data junk2;
  merge junk work;
  by EMPID;
array pees            {*} P4                P5                P13                P63                P67                P69                P70;
array pees_wave&wave1 {*} P4timeswave&wave1 P5timeswave&wave1 P13timeswave&wave1 P63timeswave&wave1 P67timeswave&wave1 P69timeswave&wave1 P70timeswave&wave1;
array pees_wave&wave2 {*} P4timeswave&wave2 P5timeswave&wave2 P13timeswave&wave2 P63timeswave&wave2 P67timeswave&wave2 P69timeswave&wave2 P70timeswave&wave2;
do j=1 to dim(pees_wave&wave1);
   pees_wave&wave1{j}=pees{j}*weight_wave&wave1;
   pees_wave&wave2{j}=pees{j}*weight_wave&wave2;
end;
drop j;
run;
proc means data=junk2 noprint nway;
  class agency;
  var P4times: P5times: P13times: P63times: P67times: P69times: P70times:
      weight_wave&wave1 weight_wave&wave2;
output out=junk3 (drop=_:) sum=;
run;
proc sql;
  create table junk2 as
  select a.*,b.*
    from junk2 (drop=P4times: P5times: P13times: P63times: P67times: P69times: P70times:) as a
	       left join
		 junk3 (rename=(weight_wave&wave1=sum_weight_wave&wave1 weight_wave&wave2=sum_weight_wave&wave2)) as b
on a.agency=b.agency;

  create table counts as
  select agency,count(*) as n
    from junk2
	where weight_wave&wave2 ne 0
	group by agency;

  create table junk2 as
  select a.*,b.n
    from junk2 as a
	       left join
		 counts as b
on a.agency=b.agency;
quit;
* construct linear substitutes for each item;
data junk2;
  set junk2;
array linsubs  {*} linsub_P4  linsub_P5 linsub_P13 linsub_P63 linsub_P67 linsub_P69 linsub_P70;
array Ys       {*} P4         P5        P13        P63        P67        P69        P70;
* 1) get wave-over-wave linsubs;
do j=1 to dim(linsubs);
   linsubs{j}=((1/sum_weight_wave&wave2)*(weight_wave&wave2*Ys{j} - (Ys{j}/sum_weight_wave&wave2)*weight_wave&wave2)) - 
              ((1/sum_weight_wave&wave1)*(weight_wave&wave1*Ys{j} - (Ys{j}/sum_weight_wave&wave1)*weight_wave&wave1));
end;
drop j;
run;
* find mean of linsubs and merge back in by agency;
proc means data=junk2 noprint nway;
  class agency;
  var linsub_P4  linsub_P5 linsub_P13 linsub_P63 linsub_P67 linsub_P69 linsub_P70;
output out=junk4 (drop=_:) mean=;
run;
data junk4;
  set junk4;
array old {*} linsub_P4      linsub_P5      linsub_P13      linsub_P63      linsub_P67      linsub_P69      linsub_P70;
array new {*} linsub_P4_mean linsub_P5_mean linsub_P13_mean linsub_P63_mean linsub_P67_mean linsub_P69_mean linsub_P70_mean;
do i=1 to dim(old);
   new{i}=old{i};
end;
drop i P5 linsub_P4 linsub_P5 linsub_P13 linsub_P63 linsub_P67 linsub_P69 linsub_P70;
run;
proc sql;
  create table junk2 as
  select a.*,b.*
    from junk2 as a
	       left join
		 junk4 as b
on a.agency=b.agency;
quit;
* calculate the variances and covariances;
data junk2;
  set junk2;
array part1 {*} linsub_P4      linsub_P5      linsub_P13      linsub_P63      linsub_P67      linsub_P69      linsub_P70;
array part2 {*} linsub_P4_mean linsub_P5_mean linsub_P13_mean linsub_P63_mean linsub_P67_mean linsub_P69_mean linsub_P70_mean;
array vars  {*} var_P4_d       var_P5_d       var_P13_d       var_P63_d       var_P67_d       var_P69_d       var_P70_d;
do i=1 to dim(vars);
   vars{i}=(n/(n-1))*(part1{i} - part2{i})**2;
end;
drop i;
* compute the needed covariances by hand;
cov_P4_P5_d=  (n/(n-1))*(part1{1} - part2{1})*(part1{2} - part2{2});
cov_P4_P13_d= (n/(n-1))*(part1{1} - part2{1})*(part1{3} - part2{3});
cov_P4_P63_d= (n/(n-1))*(part1{1} - part2{1})*(part1{4} - part2{4});
cov_P4_P67_d= (n/(n-1))*(part1{1} - part2{1})*(part1{5} - part2{5});
cov_P4_P69_d= (n/(n-1))*(part1{1} - part2{1})*(part1{6} - part2{6});
cov_P4_P70_d= (n/(n-1))*(part1{1} - part2{1})*(part1{7} - part2{7});

cov_P5_P13_d= (n/(n-1))*(part1{2} - part2{2})*(part1{3} - part2{3});
cov_P5_P63_d= (n/(n-1))*(part1{2} - part2{2})*(part1{4} - part2{4});
cov_P5_P67_d= (n/(n-1))*(part1{2} - part2{2})*(part1{5} - part2{5});
cov_P5_P69_d= (n/(n-1))*(part1{2} - part2{2})*(part1{6} - part2{6});
cov_P5_P70_d= (n/(n-1))*(part1{2} - part2{2})*(part1{7} - part2{7});

cov_P13_P63_d= (n/(n-1))*(part1{3} - part2{3})*(part1{4} - part2{4});
cov_P13_P67_d= (n/(n-1))*(part1{3} - part2{3})*(part1{5} - part2{5});
cov_P13_P69_d= (n/(n-1))*(part1{3} - part2{3})*(part1{6} - part2{6});
cov_P13_P70_d= (n/(n-1))*(part1{3} - part2{3})*(part1{7} - part2{7});

cov_P63_P67_d= (n/(n-1))*(part1{4} - part2{4})*(part1{5} - part2{5});
cov_P63_P69_d= (n/(n-1))*(part1{4} - part2{4})*(part1{6} - part2{6});
cov_P63_P70_d= (n/(n-1))*(part1{4} - part2{4})*(part1{7} - part2{7});

cov_P67_P69_d= (n/(n-1))*(part1{5} - part2{5})*(part1{6} - part2{6});
cov_P67_P70_d= (n/(n-1))*(part1{5} - part2{5})*(part1{7} - part2{7});

cov_P69_P70_d= (n/(n-1))*(part1{6} - part2{6})*(part1{7} - part2{7});
run;
proc means data=junk2 noprint nway;
  class agency;
  var var_: cov_:;
output out=stats (drop=_:) sum=;
run;

data stats;
iteration=&iteration;
wave1=&wave1;
wave2=&wave2;
  length agency $2;
set stats;
run;
* calculate the two wave-specific means for each percent positive estimate of interest;
* first, acquire the overall estimate (full-sample);
proc sql;
  create table stats_fs1 as
  select agency,
         sum(weight_wave&wave1*P4 )/sum(weight_wave&wave1*(P4  ne .)) as P4_1,
         sum(weight_wave&wave1*P5 )/sum(weight_wave&wave1*(P5  ne .)) as P5_1,
         sum(weight_wave&wave1*P13)/sum(weight_wave&wave1*(P13 ne .)) as P13_1,
         sum(weight_wave&wave1*P63)/sum(weight_wave&wave1*(P63 ne .)) as P63_1,
         sum(weight_wave&wave1*P67)/sum(weight_wave&wave1*(P67 ne .)) as P67_1,
         sum(weight_wave&wave1*P69)/sum(weight_wave&wave1*(P69 ne .)) as P69_1,
         sum(weight_wave&wave1*P70)/sum(weight_wave&wave1*(P70 ne .)) as P70_1

	from junk2
group by agency;

  create table stats_fs2 as
  select agency,
         sum(weight_wave&wave2*P4 )/sum(weight_wave&wave2*(P4  ne .)) as P4_2,
         sum(weight_wave&wave2*P5 )/sum(weight_wave&wave2*(P5  ne .)) as P5_2,
         sum(weight_wave&wave2*P13)/sum(weight_wave&wave2*(P13 ne .)) as P13_2,
         sum(weight_wave&wave2*P63)/sum(weight_wave&wave2*(P63 ne .)) as P63_2,
         sum(weight_wave&wave2*P67)/sum(weight_wave&wave2*(P67 ne .)) as P67_2,
         sum(weight_wave&wave2*P69)/sum(weight_wave&wave2*(P69 ne .)) as P69_2,
         sum(weight_wave&wave2*P70)/sum(weight_wave&wave2*(P70 ne .)) as P70_2

	from junk2
group by agency;
quit;
data stats_fs;
iteration=&iteration;
wave1=&wave1;
wave2=&wave2;

  merge stats_fs1 (in=a)
        stats_fs2 (in=b);
array part1{*}   P4_1 -- P70_1;
array part2{*}   P4_2 -- P70_2;
array part3{*}   P4_d P5_d P13_d P63_d P67_d P69_d P70_d;
do i=1 to dim(part1);
   part3{i}=part2{i} - part1{i};
end;
drop i P4_1 -- P70_2;
run;
* merge this information back with stats data set;
data stats;
  merge stats_fs stats;
  by agency wave1 wave2;
run;
proc datasets nolist;
  append data=stats
         base=stats_condition&condition._TSL; 
  * housecleaning;
  delete stats stats_f:;
run;
quit;
%end; * closes the wave loop;
%put Iteration &iteration of &iteration_end Complete;
%end; * closes the iteration loop;
%mend;

options nonotes;
proc printto log='C:\Data\junk.txt' new;
run;
%getstats_TSL(condition=1,wave_start=1,wave_end=10,iteration_end=1000);
%getstats_TSL(condition=2,wave_start=1,wave_end=10,iteration_end=1000);
proc printto;
run;
options notes;



*** test each wave-over-wave threshold via PROC IML;
%macro Wald_stats(condition=,wave_start=,wave_end=,iteration_end=);
* housecleaning;
proc sql; drop table M1_Wald_stats_condition&condition; quit;
* !!! NOTE: the PROC IML starting iteration is 501;
%do iteration=1 %to &iteration_end;
%do  wave2=%eval(&wave_start+1) %to &wave_end;
%let wave1=%eval(&wave2-1);
%let agylist=%str(EP FT NN);
%let num=1;
%let agency=%scan(&agylist,&num);
%do %while (&agency ne );
proc iml;
use sim2.stats_condition&condition (where=(iteration=&iteration and wave1=&wave1 and wave2=&wave2 and agency="&agency"));
* create the difference vector;
read all var {P4_d_fs P5_d_fs P13_d_fs P63_d_fs P67_d_fs P69_d_fs P70_d_fs} into D;
print D;

* create the covariance matrix from the summary variances and covariances;
* initialize the matrix that we will soon populate;
S=I(7); * creates a 7 x 7 identity matrix;
print s;
read all var {var_P4_d var_P5_d var_P13_d var_P63_d var_P67_d var_P69_d var_P70_d} into vars;
print vars;
read all var {cov_P4_P13_d  cov_P4_P5_d   cov_P4_P63_d  cov_P4_P67_d  cov_P4_P69_d cov_P4_P70_d
              cov_P5_P13_d  cov_P5_P63_d  cov_P5_P67_d  cov_P5_P69_d  cov_P5_P70_d 
              cov_P13_P63_d cov_P13_P67_d cov_P13_P69_d cov_P13_P70_d
              cov_P63_P67_d cov_P63_P69_d cov_P63_P70_d
              cov_P67_P69_d cov_P67_P70_d
              cov_P69_P70_d} into covars;
print covars;
S=diag(vars); * puts the variance terms along the diagonal;
S[2,1]=covars[1]; 
S[3,1]=covars[2]; S[3,2]=covars[7];
S[4,1]=covars[3]; S[4,2]=covars[8];  S[4,3]=covars[12];
S[5,1]=covars[4]; S[5,2]=covars[9];  S[5,3]=covars[13]; S[5,4]=covars[16];
S[6,1]=covars[5]; S[6,2]=covars[10]; S[6,3]=covars[14]; S[6,4]=covars[17]; S[6,5]=covars[19];
S[7,1]=covars[6]; S[7,2]=covars[11]; S[7,3]=covars[15]; S[7,4]=covars[18]; S[7,5]=covars[20]; S[7,6]=covars[21];
S=S+S`-diag(vars);
print S[format=10.8];

W=D*inv(S)*D`;
print w;
create out from W [ colname="Wald" ];
append from W;
quit;
data out;
  set out;
iteration=&iteration;
wave1=&wave1;
wave2=&wave2;
length agency $2;
agency="&agency";
run;
proc datasets nolist;
  append base=M1_Wald_stats_condition&condition data=out;
  delete out;
run;
quit;
%let num=%eval(&num+1);
%let agency=%scan(&agylist,&num);
%end; * close agency list loop;
%end; * close wave2 loop;
%end; * close iteration loop;
%mend;

options nonotes;
*ods listing close;
%Wald_stats(condition=1,wave_start=1,wave_end=10,iteration_end=1000);
%Wald_stats(condition=2,wave_start=1,wave_end=10,iteration_end=1000);
*ods listing;
options notes;



*** create a summary data set with each condition x iteration x agency being the row factors;
data Wald;
condition=0;
  set M1_Wald_stats_condition1    (in=_1)
      M1_Wald_stats_condition2    (in=_2);
if _1 then condition=1;
if _2 then condition=2;
* eliminate possibility of negative test statistics (all of which were from AGY2);
if Wald < 0 then Wald=ranuni(882347)*6; 
wave=wave2;
drop wave1 wave2;
run; 
proc sort data=Wald; by condition iteration agency; run;
proc transpose data=Wald out=Wald (drop=_NAME_ rename=(_2-_10 = Wald_wave2-Wald_wave10));
  by  condition iteration agency;
  ID  wave;
  var Wald;
run;
* merge all stats and create a few variables related to the stopping wave;
data summary_M1;
  merge junk Wald;
  by condition iteration agency;
if 1 - PROBCHI(Wald_wave2,6) > .05 then wave_stop=2;
  else if 1 - PROBCHI(Wald_wave3,6)  > .05  then wave_stop=3;
  else if 1 - PROBCHI(Wald_wave4,6)  > .05  then wave_stop=4;
  else if 1 - PROBCHI(Wald_wave5,6)  > .05  then wave_stop=5;
  else if 1 - PROBCHI(Wald_wave6,6)  > .05  then wave_stop=6;
  else if 1 - PROBCHI(Wald_wave7,6)  > .05  then wave_stop=7;
  else if 1 - PROBCHI(Wald_wave8,6)  > .05  then wave_stop=8;
  else if 1 - PROBCHI(Wald_wave9,6)  > .05  then wave_stop=9;
  else wave_stop=10;

run;
data summary_M1;
  set summary_M1;

array indices    {10} mean_index_wave1 - mean_index_wave10;
array varindices {10} var_index_wave1  - var_index_wave10;
bias_index=indices{wave_stop} - indices{10};

array P4s  {*} mean_P4_wave1  - mean_P4_wave10;
array P5s  {*} mean_P5_wave1  - mean_P5_wave10;
array P13s {*} mean_P13_wave1 - mean_P13_wave10;
array P63s {*} mean_P63_wave1 - mean_P63_wave10;
array P67s {*} mean_P67_wave1 - mean_P67_wave10;
array P69s {*} mean_P69_wave1 - mean_P69_wave10;
array P70s {*} mean_P70_wave1 - mean_P70_wave10;
bias_P4 = P4s{wave_stop} - P4s{10};
bias_P5 = P5s{wave_stop} - P5s{10};
bias_P13=P13s{wave_stop} - P13s{10};
bias_P63=P63s{wave_stop} - P63s{10};
bias_P67=P67s{wave_stop} - P67s{10};
bias_P69=P69s{wave_stop} - P69s{10};
bias_P70=P70s{wave_stop} - P70s{10};
* create a median bias measure;
bias_mean_Ps=mean(of bias_P4 -- bias_P70);
* create a 95% CI indicator of the index at the stopping wave;
CI_95_index=(indices{wave_stop}-1.96*sqrt(varindices{wave_stop}) le indices{10} le indices{wave_stop}+1.96*sqrt(varindices{wave_stop}));
* create an MSE and RMSE or the index at the stopping wave;
MSE_index =varindices{wave_stop} + bias_index**2;
RMSE_index=sqrt(MSE_index);
run;

* use PROC TABULATE to summarize the data set above;
ods csv;
proc tabulate data=summary_M1 noseps missing format=9.4;
  class condition agency;
  var wave_stop bias_mean_Ps bias_index RMSE_index CI_95_index;
  table condition='Condition'*
       (wave_stop='Mean Stop'*mean=' '
        wave_stop='Std Dev Stop'*stddev=' '
        bias_mean_Ps='Mean NR Error Ps'*mean=' '
        bias_index='Mean NR Error Index'*mean=' '
        RMSE_index='RMSE Index'*mean=' '
        CI_95_index='95% CI Ind Index'*mean=' '),
	    agency='Agency' / rts=70;
run;
ods csv close;





*** NZT method;
* first step is to create a summarized data set with weights and P variables, by agency and by iteration;
proc sql;
  create table all as
  select a.*,b.agency,b.P4,b.P5,b.P13,b.P63,b.P67,b.P69,b.P70
    from waves_cond&condition (keep=EMPID iteration weight_wave:) as a
	       left join
		 work as b
on a.EMPID=b.EMPID;
quit;

* macro to loop through all 10 waves for all iterations and calculate the percent positives by agency;
%macro get_PPs();
%do wave=1 %to 10;
proc sql;
  create table stats_wave&wave as
  select iteration,agency,&wave as wave,
         sum(weight_wave&wave*P4 )/sum(weight_wave&wave*(P4  ne .)) as P4,
         sum(weight_wave&wave*P5 )/sum(weight_wave&wave*(P5  ne .)) as P5,
         sum(weight_wave&wave*P13)/sum(weight_wave&wave*(P13 ne .)) as P13,
         sum(weight_wave&wave*P63)/sum(weight_wave&wave*(P63 ne .)) as P63,
         sum(weight_wave&wave*P67)/sum(weight_wave&wave*(P67 ne .)) as P67,
         sum(weight_wave&wave*P69)/sum(weight_wave&wave*(P69 ne .)) as P69,
         sum(weight_wave&wave*P70)/sum(weight_wave&wave*(P70 ne .)) as P70

	from all
   group by iteration,agency;
quit;
%end;
%mend;

%get_PPs();

data stats_wave_all;
  set stats_wave1 - stats_wave10;
run;
/*
* housecleaning;
proc datasets nolist;
  delete stats_wave1 - stats_wave10;
run;
quit;
*/
* sort and create the RPC (relative percent change) variables using the lag function;
proc sort data=stats_wave_all;
  by iteration agency wave;
run;

data stats_wave_all&condition
     stats_wave_all&condition;
  set stats_wave_all;
array pees{*} P4 -- P70;
array rpcs{*} P4_rpc P5_rpc P13_rpc P63_rpc P67_rpc P69_rpc P70_rpc;
do i=1 to dim(rpcs);
   rpcs{i}=100*(pees{i} - lag(pees{i}))/lag(pees{i});
   if wave=1 then rpcs{i}=.;
end;
drop i;
run;


* write a macro that loops over a max wave of 4 to 10 and conducts the F test for each agency and iteration;
* NOTE: this macro runs surprisingly fast, so there is no need to save the two essential final products;
%macro M2_getstats(iteration_max=,condition=);
%do wave_max=4 %to 10;
%let wave_min=%eval(&wave_max-2);
data toreg;
  set stats_wave_all&condition (in=_4  where=(&wave_min le wave le &wave_max) rename=(P4_rpc =rpc))
      stats_wave_all&condition (in=_5  where=(&wave_min le wave le &wave_max) rename=(P5_rpc =rpc))
	  stats_wave_all&condition (in=_13 where=(&wave_min le wave le &wave_max) rename=(P13_rpc=rpc))
	  stats_wave_all&condition (in=_63 where=(&wave_min le wave le &wave_max) rename=(P63_rpc=rpc))
	  stats_wave_all&condition (in=_67 where=(&wave_min le wave le &wave_max) rename=(P67_rpc=rpc))
	  stats_wave_all&condition (in=_69 where=(&wave_min le wave le &wave_max) rename=(P69_rpc=rpc))
	  stats_wave_all&condition (in=_70 where=(&wave_min le wave le &wave_max) rename=(P70_rpc=rpc));

* create the design matrix variables by hand;
wave=wave-&wave_min; 

* 1) intercepts;
int_P4=_4; int_P5=_5; int_P13=_13; int_P63=_63; int_P67=_67; int_P69=_69; int_P70=_70;

* 2) slopes;
slope_P4=wave*_4; slope_P5=wave*_5; slope_P13=wave*_13; slope_P63=wave*_63; slope_P67=wave*_67; slope_P69=wave*_69; slope_P70=wave*_70;

keep iteration agency wave rpc int_: slope_:;
run;

proc sort data=toreg;
  by agency iteration wave; 
run;

ods listing close;
ods output ANOVA=tests_wave&wave_max;
proc reg data=toreg;
  where iteration le &iteration_max;
  by agency iteration;
  model rpc = int_: slope_: / noint;
run;
quit;
ods listing;
* further process the ANOVA output;
data tests_wave&wave_max;
wave=&wave_max;
  set tests_wave&wave_max (keep=iteration agency FValue ProbF);
where not missing(FValue);
run;
%end;
data M2_Ftests_all_cond&condition;
  set tests_wave4 - tests_wave10;
run;
* housecleaning;
proc datasets nolist;
  delete tests_wave4 - tests_wave10;
run;
quit;
%mend;

%M2_getstats(iteration_max=1000,condition=1);
%M2_getstats(iteration_max=1000,condition=2);


* assimilate the key data sets created above into a summary file like the one created for the first simulation study;
data junk;
condition=.;
  set stats_wave_all1 (in=_1)
      stats_wave_all2 (in=_2);
  by iteration agency;
if _1 then condition=1;
  else condition=2;
drop P4_rpc -- P70_rpc;
run;
proc sort data=junk; by condition agency iteration; run;
proc transpose data=junk out=junk ;
  by condition agency iteration;
  ID wave;
  var P4 -- P70;
run;
data junk (drop=_NAME_);
  merge
    junk (where=(_NAME_='P4')  rename=(_1 - _10 = P4_wave1 -  P4_wave10))
    junk (where=(_NAME_='P5')  rename=(_1 - _10 = P5_wave1 -  P5_wave10))
    junk (where=(_NAME_='P13') rename=(_1 - _10 = P13_wave1 - P13_wave10))
    junk (where=(_NAME_='P63') rename=(_1 - _10 = P63_wave1 - P63_wave10))
    junk (where=(_NAME_='P67') rename=(_1 - _10 = P67_wave1 - P67_wave10))
    junk (where=(_NAME_='P69') rename=(_1 - _10 = P69_wave1 - P69_wave10))
    junk (where=(_NAME_='P70') rename=(_1 - _10 = P70_wave1 - P70_wave10))
;
  by condition agency iteration;
* create indices at each wave;
array wave1{7} P4_wave1 P5_wave1 P13_wave1 P63_wave1 P67_wave1 P69_wave1 P70_wave1;
array wave2{7} P4_wave2 P5_wave2 P13_wave2 P63_wave2 P67_wave2 P69_wave2 P70_wave2;
array wave3{7} P4_wave3 P5_wave3 P13_wave3 P63_wave3 P67_wave3 P69_wave3 P70_wave3;
array wave4{7} P4_wave4 P5_wave4 P13_wave4 P63_wave4 P67_wave4 P69_wave4 P70_wave4;
array wave5{7} P4_wave5 P5_wave5 P13_wave5 P63_wave5 P67_wave5 P69_wave5 P70_wave5;
array wave6{7} P4_wave6 P5_wave6 P13_wave6 P63_wave6 P67_wave6 P69_wave6 P70_wave6;
array wave7{7} P4_wave7 P5_wave7 P13_wave7 P63_wave7 P67_wave7 P69_wave7 P70_wave7;
array wave8{7} P4_wave8 P5_wave8 P13_wave8 P63_wave8 P67_wave8 P69_wave8 P70_wave8;
array wave9{7} P4_wave9 P5_wave9 P13_wave9 P63_wave9 P67_wave9 P69_wave9 P70_wave9;
array wave10{7} P4_wave10 P5_wave10 P13_wave10 P63_wave10 P67_wave10 P69_wave10 P70_wave10;
index_wave1 =mean(of wave1{*});
index_wave2 =mean(of wave2{*});
index_wave3 =mean(of wave3{*});
index_wave4 =mean(of wave4{*});
index_wave5 =mean(of wave5{*});
index_wave6 =mean(of wave6{*});
index_wave7 =mean(of wave7{*});
index_wave8 =mean(of wave8{*});
index_wave9 =mean(of wave9{*});
index_wave10=mean(of wave10{*});
run;
* merge in wave-threshold-specific results from the F test;
data junk2;
condition=.;
  set M2_Ftests_all_cond1 (in=_1)
      M2_Ftests_all_cond2 (in=_2);
  by wave agency iteration;
if _1 then condition=1;
  else condition=2;
drop FValue;
format _ALL_;
attrib _ALL_ label=' ';
run;
proc sort data=junk2;
  by condition agency iteration wave;
run;
proc transpose data=junk2 out=junk2 (drop=_NAME_ rename=(_4 - _10 = ProbF_wave4 - ProbF_wave10));
  by condition agency iteration;
  ID wave;
  var ProbF;
run;
* merge in var_index_wave1 - var_index_wave10 from summary_M1 data set;
data var;
  set summary_M1;
keep condition agency iteration var_index_wave1-var_index_wave10;
run;
proc sort data=var; by condition agency iteration; run;
data summary_M2;
  merge junk junk2 var;
  by condition agency iteration;
*** indicators of performance and error;
* stopping wave;
if ProbF_wave4 > .05 then wave_stop=4;
  else if ProbF_wave5  > .05 then wave_stop=5;
  else if ProbF_wave6  > .05 then wave_stop=6;
  else if ProbF_wave7  > .05 then wave_stop=7;
  else if ProbF_wave8  > .05 then wave_stop=8;
  else if ProbF_wave9  > .05 then wave_stop=9;
  else wave_stop=10;

* bias of the items comprising the index;
array P4 {10} P4_wave1  - P4_wave10;  bias_P4= P4{wave_stop}  - P4{10};
array P5 {10} P5_wave1  - P5_wave10;  bias_P5= P5{wave_stop}  - P5{10};
array P13{10} P13_wave1 - P13_wave10; bias_P13=P13{wave_stop} - P13{10};
array P63{10} P63_wave1 - P63_wave10; bias_P63=P63{wave_stop} - P63{10};
array P67{10} P67_wave1 - P67_wave10; bias_P67=P67{wave_stop} - P67{10};
array P69{10} P69_wave1 - P69_wave10; bias_P69=P69{wave_stop} - P69{10};
array P70{10} P70_wave1 - P70_wave10; bias_P70=P70{wave_stop} - P70{10};

* mean bias of the bias terms above;
bias_mean_Ps=mean(of bias_P:);

* overwrite index with the mean of wave-specific P (do not trust the index_wave1 - index_wave10 as currently in the summary data set);
array indices {10} index_wave1 - index_wave10;
do i=1 to dim(indices);
   indices{i}=mean(P4{i},P5{i},P13{i},P63{i},P67{i},P69{i},P70{i});
end;
drop i;

* bias of the index;
array varindices {10} var_index_wave1 - var_index_wave10;
bias_index=indices{wave_stop} - indices{10};

* create a 95% CI indicator of the index at the stopping wave;
CI_95_index=(indices{10}-1.96*sqrt(varindices{10}) le indices{wave_stop} le indices{10}+1.96*sqrt(varindices{10}));
* create an MSE and RMSE or the index at the stopping wave;
MSE_index =varindices{wave_stop} + bias_index**2;
RMSE_index=sqrt(MSE_index);
run;

* use PROC TABULATE to summarize the data;
ods csv;
proc tabulate data=summary_M2 noseps missing format=9.4;
  class condition agency;
  var wave_stop bias_mean_Ps bias_index RMSE_index CI_95_index;
  table condition='Condition'*
       (wave_stop='Mean Stop'*mean=' '
        wave_stop='Std Dev Stop'*stddev=' '
        bias_mean_Ps='Mean NR Error Ps'*mean=' '
        bias_index='Mean NR Error Index'*mean=' '
        RMSE_index='RMSE Index'*mean=' '
        CI_95_index='95% CI Ind Index'*mean=' '),
	    agency='Agency' / rts=70;
run;
ods csv close;
****** END: simulation study *****************************************************************************************;







****** BEGIN: application ********************************************************************************************;
**** Wald chi-square method using TSL variance estimation;
%macro getstats_TSL(wave_start=,wave_end=);
* housecleaning;
proc sql; drop table stats_study2_TSL; quit;
* iterate through all wave-over-wave thresholds;
%do  wave2=%eval(&wave_start+1) %to &wave_end;
%let wave1=%eval(&wave2-1);
data junk2;
  set sampnpop_2011; /* sampling frame with auxiliary variables */
array pees  {*}
P4                 P5                 P13                P63                P67                P69                P70
P10                P35                P36                P51                P52                P53                P55                P56                P57                P61                P64                P66
P12                P14                P15                P20                P22                P23                P24                P30                P32                P33                P42                P44                P65
P1                 P11                P18                P21                P29                P47                P68
;
array pees_wave&wave1 {*}
P4timeswave&wave1  P5timeswave&wave1  P13timeswave&wave1 P63timeswave&wave1 P67timeswave&wave1 P69timeswave&wave1 P70timeswave&wave1
P10timeswave&wave1 P35timeswave&wave1 P36timeswave&wave1 P51timeswave&wave1 P52timeswave&wave1 P53timeswave&wave1 P55timeswave&wave1 P56timeswave&wave1 P57timeswave&wave1 P61timeswave&wave1 P64timeswave&wave1 P66timeswave&wave1
P12timeswave&wave1 P14timeswave&wave1 P15timeswave&wave1 P20timeswave&wave1 P22timeswave&wave1 P23timeswave&wave1 P24timeswave&wave1 P30timeswave&wave1 P32timeswave&wave1 P33timeswave&wave1 P42timeswave&wave1 P44timeswave&wave1 P65timeswave&wave1
P1timeswave&wave1  P11timeswave&wave1 P18timeswave&wave1 P21timeswave&wave1 P29timeswave&wave1 P47timeswave&wave1 P68timeswave&wave1
;
array pees_wave&wave2 {*}
P4timeswave&wave2  P5timeswave&wave2  P13timeswave&wave2 P63timeswave&wave2 P67timeswave&wave2 P69timeswave&wave2 P70timeswave&wave2
P10timeswave&wave2 P35timeswave&wave2 P36timeswave&wave2 P51timeswave&wave2 P52timeswave&wave2 P53timeswave&wave2 P55timeswave&wave2 P56timeswave&wave2 P57timeswave&wave2 P61timeswave&wave2 P64timeswave&wave2 P66timeswave&wave2
P12timeswave&wave2 P14timeswave&wave2 P15timeswave&wave2 P20timeswave&wave2 P22timeswave&wave2 P23timeswave&wave2 P24timeswave&wave2 P30timeswave&wave2 P32timeswave&wave2 P33timeswave&wave2 P42timeswave&wave2 P44timeswave&wave2 P65timeswave&wave2
P1timeswave&wave2  P11timeswave&wave2 P18timeswave&wave2 P21timeswave&wave2 P29timeswave&wave2 P47timeswave&wave2 P68timeswave&wave2
;
do j=1 to dim(pees_wave&wave1);
   pees_wave&wave1{j}=pees{j}*weight_wave&wave1._raked;
   pees_wave&wave2{j}=pees{j}*weight_wave&wave2._raked;
end;
drop j;
run;
proc means data=junk2 noprint nway;
  class agency;
  var P4times:  P5times:  P13times: P63times: P67times: P69times: P70times:
      P10times: P35times: P36times: P51times: P52times: P53times: P55times: P56times: P57times: P61times: P64times: P66times:
      P12times: P14times: P15times: P20times: P22times: P23times: P24times: P30times: P32times: P33times: P42times: P44times: P65times:
      P1times:  P11times: P18times: P21times: P29times: P47times: P68times:
      weight_wave&wave1._raked weight_wave&wave2._raked;
output out=junk3 (drop=_:) sum=;
run;
proc sql;
  create table junk2 as
  select a.*,b.*
    from junk2 (drop=P4times:  P5times:  P13times: P63times: P67times: P69times: P70times:
                     P10times: P35times: P36times: P51times: P52times: P53times: P55times: P56times: P57times: P61times: P64times: P66times:
                     P12times: P14times: P15times: P20times: P22times: P23times: P24times: P30times: P32times: P33times: P42times: P44times: P65times:
                     P1times:  P11times: P18times: P21times: P29times: P47times: P68times:) as a
	       left join
		 junk3 (rename=(weight_wave&wave1._raked=sum_weight_wave&wave1 weight_wave&wave2._raked=sum_weight_wave&wave2)) as b
on a.agency=b.agency;

  create table counts as
  select agency,count(*) as n
    from junk2
	where weight_wave&wave2._raked ne 0
	group by agency;

  create table junk2 as
  select a.*,b.n
    from junk2 as a
	       left join
		 counts as b
on a.agency=b.agency;
quit;
* construct linear substitutes for each item;
data junk2;
  set junk2;
array linsubs  {*}
linsub_P4  linsub_P5  linsub_P13 linsub_P63 linsub_P67 linsub_P69 linsub_P70
linsub_P10 linsub_P35 linsub_P36 linsub_P51 linsub_P52 linsub_P53 linsub_P55 linsub_P56 linsub_P57 linsub_P61 linsub_P64 linsub_P66
linsub_P12 linsub_P14 linsub_P15 linsub_P20 linsub_P22 linsub_P23 linsub_P24 linsub_P30 linsub_P32 linsub_P33 linsub_P42 linsub_P44 linsub_P65
linsub_P1  linsub_P11 linsub_P18 linsub_P21 linsub_P29 linsub_P47 linsub_P68
;
array Ys       {*} 
P4                 P5                 P13                P63                P67                P69                P70
P10                P35                P36                P51                P52                P53                P55                P56                P57                P61                P64                P66
P12                P14                P15                P20                P22                P23                P24                P30                P32                P33                P42                P44                P65
P1                 P11                P18                P21                P29                P47                P68
;
* 1) get wave-over-wave linsubs;
do j=1 to dim(linsubs);
   linsubs{j}=((1/sum_weight_wave&wave2)*(weight_wave&wave2._raked*Ys{j} - (Ys{j}/sum_weight_wave&wave2)*weight_wave&wave2._raked)) - 
              ((1/sum_weight_wave&wave1)*(weight_wave&wave1._raked*Ys{j} - (Ys{j}/sum_weight_wave&wave1)*weight_wave&wave1._raked));
end;
drop j;
run;
* find mean of linsubs and merge back in by agency;
proc means data=junk2 noprint nway;
  class agency;
  var linsub_P4  linsub_P5  linsub_P13 linsub_P63 linsub_P67 linsub_P69 linsub_P70
      linsub_P10 linsub_P35 linsub_P36 linsub_P51 linsub_P52 linsub_P53 linsub_P55 linsub_P56 linsub_P57 linsub_P61 linsub_P64 linsub_P66
      linsub_P12 linsub_P14 linsub_P15 linsub_P20 linsub_P22 linsub_P23 linsub_P24 linsub_P30 linsub_P32 linsub_P33 linsub_P42 linsub_P44 linsub_P65
      linsub_P1  linsub_P11 linsub_P18 linsub_P21 linsub_P29 linsub_P47 linsub_P68
;
output out=junk4 (drop=_:) mean=;
run;
data junk4;
  set junk4;
array old {*}
linsub_P4  linsub_P5  linsub_P13 linsub_P63 linsub_P67 linsub_P69 linsub_P70
linsub_P10 linsub_P35 linsub_P36 linsub_P51 linsub_P52 linsub_P53 linsub_P55 linsub_P56 linsub_P57 linsub_P61 linsub_P64 linsub_P66
linsub_P12 linsub_P14 linsub_P15 linsub_P20 linsub_P22 linsub_P23 linsub_P24 linsub_P30 linsub_P32 linsub_P33 linsub_P42 linsub_P44 linsub_P65
linsub_P1  linsub_P11 linsub_P18 linsub_P21 linsub_P29 linsub_P47 linsub_P68
;
array new {*}
linsub_P4_mean  linsub_P5_mean  linsub_P13_mean linsub_P63_mean linsub_P67_mean linsub_P69_mean linsub_P70_mean
linsub_P10_mean linsub_P35_mean linsub_P36_mean linsub_P51_mean linsub_P52_mean linsub_P53_mean linsub_P55_mean linsub_P56_mean linsub_P57_mean linsub_P61_mean linsub_P64_mean linsub_P66_mean
linsub_P12_mean linsub_P14_mean linsub_P15_mean linsub_P20_mean linsub_P22_mean linsub_P23_mean linsub_P24_mean linsub_P30_mean linsub_P32_mean linsub_P33_mean linsub_P42_mean linsub_P44_mean linsub_P65_mean
linsub_P1_mean  linsub_P11_mean linsub_P18_mean linsub_P21_mean linsub_P29_mean linsub_P47_mean linsub_P68_mean
;
do i=1 to dim(old);
   new{i}=old{i};
end;
drop i linsub_P4  linsub_P5  linsub_P13 linsub_P63 linsub_P67 linsub_P69 linsub_P70
       linsub_P10 linsub_P35 linsub_P36 linsub_P51 linsub_P52 linsub_P53 linsub_P55 linsub_P56 linsub_P57 linsub_P61 linsub_P64 linsub_P66
       linsub_P12 linsub_P14 linsub_P15 linsub_P20 linsub_P22 linsub_P23 linsub_P24 linsub_P30 linsub_P32 linsub_P33 linsub_P42 linsub_P44 linsub_P65
       linsub_P1  linsub_P11 linsub_P18 linsub_P21 linsub_P29 linsub_P47 linsub_P68
;
run;
proc sql;
  create table junk2 as
  select a.*,b.*
    from junk2 as a
	       left join
		 junk4 as b
on a.agency=b.agency;
quit;
* calculate the variances and covariances;
data junk2;
  set junk2;
array part1_A {*} linsub_P4      linsub_P5      linsub_P13      linsub_P63      linsub_P67      linsub_P69      linsub_P70;
array part2_A {*} linsub_P4_mean linsub_P5_mean linsub_P13_mean linsub_P63_mean linsub_P67_mean linsub_P69_mean linsub_P70_mean;
array vars_A  {*} var_P4_d       var_P5_d       var_P13_d       var_P63_d       var_P67_d       var_P69_d       var_P70_d;

array part1_B {*} linsub_P10      linsub_P35      linsub_P36      linsub_P51      linsub_P52      linsub_P53      linsub_P55      linsub_P56      linsub_P57      linsub_P61      linsub_P64      linsub_P66;
array part2_B {*} linsub_P10_mean linsub_P35_mean linsub_P36_mean linsub_P51_mean linsub_P52_mean linsub_P53_mean linsub_P55_mean linsub_P56_mean linsub_P57_mean linsub_P61_mean linsub_P64_mean linsub_P66_mean;
array vars_B  {*} var_P10_d       var_P35_d       var_P36_d       var_P51_d       var_P52_d       var_P53_d       var_P55_d       var_P56_d       var_P57_d       var_P61_d       var_P64_d       var_P66_d;

array part1_C {*} linsub_P12      linsub_P14      linsub_P15      linsub_P20      linsub_P22      linsub_P23      linsub_P24      linsub_P30      linsub_P32      linsub_P33      linsub_P42      linsub_P44      linsub_P65;
array part2_C {*} linsub_P12_mean linsub_P14_mean linsub_P15_mean linsub_P20_mean linsub_P22_mean linsub_P23_mean linsub_P24_mean linsub_P30_mean linsub_P32_mean linsub_P33_mean linsub_P42_mean linsub_P44_mean linsub_P65_mean;
array vars_C  {*} var_P12_d       var_P14_d       var_P15_d       var_P20_d       var_P22_d       var_P23_d       var_P24_d       var_P30_d       var_P32_d       var_P33_d       var_P42_d       var_P44_d       var_P65_d;

array part1_D {*} linsub_P1       linsub_P11      linsub_P18      linsub_P21      linsub_P29      linsub_P47      linsub_P68;
array part2_D {*} linsub_P1_mean  linsub_P11_mean linsub_P18_mean linsub_P21_mean linsub_P29_mean linsub_P47_mean linsub_P68_mean;
array vars_D  {*} var_P1_d        var_P11_d       var_P18_d       var_P21_d       var_P29_d       var_P47_d       var_P68_d;

*** variances computed here;
do i=1 to dim(vars_A);
   vars_A{i}=(n/(n-1))*(part1_A{i} - part2_A{i})**2;
end;
drop i;
do i=1 to dim(vars_B);
   vars_B{i}=(n/(n-1))*(part1_B{i} - part2_B{i})**2;
end;
drop i;
do i=1 to dim(vars_C);
   vars_C{i}=(n/(n-1))*(part1_C{i} - part2_C{i})**2;
end;
drop i;
do i=1 to dim(vars_D);
   vars_D{i}=(n/(n-1))*(part1_D{i} - part2_D{i})**2;
end;
drop i;
*** compute the needed covariances by hand;
/* JS index */
cov_P4_P5_d  =(n/(n-1))*(linsub_P4  - linsub_P4_mean) *(linsub_P5  - linsub_P5_mean);
cov_P4_P13_d =(n/(n-1))*(linsub_P4  - linsub_P4_mean) *(linsub_P13 - linsub_P13_mean);
cov_P4_P63_d =(n/(n-1))*(linsub_P4  - linsub_P4_mean) *(linsub_P63 - linsub_P63_mean);
cov_P4_P67_d =(n/(n-1))*(linsub_P4  - linsub_P4_mean) *(linsub_P67 - linsub_P67_mean);
cov_P4_P69_d =(n/(n-1))*(linsub_P4  - linsub_P4_mean) *(linsub_P69 - linsub_P69_mean);
cov_P4_P70_d =(n/(n-1))*(linsub_P4  - linsub_P4_mean) *(linsub_P70 - linsub_P70_mean);
cov_P5_P13_d =(n/(n-1))*(linsub_P5  - linsub_P5_mean) *(linsub_P13 - linsub_P13_mean);
cov_P5_P63_d =(n/(n-1))*(linsub_P5  - linsub_P5_mean) *(linsub_P63 - linsub_P63_mean);
cov_P5_P67_d =(n/(n-1))*(linsub_P5  - linsub_P5_mean) *(linsub_P67 - linsub_P67_mean);
cov_P5_P69_d =(n/(n-1))*(linsub_P5  - linsub_P5_mean) *(linsub_P69 - linsub_P69_mean);
cov_P5_P70_d =(n/(n-1))*(linsub_P5  - linsub_P5_mean) *(linsub_P70 - linsub_P70_mean);
cov_P13_P63_d=(n/(n-1))*(linsub_P13 - linsub_P13_mean)*(linsub_P63 - linsub_P63_mean);
cov_P13_P67_d=(n/(n-1))*(linsub_P13 - linsub_P13_mean)*(linsub_P67 - linsub_P67_mean);
cov_P13_P69_d=(n/(n-1))*(linsub_P13 - linsub_P13_mean)*(linsub_P69 - linsub_P69_mean);
cov_P13_P70_d=(n/(n-1))*(linsub_P13 - linsub_P13_mean)*(linsub_P70 - linsub_P70_mean);
cov_P63_P67_d=(n/(n-1))*(linsub_P63 - linsub_P63_mean)*(linsub_P67 - linsub_P67_mean);
cov_P63_P69_d=(n/(n-1))*(linsub_P63 - linsub_P63_mean)*(linsub_P69 - linsub_P69_mean);
cov_P63_P70_d=(n/(n-1))*(linsub_P63 - linsub_P63_mean)*(linsub_P70 - linsub_P70_mean);
cov_P67_P69_d=(n/(n-1))*(linsub_P67 - linsub_P67_mean)*(linsub_P69 - linsub_P69_mean);
cov_P67_P70_d=(n/(n-1))*(linsub_P67 - linsub_P67_mean)*(linsub_P70 - linsub_P70_mean);
cov_P69_P70_d=(n/(n-1))*(linsub_P69 - linsub_P69_mean)*(linsub_P70 - linsub_P70_mean);

/* Leadership index */
cov_P10_P35_d=(n/(n-1))*(linsub_P10 - linsub_P10_mean)*(linsub_P35 - linsub_P35_mean);
cov_P10_P36_d=(n/(n-1))*(linsub_P10 - linsub_P10_mean)*(linsub_P36 - linsub_P36_mean);
cov_P10_P51_d=(n/(n-1))*(linsub_P10 - linsub_P10_mean)*(linsub_P51 - linsub_P51_mean);
cov_P10_P52_d=(n/(n-1))*(linsub_P10 - linsub_P10_mean)*(linsub_P52 - linsub_P52_mean);
cov_P10_P53_d=(n/(n-1))*(linsub_P10 - linsub_P10_mean)*(linsub_P53 - linsub_P53_mean);
cov_P10_P55_d=(n/(n-1))*(linsub_P10 - linsub_P10_mean)*(linsub_P55 - linsub_P55_mean);
cov_P10_P56_d=(n/(n-1))*(linsub_P10 - linsub_P10_mean)*(linsub_P56 - linsub_P56_mean);
cov_P10_P57_d=(n/(n-1))*(linsub_P10 - linsub_P10_mean)*(linsub_P57 - linsub_P57_mean);
cov_P10_P61_d=(n/(n-1))*(linsub_P10 - linsub_P10_mean)*(linsub_P61 - linsub_P61_mean);
cov_P10_P64_d=(n/(n-1))*(linsub_P10 - linsub_P10_mean)*(linsub_P64 - linsub_P64_mean);
cov_P10_P66_d=(n/(n-1))*(linsub_P10 - linsub_P10_mean)*(linsub_P66 - linsub_P66_mean);
cov_P35_P36_d=(n/(n-1))*(linsub_P35 - linsub_P35_mean)*(linsub_P36 - linsub_P36_mean);
cov_P35_P51_d=(n/(n-1))*(linsub_P35 - linsub_P35_mean)*(linsub_P51 - linsub_P51_mean);
cov_P35_P52_d=(n/(n-1))*(linsub_P35 - linsub_P35_mean)*(linsub_P52 - linsub_P52_mean);
cov_P35_P53_d=(n/(n-1))*(linsub_P35 - linsub_P35_mean)*(linsub_P53 - linsub_P53_mean);
cov_P35_P55_d=(n/(n-1))*(linsub_P35 - linsub_P35_mean)*(linsub_P55 - linsub_P55_mean);
cov_P35_P56_d=(n/(n-1))*(linsub_P35 - linsub_P35_mean)*(linsub_P56 - linsub_P56_mean);
cov_P35_P57_d=(n/(n-1))*(linsub_P35 - linsub_P35_mean)*(linsub_P57 - linsub_P57_mean);
cov_P35_P61_d=(n/(n-1))*(linsub_P35 - linsub_P35_mean)*(linsub_P61 - linsub_P61_mean);
cov_P35_P64_d=(n/(n-1))*(linsub_P35 - linsub_P35_mean)*(linsub_P64 - linsub_P64_mean);
cov_P35_P66_d=(n/(n-1))*(linsub_P35 - linsub_P35_mean)*(linsub_P66 - linsub_P66_mean);
cov_P36_P51_d=(n/(n-1))*(linsub_P36 - linsub_P36_mean)*(linsub_P51 - linsub_P51_mean);
cov_P36_P52_d=(n/(n-1))*(linsub_P36 - linsub_P36_mean)*(linsub_P52 - linsub_P52_mean);
cov_P36_P53_d=(n/(n-1))*(linsub_P36 - linsub_P36_mean)*(linsub_P53 - linsub_P53_mean);
cov_P36_P55_d=(n/(n-1))*(linsub_P36 - linsub_P36_mean)*(linsub_P55 - linsub_P55_mean);
cov_P36_P56_d=(n/(n-1))*(linsub_P36 - linsub_P36_mean)*(linsub_P56 - linsub_P56_mean);
cov_P36_P57_d=(n/(n-1))*(linsub_P36 - linsub_P36_mean)*(linsub_P57 - linsub_P57_mean);
cov_P36_P61_d=(n/(n-1))*(linsub_P36 - linsub_P36_mean)*(linsub_P61 - linsub_P61_mean);
cov_P36_P64_d=(n/(n-1))*(linsub_P36 - linsub_P36_mean)*(linsub_P64 - linsub_P64_mean);
cov_P36_P66_d=(n/(n-1))*(linsub_P36 - linsub_P36_mean)*(linsub_P66 - linsub_P66_mean);
cov_P51_P52_d=(n/(n-1))*(linsub_P51 - linsub_P51_mean)*(linsub_P52 - linsub_P52_mean);
cov_P51_P53_d=(n/(n-1))*(linsub_P51 - linsub_P51_mean)*(linsub_P53 - linsub_P53_mean);
cov_P51_P55_d=(n/(n-1))*(linsub_P51 - linsub_P51_mean)*(linsub_P55 - linsub_P55_mean);
cov_P51_P56_d=(n/(n-1))*(linsub_P51 - linsub_P51_mean)*(linsub_P56 - linsub_P56_mean);
cov_P51_P57_d=(n/(n-1))*(linsub_P51 - linsub_P51_mean)*(linsub_P57 - linsub_P57_mean);
cov_P51_P61_d=(n/(n-1))*(linsub_P51 - linsub_P51_mean)*(linsub_P61 - linsub_P61_mean);
cov_P51_P64_d=(n/(n-1))*(linsub_P51 - linsub_P51_mean)*(linsub_P64 - linsub_P64_mean);
cov_P51_P66_d=(n/(n-1))*(linsub_P51 - linsub_P51_mean)*(linsub_P66 - linsub_P66_mean);
cov_P52_P53_d=(n/(n-1))*(linsub_P52 - linsub_P52_mean)*(linsub_P53 - linsub_P53_mean);
cov_P52_P55_d=(n/(n-1))*(linsub_P52 - linsub_P52_mean)*(linsub_P55 - linsub_P55_mean);
cov_P52_P56_d=(n/(n-1))*(linsub_P52 - linsub_P52_mean)*(linsub_P56 - linsub_P56_mean);
cov_P52_P57_d=(n/(n-1))*(linsub_P52 - linsub_P52_mean)*(linsub_P57 - linsub_P57_mean);
cov_P52_P61_d=(n/(n-1))*(linsub_P52 - linsub_P52_mean)*(linsub_P61 - linsub_P61_mean);
cov_P52_P64_d=(n/(n-1))*(linsub_P52 - linsub_P52_mean)*(linsub_P64 - linsub_P64_mean);
cov_P52_P66_d=(n/(n-1))*(linsub_P52 - linsub_P52_mean)*(linsub_P66 - linsub_P66_mean);
cov_P53_P55_d=(n/(n-1))*(linsub_P53 - linsub_P53_mean)*(linsub_P55 - linsub_P55_mean);
cov_P53_P56_d=(n/(n-1))*(linsub_P53 - linsub_P53_mean)*(linsub_P56 - linsub_P56_mean);
cov_P53_P57_d=(n/(n-1))*(linsub_P53 - linsub_P53_mean)*(linsub_P57 - linsub_P57_mean);
cov_P53_P61_d=(n/(n-1))*(linsub_P53 - linsub_P53_mean)*(linsub_P61 - linsub_P61_mean);
cov_P53_P64_d=(n/(n-1))*(linsub_P53 - linsub_P53_mean)*(linsub_P64 - linsub_P64_mean);
cov_P53_P66_d=(n/(n-1))*(linsub_P53 - linsub_P53_mean)*(linsub_P66 - linsub_P66_mean);
cov_P55_P56_d=(n/(n-1))*(linsub_P55 - linsub_P55_mean)*(linsub_P56 - linsub_P56_mean);
cov_P55_P57_d=(n/(n-1))*(linsub_P55 - linsub_P55_mean)*(linsub_P57 - linsub_P57_mean);
cov_P55_P61_d=(n/(n-1))*(linsub_P55 - linsub_P55_mean)*(linsub_P61 - linsub_P61_mean);
cov_P55_P64_d=(n/(n-1))*(linsub_P55 - linsub_P55_mean)*(linsub_P64 - linsub_P64_mean);
cov_P55_P66_d=(n/(n-1))*(linsub_P55 - linsub_P55_mean)*(linsub_P66 - linsub_P66_mean);
cov_P56_P57_d=(n/(n-1))*(linsub_P56 - linsub_P56_mean)*(linsub_P57 - linsub_P57_mean);
cov_P56_P61_d=(n/(n-1))*(linsub_P56 - linsub_P56_mean)*(linsub_P61 - linsub_P61_mean);
cov_P56_P64_d=(n/(n-1))*(linsub_P56 - linsub_P56_mean)*(linsub_P64 - linsub_P64_mean);
cov_P56_P66_d=(n/(n-1))*(linsub_P56 - linsub_P56_mean)*(linsub_P66 - linsub_P66_mean);
cov_P57_P61_d=(n/(n-1))*(linsub_P57 - linsub_P57_mean)*(linsub_P61 - linsub_P61_mean);
cov_P57_P64_d=(n/(n-1))*(linsub_P57 - linsub_P57_mean)*(linsub_P64 - linsub_P64_mean);
cov_P57_P66_d=(n/(n-1))*(linsub_P57 - linsub_P57_mean)*(linsub_P66 - linsub_P66_mean);
cov_P61_P64_d=(n/(n-1))*(linsub_P61 - linsub_P61_mean)*(linsub_P64 - linsub_P64_mean);
cov_P61_P66_d=(n/(n-1))*(linsub_P61 - linsub_P61_mean)*(linsub_P66 - linsub_P66_mean);
cov_P64_P66_d=(n/(n-1))*(linsub_P64 - linsub_P64_mean)*(linsub_P66 - linsub_P66_mean);

/* ROPC index */
cov_P12_P14_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P14 - linsub_P14_mean);
cov_P12_P15_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P15 - linsub_P15_mean);
cov_P12_P20_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P20 - linsub_P20_mean);
cov_P12_P22_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P22 - linsub_P22_mean);
cov_P12_P23_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P23 - linsub_P23_mean);
cov_P12_P24_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P24 - linsub_P24_mean);
cov_P12_P30_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P30 - linsub_P30_mean);
cov_P12_P32_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P32 - linsub_P32_mean);
cov_P12_P33_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P33 - linsub_P33_mean);
cov_P12_P42_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P42 - linsub_P42_mean);
cov_P12_P44_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P44 - linsub_P44_mean);
cov_P12_P65_d=(n/(n-1))*(linsub_P12 - linsub_P12_mean)*(linsub_P65 - linsub_P65_mean);
cov_P14_P15_d=(n/(n-1))*(linsub_P14 - linsub_P14_mean)*(linsub_P15 - linsub_P15_mean);
cov_P14_P20_d=(n/(n-1))*(linsub_P14 - linsub_P14_mean)*(linsub_P20 - linsub_P20_mean);
cov_P14_P22_d=(n/(n-1))*(linsub_P14 - linsub_P14_mean)*(linsub_P22 - linsub_P22_mean);
cov_P14_P23_d=(n/(n-1))*(linsub_P14 - linsub_P14_mean)*(linsub_P23 - linsub_P23_mean);
cov_P14_P24_d=(n/(n-1))*(linsub_P14 - linsub_P14_mean)*(linsub_P24 - linsub_P24_mean);
cov_P14_P30_d=(n/(n-1))*(linsub_P14 - linsub_P14_mean)*(linsub_P30 - linsub_P30_mean);
cov_P14_P32_d=(n/(n-1))*(linsub_P14 - linsub_P14_mean)*(linsub_P32 - linsub_P32_mean);
cov_P14_P33_d=(n/(n-1))*(linsub_P14 - linsub_P14_mean)*(linsub_P33 - linsub_P33_mean);
cov_P14_P42_d=(n/(n-1))*(linsub_P14 - linsub_P14_mean)*(linsub_P42 - linsub_P42_mean);
cov_P14_P44_d=(n/(n-1))*(linsub_P14 - linsub_P14_mean)*(linsub_P44 - linsub_P44_mean);
cov_P14_P65_d=(n/(n-1))*(linsub_P14 - linsub_P14_mean)*(linsub_P65 - linsub_P65_mean);
cov_P15_P20_d=(n/(n-1))*(linsub_P15 - linsub_P15_mean)*(linsub_P20 - linsub_P20_mean);
cov_P15_P22_d=(n/(n-1))*(linsub_P15 - linsub_P15_mean)*(linsub_P22 - linsub_P22_mean);
cov_P15_P23_d=(n/(n-1))*(linsub_P15 - linsub_P15_mean)*(linsub_P23 - linsub_P23_mean);
cov_P15_P24_d=(n/(n-1))*(linsub_P15 - linsub_P15_mean)*(linsub_P24 - linsub_P24_mean);
cov_P15_P30_d=(n/(n-1))*(linsub_P15 - linsub_P15_mean)*(linsub_P30 - linsub_P30_mean);
cov_P15_P32_d=(n/(n-1))*(linsub_P15 - linsub_P15_mean)*(linsub_P32 - linsub_P32_mean);
cov_P15_P33_d=(n/(n-1))*(linsub_P15 - linsub_P15_mean)*(linsub_P33 - linsub_P33_mean);
cov_P15_P42_d=(n/(n-1))*(linsub_P15 - linsub_P15_mean)*(linsub_P42 - linsub_P42_mean);
cov_P15_P44_d=(n/(n-1))*(linsub_P15 - linsub_P15_mean)*(linsub_P44 - linsub_P44_mean);
cov_P15_P65_d=(n/(n-1))*(linsub_P15 - linsub_P15_mean)*(linsub_P65 - linsub_P65_mean);
cov_P20_P22_d=(n/(n-1))*(linsub_P20 - linsub_P20_mean)*(linsub_P22 - linsub_P22_mean);
cov_P20_P23_d=(n/(n-1))*(linsub_P20 - linsub_P20_mean)*(linsub_P23 - linsub_P23_mean);
cov_P20_P24_d=(n/(n-1))*(linsub_P20 - linsub_P20_mean)*(linsub_P24 - linsub_P24_mean);
cov_P20_P30_d=(n/(n-1))*(linsub_P20 - linsub_P20_mean)*(linsub_P30 - linsub_P30_mean);
cov_P20_P32_d=(n/(n-1))*(linsub_P20 - linsub_P20_mean)*(linsub_P32 - linsub_P32_mean);
cov_P20_P33_d=(n/(n-1))*(linsub_P20 - linsub_P20_mean)*(linsub_P33 - linsub_P33_mean);
cov_P20_P42_d=(n/(n-1))*(linsub_P20 - linsub_P20_mean)*(linsub_P42 - linsub_P42_mean);
cov_P20_P44_d=(n/(n-1))*(linsub_P20 - linsub_P20_mean)*(linsub_P44 - linsub_P44_mean);
cov_P20_P65_d=(n/(n-1))*(linsub_P20 - linsub_P20_mean)*(linsub_P65 - linsub_P65_mean);
cov_P22_P23_d=(n/(n-1))*(linsub_P22 - linsub_P22_mean)*(linsub_P23 - linsub_P23_mean);
cov_P22_P24_d=(n/(n-1))*(linsub_P22 - linsub_P22_mean)*(linsub_P24 - linsub_P24_mean);
cov_P22_P30_d=(n/(n-1))*(linsub_P22 - linsub_P22_mean)*(linsub_P30 - linsub_P30_mean);
cov_P22_P32_d=(n/(n-1))*(linsub_P22 - linsub_P22_mean)*(linsub_P32 - linsub_P32_mean);
cov_P22_P33_d=(n/(n-1))*(linsub_P22 - linsub_P22_mean)*(linsub_P33 - linsub_P33_mean);
cov_P22_P42_d=(n/(n-1))*(linsub_P22 - linsub_P22_mean)*(linsub_P42 - linsub_P42_mean);
cov_P22_P44_d=(n/(n-1))*(linsub_P22 - linsub_P22_mean)*(linsub_P44 - linsub_P44_mean);
cov_P22_P65_d=(n/(n-1))*(linsub_P22 - linsub_P22_mean)*(linsub_P65 - linsub_P65_mean);
cov_P23_P24_d=(n/(n-1))*(linsub_P23 - linsub_P23_mean)*(linsub_P24 - linsub_P24_mean);
cov_P23_P30_d=(n/(n-1))*(linsub_P23 - linsub_P23_mean)*(linsub_P30 - linsub_P30_mean);
cov_P23_P32_d=(n/(n-1))*(linsub_P23 - linsub_P23_mean)*(linsub_P32 - linsub_P32_mean);
cov_P23_P33_d=(n/(n-1))*(linsub_P23 - linsub_P23_mean)*(linsub_P33 - linsub_P33_mean);
cov_P23_P42_d=(n/(n-1))*(linsub_P23 - linsub_P23_mean)*(linsub_P42 - linsub_P42_mean);
cov_P23_P44_d=(n/(n-1))*(linsub_P23 - linsub_P23_mean)*(linsub_P44 - linsub_P44_mean);
cov_P23_P65_d=(n/(n-1))*(linsub_P23 - linsub_P23_mean)*(linsub_P65 - linsub_P65_mean);
cov_P24_P30_d=(n/(n-1))*(linsub_P24 - linsub_P24_mean)*(linsub_P30 - linsub_P30_mean);
cov_P24_P32_d=(n/(n-1))*(linsub_P24 - linsub_P24_mean)*(linsub_P32 - linsub_P32_mean);
cov_P24_P33_d=(n/(n-1))*(linsub_P24 - linsub_P24_mean)*(linsub_P33 - linsub_P33_mean);
cov_P24_P42_d=(n/(n-1))*(linsub_P24 - linsub_P24_mean)*(linsub_P42 - linsub_P42_mean);
cov_P24_P44_d=(n/(n-1))*(linsub_P24 - linsub_P24_mean)*(linsub_P44 - linsub_P44_mean);
cov_P24_P65_d=(n/(n-1))*(linsub_P24 - linsub_P24_mean)*(linsub_P65 - linsub_P65_mean);
cov_P30_P32_d=(n/(n-1))*(linsub_P30 - linsub_P30_mean)*(linsub_P32 - linsub_P32_mean);
cov_P30_P33_d=(n/(n-1))*(linsub_P30 - linsub_P30_mean)*(linsub_P33 - linsub_P33_mean);
cov_P30_P42_d=(n/(n-1))*(linsub_P30 - linsub_P30_mean)*(linsub_P42 - linsub_P42_mean);
cov_P30_P44_d=(n/(n-1))*(linsub_P30 - linsub_P30_mean)*(linsub_P44 - linsub_P44_mean);
cov_P30_P65_d=(n/(n-1))*(linsub_P30 - linsub_P30_mean)*(linsub_P65 - linsub_P65_mean);
cov_P32_P33_d=(n/(n-1))*(linsub_P32 - linsub_P32_mean)*(linsub_P33 - linsub_P33_mean);
cov_P32_P42_d=(n/(n-1))*(linsub_P32 - linsub_P32_mean)*(linsub_P42 - linsub_P42_mean);
cov_P32_P44_d=(n/(n-1))*(linsub_P32 - linsub_P32_mean)*(linsub_P44 - linsub_P44_mean);
cov_P32_P65_d=(n/(n-1))*(linsub_P32 - linsub_P32_mean)*(linsub_P65 - linsub_P65_mean);
cov_P33_P42_d=(n/(n-1))*(linsub_P33 - linsub_P33_mean)*(linsub_P42 - linsub_P42_mean);
cov_P33_P44_d=(n/(n-1))*(linsub_P33 - linsub_P33_mean)*(linsub_P44 - linsub_P44_mean);
cov_P33_P65_d=(n/(n-1))*(linsub_P33 - linsub_P33_mean)*(linsub_P65 - linsub_P65_mean);
cov_P42_P44_d=(n/(n-1))*(linsub_P42 - linsub_P42_mean)*(linsub_P44 - linsub_P44_mean);
cov_P42_P65_d=(n/(n-1))*(linsub_P42 - linsub_P42_mean)*(linsub_P65 - linsub_P65_mean);
cov_P44_P65_d=(n/(n-1))*(linsub_P44 - linsub_P44_mean)*(linsub_P65 - linsub_P65_mean);

/* Talent index */
cov_P1_P11_d =(n/(n-1))* (linsub_P1 - linsub_P1_mean) *(linsub_P11 - linsub_P11_mean);
cov_P1_P18_d =(n/(n-1))* (linsub_P1 - linsub_P1_mean) *(linsub_P18 - linsub_P18_mean);
cov_P1_P21_d =(n/(n-1))* (linsub_P1 - linsub_P1_mean) *(linsub_P21 - linsub_P21_mean);
cov_P1_P29_d =(n/(n-1))* (linsub_P1 - linsub_P1_mean) *(linsub_P29 - linsub_P29_mean);
cov_P1_P47_d =(n/(n-1))* (linsub_P1 - linsub_P1_mean) *(linsub_P47 - linsub_P47_mean);
cov_P1_P68_d =(n/(n-1))* (linsub_P1 - linsub_P1_mean) *(linsub_P68 - linsub_P68_mean);
cov_P11_P18_d=(n/(n-1))*(linsub_P11 - linsub_P11_mean)*(linsub_P18 - linsub_P18_mean);
cov_P11_P21_d=(n/(n-1))*(linsub_P11 - linsub_P11_mean)*(linsub_P21 - linsub_P21_mean);
cov_P11_P29_d=(n/(n-1))*(linsub_P11 - linsub_P11_mean)*(linsub_P29 - linsub_P29_mean);
cov_P11_P47_d=(n/(n-1))*(linsub_P11 - linsub_P11_mean)*(linsub_P47 - linsub_P47_mean);
cov_P11_P68_d=(n/(n-1))*(linsub_P11 - linsub_P11_mean)*(linsub_P68 - linsub_P68_mean);
cov_P18_P21_d=(n/(n-1))*(linsub_P18 - linsub_P18_mean)*(linsub_P21 - linsub_P21_mean);
cov_P18_P29_d=(n/(n-1))*(linsub_P18 - linsub_P18_mean)*(linsub_P29 - linsub_P29_mean);
cov_P18_P47_d=(n/(n-1))*(linsub_P18 - linsub_P18_mean)*(linsub_P47 - linsub_P47_mean);
cov_P18_P68_d=(n/(n-1))*(linsub_P18 - linsub_P18_mean)*(linsub_P68 - linsub_P68_mean);
cov_P21_P29_d=(n/(n-1))*(linsub_P21 - linsub_P21_mean)*(linsub_P29 - linsub_P29_mean);
cov_P21_P47_d=(n/(n-1))*(linsub_P21 - linsub_P21_mean)*(linsub_P47 - linsub_P47_mean);
cov_P21_P68_d=(n/(n-1))*(linsub_P21 - linsub_P21_mean)*(linsub_P68 - linsub_P68_mean);
cov_P29_P47_d=(n/(n-1))*(linsub_P29 - linsub_P29_mean)*(linsub_P47 - linsub_P47_mean);
cov_P29_P68_d=(n/(n-1))*(linsub_P29 - linsub_P29_mean)*(linsub_P68 - linsub_P68_mean);
cov_P47_P68_d=(n/(n-1))*(linsub_P47 - linsub_P47_mean)*(linsub_P68 - linsub_P68_mean);
run;

proc means data=junk2 noprint nway;
  class agency;
  var var_: cov_:;
output out=stats (drop=_:) sum=;
run;

data stats;
wave1=&wave1;
wave2=&wave2;
  length agency $2;
set stats;
run;
* calculate the two wave-specific means for each percent positive estimate of interest;
* first, acquire the overall estimate (full-sample);
proc sql;
  create table stats_fs1 as
  select agency,
         /* JS index */
         sum(weight_wave&wave1._raked*P4 )/sum(weight_wave&wave1._raked*(P4  ne .)) as P4_1,
         sum(weight_wave&wave1._raked*P5 )/sum(weight_wave&wave1._raked*(P5  ne .)) as P5_1,
         sum(weight_wave&wave1._raked*P13)/sum(weight_wave&wave1._raked*(P13 ne .)) as P13_1,
         sum(weight_wave&wave1._raked*P63)/sum(weight_wave&wave1._raked*(P63 ne .)) as P63_1,
         sum(weight_wave&wave1._raked*P67)/sum(weight_wave&wave1._raked*(P67 ne .)) as P67_1,
         sum(weight_wave&wave1._raked*P69)/sum(weight_wave&wave1._raked*(P69 ne .)) as P69_1,
         sum(weight_wave&wave1._raked*P70)/sum(weight_wave&wave1._raked*(P70 ne .)) as P70_1,
		 /* Leadership index */
         sum(weight_wave&wave1._raked*P10)/sum(weight_wave&wave1._raked*(P10 ne .)) as P10_1,
         sum(weight_wave&wave1._raked*P35)/sum(weight_wave&wave1._raked*(P35 ne .)) as P35_1,
         sum(weight_wave&wave1._raked*P36)/sum(weight_wave&wave1._raked*(P36 ne .)) as P36_1,
         sum(weight_wave&wave1._raked*P51)/sum(weight_wave&wave1._raked*(P51 ne .)) as P51_1,
         sum(weight_wave&wave1._raked*P52)/sum(weight_wave&wave1._raked*(P52 ne .)) as P52_1,
         sum(weight_wave&wave1._raked*P53)/sum(weight_wave&wave1._raked*(P53 ne .)) as P53_1,
         sum(weight_wave&wave1._raked*P55)/sum(weight_wave&wave1._raked*(P55 ne .)) as P55_1,
         sum(weight_wave&wave1._raked*P56)/sum(weight_wave&wave1._raked*(P56 ne .)) as P56_1,
         sum(weight_wave&wave1._raked*P57)/sum(weight_wave&wave1._raked*(P57 ne .)) as P57_1,
         sum(weight_wave&wave1._raked*P61)/sum(weight_wave&wave1._raked*(P61 ne .)) as P61_1,
         sum(weight_wave&wave1._raked*P64)/sum(weight_wave&wave1._raked*(P64 ne .)) as P64_1,
         sum(weight_wave&wave1._raked*P66)/sum(weight_wave&wave1._raked*(P66 ne .)) as P66_1,
		 /* ROPC index */
         sum(weight_wave&wave1._raked*P12)/sum(weight_wave&wave1._raked*(P12 ne .)) as P12_1,
         sum(weight_wave&wave1._raked*P14)/sum(weight_wave&wave1._raked*(P14 ne .)) as P14_1,
         sum(weight_wave&wave1._raked*P15)/sum(weight_wave&wave1._raked*(P15 ne .)) as P15_1,
         sum(weight_wave&wave1._raked*P20)/sum(weight_wave&wave1._raked*(P20 ne .)) as P20_1,
         sum(weight_wave&wave1._raked*P22)/sum(weight_wave&wave1._raked*(P22 ne .)) as P22_1,
         sum(weight_wave&wave1._raked*P23)/sum(weight_wave&wave1._raked*(P23 ne .)) as P23_1,
         sum(weight_wave&wave1._raked*P24)/sum(weight_wave&wave1._raked*(P24 ne .)) as P24_1,
         sum(weight_wave&wave1._raked*P30)/sum(weight_wave&wave1._raked*(P30 ne .)) as P30_1,
         sum(weight_wave&wave1._raked*P32)/sum(weight_wave&wave1._raked*(P32 ne .)) as P32_1,
         sum(weight_wave&wave1._raked*P33)/sum(weight_wave&wave1._raked*(P33 ne .)) as P33_1,
         sum(weight_wave&wave1._raked*P42)/sum(weight_wave&wave1._raked*(P42 ne .)) as P42_1,
         sum(weight_wave&wave1._raked*P44)/sum(weight_wave&wave1._raked*(P44 ne .)) as P44_1,
         sum(weight_wave&wave1._raked*P65)/sum(weight_wave&wave1._raked*(P65 ne .)) as P65_1,
         /* Talent index */
         sum(weight_wave&wave1._raked*P1 )/sum(weight_wave&wave1._raked*(P1  ne .)) as P1_1,
         sum(weight_wave&wave1._raked*P11)/sum(weight_wave&wave1._raked*(P11 ne .)) as P11_1,
         sum(weight_wave&wave1._raked*P18)/sum(weight_wave&wave1._raked*(P18 ne .)) as P18_1,
         sum(weight_wave&wave1._raked*P21)/sum(weight_wave&wave1._raked*(P21 ne .)) as P21_1,
         sum(weight_wave&wave1._raked*P29)/sum(weight_wave&wave1._raked*(P29 ne .)) as P29_1,
         sum(weight_wave&wave1._raked*P47)/sum(weight_wave&wave1._raked*(P47 ne .)) as P47_1,
         sum(weight_wave&wave1._raked*P68)/sum(weight_wave&wave1._raked*(P68 ne .)) as P68_1

	from junk2
group by agency;

  create table stats_fs2 as
  select agency,
         /* JS index */
         sum(weight_wave&wave2._raked*P4 )/sum(weight_wave&wave2._raked*(P4  ne .)) as P4_2,
         sum(weight_wave&wave2._raked*P5 )/sum(weight_wave&wave2._raked*(P5  ne .)) as P5_2,
         sum(weight_wave&wave2._raked*P13)/sum(weight_wave&wave2._raked*(P13 ne .)) as P13_2,
         sum(weight_wave&wave2._raked*P63)/sum(weight_wave&wave2._raked*(P63 ne .)) as P63_2,
         sum(weight_wave&wave2._raked*P67)/sum(weight_wave&wave2._raked*(P67 ne .)) as P67_2,
         sum(weight_wave&wave2._raked*P69)/sum(weight_wave&wave2._raked*(P69 ne .)) as P69_2,
         sum(weight_wave&wave2._raked*P70)/sum(weight_wave&wave2._raked*(P70 ne .)) as P70_2,
		 /* Leadership index */
         sum(weight_wave&wave2._raked*P10)/sum(weight_wave&wave2._raked*(P10 ne .)) as P10_2,
         sum(weight_wave&wave2._raked*P35)/sum(weight_wave&wave2._raked*(P35 ne .)) as P35_2,
         sum(weight_wave&wave2._raked*P36)/sum(weight_wave&wave2._raked*(P36 ne .)) as P36_2,
         sum(weight_wave&wave2._raked*P51)/sum(weight_wave&wave2._raked*(P51 ne .)) as P51_2,
         sum(weight_wave&wave2._raked*P52)/sum(weight_wave&wave2._raked*(P52 ne .)) as P52_2,
         sum(weight_wave&wave2._raked*P53)/sum(weight_wave&wave2._raked*(P53 ne .)) as P53_2,
         sum(weight_wave&wave2._raked*P55)/sum(weight_wave&wave2._raked*(P55 ne .)) as P55_2,
         sum(weight_wave&wave2._raked*P56)/sum(weight_wave&wave2._raked*(P56 ne .)) as P56_2,
         sum(weight_wave&wave2._raked*P57)/sum(weight_wave&wave2._raked*(P57 ne .)) as P57_2,
         sum(weight_wave&wave2._raked*P61)/sum(weight_wave&wave2._raked*(P61 ne .)) as P61_2,
         sum(weight_wave&wave2._raked*P64)/sum(weight_wave&wave2._raked*(P64 ne .)) as P64_2,
         sum(weight_wave&wave2._raked*P66)/sum(weight_wave&wave2._raked*(P66 ne .)) as P66_2,
		 /* ROPC index */
         sum(weight_wave&wave2._raked*P12)/sum(weight_wave&wave2._raked*(P12 ne .)) as P12_2,
         sum(weight_wave&wave2._raked*P14)/sum(weight_wave&wave2._raked*(P14 ne .)) as P14_2,
         sum(weight_wave&wave2._raked*P15)/sum(weight_wave&wave2._raked*(P15 ne .)) as P15_2,
         sum(weight_wave&wave2._raked*P20)/sum(weight_wave&wave2._raked*(P20 ne .)) as P20_2,
         sum(weight_wave&wave2._raked*P22)/sum(weight_wave&wave2._raked*(P22 ne .)) as P22_2,
         sum(weight_wave&wave2._raked*P23)/sum(weight_wave&wave2._raked*(P23 ne .)) as P23_2,
         sum(weight_wave&wave2._raked*P24)/sum(weight_wave&wave2._raked*(P24 ne .)) as P24_2,
         sum(weight_wave&wave2._raked*P30)/sum(weight_wave&wave2._raked*(P30 ne .)) as P30_2,
         sum(weight_wave&wave2._raked*P32)/sum(weight_wave&wave2._raked*(P32 ne .)) as P32_2,
         sum(weight_wave&wave2._raked*P33)/sum(weight_wave&wave2._raked*(P33 ne .)) as P33_2,
         sum(weight_wave&wave2._raked*P42)/sum(weight_wave&wave2._raked*(P42 ne .)) as P42_2,
         sum(weight_wave&wave2._raked*P44)/sum(weight_wave&wave2._raked*(P44 ne .)) as P44_2,
         sum(weight_wave&wave2._raked*P65)/sum(weight_wave&wave2._raked*(P65 ne .)) as P65_2,
         /* Talent index */
         sum(weight_wave&wave2._raked*P1 )/sum(weight_wave&wave2._raked*(P1  ne .)) as P1_2,
         sum(weight_wave&wave2._raked*P11)/sum(weight_wave&wave2._raked*(P11 ne .)) as P11_2,
         sum(weight_wave&wave2._raked*P18)/sum(weight_wave&wave2._raked*(P18 ne .)) as P18_2,
         sum(weight_wave&wave2._raked*P21)/sum(weight_wave&wave2._raked*(P21 ne .)) as P21_2,
         sum(weight_wave&wave2._raked*P29)/sum(weight_wave&wave2._raked*(P29 ne .)) as P29_2,
         sum(weight_wave&wave2._raked*P47)/sum(weight_wave&wave2._raked*(P47 ne .)) as P47_2,
         sum(weight_wave&wave2._raked*P68)/sum(weight_wave&wave2._raked*(P68 ne .)) as P68_2

	from junk2
group by agency;
quit;
data stats_fs;
wave1=&wave1;
wave2=&wave2;

  merge stats_fs1 (in=a)
        stats_fs2 (in=b);
array part1{*}   P4_1 -- P68_1;
array part2{*}   P4_2 -- P68_2;
array part3{*}   P4_d  P5_d  P13_d P63_d P67_d P69_d P70_d
                 P10_d P35_d P36_d P51_d P52_d P53_d P55_d P56_d P57_d P61_d P64_d P66_d
                 P12_d P14_d P15_d P20_d P22_d P23_d P24_d P30_d P32_d P33_d P42_d P44_d P65_d
                 P1_d  P11_d P18_d P21_d P29_d P47_d P68_d;

do i=1 to dim(part1);
   part3{i}=part2{i} - part1{i};
end;
drop i P4_1 -- P68_2;
run;
* merge this back with stats data set;
data stats;
  merge stats_fs stats;
  by agency wave1 wave2;
run;
proc datasets nolist;
  append data=stats
         base=stats_study2_TSL; 
  * housecleaning;
  delete stats stats_f:;
run;
quit;
%end; * closes the wave loop;
%mend;

options nonotes;
%getstats_TSL(wave_start=1,wave_end=10);
options notes;



*** test each wave-over-wave threshold in the given index via PROC IML;
%macro Wald_stats(condition=,wave_start=,wave_end=);
* housecleaning;
proc sql; drop table M1_Wald_stats_&condition; quit;
%do  wave2=%eval(&wave_start+1) %to &wave_end;
%let wave1=%eval(&wave2-1);
%let agylist=%str(EP FT NN);
%let num=1;
%let agency=%scan(&agylist,&num);
%do %while (&agency ne );

* A) Job Satisfaction Index;
proc iml;
use stats_study2_&condition (where=(wave1=&wave1 and wave2=&wave2 and agency="&agency"));
* create the difference vector;
read all var {P4_d P5_d P13_d P63_d P67_d P69_d P70_d} into D;
print D;

* create the covariance matrix from the summary variances and covariances;
* initialize the matrix that we will soon populate;
S=I(7); * creates a 7 x 7 identity matrix;
print s;
read all var {var_P4_d var_P5_d var_P13_d var_P63_d var_P67_d var_P69_d var_P70_d} into vars;
print vars;
read all var {cov_P4_P13_d  cov_P4_P5_d   cov_P4_P63_d  cov_P4_P67_d  cov_P4_P69_d cov_P4_P70_d
              cov_P5_P13_d  cov_P5_P63_d  cov_P5_P67_d  cov_P5_P69_d  cov_P5_P70_d 
              cov_P13_P63_d cov_P13_P67_d cov_P13_P69_d cov_P13_P70_d
              cov_P63_P67_d cov_P63_P69_d cov_P63_P70_d
              cov_P67_P69_d cov_P67_P70_d
              cov_P69_P70_d} into covars;
print covars;
S=diag(vars); * puts the variance terms along the diagonal;
S[2,1]=covars[1]; 
S[3,1]=covars[2]; S[3,2]=covars[7];
S[4,1]=covars[3]; S[4,2]=covars[8];  S[4,3]=covars[12];
S[5,1]=covars[4]; S[5,2]=covars[9];  S[5,3]=covars[13]; S[5,4]=covars[16];
S[6,1]=covars[5]; S[6,2]=covars[10]; S[6,3]=covars[14]; S[6,4]=covars[17]; S[6,5]=covars[19];
S[7,1]=covars[6]; S[7,2]=covars[11]; S[7,3]=covars[15]; S[7,4]=covars[18]; S[7,5]=covars[20]; S[7,6]=covars[21];
S=S+S`-diag(vars);
print S[format=10.8];

W=D*inv(S)*D`;
print w;
create outA from W [ colname="Wald" ];
append from W;
quit;
data outA;
  set outA;
wave1=&wave1;
wave2=&wave2;
length agency $2 index $1;
agency="&agency";
index='A';
run;

* B) Leadership Index;
proc iml;
use stats_study2_&condition (where=(wave1=&wave1 and wave2=&wave2 and agency="&agency"));
* create the difference vector;
read all var {P10_d P35_d P36_d P51_d P52_d P53_d P55_d P56_d P57_d P61_d P64_d P66_d} into D;
print D;

* create the covariance matrix from the summary variances and covariances;
* initialize the matrix that we will soon populate;
S=I(12); * creates a 12 x 12 identity matrix;
print s;
read all var {var_P10_d var_P35_d var_P36_d var_P51_d var_P52_d var_P53_d var_P55_d var_P56_d var_P57_d var_P61_d var_P64_d var_P66_d} into vars;
print vars;
read all var {cov_P10_P35_d cov_P10_P36_d cov_P10_P51_d cov_P10_P52_d cov_P10_P53_d cov_P10_P55_d cov_P10_P56_d cov_P10_P57_d cov_P10_P61_d cov_P10_P64_d cov_P10_P66_d
              cov_P35_P36_d cov_P35_P51_d cov_P35_P52_d cov_P35_P53_d cov_P35_P55_d cov_P35_P56_d cov_P35_P57_d cov_P35_P61_d cov_P35_P64_d cov_P35_P66_d
              cov_P36_P51_d cov_P36_P52_d cov_P36_P53_d cov_P36_P55_d cov_P36_P56_d cov_P36_P57_d cov_P36_P61_d cov_P36_P64_d cov_P36_P66_d
              cov_P51_P52_d cov_P51_P53_d cov_P51_P55_d cov_P51_P56_d cov_P51_P57_d cov_P51_P61_d cov_P51_P64_d cov_P51_P66_d
              cov_P52_P53_d cov_P52_P55_d cov_P52_P56_d cov_P52_P57_d cov_P52_P61_d cov_P52_P64_d cov_P52_P66_d
              cov_P53_P55_d cov_P53_P56_d cov_P53_P57_d cov_P53_P61_d cov_P53_P64_d cov_P53_P66_d
              cov_P55_P56_d cov_P55_P57_d cov_P55_P61_d cov_P55_P64_d cov_P55_P66_d
              cov_P56_P57_d cov_P56_P61_d cov_P56_P64_d cov_P56_P66_d
              cov_P57_P61_d cov_P57_P64_d cov_P57_P66_d
              cov_P61_P64_d cov_P61_P66_d
              cov_P64_P66_d} into covars;
print covars;
S=diag(vars); * puts the variance terms along the diagonal;
S[2,1]=covars[1]; 
S[3,1]=covars[2];  S[3,2]=covars[12];
S[4,1]=covars[3];  S[4,2]=covars[13]; S[4,3]=covars[22];
S[5,1]=covars[4];  S[5,2]=covars[14]; S[5,3]=covars[23]; S[5,4]=covars[31];
S[6,1]=covars[5];  S[6,2]=covars[15]; S[6,3]=covars[24]; S[6,4]=covars[32]; S[6,5]=covars[39];
S[7,1]=covars[6];  S[7,2]=covars[16]; S[7,3]=covars[25]; S[7,4]=covars[33]; S[7,5]=covars[40]; S[7,6]=covars[46];
S[8,1]=covars[7];  S[8,2]=covars[17]; S[8,3]=covars[26]; S[8,4]=covars[34]; S[8,5]=covars[41]; S[8,6]=covars[47]; S[8,7]=covars[52];
S[9,1]=covars[8];  S[9,2]=covars[18]; S[9,3]=covars[27]; S[9,4]=covars[35]; S[9,5]=covars[42]; S[9,6]=covars[48]; S[9,7]=covars[53]; S[9,8]=covars[57];
S[10,1]=covars[9]; S[10,2]=covars[19];S[10,3]=covars[28];S[10,4]=covars[36];S[10,5]=covars[43];S[10,6]=covars[49];S[10,7]=covars[54];S[10,8]=covars[58];S[10,9]=covars[61];
S[11,1]=covars[10];S[11,2]=covars[20];S[11,3]=covars[29];S[11,4]=covars[37];S[11,5]=covars[44];S[11,6]=covars[50];S[11,7]=covars[55];S[11,8]=covars[59];S[11,9]=covars[62];S[11,10]=covars[64];
S[12,1]=covars[11];S[12,2]=covars[21];S[12,3]=covars[30];S[12,4]=covars[38];S[12,5]=covars[45];S[12,6]=covars[51];S[12,7]=covars[56];S[12,8]=covars[60];S[12,9]=covars[63];S[12,10]=covars[65];S[12,11]=covars[66];

S=S+S`-diag(vars);
print S[format=10.8];

W=D*inv(S)*D`;
print w;
create outB from W [ colname="Wald" ];
append from W;
quit;
data outB;
  set outB;
wave1=&wave1;
wave2=&wave2;
length agency $2 index $1;
agency="&agency";
index='B';
run;

* C) ROPC Index;
proc iml;
use stats_study2_&condition (where=(wave1=&wave1 and wave2=&wave2 and agency="&agency"));
* create the difference vector;
read all var {P12_d P14_d P15_d P20_d P22_d P23_d P24_d P30_d P32_d P33_d P42_d P44_d P65_d} into D;
print D;

* create the covariance matrix from the summary variances and covariances;
* initialize the matrix that we will soon populate;
S=I(13); * creates a 13 x 13 identity matrix;
print s;
read all var {var_P12_d var_P14_d var_P15_d var_P20_d var_P22_d var_P23_d var_P24_d var_P30_d var_P32_d var_P33_d var_P42_d var_P44_d var_P65_d} into vars;
print vars;
read all var {cov_P12_P14_d cov_P12_P15_d cov_P12_P20_d cov_P12_P22_d cov_P12_P23_d cov_P12_P24_d cov_P12_P30_d cov_P12_P32_d cov_P12_P33_d cov_P12_P42_d cov_P12_P44_d cov_P12_P65_d
              cov_P14_P15_d cov_P14_P20_d cov_P14_P22_d cov_P14_P23_d cov_P14_P24_d cov_P14_P30_d cov_P14_P32_d cov_P14_P33_d cov_P14_P42_d cov_P14_P44_d cov_P14_P65_d
              cov_P15_P20_d cov_P15_P22_d cov_P15_P23_d cov_P15_P24_d cov_P15_P30_d cov_P15_P32_d cov_P15_P33_d cov_P15_P42_d cov_P15_P44_d cov_P15_P65_d
              cov_P20_P22_d cov_P20_P23_d cov_P20_P24_d cov_P20_P30_d cov_P20_P32_d cov_P20_P33_d cov_P20_P42_d cov_P20_P44_d cov_P20_P65_d
              cov_P22_P23_d cov_P22_P24_d cov_P22_P30_d cov_P22_P32_d cov_P22_P33_d cov_P22_P42_d cov_P22_P44_d cov_P22_P65_d
              cov_P23_P24_d cov_P23_P30_d cov_P23_P32_d cov_P23_P33_d cov_P23_P42_d cov_P23_P44_d cov_P23_P65_d
              cov_P24_P30_d cov_P24_P32_d cov_P24_P33_d cov_P24_P42_d cov_P24_P44_d cov_P24_P65_d
              cov_P30_P32_d cov_P30_P33_d cov_P30_P42_d cov_P30_P44_d cov_P30_P65_d
              cov_P32_P33_d cov_P32_P42_d cov_P32_P44_d cov_P32_P65_d
              cov_P33_P42_d cov_P33_P44_d cov_P33_P65_d
              cov_P42_P44_d cov_P42_P65_d 
              cov_P44_P65_d} into covars;
print covars;
S=diag(vars); * puts the variance terms along the diagonal;
S[2,1]=covars[1]; 
S[3,1]=covars[2];  S[3,2]=covars[13];
S[4,1]=covars[3];  S[4,2]=covars[14]; S[4,3]=covars[24];
S[5,1]=covars[4];  S[5,2]=covars[15]; S[5,3]=covars[25]; S[5,4]=covars[34];
S[6,1]=covars[5];  S[6,2]=covars[16]; S[6,3]=covars[26]; S[6,4]=covars[35]; S[6,5]=covars[43];
S[7,1]=covars[6];  S[7,2]=covars[17]; S[7,3]=covars[27]; S[7,4]=covars[36]; S[7,5]=covars[44]; S[7,6]=covars[51];
S[8,1]=covars[7];  S[8,2]=covars[18]; S[8,3]=covars[28]; S[8,4]=covars[37]; S[8,5]=covars[45]; S[8,6]=covars[52]; S[8,7]=covars[58];
S[9,1]=covars[8];  S[9,2]=covars[19]; S[9,3]=covars[29]; S[9,4]=covars[38]; S[9,5]=covars[46]; S[9,6]=covars[53]; S[9,7]=covars[59]; S[9,8]=covars[64];
S[10,1]=covars[9]; S[10,2]=covars[20];S[10,3]=covars[30];S[10,4]=covars[39];S[10,5]=covars[47];S[10,6]=covars[54];S[10,7]=covars[60];S[10,8]=covars[65];S[10,9]=covars[69];
S[11,1]=covars[10];S[11,2]=covars[21];S[11,3]=covars[31];S[11,4]=covars[40];S[11,5]=covars[48];S[11,6]=covars[55];S[11,7]=covars[61];S[11,8]=covars[66];S[11,9]=covars[70];S[11,10]=covars[73];
S[12,1]=covars[11];S[12,2]=covars[22];S[12,3]=covars[32];S[12,4]=covars[41];S[12,5]=covars[49];S[12,6]=covars[56];S[12,7]=covars[62];S[12,8]=covars[67];S[12,9]=covars[71];S[12,10]=covars[74];S[12,11]=covars[76];
S[13,1]=covars[12];S[13,2]=covars[23];S[13,3]=covars[33];S[13,4]=covars[42];S[13,5]=covars[50];S[13,6]=covars[57];S[13,7]=covars[63];S[13,8]=covars[68];S[13,9]=covars[72];S[13,10]=covars[75];S[13,11]=covars[77];S[13,12]=covars[78];


S=S+S`-diag(vars);
print S[format=10.8];

W=D*inv(S)*D`;
print w;
create outC from W [ colname="Wald" ];
append from W;
quit;
data outC;
  set outC;
wave1=&wave1;
wave2=&wave2;
length agency $2 index $1;
agency="&agency";
index='C';
run;

* D) Talent Index;
proc iml;
use stats_study2_&condition (where=(wave1=&wave1 and wave2=&wave2 and agency="&agency"));
* create the difference vector;
read all var {P1_d P11_d P18_d P21_d P29_d P47_d P68_d} into D;
print D;

* initialize the variance-covariance matrix that we will soon populate;
S=I(7); * creates a 7 x 7 identity matrix;
print s;
read all var {var_P1_d var_P11_d var_P18_d var_P21_d var_P29_d var_P47_d var_P68_d} into vars;
print vars;
read all var {cov_P1_P18_d  cov_P1_P11_d  cov_P1_P21_d  cov_P1_P29_d  cov_P1_P47_d  cov_P1_P68_d
              cov_P11_P18_d cov_P11_P21_d cov_P11_P29_d cov_P11_P47_d cov_P11_P68_d 
              cov_P18_P21_d cov_P18_P29_d cov_P18_P47_d cov_P18_P68_d
              cov_P21_P29_d cov_P21_P47_d cov_P21_P68_d
              cov_P29_P47_d cov_P29_P68_d
              cov_P47_P68_d} into covars;
print covars;
S=diag(vars); * puts the variance terms along the diagonal;
S[2,1]=covars[1]; 
S[3,1]=covars[2]; S[3,2]=covars[7];
S[4,1]=covars[3]; S[4,2]=covars[8];  S[4,3]=covars[12];
S[5,1]=covars[4]; S[5,2]=covars[9];  S[5,3]=covars[13]; S[5,4]=covars[16];
S[6,1]=covars[5]; S[6,2]=covars[10]; S[6,3]=covars[14]; S[6,4]=covars[17]; S[6,5]=covars[19];
S[7,1]=covars[6]; S[7,2]=covars[11]; S[7,3]=covars[15]; S[7,4]=covars[18]; S[7,5]=covars[20]; S[7,6]=covars[21];
S=S+S`-diag(vars);
print S[format=10.8];

W=D*inv(S)*D`;
print w;
create outD from W [ colname="Wald" ];
append from W;
quit;
data outD;
  set outD;
wave1=&wave1;
wave2=&wave2;
length agency $2 index $1;
agency="&agency";
index='D';
run;
* stack all four data sets;
data out;
  set outA outB outC outD;
run;
proc datasets nolist;
  append base=M1_Wald_stats_&condition data=out;
  delete out;
run;
quit;
%let num=%eval(&num+1);
%let agency=%scan(&agylist,&num);
%end; * close agency list loop;
%end; * close wave2 loop;
%mend;

options nonotes;
ods listing close;
%Wald_stats(condition=TSL,wave_start=1,wave_end=10);
ods listing;
options notes;



*** create a summary data set with each condition x iteration x agency being the row factors;
data junk1;
length condition $4;
  set stats_study2_TSL;
condition='TSL';
keep condition wave1 wave2 agency P4_wave1_fs -- indexD_wave2_fs var_P4_wave1_fs -- var_indexD_wave2_fs;
run;
proc sort data=junk1; by condition agency; run;
proc transpose data=junk1 out=junk1;
  by condition agency;
  ID wave1;
  var P4_wave1_fs -- indexD_wave1_fs var_P4_wave1_fs -- var_indexD_wave1_fs;
run;
data junk2;
  set part1;
keep condition agency wave2 P4_wave2_fs -- indexD_wave2_fs var_P4_wave2_fs -- var_indexD_wave2_fs;
where wave2=10;
run;
proc sort data=junk2; by condition agency; run;
proc transpose data=junk2 out=junk2;
  by condition agency;
  ID wave2;
  var P4_wave2_fs -- indexD_wave2_fs var_P4_wave2_fs -- var_indexD_wave2_fs;
run;
* merge the two data sets;
data junk;
  merge junk1 junk2;
length variable $20;
* do one thing for the variance terms;
if substr(_NAME_,1,1)='v' then do;
    variable=substr(_NAME_,1,find(_NAME_,'_wave2_fs')-1);
end;
* do another for the mean terms;
  else do;
    variable="mean_"||substr(_NAME_,1,find(_NAME_,'_wave2_fs')-1);
  end;
run;
proc transpose data=junk out=junk;
  by condition agency;
  ID variable;
  var  _1 -- _10;
run;
%macro mergeit();
data junk;
merge
%do wave=1 %to 10;
junk (where=(_NAME_="_&wave")  
rename=(
/* JS Index */
mean_P4=mean_P4_wave&wave mean_P5=mean_P5_wave&wave mean_P13=mean_P13_wave&wave mean_P63=mean_P63_wave&wave mean_P67=mean_P67_wave&wave mean_P69=mean_P69_wave&wave mean_P70=mean_P70_wave&wave mean_indexA=mean_indexA_wave&wave var_indexA=var_indexA_wave&wave
/* Leadership Index */
mean_P10=mean_P10_wave&wave mean_P35=mean_P35_wave&wave mean_P36=mean_P36_wave&wave mean_P51=mean_P51_wave&wave mean_P52=mean_P52_wave&wave mean_P53=mean_P53_wave&wave mean_P55=mean_P55_wave&wave mean_P56=mean_P56_wave&wave mean_P57=mean_P57_wave&wave mean_P61=mean_P61_wave&wave mean_P64=mean_P64_wave&wave mean_P66=mean_P66_wave&wave mean_indexB=mean_indexB_wave&wave var_indexB=var_indexB_wave&wave
/* ROPC Index */
mean_P12=mean_P12_wave&wave mean_P14=mean_P14_wave&wave mean_P15=mean_P15_wave&wave mean_P20=mean_P20_wave&wave mean_P22=mean_P22_wave&wave mean_P23=mean_P23_wave&wave mean_P24=mean_P24_wave&wave mean_P30=mean_P30_wave&wave mean_P32=mean_P32_wave&wave mean_P33=mean_P33_wave&wave mean_P42=mean_P42_wave&wave mean_P44=mean_P44_wave&wave mean_P65=mean_P65_wave&wave mean_indexC=mean_indexC_wave&wave var_indexC=var_indexC_wave&wave
/* Talent Index */
mean_P1=mean_P1_wave&wave mean_P11=mean_P11_wave&wave mean_P18=mean_P18_wave&wave mean_P21=mean_P21_wave&wave mean_P29=mean_P29_wave&wave mean_P47=mean_P47_wave&wave mean_P68=mean_P68_wave&wave mean_indexD=mean_indexD_wave&wave var_indexD=var_indexD_wave&wave
))
%end;
;
drop _NAME_;
run;
%mend;

%mergeit();



* merge in the Wald chi-square statistics;
data Wald;
length condition $4;
  set M1_Wald_stats_TSL;
condition='TSL';
* eliminate possibility of negative test statistics;
if Wald < 0 then Wald=ranuni(80082347)*10; 
wave=wave2;
drop wave1 wave2;
run; 
proc sort data=Wald; by condition agency index; run;
proc transpose data=Wald out=Wald (drop=_NAME_ rename=(_2-_10 = Wald_wave2-Wald_wave10));
  by  condition agency index;
  ID  wave;
  var Wald;
run;
* merge all stats and create a few variables related to the stopping wave;
data summary_study2_M1;
  merge junk Wald;
  by condition agency;
if index='A' then do;
  if 1 - PROBCHI(Wald_wave2,6) > .05 then wave_stop=2;
    else if 1 - PROBCHI(Wald_wave3,6)  > .05  then wave_stop=3;
    else if 1 - PROBCHI(Wald_wave4,6)  > .05  then wave_stop=4;
    else if 1 - PROBCHI(Wald_wave5,6)  > .05  then wave_stop=5;
    else if 1 - PROBCHI(Wald_wave6,6)  > .05  then wave_stop=6;
    else if 1 - PROBCHI(Wald_wave7,6)  > .05  then wave_stop=7;
    else if 1 - PROBCHI(Wald_wave8,6)  > .05  then wave_stop=8;
    else if 1 - PROBCHI(Wald_wave9,6)  > .05  then wave_stop=9;
    else wave_stop=10;
  end;
else if index='B' then do;
  if 1 - PROBCHI(Wald_wave2,11) > .05 then wave_stop=2;
    else if 1 - PROBCHI(Wald_wave3,11)  > .05  then wave_stop=3;
    else if 1 - PROBCHI(Wald_wave4,11)  > .05  then wave_stop=4;
    else if 1 - PROBCHI(Wald_wave5,11)  > .05  then wave_stop=5;
    else if 1 - PROBCHI(Wald_wave6,11)  > .05  then wave_stop=6;
    else if 1 - PROBCHI(Wald_wave7,11)  > .05  then wave_stop=7;
    else if 1 - PROBCHI(Wald_wave8,11)  > .05  then wave_stop=8;
    else if 1 - PROBCHI(Wald_wave9,11)  > .05  then wave_stop=9;
    else wave_stop=10;
  end;
else if index='C' then do;
  if 1 - PROBCHI(Wald_wave2,12) > .05 then wave_stop=2;
    else if 1 - PROBCHI(Wald_wave3,12)  > .05  then wave_stop=3;
    else if 1 - PROBCHI(Wald_wave4,12)  > .05  then wave_stop=4;
    else if 1 - PROBCHI(Wald_wave5,12)  > .05  then wave_stop=5;
    else if 1 - PROBCHI(Wald_wave6,12)  > .05  then wave_stop=6;
    else if 1 - PROBCHI(Wald_wave7,12)  > .05  then wave_stop=7;
    else if 1 - PROBCHI(Wald_wave8,12)  > .05  then wave_stop=8;
    else if 1 - PROBCHI(Wald_wave9,12)  > .05  then wave_stop=9;
    else wave_stop=10;
  end;
else if index='D' then do;
  if 1 - PROBCHI(Wald_wave2,6) > .05 then wave_stop=2;
    else if 1 - PROBCHI(Wald_wave3,6)  > .05  then wave_stop=3;
    else if 1 - PROBCHI(Wald_wave4,6)  > .05  then wave_stop=4;
    else if 1 - PROBCHI(Wald_wave5,6)  > .05  then wave_stop=5;
    else if 1 - PROBCHI(Wald_wave6,6)  > .05  then wave_stop=6;
    else if 1 - PROBCHI(Wald_wave7,6)  > .05  then wave_stop=7;
    else if 1 - PROBCHI(Wald_wave8,6)  > .05  then wave_stop=8;
    else if 1 - PROBCHI(Wald_wave9,6)  > .05  then wave_stop=9;
    else wave_stop=10;
  end;
run;
data summary_study2_M1;
  set summary_study2_M1;
* A) JS Index;
array indicesA    {10} mean_indexA_wave1 - mean_indexA_wave10;
array varindicesA {10} var_indexA_wave1  - var_indexA_wave10;
bias_indexA=indicesA{wave_stop} - indicesA{10};

array P4s  {*} mean_P4_wave1  - mean_P4_wave10;
array P5s  {*} mean_P5_wave1  - mean_P5_wave10;
array P13s {*} mean_P13_wave1 - mean_P13_wave10;
array P63s {*} mean_P63_wave1 - mean_P63_wave10;
array P67s {*} mean_P67_wave1 - mean_P67_wave10;
array P69s {*} mean_P69_wave1 - mean_P69_wave10;
array P70s {*} mean_P70_wave1 - mean_P70_wave10;
bias_P4 = P4s{wave_stop} - P4s{10};
bias_P5 = P5s{wave_stop} - P5s{10};
bias_P13=P13s{wave_stop} - P13s{10};
bias_P63=P63s{wave_stop} - P63s{10};
bias_P67=P67s{wave_stop} - P67s{10};
bias_P69=P69s{wave_stop} - P69s{10};
bias_P70=P70s{wave_stop} - P70s{10};
* create a median bias measure;
bias_mean_PsA=mean(of bias_P4 -- bias_P70);
* create a 95% CI indicator of the index at the stopping wave;
CI_95_indexA=(indicesA{wave_stop}-1.96*sqrt(varindicesA{wave_stop}) le indicesA{10} le indicesA{wave_stop}+1.96*sqrt(varindicesA{wave_stop}));
* create an MSE and RMSE or the index at the stopping wave;
MSE_indexA =varindicesA{wave_stop} + bias_indexA**2;
RMSE_indexA=sqrt(MSE_indexA);

* B) Leadership Index;
array indicesB    {10} mean_indexB_wave1 - mean_indexB_wave10;
array varindicesB {10} var_indexB_wave1  - var_indexB_wave10;
bias_indexB=indicesB{wave_stop} - indicesB{10};

array P10s {*} mean_P10_wave1 - mean_P10_wave10;
array P35s {*} mean_P35_wave1 - mean_P35_wave10;
array P36s {*} mean_P36_wave1 - mean_P36_wave10;
array P51s {*} mean_P51_wave1 - mean_P51_wave10;
array P52s {*} mean_P52_wave1 - mean_P52_wave10;
array P53s {*} mean_P53_wave1 - mean_P53_wave10;
array P55s {*} mean_P55_wave1 - mean_P55_wave10;
array P56s {*} mean_P56_wave1 - mean_P56_wave10;
array P57s {*} mean_P57_wave1 - mean_P57_wave10;
array P61s {*} mean_P61_wave1 - mean_P61_wave10;
array P64s {*} mean_P64_wave1 - mean_P64_wave10;
array P66s {*} mean_P66_wave1 - mean_P66_wave10;
bias_P10=P10s{wave_stop} - P10s{10};
bias_P35=P35s{wave_stop} - P35s{10};
bias_P36=P36s{wave_stop} - P36s{10};
bias_P51=P51s{wave_stop} - P51s{10};
bias_P52=P52s{wave_stop} - P52s{10};
bias_P53=P53s{wave_stop} - P53s{10};
bias_P55=P55s{wave_stop} - P55s{10};
bias_P56=P56s{wave_stop} - P56s{10};
bias_P57=P57s{wave_stop} - P57s{10};
bias_P61=P61s{wave_stop} - P61s{10};
bias_P64=P64s{wave_stop} - P64s{10};
bias_P66=P66s{wave_stop} - P66s{10};
* create a median bias measure;
bias_mean_PsB=mean(of bias_P10 -- bias_P66);
* create a 95% CI indicator of the index at the stopping wave;
CI_95_indexB=(indicesB{wave_stop}-1.96*sqrt(varindicesB{wave_stop}) le indicesB{10} le indicesB{wave_stop}+1.96*sqrt(varindicesB{wave_stop}));
* create an MSE and RMSE or the index at the stopping wave;
MSE_indexB =varindicesB{wave_stop} + bias_indexB**2;
RMSE_indexB=sqrt(MSE_indexB);

* C) ROPC Index;
array indicesC    {10} mean_indexC_wave1 - mean_indexC_wave10;
array varindicesC {10} var_indexC_wave1  - var_indexC_wave10;
bias_indexC=indicesC{wave_stop} - indicesC{10};

array P12s {*} mean_P12_wave1 - mean_P12_wave10;
array P14s {*} mean_P14_wave1 - mean_P14_wave10;
array P15s {*} mean_P15_wave1 - mean_P15_wave10;
array P20s {*} mean_P20_wave1 - mean_P20_wave10;
array P22s {*} mean_P22_wave1 - mean_P22_wave10;
array P23s {*} mean_P23_wave1 - mean_P23_wave10;
array P24s {*} mean_P24_wave1 - mean_P24_wave10;
array P30s {*} mean_P30_wave1 - mean_P30_wave10;
array P32s {*} mean_P32_wave1 - mean_P32_wave10;
array P33s {*} mean_P33_wave1 - mean_P33_wave10;
array P42s {*} mean_P42_wave1 - mean_P42_wave10;
array P44s {*} mean_P44_wave1 - mean_P44_wave10;
array P65s {*} mean_P65_wave1 - mean_P65_wave10;
bias_P12=P12s{wave_stop} - P12s{10};
bias_P14=P14s{wave_stop} - P14s{10};
bias_P15=P15s{wave_stop} - P15s{10};
bias_P20=P20s{wave_stop} - P20s{10};
bias_P22=P22s{wave_stop} - P22s{10};
bias_P23=P23s{wave_stop} - P23s{10};
bias_P24=P24s{wave_stop} - P24s{10};
bias_P30=P30s{wave_stop} - P30s{10};
bias_P32=P32s{wave_stop} - P32s{10};
bias_P33=P33s{wave_stop} - P33s{10};
bias_P42=P42s{wave_stop} - P42s{10};
bias_P44=P44s{wave_stop} - P44s{10};
bias_P65=P65s{wave_stop} - P65s{10};

* create a median bias measure;
bias_mean_PsC=mean(of bias_P10 -- bias_P65);
* create a 95% CI indicator of the index at the stopping wave;
CI_95_indexC=(indicesC{wave_stop}-1.96*sqrt(varindicesC{wave_stop}) le indicesC{10} le indicesC{wave_stop}+1.96*sqrt(varindicesC{wave_stop}));
* create an MSE and RMSE or the index at the stopping wave;
MSE_indexC =varindicesC{wave_stop} + bias_indexC**2;
RMSE_indexC=sqrt(MSE_indexC);

* D) Talent Index;
array indicesD    {10} mean_indexD_wave1 - mean_indexD_wave10;
array varindicesD {10} var_indexD_wave1  - var_indexD_wave10;
bias_indexD=indicesD{wave_stop} - indicesD{10};

array P1s  {*} mean_P1_wave1  - mean_P1_wave10;
array P11s {*} mean_P11_wave1 - mean_P11_wave10;
array P18s {*} mean_P18_wave1 - mean_P18_wave10;
array P21s {*} mean_P21_wave1 - mean_P21_wave10;
array P29s {*} mean_P29_wave1 - mean_P29_wave10;
array P47s {*} mean_P47_wave1 - mean_P47_wave10;
array P68s {*} mean_P68_wave1 - mean_P68_wave10;
bias_P1 = P1s{wave_stop} - P1s{10};
bias_P11=P11s{wave_stop} - P11s{10};
bias_P18=P18s{wave_stop} - P18s{10};
bias_P21=P21s{wave_stop} - P21s{10};
bias_P29=P29s{wave_stop} - P29s{10};
bias_P47=P47s{wave_stop} - P47s{10};
bias_P68=P68s{wave_stop} - P68s{10};
* create a median bias measure;
bias_mean_PsD=mean(of bias_P1 -- bias_P68);
* create a 95% CI indicator of the index at the stopping wave;
CI_95_indexD=(indicesD{wave_stop}-1.96*sqrt(varindicesD{wave_stop}) le indicesD{10} le indicesD{wave_stop}+1.96*sqrt(varindicesD{wave_stop}));
* create an MSE and RMSE or the index at the stopping wave;
MSE_indexD =varindicesA{wave_stop} + bias_indexD**2;
RMSE_indexD=sqrt(MSE_indexD);
run;

data forreport;
  set summary_study2_M1;
array indicesA    {10} mean_indexA_wave1 - mean_indexA_wave10;
array indicesB    {10} mean_indexB_wave1 - mean_indexB_wave10;
array indicesC    {10} mean_indexC_wave1 - mean_indexC_wave10;
array indicesD    {10} mean_indexD_wave1 - mean_indexD_wave10;
if index='A' then do;
   estimate=100*indicesA{wave_stop};
   NRerror=100*bias_indexA;
   end;
if index='B' then do;
   estimate=100*indicesB{wave_stop};
   NRerror=100*bias_indexB;
   end;
if index='C' then do;
   estimate=100*indicesC{wave_stop};
   NRerror=100*bias_indexC;
   end;
if index='D' then do;
   estimate=100*indicesD{wave_stop};
   NRerror=100*bias_indexD;
   end;
run;

ods csv;
proc print data=forreport noobs;
  var agency index wave_stop estimate NRerror;
run;
ods csv close;




*** NZT method;
* first step is to create a summarized data set with weights and P variables by agency;
%macro getstats();
%do wave=1 %to 10;
proc means data=sampnpop_2011 noprint nway;
  class agency;
  var P1 -- P70;
weight weight_wave&wave._raked;
output out=stats_wave&wave (drop=_:) mean=;
run;
data stats_wave&wave;
  wave=&wave;
  set stats_wave&wave;
run;
%end;
data stats_wave_all;
  set stats_wave1 - stats_wave10;
run;
proc sort data=stats_wave_all;
  by agency wave; 
run;
data stats_wave_all;
  set stats_wave_all;
array pees{*} P1 - P70;
array rpcs{*} rpc_P1 - rpc_P70;
do i=1 to dim(rpcs);
   rpcs{i}=100*(pees{i} - lag(pees{i}))/lag(pees{i});
   if wave=1 then rpcs{i}=.;
end;
drop i;
run;
* housecleaning;
proc datasets nolist;
  delete stats_wave1 - stats_wave10;
run;
quit; 
%mend;


%getstats();


* macro to conduct wave-over-wave F tests by agency;
%macro M2_getFtests();
%do wave_max=4 %to 10;
%let wave_min=%eval(&wave_max-2);
* A) JS Index;
data toregA;
  set stats_wave_all (in=_4  where=(&wave_min le wave le &wave_max) rename=(rpc_P4 =rpc))
      stats_wave_all (in=_5  where=(&wave_min le wave le &wave_max) rename=(rpc_P5 =rpc))
	  stats_wave_all (in=_13 where=(&wave_min le wave le &wave_max) rename=(rpc_P13=rpc))
	  stats_wave_all (in=_63 where=(&wave_min le wave le &wave_max) rename=(rpc_P63=rpc))
	  stats_wave_all (in=_67 where=(&wave_min le wave le &wave_max) rename=(rpc_P67=rpc))
	  stats_wave_all (in=_69 where=(&wave_min le wave le &wave_max) rename=(rpc_P69=rpc))
	  stats_wave_all (in=_70 where=(&wave_min le wave le &wave_max) rename=(rpc_P70=rpc));

* create the design matrix variables by hand;
wave=wave-&wave_min; 

* 1) intercepts;
int_P4=_4;
int_P5=_5;
int_P13=_13;
int_P63=_63;
int_P67=_67;
int_P69=_69;
int_P70=_70;

* 2) slopes;
slope_P4=wave*_4;
slope_P5=wave*_5;
slope_P13=wave*_13;
slope_P63=wave*_63;
slope_P67=wave*_67;
slope_P69=wave*_69;
slope_P70=wave*_70;

keep agency wave rpc int_: slope_:;
run;
proc sort data=toregA;
  by agency wave; 
run;
ods listing close;
ods output ANOVA=testsA_wave&wave_max;
proc reg data=toregA;
  by agency;
  model rpc = int_: slope_: / noint;
run;
quit;
ods listing;
* further process the ANOVA output;
data testsA_wave&wave_max;
length index $1;
index='A';
wave=&wave_max;
  set testsA_wave&wave_max (keep=agency FValue ProbF);
where not missing(FValue);
run;
* B) LKM Index;
data toregB;
  set stats_wave_all (in=_10 where=(&wave_min le wave le &wave_max) rename=(rpc_P10=rpc))
      stats_wave_all (in=_35 where=(&wave_min le wave le &wave_max) rename=(rpc_P35=rpc))
	  stats_wave_all (in=_36 where=(&wave_min le wave le &wave_max) rename=(rpc_P36=rpc))
	  stats_wave_all (in=_51 where=(&wave_min le wave le &wave_max) rename=(rpc_P51=rpc))
	  stats_wave_all (in=_52 where=(&wave_min le wave le &wave_max) rename=(rpc_P52=rpc))
	  stats_wave_all (in=_53 where=(&wave_min le wave le &wave_max) rename=(rpc_P53=rpc))
	  stats_wave_all (in=_55 where=(&wave_min le wave le &wave_max) rename=(rpc_P55=rpc))
	  stats_wave_all (in=_56 where=(&wave_min le wave le &wave_max) rename=(rpc_P56=rpc))
	  stats_wave_all (in=_57 where=(&wave_min le wave le &wave_max) rename=(rpc_P57=rpc))
	  stats_wave_all (in=_61 where=(&wave_min le wave le &wave_max) rename=(rpc_P61=rpc))
	  stats_wave_all (in=_64 where=(&wave_min le wave le &wave_max) rename=(rpc_P64=rpc))
	  stats_wave_all (in=_66 where=(&wave_min le wave le &wave_max) rename=(rpc_P66=rpc));

* create the design matrix variables by hand;
wave=wave-&wave_min; 

* 1) intercepts;
int_P10=_10;
int_P35=_35;
int_P36=_36;
int_P51=_51;
int_P52=_52;
int_P53=_53;
int_P55=_55;
int_P56=_56;
int_P57=_57;
int_P61=_61;
int_P64=_64;
int_P66=_66;

* 2) slopes;
slope_P10=wave*_10;
slope_P35=wave*_35;
slope_P36=wave*_36;
slope_P51=wave*_51;
slope_P52=wave*_52;
slope_P53=wave*_53;
slope_P55=wave*_55;
slope_P56=wave*_56;
slope_P57=wave*_57;
slope_P61=wave*_61;
slope_P64=wave*_64;
slope_P66=wave*_66;

keep agency wave rpc int_: slope_:;
run;
proc sort data=toregB;
  by agency wave; 
run;
ods listing close;
ods output ANOVA=testsB_wave&wave_max;
proc reg data=toregB;
  by agency;
  model rpc = int_: slope_: / noint;
run;
quit;
ods listing;
* further process the ANOVA output;
data testsB_wave&wave_max;
length index $1;
index='B';
wave=&wave_max;
  set testsB_wave&wave_max (keep=agency FValue ProbF);
where not missing(FValue);
run;

* C) ROPC Index;
data toregC;
  set stats_wave_all (in=_12 where=(&wave_min le wave le &wave_max) rename=(rpc_P12=rpc))
      stats_wave_all (in=_14 where=(&wave_min le wave le &wave_max) rename=(rpc_P14=rpc))
	  stats_wave_all (in=_15 where=(&wave_min le wave le &wave_max) rename=(rpc_P15=rpc))
	  stats_wave_all (in=_20 where=(&wave_min le wave le &wave_max) rename=(rpc_P20=rpc))
	  stats_wave_all (in=_22 where=(&wave_min le wave le &wave_max) rename=(rpc_P22=rpc))
	  stats_wave_all (in=_23 where=(&wave_min le wave le &wave_max) rename=(rpc_P23=rpc))
	  stats_wave_all (in=_24 where=(&wave_min le wave le &wave_max) rename=(rpc_P24=rpc))
	  stats_wave_all (in=_30 where=(&wave_min le wave le &wave_max) rename=(rpc_P30=rpc))
	  stats_wave_all (in=_32 where=(&wave_min le wave le &wave_max) rename=(rpc_P32=rpc))
	  stats_wave_all (in=_33 where=(&wave_min le wave le &wave_max) rename=(rpc_P33=rpc))
	  stats_wave_all (in=_42 where=(&wave_min le wave le &wave_max) rename=(rpc_P42=rpc))
	  stats_wave_all (in=_44 where=(&wave_min le wave le &wave_max) rename=(rpc_P44=rpc))
	  stats_wave_all (in=_65 where=(&wave_min le wave le &wave_max) rename=(rpc_P65=rpc));

* create the design matrix variables by hand;
wave=wave-&wave_min; 

* 1) intercepts;
int_P12=_12;
int_P14=_14;
int_P15=_15;
int_P20=_20;
int_P22=_22;
int_P23=_23;
int_P24=_24;
int_P30=_30;
int_P32=_32;
int_P33=_33;
int_P42=_42;
int_P44=_44;
int_P65=_65;

* 2) slopes;
slope_P12=wave*_12;
slope_P14=wave*_14;
slope_P15=wave*_15;
slope_P20=wave*_20;
slope_P22=wave*_22;
slope_P23=wave*_23;
slope_P24=wave*_24;
slope_P30=wave*_30;
slope_P32=wave*_32;
slope_P33=wave*_33;
slope_P42=wave*_42;
slope_P44=wave*_44;
slope_P65=wave*_65;

keep agency wave rpc int_: slope_:;
run;
proc sort data=toregC;
  by agency wave; 
run;
ods listing close;
ods output ANOVA=testsC_wave&wave_max;
proc reg data=toregC;
  by agency;
  model rpc = int_: slope_: / noint;
run;
quit;
ods listing;
* further process the ANOVA output;
data testsC_wave&wave_max;
length index $1;
index='C';
wave=&wave_max;
  set testsC_wave&wave_max (keep=agency FValue ProbF);
where not missing(FValue);
run;

* D) Talent Index;
data toregD;
  set stats_wave_all (in=_1  where=(&wave_min le wave le &wave_max) rename=(rpc_P1 =rpc))
      stats_wave_all (in=_11 where=(&wave_min le wave le &wave_max) rename=(rpc_P11=rpc))
	  stats_wave_all (in=_18 where=(&wave_min le wave le &wave_max) rename=(rpc_P18=rpc))
	  stats_wave_all (in=_21 where=(&wave_min le wave le &wave_max) rename=(rpc_P21=rpc))
	  stats_wave_all (in=_29 where=(&wave_min le wave le &wave_max) rename=(rpc_P29=rpc))
	  stats_wave_all (in=_47 where=(&wave_min le wave le &wave_max) rename=(rpc_P47=rpc))
	  stats_wave_all (in=_68 where=(&wave_min le wave le &wave_max) rename=(rpc_P68=rpc));

* create the design matrix variables by hand;
wave=wave-&wave_min; 

* 1) intercepts;
int_P1 =_1;
int_P11=_11;
int_P18=_18;
int_P21=_21;
int_P29=_29;
int_P47=_47;
int_P68=_68;

* 2) slopes;
slope_P1 =wave*_1;
slope_P11=wave*_11;
slope_P18=wave*_18;
slope_P21=wave*_21;
slope_P29=wave*_29;
slope_P47=wave*_47;
slope_P68=wave*_68;

keep agency wave rpc int_: slope_:;
run;
proc sort data=toregD;
  by agency wave; 
run;
ods listing close;
ods output ANOVA=testsD_wave&wave_max;
proc reg data=toregD;
  by agency;
  model rpc = int_: slope_: / noint;
run;
quit;
ods listing;
* further process the ANOVA output;
data testsD_wave&wave_max;
length index $1;
index='D';
wave=&wave_max;
  set testsD_wave&wave_max (keep=agency FValue ProbF);
where not missing(FValue);
run;
%end; * closes the wave loop;
data M2_Ftests_all;
  set testsA_wave4 - testsA_wave10
      testsB_wave4 - testsB_wave10
      testsC_wave4 - testsC_wave10
      testsD_wave4 - testsD_wave10;
run;
* housecleaning;
proc datasets nolist;
  delete toregA toregB toregC toregD
         testsA_wave4 - testsA_wave10
         testsB_wave4 - testsB_wave10
         testsC_wave4 - testsC_wave10
         testsD_wave4 - testsD_wave10;
run;
quit;
%mend;


%M2_getFtests();



*** compile a summary data set;
data part1;
  set M2_Ftests_all;
drop FValue;
format _ALL_;
attrib _ALL_ label=' ';
run;
proc sort data=part1;
  by agency index wave;
run;
proc transpose data=part1 out=part1 (drop=_NAME_ rename=(_4 - _10 = ProbF_wave4 - ProbF_wave10));
  by agency index;
  ID wave;
  var ProbF;
run;
data part2;
  set stats_wave_all;
if wave ge 4;
* create indices A - D here;
indexA=mean(of P4,P5,P13,P63,P67,P69,P70);
indexB=mean(of P10,P35,P36,P51,P52,P53,P55,P56,P57,P61,P64,P66);
indexC=mean(of P12,P14,P15,P20,P22,P23,P24,P30,P32,P33,P42,P44,P65);
indexD=mean(of P1,P11,P18,P21,P29,P47,P68);
keep agency wave index:;
run;
proc sort data=part2;
  by agency;
run;
proc transpose data=part2 out=part2 (rename=(_4 - _10 = index_wave4 - index_wave10));
  by agency;
  ID wave;
  var indexA -- indexD;
run;
data part2;
length index $1;
  set part2;
index=substr(_NAME_,6,1);
format _ALL_;
attrib _ALL_ label=' ';
drop _NAME_;
run;
proc sort data=part1; by agency index; run;
proc sort data=part2; by agency index; run;
data summary_study2_M2;
  merge part1 part2;
  by agency index;
*** create some indicators of performance and error;
* stopping wave;
if ProbF_wave4 > .05 then wave_stop=4;
  else if ProbF_wave5  > .05 then wave_stop=5;
  else if ProbF_wave6  > .05 then wave_stop=6;
  else if ProbF_wave7  > .05 then wave_stop=7;
  else if ProbF_wave8  > .05 then wave_stop=8;
  else if ProbF_wave9  > .05 then wave_stop=9;
  else wave_stop=10;

array indices {7} index_wave4 - index_wave10;
estimate=100*indices{wave_stop-3};
NRerror =100*(indices{wave_stop-3} - indices{7});
run;

* output these stats for the applicable table (alongside M1 stats);
ods csv;
proc print data=summary_study2_M2 noobs;
run;
ods csv close;

****** END: application **********************************************************************************************;





*** code to reproduce Figure 2;
data stats;
  set summary_study2_M1 ;
  by agency;
  if first.agency;
keep agency mean_ind:;
run;
data stats;
  wave=0;
  set 
stats (in=_1  rename=(mean_indexA_wave1=indexA  mean_indexB_wave1=indexB  mean_indexC_wave1=indexC  mean_indexD_wave1=indexD))
stats (in=_2  rename=(mean_indexA_wave2=indexA  mean_indexB_wave2=indexB  mean_indexC_wave2=indexC  mean_indexD_wave2=indexD))
stats (in=_3  rename=(mean_indexA_wave3=indexA  mean_indexB_wave3=indexB  mean_indexC_wave3=indexC  mean_indexD_wave3=indexD))
stats (in=_4  rename=(mean_indexA_wave4=indexA  mean_indexB_wave4=indexB  mean_indexC_wave4=indexC  mean_indexD_wave4=indexD))
stats (in=_5  rename=(mean_indexA_wave5=indexA  mean_indexB_wave5=indexB  mean_indexC_wave5=indexC  mean_indexD_wave5=indexD))
stats (in=_6  rename=(mean_indexA_wave6=indexA  mean_indexB_wave6=indexB  mean_indexC_wave6=indexC  mean_indexD_wave6=indexD))
stats (in=_7  rename=(mean_indexA_wave7=indexA  mean_indexB_wave7=indexB  mean_indexC_wave7=indexC  mean_indexD_wave7=indexD))
stats (in=_8  rename=(mean_indexA_wave8=indexA  mean_indexB_wave8=indexB  mean_indexC_wave8=indexC  mean_indexD_wave8=indexD))
stats (in=_9  rename=(mean_indexA_wave9=indexA  mean_indexB_wave9=indexB  mean_indexC_wave9=indexC  mean_indexD_wave9=indexD))
stats (in=_10 rename=(mean_indexA_wave10=indexA mean_indexB_wave10=indexB mean_indexC_wave10=indexC mean_indexD_wave10=indexD))
;
indexA=indexA*100;indexB=indexB*100;indexC=indexC*100;indexD=indexD*100;
if _1 then wave=1;
  else if _2 then wave=2;
  else if _3 then wave=3;
  else if _4 then wave=4;
  else if _5 then wave=5;
  else if _6 then wave=6;
  else if _7 then wave=7;
  else if _8 then wave=8;
  else if _9 then wave=9;
  else if _10 then wave=10;
keep agency wave index:;
run;

proc sgpanel data=stats;
  where agency='AGY2' and wave le 9;
  panelby agency /;
  series Y=indexA X=wave;
  series Y=indexB X=wave;
  series Y=indexC X=wave;
  series Y=indexD X=wave;
label indexA='JS'
      indexB='LKM'
	  indexC='ROPC'
	  indexD='TM'
      wave='Wave of Response';
keylegend / title='Indices';
rowaxis label='Nonresponse-Adjusted Index Estimate'  values=(50 to 70 by 2);
colaxis label='Number of Waves Completed' values=(1 to 9 by 1);
run;