new;

/* SUBSETR - program for computing Akaike AIC, Schwarz BIC and Rissanen RIC (Sclove, 1987) statistics

for ordered subsets of means based on repeated observations. Mean are organized by repeated measures

time period in one flat flat text file and the variance-covariance matrix in a second flat test file.

Reference: C. Dayton (1998), Information Criteria for the Paired-Comparisons Problem, American

Statistician, 52, 144-151. */

/* All screen output is also directed to the file subset.out in the current directory */

print "Default file for all printed output is Subset.out in the current directory";

output file = subset.out reset;

print "Program SUBSETR for Ordered Subsets of Means or Proportions";

print "Prepared by: C. Mitchell Dayton";

print "Department of Measurement & Statistics";

print "University of Maryland";

print "E-Mail: CD4@UMAIL.UMD.EDU";

/* Determine number of patterns to display on output */

print "Number of AIC/BIC values to display (5 is recommended; for 3 groups, use 4) ";; nsave = con(1,1);

print "Sample size = ";n=con(1,1);

/* Read means */

print "Name of input file for means? ";; fmt = cons;

loadm ydat[] = ^fmt;

k = rows(ydat); /* k is number of repeated measures */

print "Name of input file for variance-covariance matrix? ";; fmt1 = cons;

loadm cov[] = ^fmt1;

cov = reshape(cov,k,k);

aa = (n-1)*cov;

a = (k*n/2)*ln(2*pi);

xx=cumsumc(ones(k,1));

base=2;

aic1=zeros(nsave,k+1); bic1=aic1; ric1=aic1;

x=1;

do until x > 2^(k-1); x1=x-1;

/* Generate v with binary equivalent of integer = x-1 */

if x1 == 0; v = zeros(1,k);

elseif x1 > 0;

nbase = maxc(log(x1)./log(base))+1;

zbase=reshape(base,nbase,rows(x1));

vv = rev(recserrc(x1,zbase))';

v = zeros(1,k-rows(vv'))~vv;

endif;

/* Compute u that contains ordered subsets based on changes in v.

For example, if v = {0010}, then u = {1123} */

u = zeros(1,k);

i=2;

u[1]=1;

do until i > k;

if v[i] == v[i-1]; u[i] = u[i-1];

elseif v[i] /= v[i-1]; u[i] = u[i-1]+1; endif;

i = i+1; endo;

/* Create k-mean vector, mms, with equality restrictions imposed based on patterns in u vector

*/

mns = zeros(1,k); nns=mns;

xb = mns;

nn = ones(1,k);

i=1;

do until i > maxc(u');

mns[i]= sumc(submat(ydat,indexcat(u',i),1));

nns[i]= sumc(submat(nn',indexcat(u',i),1));

mns[i]=mns[i]/rows(mns[i]);

i=i+1; endo;

mms = (subscat(u',xx,mns'))';

nns = (subscat(u',xx,nns'))';

xb = mms./nns;

/* Compute log-likelihoods and AIC, BIC values*/

b = (n/2)*ln(det((aa+n*(ydat-xb')*(ydat-xb')')/n));

lik = -.5*sumc(diag(inv((aa+n*(ydat-xb')*(ydat-xb')')/n)*(aa+n*(ydat-xb')*(ydat-xb')'))) -a - b;

mod=k!/((k-2)!*2)+k;

aic = -2*lik + 2*(mod+maxc(u'));

bic = -2*lik + ln(n)*(mod+maxc(u'));

ric = -2*lik + ln((n+2)/24)*(mod+maxc(u'));

/* print aic~u~v; used as diagnostic only */

/* Save AIC, BIC if smaller than any currently saved case */

if x < nsave+1; aic1[x,.] = aic~u; bic1[x,.] = bic~u; ric1[x,.] = ric~u; goto next1; endif;

if maxc(aic1[.,1]) > aic; aic1[maxindc(aic1[.,1]),.] = aic~u; endif;

if maxc(bic1[.,1]) > bic; bic1[maxindc(bic1[.,1]),.] = bic~u; endif;

if maxc(ric1[.,1]) > ric; ric1[maxindc(ric1[.,1]),.] = ric~u; endif;

next1:;

x = x+1;

endo; /* End of loop for generating mean patterns */

/* Sort saved AIC, BIC values in ascending order */

aic1 = sortc(aic1,1); bic1 = sortc(bic1,1); ric1 = sortc(ric1,1);

/* Print results for case of means */

format /mat /rd 10, 3;

ii=1;

print " "; print "Ordered means are: "; print ydat';

print " "; print "Best models";

print " "; print "Smallest Akaike AIC values and ordered subsets: ";

print aic1[ii:nsave,.];

print " "; print "Smallest Schwarz BIC values and ordered subsets: ";

print bic1[ii:nsave,.];

print " "; print "Smallest Rissanen RIC values and ordered subsets: ";

print ric1[ii:nsave,.];

format /mat /rd 15, 4;

output off;

end