clear;tic;format short;

% Specify the data to be analyzed. 
% The data matrix must be named "scores".

% load esteem2.dat; scores = esteem2;

%***************  START OF TRIAL-RUN DATA  ************************************

% The following data are from Brennan, 2001, p 73, Table 3.2, for a
% p x i x o design, with 10 persons, 4 items, and 3 occasions.
% To analyze your own data, the trial-run data matrix must be removed/deleted.

% Brennan (2001, p. 28, Table 2.2, & p. 36);
% Type 1, one-facet crossed design, as in p x i;
% 10 persons, 12 items.
scores = [
1   0   1   0   0   0   0   0   0   0   0   0
1   1   1   0   0   1   0   0   0   0   0   0
1   1   1   1   1   0   0   0   0   0   0   0
1   1   0   1   1   0   0   1   0   0   0   0
1   1   1   1   1   0   1   0   0   0   0   0
1   1   1   0   1   1   1   0   0   0   0   0
1   1   1   1   1   1   1   0   0   0   0   0
1   1   1   1   0   1   1   1   1   1   0   0
1   1   1   1   1   1   1   1   1   1   1   0
1   1   1   1   1   1   1   1   1   1   1   1 ];

% Brennan (2001, p. 43, Table 2.8, & p. 44);
% Type 2, one-facet nested design, as in i : p ;
% 10 persons, 8 items.  
scores = [
   2   6   7   5   2   5   5   5
   4   5   6   7   6   7   5   7
   5   5   4   6   5   4   5   5
   5   9   8   6   5   7   7   6
   4   3   5   6   4   5   6   4
   4   4   4   7   6   4   7   8
   2   6   6   5   2   7   7   5
   3   4   4   5   6   6   6   4
   0   5   4   5   5   5   5   3
   6   8   7   6   6   8   8   6  ];

% Brennan (2001, p. 72, Table 3.1, & p. 111;
% Type 3, two-facet crossed design, as in  p x i x o; 
% 10 persons, 4 items, 2 occasions. 
scores = [
2 6 7 5  2 5 5 5
4 5 6 7  6 7 5 7
5 5 4 6  5 4 5 5
5 9 8 6  5 7 7 6
4 3 5 6  4 5 6 4
4 4 4 7  6 4 7 8
2 6 6 5  2 7 7 5
3 4 4 5  6 6 6 4
0 5 4 5  5 5 5 3
6 8 7 6  6 8 8 6 ];

% Brennan (2001, p. 73, Table 3.2, & p. 116;
% Type 4, two-facet nested design, as in  p x (r : t); 
% 10 persons, 4 raters, 3 tasks. 
scores = [
5 6 5 5  5 3 4 5  6 7 3 3
9 3 7 7  7 5 5 5  7 7 5 2
3 4 3 3  5 3 3 5  6 5 1 6
7 5 5 3  3 1 4 3  5 3 3 5
9 2 9 7  7 7 3 7  2 7 5 3
3 4 3 5  3 3 6 3  4 5 1 2
7 3 7 7  7 5 5 7  5 5 5 4
5 8 5 7  7 5 5 4  3 2 1 1
9 9 8 8  6 6 6 5  5 8 1 1
4 4 4 3  3 5 6 5  5 7 1 1 ];


%***************  END OF TRIAL-RUN DATA  ************************************


% Enter the number of levels/conditions of Facet 1 (e.g., # of items)
nfacet1 = 4;

% Enter the number of levels/conditions of Facet 2 (e.g., # of occasions)
% You can ignore this step for single-facet designs
nfacet2 = 3;

% For two-facet designs, Facet 1 is the facet with the fastest-changing conditions
% in the columns of your data matrix.  For example, if the first 10 columns/variables
% contained the data for 10 different items measured on occasion 1, and if the 
% next 10 columns/variables contained the data for the same 10 items measured on
% occasion 2, then items would be the fastest-changing facet.  As you slide from
% one column to the next across the data matrix, it is the item levels
% that change most quickly.  You would therefore enter a value of "10" for
% NFACET1 and a value of "2" for NFACET2 on the above statements;

% Enter the design of your data on the "COMPUTE TYPE =" statement below:
%    enter "1" for a single-facet fully-crossed design, as in P * F1
%    enter "2" for a single-facet nested nested design, as in F1 : P
%    enter "3" for a two-facet fully-crossed design, as in P * F1 * F2
%    enter "4" for a two-facet nested design, as in P * (F1 : F2)
%    enter "5" for a two-facet nested design, as in (F1 : P) * F2
%    enter "6" for a two-facet nested design, as in F1 : (P * F2)
%    enter "7" for a two-facet nested design, as in (F1 * F2) : P
%    enter "8" for a two-facet nested design, as in F1 : F2 : P;
type = 4;

% Enter D-study values for Facet 1; enter the values inside square
% brackets, and place a comma between the values.
dfacet1 = [12,6,4,3,2,2];
dfacet1 = [5,10,15,20];
dfacet1 = [4,8,16];

% Enter D-study values for Facet 2; enter the values inside square
% brackets, and place a comma between the values.  You can ignore this
% step for single-facet designs.
dfacet2 = [1,2,3,4,5,6];
dfacet2 = [1,2,3,4,5];
dfacet2 = [2,4];

% At the bottom of this file are commands for plotting
% the results for the D-study values that you specified above.
% Specify the data that you would like to plot by entering the
% appropriate number on the GRAPHDAT statement:
%   enter "1" for relative error variances
%   enter "2" for absolute error variances
%   enter "3" for G-coefficients
%   enter "4" for phi coefficients
graphdat = 3;

% End of user specifications
%**************************************************************************

disp([' ']);disp([' ']);disp(['GENERALIZABILITY THEORY ANALYSES:']);disp([' ']);

if (type == 1); disp(['Design Type 1: single-facet fully-crossed design, as in P * F1']);end
if (type == 2); disp(['Design Type 2: single-facet nested nested design, as in F1 : P']);end
if (type == 3); disp(['Design Type 3: two-facet fully-crossed design, as in P * F1 * F2']);end
if (type == 4); disp(['Design Type 4: two-facet nested design, as in P * (F1 : F2)']);end
if (type == 5); disp(['Design Type 5: two-facet nested design, as in (F1 : P) * F2']);end
if (type == 6); disp(['Design Type 6: two-facet nested design, as in F1 : (P * F2)']);end
if (type == 7); disp(['Design Type 7: two-facet nested design, as in (F1 * F2) : P']);end
if (type == 8); disp(['Design Type 8: two-facet nested design, as in F1 : F2 : P']);end
disp([' ']);

p = size(scores,1);
i = nfacet1;
h = nfacet2;

disp(['Number of persons/objects ("P") = ' num2str(p)]);disp([' ']);
disp(['Number of levels for Facet 1 ("F1") = ' num2str(i)]);disp([' ']);
if type > 2; disp(['Number of levels for Facet 2 ("F2") = ' num2str(h)]);end;disp([' ']);


% single-facet fully-crossed design, as in P * F1
if type == 1;

% degrees of freedom
dfp = p-1;
dfi = i-1;
dfpi = dfp*dfi;

% raw data computations needed for later sum of squares calculations
hABCD = scores;
pmean = sum(hABCD,2)/i;
p2mean = sum((pmean.^2),1);
imean = sum(hABCD,1)/p;
i2mean = sum((imean.^2),2);
pi2tot = sum(sum(hABCD.^2));
gmean = sum(imean,2)/i;

Tpers = i*p2mean;
Titem = p*i2mean;
Tpi   = pi2tot;
Tu    = p*i*(gmean*gmean);

% sum of squares
pss  = Tpers - Tu;
iss  = Titem - Tu;
piss = Tpi - Tpers - Titem + Tu;

% mean squares
pms = pss/dfp;
ims = iss/dfi;
pims = piss/dfpi;

% variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = (inv([ i 0 1; 0 p 1; 0 0 1]) *  [ pms; ims; pims ])';
% setting negative variances to zero
for luper = 1:size(vcs,2); if vcs(1,luper) < 0; 
vcs(1,luper) = 0;
disp(['One or more negative variance estimates have been set to zero']);disp([' ']);
end;end
pvar  = vcs(1,1);
ivar  = vcs(1,2);
pivar = vcs(1,3);
totvar = sum(vcs);

% relative and absolute error variances
relv = pivar/i;
absv = (ivar/i) + relv;

% matrices for D-studies
persons  = [ 1, 0, 1 ];
vfacet1  = [ 0, i, i ];

disp(['ANOVA Table']);disp([' ']);
disp (['Term:              df          SS          MS     Variance        %']);
disp (['P        ' sprintf('%12.3f',[dfp,pss,pms,pvar,(pvar/totvar)]) ]);
disp (['F1       ' sprintf('%12.3f',[dfi,iss,ims,ivar,(ivar/totvar)]) ]);
disp (['P*F1     ' sprintf('%12.3f',[dfpi,piss,pims,pivar,(pivar/totvar)]) ]);disp([' ']);

end;


% single-facet nested nested design, as in F1 : P
if type == 2;

% degrees of freedom
dfp = p-1;
dfip = p*(i-1);

% raw data computations needed for later sum of squares calculations
hABCD = scores;
pmean = sum(hABCD,2)/i;
p2mean = sum(pmean.^2);
pi2tot = sum(sum(hABCD.^2));
gmean = sum(sum(pmean))/p;

Tpers = i*p2mean;
Tip   = pi2tot;
Tu    = p*i*(gmean*gmean);

% sum of squares
pss = Tpers - Tu;
ipss = Tip - Tpers;

% mean squares
pms = pss/dfp;
ipms = ipss/dfip;

% variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = (inv([ i 1; 0 1 ]) *  [ pms; ipms ])';
% setting negative variances to zero
for luper = 1:size(vcs,2); if vcs(1,luper) < 0; 
vcs(1,luper) = 0;
disp(['One or more negative variance estimates have been set to zero']);disp([' ']);
end;end
pvar   = vcs(1,1);
ipvar  = vcs(1,2);
totvar = sum(vcs);

% relative and absolute error variances
relv = ipvar/i;
absv = relv;

% matrices for D-studies
persons  = [ 1, 1 ];
vfacet1  = [ 0, i ];

disp(['ANOVA Table']);disp([' ']);
disp (['Term:              df          SS          MS     Variance        %']);
disp (['P        ' sprintf('%12.3f',[dfp,pss,pms,pvar,(pvar/totvar)]) ]);
disp (['F1:P     ' sprintf('%12.3f',[dfip,ipss,ipms,ipvar,(ipvar/totvar)]) ]);disp([' ']);

end;


% two-facet fully-crossed design, as in P * F1 * F2
if type == 3;

% degrees of freedom
dfp = (size(scores,1) - 1);
dfo = (h - 1);
dfi = ((size(scores,2)/(h)) - 1);
dfpo = (dfp*dfo);
dfpi = (dfp*dfi);
dfio = (dfi*dfo);
dfpoi = (dfp*dfo*dfi);

% pointers/counters
cntfirst = 1;
cntlast = i;

% raw data computations needed for later sum of squares calculations
pimean = repmat(0,p,i);
pmean = repmat(0,p,1);
omean = repmat(-9999,h,1);
pomean = repmat(-9999,p,h);
iomean = repmat(-9999,h,i);
for cols = 1:h;
hABCD = scores(:, cntfirst:cntlast);
iomean(cols,:) = sum(hABCD,1)/p;
pomean(:,cols) = sum(hABCD,2)/i;
pimean = pimean+hABCD;
cntfirst = cntfirst + i;
cntlast = cntlast + i;
end;

pmean = sum(pomean,2)/h;
omean = sum(pomean,1)/p;
pimean = pimean/h;
imean = sum(pimean,1)/p;
gmean = sum(pmean,1)/p;

Tpers = (i*h)*sum(pmean.^2);
Titem = (p*h)*sum((imean.^2),2);
Tocca = (p*i)*sum(omean.^2);
Tu  = p*(i)*(h)*((gmean)*(gmean));
Tpi = h*sum(sum(pimean.^2));
Tpo = sum((pomean.^2),1);
Tpo = i*sum(Tpo,2);
Tio = p*sum(sum(iomean.^2));
Tpoi = sum(sum(scores.^2));

% sum of squares
pss = (Tpers-Tu);
iss = (Titem-Tu);
oss = (Tocca-Tu);
piss = (Tpi-Tpers-Titem+Tu);
poss = (Tpo-Tpers-Tocca+Tu);
ioss = (Tio-Titem-Tocca+Tu);
poiss = (Tpoi-Tpi-Tpo-Tio+Tpers+Titem+Tocca-Tu);

% mean squares
poims = poiss/dfpoi;
ioms = ioss/dfio;
poms = poss/dfpo;
pims = piss/dfpi;
oms = oss/dfo;
ims = iss/dfi;
pms = pss/dfp;

% variance components: Brennan, 2001, p 79, 83, 85, 442
bigP = ([ i*h,0,0,h,i,0,1; 0,p*h,0,h,0,p,1; 0,0,p*i,0,i,p,1; 0,0,0,h,0,0,1;
          0,0,0,0,i,0,1;   0,0,0,0,0,p,1;   0,0,0,0,0,0,1  ]);
vcs = (inv(bigP) * [ pms; ims; oms; pims; poms; ioms; poims ])';
% setting negative variances to zero
for luper = 1:size(vcs,2); if vcs(1,luper) < 0; 
vcs(1,luper) = 0;
disp(['One or more negative variance estimates have been set to zero']);disp([' ']);
end;end
pvar  = vcs(1,1);
ivar  = vcs(1,2);
ovar  = vcs(1,3);
pivar = vcs(1,4);
povar = vcs(1,5);
iovar = vcs(1,6);
poivar = vcs(1,7);
totvar = sum(vcs);

% relative and absolute error variances
relv = (pivar/i)+(povar/h)+(poivar/(i*h));
absv = (ivar/i)+(ovar/h)+(iovar/(i*h))+relv;

% matrices for D-studies
persons  = [ 1, 0, 0, 1, 1, 0, 1 ];
vfacet1  = [ 0, i, 0, i, 0, i, i ];
vfacet2  = [ 0, 0, h, 0, h, h, h ];

disp(['ANOVA Table']);disp([' ']);
disp (['Term:              df          SS          MS     Variance        %']);
disp (['P        ' sprintf('%12.3f',[dfp,pss,pms,pvar,(pvar/totvar)]) ]);
disp (['F1       ' sprintf('%12.3f',[dfi,iss,ims,ivar,(ivar/totvar)]) ]);
disp (['F2       ' sprintf('%12.3f',[dfo,oss,oms,ovar,(ovar/totvar)]) ]);
disp (['P*F1     ' sprintf('%12.3f',[dfpi,piss,pims,pivar,(pivar/totvar)]) ]);
disp (['P*F2     ' sprintf('%12.3f',[dfpo,poss,poms,povar,(povar/totvar)]) ]);
disp (['F1*F2    ' sprintf('%12.3f',[dfio,ioss,ioms,iovar,(iovar/totvar)]) ]);
disp (['P*F1*F2  ' sprintf('%12.3f',[dfpoi,poiss,poims,poivar,(poivar/totvar)]) ]);disp([' ']);

end;


% two-facet nested design, as in P * (F1 : F2)
if type == 4;

% degrees of freedom
dfp = (size(scores,1) - 1);
dfh = (h-1);
dfih = (h)*(i-1);
dfph = dfp*dfh;
dfpih = h*dfp*(i-1);
dftot = dfp+dfh+dfih+dfph+dfpih;

% pointers/counters
cntfirst = 1;
cntlast = i;

% raw data computations needed for later sum of squares calculations
pmean = repmat(0,p,1);
hmean = repmat(-9999,h,1);
phmean = repmat(-9999,p,h);
ihmean = repmat(-9999,h,i);
for cols = 1:h;
hABCD = scores(:, cntfirst:cntlast);
ihmean(cols,:) = sum(hABCD,1)/p;
phmean(:,cols) = sum(hABCD,2)/i;
cntfirst = cntfirst + i;
cntlast = cntlast + i;
end;

pmean = sum(phmean,2)/h;
hmean = sum(phmean,1)/p;
gmean = sum(pmean,1)/p;

Tp = (i*h)*(sum(pmean.^2));
Th = (p*i)*sum((hmean.^2),2);
Tih = p*sum(sum(ihmean.^2));
Tph = i*sum(sum(phmean.^2));
Tpih = sum(sum(scores.^2));
Tu = p*i*h*gmean*gmean;

% sum of squares
pss = Tp-Tu;
hss = Th-Tu;
ihss = Tih-Th;
phss = Tph-Tp-Th+Tu;
pihss = Tpih-Tph-Tih+Th;
sstot = pss+hss+ihss+phss+pihss;

% mean squares
pms = pss/dfp;
hms = hss/dfh;
ihms = ihss/dfih;
phms = phss/dfph;
pihms = pihss/dfpih;

% variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = (inv([ i*h,0,0,i,1; 0,p*i,p,i,1; 0,0,p,0,1; 0,0,0,i,1; 0,0,0,0,1 ]) *  [ pms; hms; ihms; phms; pihms ])';
% setting negative variances to zero
for luper = 1:size(vcs,2); if vcs(1,luper) < 0; 
vcs(1,luper) = 0;
disp(['One or more negative variance estimates have been set to zero']);disp([' ']);
end;end
pvar   = vcs(1,1);
hvar   = vcs(1,2);
ihvar  = vcs(1,3);
phvar  = vcs(1,4);
pihvar = vcs(1,5);
totvar = sum(vcs);

% relative and absolute error variances
relv = (pihvar/(i*h))+(phvar/h);
absv = (hvar/h)+(ihvar/(i*h)) + relv;

% matrices for D-studies
persons  = [ 1, 0, 0, 1, 1 ];
vfacet1  = [ 0, 0, i, 0, i ];
vfacet2  = [ 0, h, h, h, h ];

disp(['ANOVA Table']);disp([' ']);
disp (['Term:              df          SS          MS     Variance        %']);
disp (['P        ' sprintf('%12.3f',[dfp,pss,pms,pvar,(pvar/totvar)]) ]);
disp (['F2       ' sprintf('%12.3f',[dfh,hss,hms,hvar,(hvar/totvar)]) ]);
disp (['F1:F2    ' sprintf('%12.3f',[dfih,ihss,ihms,ihvar,(ihvar/totvar)]) ]);
disp (['P*F2     ' sprintf('%12.3f',[dfph,phss,phms,phvar,(phvar/totvar)]) ]);
disp (['P*F1:F2  ' sprintf('%12.3f',[dfpih,pihss,pihms,pihvar,(pihvar/totvar)]) ]);disp([' ']);

end;


% two-facet nested design, as in (F1 : P) * F2
if type == 5;

% degrees of freedom
dfp = (size(scores,1) -1);
dfh = h-1;
dfip = p*(i-1);
dfph = (p-1)*(h-1);
dfhip = p*(i-1)*(h-1);
dftot = dfp+dfh+dfip+dfph+dfhip;

% initializing pointers/counters
cntfirst = 1;
cntlast = i;

% raw data computations needed for later sum of squares calculations
pmean = repmat(0,p,1);
hmean = repmat(-9999,h,1);
phmean = repmat(-9999,p,h);
ipmean = repmat(0,p,i);
for cols = 1:h;
hABCD = scores(:, cntfirst:cntlast);
phmean(:,cols) = sum(hABCD,2)/i;
ipmean = ipmean+ hABCD;
cntfirst = cntfirst + i;
cntlast = cntlast + i;
end;

pmean = sum(phmean,2)/h;
hmean = sum(phmean,1)/p;
gmean = sum(pmean,1)/p;
ipmean = ipmean/h;

Tp = (i*h)*(sum(pmean.^2));
Th = (p*i)*sum((hmean.^2),2);
Tip = h*sum(sum(ipmean.^2));
Tph = i*sum(sum(phmean.^2));
Thip = sum(sum(scores.^2));
Tu = p*i*h*gmean*gmean;

% sum of squares
pss = Tp-Tu;
hss = Th-Tu;
ipss = Tip-Tp;
phss = Tph-Tp-Th+Tu;
hipss = Thip-Tph-Tip+Tp;
sstot = pss+hss+ipss+phss+hipss;

% mean squares
pms = pss/dfp;
hms = hss/dfh;
ipms = ipss/dfip;
phms = phss/dfph;
hipms = hipss/dfhip;

% variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = (inv([ i*h,0,0,i,1; 0,p*i,0,i,1; 0,0,h,0,1; 0,0,0,i,1; 0,0,0,0,1 ]) *  [ pms; hms; ipms; phms; hipms ])';
% setting negative variances to zero
for luper = 1:size(vcs,2); if vcs(1,luper) < 0; 
vcs(1,luper) = 0;
disp(['One or more negative variance estimates have been set to zero']);disp([' ']);
end;end
pvar   = vcs(1,1);
hvar   = vcs(1,2);
ipvar  = vcs(1,3);
phvar  = vcs(1,4);
hipvar = vcs(1,5);
totvar = sum(vcs);

% relative and absolute error variances
relv = (ipvar/i)+(phvar/h)+(hipvar/(i*h));
absv = (hvar/h)+relv;

% matrices for D-studies
persons  = [ 1, 0, 1, 1, 1 ];
vfacet1  = [ 0, 0, i, 0, i ];
vfacet2  = [ 0, h, 0, h, h ];

disp(['ANOVA Table']);disp([' ']);
disp (['Term:              df          SS          MS     Variance        %']);
disp (['P        ' sprintf('%12.3f',[dfp,pss,pms,pvar,(pvar/totvar)]) ]);
disp (['F2       ' sprintf('%12.3f',[dfh,hss,hms,hvar,(hvar/totvar)]) ]);
disp (['F1:P     ' sprintf('%12.3f',[dfip,ipss,ipms,ipvar,(ipvar/totvar)]) ]);
disp (['P*F2     ' sprintf('%12.3f',[dfph,phss,phms,phvar,(phvar/totvar)]) ]);
disp (['F2*F1:P  ' sprintf('%12.3f',[dfhip,hipss,hipms,hipvar,(hipvar/totvar)]) ]);disp([' ']);

end;


% two-facet nested design, as in F1 : (P * F2)
if type == 6;

% degrees of freedom
dfp = (size(scores,1) -1);
dfh = h-1;
dfph = (p-1)*(h-1);
dfiph = p*h*(i-1);
dftot = dfp+dfh+dfph+dfiph;

% pointers/counters
cntfirst = 1;
cntlast = i;

% raw data computations needed for later sum of squares calculations
pmean = repmat(0,p,1);
hmean = repmat(-9999,h,1);
phmean = repmat(-9999,p,h);
for cols = 1:h;
hABCD = scores(:, cntfirst:cntlast);
phmean(:,cols) = sum(hABCD,2)/i;
cntfirst = cntfirst + i;
cntlast = cntlast + i;
end;

pmean = sum(phmean,2)/h;
hmean = sum(phmean,1)/p;
gmean = sum(pmean,1)/p;

Tp = (i*h)*(sum(pmean.^2));
Th = (p*i)*sum((hmean.^2),2);
Tph = i*sum(sum(phmean.^2));
Tiph = sum(sum(scores.^2));
Tu = p*i*h*gmean*gmean;

% sum of squares
pss = Tp-Tu;
hss = Th-Tu;
phss = Tph-Tp-Th+Tu;
iphss = Tiph-Tph;
sstot = pss+hss+phss+iphss;

% mean squares
pms = pss/dfp;
hms = hss/dfh;
phms = phss/dfph;
iphms = iphss/dfiph;

% variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = (inv([ i*h,0,i,1; 0,p*i,i,1; 0,0,i,1; 0,0,0,1 ]) * [ pms; hms; phms; iphms ])';
% setting negative variances to zero
for luper = 1:size(vcs,2); if vcs(1,luper) < 0; 
vcs(1,luper) = 0;
disp(['One or more negative variance estimates have been set to zero']);disp([' ']);
end;end
pvar   = vcs(1,1);
hvar   = vcs(1,2);
phvar  = vcs(1,3);
iphvar = vcs(1,4);
totvar = sum(vcs);

% relative and absolute error variances & G & phi coefficients
relv = (phvar/h)+(iphvar/(i*h));
absv = (hvar/h)+relv;

% matrices for D-studies
persons  = [ 1, 0, 1, 1 ];
vfacet1  = [ 0, 0, 0, i ];
vfacet2  = [ 0, h, h, h ];

disp(['ANOVA Table']);disp([' ']);
disp (['Term:              df          SS          MS     Variance        %']);
disp (['P        ' sprintf('%12.3f',[dfp,pss,pms,pvar,(pvar/totvar)]) ]);
disp (['F2       ' sprintf('%12.3f',[dfh,hss,hms,hvar,(hvar/totvar)]) ]);
disp (['P*F2     ' sprintf('%12.3f',[dfph,phss,phms,phvar,(phvar/totvar)]) ]); 
disp (['F1:P*F2  ' sprintf('%12.3f',[dfiph,iphss,iphms,iphvar,(iphvar/totvar)]) ]);disp([' ']);
end;


% two-facet nested design, as in (F1 * F2) : P
if type == 7;

% degrees of freedom
dfp = (size(scores,1) -1);
dfip = p*(i-1);
dfhp = p*(h-1);
dfihp = p*(i-1)*(h-1);
dftot = dfp+dfip+dfhp+dfihp;

% pointers/counters
cntfirst = 1;
cntlast = i;

% raw data computations needed for later sum of squares calculations
pmean = repmat(0,p,1);
phmean = repmat(-9999,p,h);
ipmean = repmat(0,p,i);
hpmean = repmat(-9999,p,h);
for cols = 1:h;
hABCD = scores(:, cntfirst:cntlast);
phmean(:,cols) = sum(hABCD,2)/i;
ipmean = ipmean+ hABCD;
hpmean(:,cols) = sum(hABCD,2)/i;
cntfirst = cntfirst + i;
cntlast = cntlast + i;
end;

pmean = sum(hpmean,2)/h;
hmean = sum(phmean,1)/p;
gmean = sum(pmean,1)/p;
ipmean = ipmean/h;

Tp = (i*h)*(sum(pmean.^2));
Th = (p*i)*sum(sum(hmean.^2));
Tip = h*sum(sum(ipmean.^2));
Thp = i*sum(sum(hpmean.^2));
Tihp = sum(sum(scores.^2));
Tu = p*i*h*gmean*gmean;

% sum of squares
pss = Tp-Tu;
ipss = Tip-Tp;
hpss = Thp-Tp;
ihpss = Tihp-Tip-Thp+Tp;
sstot = pss+ipss+hpss+ihpss;

% mean squares
pms = pss/dfp;
ipms = ipss/dfip;
hpms = hpss/dfhp;
ihpms = ihpss/dfihp;

% variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = (inv([ i*h,h,i,1; 0,h,0,1; 0,0,i,1; 0,0,0,1 ]) *  [ pms; ipms; hpms; ihpms ])';
% setting negative variances to zero
for luper = 1:size(vcs,2); if vcs(1,luper) < 0; 
vcs(1,luper) = 0;
disp(['One or more negative variance estimates have been set to zero']);disp([' ']);
end;end
pvar   = vcs(1,1);
ipvar  = vcs(1,2);
hpvar  = vcs(1,3);
ihpvar = vcs(1,4);
totvar = sum(vcs);

% relative and absolute error variances
relv = (ipvar/i)+(hpvar/h)+(ihpvar/(i*h));
absv = relv;

% matrices for D-studies
persons  = [ 1, 1, 1, 1 ];
vfacet1  = [ 0, i, 0, i ];
vfacet2  = [ 0, 0, h, h ];

disp(['ANOVA Table']);disp([' ']);
disp (['Term:              df          SS          MS     Variance        %']);
disp (['P        ' sprintf('%12.3f',[dfp,pss,pms,pvar,(pvar/totvar)]) ]);
disp (['F1*P     ' sprintf('%12.3f',[dfip,ipss,ipms,ipvar,(ipvar/totvar)]) ]); 
disp (['F2*P     ' sprintf('%12.3f',[dfhp,hpss,hpms,hpvar,(hpvar/totvar)]) ]);
disp (['F1*F2:P  ' sprintf('%12.3f',[dfihp,ihpss,ihpms,ihpvar,(ihpvar/totvar)]) ]);disp([' ']);

end;


% two-facet nested design, as in F1 : F2 : P
if type == 8;

% degrees of freedom
dfp = (size(scores,1) -1);
dfhp = p*(h-1);
dfihp = p*h*(i-1);
dftot = dfp+dfhp+dfihp;

% pointers/counters
cntfirst = 1;
cntlast = i;

% raw data computations needed for later sum of squares calculations
pmean = repmat(0,p,1);
phmean = repmat(-9999,p,h);
hpmean = repmat(-9999,p,h);
for cols = 1:h;
hABCD = scores(:, cntfirst:cntlast);
phmean(:,cols) = sum(hABCD,2)/i;
hpmean(:,cols) = sum(hABCD,2)/i;
cntfirst = cntfirst + i;
cntlast = cntlast + i;
end;

pmean = sum(hpmean,2)/h;
hmean = sum(phmean,1)/p;
gmean = sum(pmean,1)/p;

Tp = (i*h)*(sum(pmean.^2));
Th = (p*i)*sum((hmean.^2),2);
Thp = i*sum(sum(hpmean.^2));
Tihp = sum(sum(scores.^2));
Tu = p*i*h*gmean*gmean;

% sum of squares
pss = Tp-Tu;
hpss = Thp-Tp;
ihpss = Tihp-Thp;
sstot = pss+hpss+ihpss;

% mean squares
pms = pss/dfp;
hpms = hpss/dfhp;
ihpms = ihpss/dfihp;

% variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = (inv([ i*h, i, 1; 0, i, 1; 0, 0, 1 ]) *  [ pms; hpms; ihpms ])';
% setting negative variances to zero
for luper = 1:size(vcs,2); if vcs(1,luper) < 0; 
vcs(1,luper) = 0;
disp(['One or more negative variance estimates have been set to zero']);disp([' ']);
end;end
pvar   = vcs(1,1);
hpvar  = vcs(1,2);
ihpvar = vcs(1,3);
totvar = sum(vcs);

% relative and absolute error variances
relv =  (hpvar/h)+(ihpvar/(i*h));
absv = relv;

% matrices for D-studies
persons  = [ 1, 1, 1 ];
vfacet1  = [ 0, 0, i ];
vfacet2  = [ 0, h, h ];

disp(['ANOVA Table']);disp([' ']);
disp (['Term:              df          SS          MS     Variance        %']);
disp (['P        ' sprintf('%12.3f',[dfp,pss,pms,pvar,(pvar/totvar)]) ]);
disp (['F2:P     ' sprintf('%12.3f',[dfhp,hpss,hpms,hpvar,(hpvar/totvar)]) ]);
disp (['F1:F2:P  ' sprintf('%12.3f',[dfihp,ihpss,ihpms,ihpvar,(ihpvar/totvar)]) ]);disp([' ']);

end;


% G & phi coefficients
relG = pvar/(pvar+relv);
absG = pvar/(pvar+absv);

disp(['Relative Error Variance = ' num2str(relv) ]);disp([' ']);
disp(['Absolute Error Variance = ' num2str(absv) ]);disp([' ']);
disp(['G coefficient   = ' num2str(relG) ]);disp([' ']);
disp(['Phi coefficient = ' num2str(absG) ]);disp([' ']);


% D studies
disp([' ']);disp(['D-Study:']);disp([' ']);
disp(['Entered D-Study values for Facet 1:  ' num2str(dfacet1) ]);disp([' ']);
disp(['Entered D-Study values for Facet 2:  ' num2str(dfacet2) ]);disp([' ']);

dabsv  = dfacet1';
drelv  = dabsv;
drelG  = dabsv;
dabsG  = dabsv;
if (type <= 2); dfacet2 = 1; vfacet2=repmat(1,1,size(vcs,2));end
for f2 = 1:size(dfacet2,2);
resdum = repmat(-9999,size(dfacet1,2),4);
counter = 1;
for f1 = 1:size(dfacet1,2);

absvdum = 0;
relvdum = 0;
for lupe = 2:size(vcs,2);
f2dum = 1;
f1dum = 1;
if ( vfacet1(1,lupe)  > 1 );
f1dum = dfacet1(1,f1);
end;
if ( vfacet2(1,lupe) > 1 );
f2dum = dfacet2(1,f2);
end;
absvdum = absvdum + ( vcs(1,lupe) / (f1dum * f2dum ) );
if ( persons(1,lupe) ~= 0 );
relvdum = relvdum + ( vcs(1,lupe) / (f1dum * f2dum ) );
end;
end;
gdum = vcs(1,1) / ( vcs(1,1) + relvdum );
pdum = vcs(1,1) / ( vcs(1,1) + absvdum );

resdum(counter,:) = [ absvdum, relvdum, gdum, pdum ];
counter = counter + 1;
end;

dabsv = [ dabsv , resdum(:,1) ];
drelv = [ drelv , resdum(:,2) ];
drelG = [ drelG , resdum(:,3) ];
dabsG = [ dabsG , resdum(:,4) ];
end;

disp(['In the D-study results below, the levels of Facet 1 appear in']);
disp(['the first column, and the levels of Facet 2 appear in the first row.']);disp([' ']);

disp(['D-Study Absolute Error Variances']);
disp([0,dfacet2;dabsv]);
disp(['D-Study Relative Error Variances']);
disp([0,dfacet2;drelv]);
disp(['D-Study G Coefficients']);
disp([0,dfacet2;drelG]);
disp(['D-Study Phi Coefficients']);
disp([0,dfacet2;dabsG]);

% saving data for plots
if (graphdat == 1); datagraph = drelv;end
if (graphdat == 2); datagraph = dabsv;end
if (graphdat == 3); datagraph = drelG;end
if (graphdat == 4); datagraph = dabsG;end

% plots
plot ( datagraph(:,1),datagraph(:,2:size(datagraph,2)),'LineWidth',1.5)
axis([  min(datagraph(:,1))  max(datagraph(:,1))  0 1 ]); colordef white;
set(gca,'XTick',dfacet1);set(gca,'YTick',0:.1:1)
set(gca,'FontName','Times', 'FontSize',14 );
xlabel('Levels of Facet 1');
if (graphdat == 1); datagraph = drelv; ylabel('Absolute Error Variance');end
if (graphdat == 2); datagraph = dabsv; ylabel('Relative Error Variance');end
if (graphdat == 3); datagraph = drelG; ylabel('G Coefficient');end
if (graphdat == 4); datagraph = dabsG; ylabel('Phi Coefficient');end
if type <=2; title ('D-Study Results For Levels of Facet 1');end
if type >=3;
title ('D-Study Results For Levels of Facet 1 and Facet 2');
% legend     
textobj = findobj('type', 'text'); set(textobj, 'fontunits', 'points'); 
set(textobj, 'fontsize', 14); set(textobj, 'fontweight', 'normal');
leg = [ 'a         ' ,'b'];
for lupe = 1:size(dfacet2,2); leg = [ leg ; 'Facet 2 = ', int2str(dfacet2(1,lupe)) ];end
leg = leg(2:size(leg,1),:);
legend(leg,0);legend boxoff;
end




