options nocenter nodate number linesize=90; title; proc iml; /* There are two main ways of entering data into this program. Method 1: You can specify the data to be analyzed on USE and READ statements, as in the example below. The USE statement indicates the SAS file containing the data. On the READ statement, ALL indicates that all of the cases will be used in the analyses. To include only some of the variables from the data file in the analyses, insert "var {variable names}" before INTO. The name of the data matrix that will be used within PROC IML must be "SCORES" (as in the example). /* use libraryname.filename; read all into scores; */ /* Method 2: A data matrix can be entered directly, within the PROC IML environment, as in the trial-run data below. The data matrix must be named "scores"; the data must appear inside curly brackets; place a space between the data values for each case; and place a comma at the end of the data values for each case. */ /**************** START OF TRIAL-RUN DATA ************************************/ /* 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 }; /* The following data are from Brennan, 2001, p 72, Table 3.1, for a p x i x o design, with 10 persons, 4 items, and 2 occasions. To analyze your own data, the trial-run data matrix must be removed/deleted. */ 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 }; /**************** 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 = 2; /* 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 "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 = 3; /* Enter D-study values for Facet 1; enter the values inside curly brackets, and leave a space between the values. */ dfacet1 = {5 10 15 20}; /* Enter D-study values for Facet 2; enter the values inside curly brackets, and leave a space between the values. You can ignore this step for single-facet designs. */ dfacet2 = {1 2 3 4 5}; /* At the very bottom of this file, after the QUIT statement, is a PLOT command that can be used to plot 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. Now just run this whole file. */ /***************************************************************************/ /* row sums module */ start rsum(matname); rsums =j(nrow(matname),1); do rows = 1 to nrow(matname); dumr = matname[rows,]; rsums[rows,1]=sum(dumr); end; return(rsums); finish; /* column sums module */ start csum(matname); csums =j(1,ncol(matname)); do cols = 1 to ncol(matname); dumc = matname[,cols]; csums[1,cols]=sum(dumc); end; return(csums); finish; /* row sum of squares module */ start rssq(matname); rssqs =j(nrow(matname),1); do rows = 1 to nrow(matname); dumr = matname[rows,]; rssqs[rows,1]=ssq(dumr); end; return(rssqs); finish; /* column sum of squares module */ start cssq(matname); cssqs =j(1,ncol(matname),-9999); do cols = 1 to ncol(matname); dumc = matname[,cols]; cssqs[1,cols]=ssq(dumc); end; return(cssqs); finish; print 'GENERALIZABILITY THEORY ANALYSES:' ; if type = 1 then; print "Design Type 1: single-facet fully-crossed design, as in P * F1" ; if type = 2 then; print "Design Type 2: single-facet nested nested design, as in F1 : P"; if type = 3 then; print "Design Type 3: two-facet fully-crossed design, as in P * F1 * F2"; if type = 4 then; print "Design Type 4: two-facet nested design, as in P * (F1 : F2)"; if type = 5 then; print "Design Type 5: two-facet nested design, as in (F1 : P) * F2"; if type = 6 then; print "Design Type 6: two-facet nested design, as in F1 : (P * F2)"; if type = 7 then; print "Design Type 7: two-facet nested design, as in (F1 * F2) : P"; if type = 8 then; print "Design Type 8: two-facet nested design, as in F1 : F2 : P"; p = nrow(scores); i = nfacet1; h = nfacet2; print "Number of persons/objects ('P'):" p; print "Number of levels for Facet 1 ('F1'):" nfacet1; if type > 2 then; print "Number of levels for Facet 2 ('F2'):" nfacet2; /* single-facet fully-crossed design, as in P * F1 */ if type = 1 then do; /* 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 = rsum(hABCD)/i; p2mean = cssq(pmean); imean = csum(hABCD)/p; i2mean =rssq(imean); pi2tot = sum(hABCD##2); gmean = rsum(imean)/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 = t(inv( (i||0||1) // (0||p||1) // (0||0||1) ) * ( pms // ims // pims )); /* setting negative variances to zero */ do luper = 1 to ncol(vcs); if(vcs[1,luper] < 0) then do; vcs[1,luper] = 0; print, 'One or more negative variance estimates have been set to zero'; 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 ); c = {"df" "SS" "MS" "Variance" "Proportion"}; r = {"P" "F1" "P*F1"}; print,,, "ANOVA Table:", ( (dfp || pss || pms || pvar || (pvar/totvar)) // (dfi || iss || ims || ivar || (ivar/totvar)) // (dfpi || piss || pims || pivar || (pivar/totvar)) ) [colname=c rowname=r format=10.3]; end; /* single-facet nested nested design, as in F1 : P */ if type = 2 then do; /* degrees of freedom */ dfp = p-1; dfip = p*(i-1); /* raw data computations needed for later sum of squares calculations */ hABCD = scores; pmean = (rsum(hABCD))/i; p2mean = ssq(pmean); pi2mean = ssq(hABCD); mean = sum(pmean)/p; /* sum of squares */ pss = i*p2mean-p*i*mean*mean; ipss = (pi2mean)-(i*p2mean); /* mean squares */ pms = pss/dfp; ipms = ipss/dfip; /* variance components: Brennan, 2001, p 79, 83, 85, 442 */ vcs = t(inv( (i||1) // (0||1) ) * ( pms // ipms )); /* setting negative variances to zero */ do luper = 1 to ncol(vcs); if(vcs[1,luper] < 0) then do; vcs[1,luper] = 0; print, 'One or more negative variance estimates have been set to zero'; 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 ); c = {"df","SS","MS" ,"Variance" ,"Proportion"}; r = {"P","F1:P"}; print,,, "ANOVA Table:", ( (dfp || pss || pms || pvar || (pvar/totvar) ) // (dfip || ipss || ipms || ipvar || (ipvar/totvar) ) ) [colname=c rowname=r format=10.3]; end; /* two-facet fully-crossed design, as in P * F1 * F2 */ if type = 3 then do; /* degrees of freedom */ dfp = (nrow(scores) - 1); dfo = (h - 1); dfi = ((ncol(scores)/(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 = j(p,i,0); pmean = j(p,1,0); omean = j(h,1,-9999); pomean = j(p,h,-9999); iomean = j(h,i,-9999); do cols = 1 to h; hABCD = scores[,cntfirst:cntlast]; iomean[cols,] = csum(hABCD)/p; pomean[,cols] = rsum(hABCD)/i; pimean = pimean+hABCD; cntfirst = cntfirst + i; cntlast = cntlast + i; end; pmean = rsum(pomean)/h; omean = csum(pomean)/p; pimean = pimean/h; imean = csum(pimean)/p; mean = csum(pmean)/p; Tpers = (i*h)*(ssq(pmean)); Titem = (p*h)*(rssq(imean)); Tocca = (p*i)*(ssq(omean)); Tu = (p)*(i)*(h)*((mean)*(mean)); Tpi = (h)*(ssq(pimean)); Tpo = cssq(pomean); Tpo = rsum(Tpo)*i; Tio = (p)*(ssq(iomean)); Tpoi = (ssq(scores)); /* 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 = t(inv(bigP) * ( pms // ims // oms // pims // poms // ioms // poims )); /* setting negative variances to zero */ do luper = 1 to ncol(vcs); if(vcs[1,luper] < 0) then do; vcs[1,luper] = 0; print, 'One or more negative variance estimates have been set to zero'; 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 ); c = {"df","SS","MS" ,"Variance" ,"Proportion"}; r = {"P","F1","F2","P*F1","P*F2","F1*F2","P*F1*F2"}; print,,, "ANOVA Table:", ( (dfp || pss || pms || pvar || (pvar/totvar) ) // (dfi || iss || ims || ivar || (ivar/totvar) ) // (dfo || oss || oms || ovar || (ovar/totvar) ) // (dfpi || piss || pims || pivar || (pivar/totvar) )// (dfpo || poss || poms || povar || (povar/totvar) )// (dfio || ioss || ioms || iovar || (iovar/totvar) )// (dfpoi || poiss || poims || poivar || (poivar/totvar)) ) [colname=c rowname=r format=10.3]; end; /* two-facet nested design, as in P * (F1 : F2) */ if type = 4 then do; /* degrees of freedom */ dfp = (nrow(scores) - 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 = j(p,1,0); hmean = j(h,1,-9999); phmean = j(p,h,-9999); ihmean = j(h,i,-9999); do cols = 1 to h; hABCD = scores[,cntfirst:cntlast]; ihmean[cols,] = csum(hABCD)/p; phmean[,cols] = rsum(hABCD)/i; cntfirst = cntfirst + i; cntlast = cntlast + i; end; pmean = rsum(phmean)/h; hmean = csum(phmean)/p; mean = csum(pmean)/p; Tp = (i*h)*(ssq(pmean)); Th = (p*i)*(rssq(hmean)); Tih = p*(ssq(ihmean)); Tph = i*(ssq(phmean)); Tpih = ssq(scores); Tu = p*i*h*mean*mean; /* 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 = t(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 */ do luper = 1 to ncol(vcs); if(vcs[1,luper] < 0) then do; vcs[1,luper] = 0; print, 'One or more negative variance estimates have been set to zero'; 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 ); c = {"df","SS","MS","Variance","Proportion"}; r = {"P","F2","F1:F2","P*F2","P*F1:F2"}; print,,, "ANOVA Table:", ( (dfp || pss || pms || pvar || (pvar/totvar) ) // (dfh || hss || hms || hvar || (hvar/totvar) ) // (dfih || ihss || ihms || ihvar || (ihvar/totvar) ) // (dfph || phss || phms || phvar || (phvar/totvar) ) // (dfpih || pihss || pihms || pihvar || (pihvar/totvar)) ) [colname=c rowname=r format=10.3]; end; /* two-facet nested design, as in (F1 : P) * F2 */ if type = 5 then do; /* degrees of freedom */ dfp = (nrow(scores) -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 = j(p,1,0); hmean = j(h,1,-9999); phmean = j(p,h,-9999); ipmean = j(p,i,0); do cols = 1 to h; hABCD = scores[,cntfirst:cntlast]; phmean[,cols] = rsum(hABCD)/i; ipmean = ipmean+ hABCD; cntfirst = cntfirst + i; cntlast = cntlast + i; end; pmean = rsum(phmean)/h; hmean = csum(phmean)/p; mean = csum(pmean)/p; ipmean = ipmean/h; Tp = (i*h)*(ssq(pmean)); Th = (p*i)*(rssq(hmean)); Tip = h*(ssq(ipmean)); Tph = (i*(ssq(phmean))); Thip = ssq(scores); Tu = p*i*h*mean*mean; /* 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 = t(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 */ do luper = 1 to ncol(vcs); if(vcs[1,luper] < 0) then do; vcs[1,luper] = 0; print, 'One or more negative variance estimates have been set to zero'; 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 ); c = {"df","SS","MS","Variance","Proportion"}; r = {"P","F2","F1:P","P*F2","F2*F1:P"}; print,,, "ANOVA Table:", ( (dfp || pss || pms || pvar || (pvar/totvar) ) // (dfh || hss || hms || hvar || (hvar/totvar) ) // (dfip || ipss || ipms || ipvar || (ipvar/totvar) ) // (dfph || phss || phms || phvar || (phvar/totvar) ) // (dfhip || hipss || hipms || hipvar || (hipvar/totvar)) ) [colname=c rowname=r format=10.3]; end; /* two-facet nested design, as in F1 : (P * F2) */ if type = 6 then do; /* degrees of freedom */ dfp = (nrow(scores) -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 = j(p,1,0); hmean = j(h,1,-9999); phmean = j(p,h,-9999); do cols = 1 to h; hABCD = scores[,cntfirst:cntlast]; phmean[,cols] = rsum(hABCD)/i; cntfirst = cntfirst + i; cntlast = cntlast + i; end; pmean = rsum(phmean)/h; hmean = csum(phmean)/p; mean = csum(pmean)/p; Tp = (i*h)*(ssq(pmean)); Th = (p*i)*(rssq(hmean)); Tph = (i*(ssq(phmean))); Tiph = ssq(scores); Tu = p*i*h*mean*mean; /* 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 = t(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 */ do luper = 1 to ncol(vcs); if(vcs[1,luper] < 0) then do; vcs[1,luper] = 0; print, 'One or more negative variance estimates have been set to zero'; 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 ); c = {"df","SS","MS","Variance","Proportion"}; r = {"P","F2","P*F2","F1:P*F2"}; print,,, "ANOVA Table:", ( (dfp || pss || pms || pvar || (pvar/totvar) ) // (dfh || hss || hms || hvar || (hvar/totvar) ) // (dfph || phss || phms || phvar || (phvar/totvar) ) // (dfiph || iphss || iphms || iphvar || (iphvar/totvar)) ) [colname=c rowname=r format=10.3]; end; /* two-facet nested design, as in (F1 * F2) : P */ if type = 7 then do; /* degrees of freedom */ dfp = (nrow(scores) -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 = j(p,1,0); phmean = j(p,h,-9999); ipmean = j(p,i,0); hpmean = j(p,h,-9999); do cols = 1 to h; hABCD = scores[,cntfirst:cntlast]; phmean[,cols] = rsum(hABCD)/i; ipmean = ipmean+ hABCD; hpmean[,cols] = rsum(hABCD)/i; cntfirst = cntfirst + i; cntlast = cntlast + i; end; pmean = rsum(hpmean)/h; hmean = csum(phmean)/p; mean = csum(pmean)/p; ipmean = ipmean/h; Tp = (i*h)*(ssq(pmean)); Th = (p*i)*(ssq(hmean)); Tip = h*(ssq(ipmean)); Thp = (i*(ssq(hpmean))); Tihp = ssq(scores); Tu = p*i*h*mean*mean; /* 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 = t(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 */ do luper = 1 to ncol(vcs); if(vcs[1,luper] < 0) then do; vcs[1,luper] = 0; print, 'One or more negative variance estimates have been set to zero'; end;end; pvar = vcs[1,1]; ipvar = vcs[1,2]; hpvar = vcs[1,3]; ihpvar = vcs[1,4]; totvar = sum(vcs); /* matrices for D-studies */ persons = ( 1 || 1 || 1 || 1 ); vfacet1 = ( 0 || i || 0 || i ); vfacet2 = ( 0 || 0 || h || h ); c = {"df","SS","MS","Variance","Proportion"}; r = {"P","F1*P","F2*P","F1*F2:P"}; print,,, "ANOVA Table:", ( (dfp || pss || pms || pvar || (pvar/totvar) ) // (dfip || ipss || ipms || ipvar || (ipvar/totvar) ) // (dfhp || hpss || hpms || hpvar || (hpvar/totvar) ) // (dfihp || ihpss || ihpms || ihpvar || (ihpvar/totvar) ) ) [colname=c rowname=r format=10.3]; end; /* two-facet nested design, as in F1 : F2 : P */ if type = 8 then do; /* degrees of freedom */ dfp = (nrow(scores) -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 = j(p,1,0); phmean = j(p,h,-9999); hpmean = j(p,h,-9999); do cols = 1 to h; hABCD = scores[,cntfirst:cntlast]; phmean[,cols] = rsum(hABCD)/i; hpmean[,cols] = rsum(hABCD)/i; cntfirst = cntfirst + i; cntlast = cntlast + i; end; pmean = rsum(hpmean)/h; hmean = csum(phmean)/p; mean = csum(pmean)/p; Tp = (i*h)*(ssq(pmean)); Th = (p*i)*(rssq(hmean)); Thp = (i*(ssq(hpmean))); Tihp = ssq(scores); Tu = p*i*h*mean*mean; /* 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 = t(inv( (i*h||i||1) // (0||i||1) // (0||0||1 ) ) * ( pms // hpms // ihpms )); /* setting negative variances to zero */ do luper = 1 to ncol(vcs); if(vcs[1,luper] < 0) then do; vcs[1,luper] = 0; print, 'One or more negative variance estimates have been set to zero'; 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 ); c = {"df","SS","MS","Variance","Proportion"}; r = {"P","F2:P","F1:F2:P"}; print,,, "ANOVA Table:", ( (dfp || pss || pms || pvar || (pvar/totvar) ) // (dfhp || hpss || hpms || hpvar || (hpvar/totvar) ) // (dfihp || ihpss || ihpms || ihpvar || (ihpvar/totvar)) ) [colname=c rowname=r format=10.3]; end; /* G & phi coefficients */ relG = pvar/(pvar+relv); absG = pvar/(pvar+absv); print,,"Error Variances:", ( relv || absv ) [colname={"Relative","Absolute"} format=10.3]; print,,"G-coefficients:", ( relG || absG ) [colname={"G", "Phi"} format=10.3]; /* D study */ print,,'D-Study:'; print,'Entered D-Study values for Facet 1:', dfacet1; print,'Entered D-Study values for Facet 2:', dfacet2; dabsv = t(dfacet1); drelv = dabsv; drelG = dabsv; dabsG = dabsv; dabsv2 = { -9999 -9999 -9999 }; drelv2 = dabsv2; drelG2 = dabsv2; dabsG2 = dabsv2; if type <= 2 then do; dfacet2 = 1; vfacet2 = j(1,ncol(vcs),1); end; do f2 = 1 to ncol(dfacet2); resdum = j(ncol(dfacet1),4,-9999); counter = 1; do f1 = 1 to ncol(dfacet1); absvdum = 0; relvdum = 0; do lupe = 2 to ncol(vcs); f2dum = 1; f1dum = 1; if vfacet1[1,lupe] > 1 then; f1dum = dfacet1[1,f1]; if vfacet2[1,lupe] > 1 then; f2dum = dfacet2[1,f2]; absvdum = absvdum + ( vcs[1,lupe] / (f1dum * f2dum ) ); if persons[1,lupe] ^= 0 then; relvdum = relvdum + ( vcs[1,lupe] / (f1dum * f2dum ) ); end; gdum = vcs[1,1] / ( vcs[1,1] + relvdum ); pdum = vcs[1,1] / ( vcs[1,1] + absvdum ); resdum[counter,] = ( absvdum || relvdum || gdum || pdum ); dabsv2 = ( dabsv2 // ( dfacet2[1,f2] || dfacet1[1,f1] || absvdum) ); drelv2 = ( drelv2 // ( dfacet2[1,f2] || dfacet1[1,f1] || relvdum) ); drelG2 = ( drelG2 // ( dfacet2[1,f2] || dfacet1[1,f1] || gdum) ); dabsG2 = ( dabsG2 // ( dfacet2[1,f2] || dfacet1[1,f1] || pdum) ); counter = counter + 1; end; dabsv = ( dabsv || resdum[,1] ); drelv = ( drelv || resdum[,2] ); drelG = ( drelG || resdum[,3] ); dabsG = ( dabsG || resdum[,4] ); end; dabsv2 = dabsv2[2:nrow(dabsv2),]; drelv2 = drelv2[2:nrow(drelv2),]; drelG2 = drelG2[2:nrow(drelG2),]; dabsG2 = dabsG2[2:nrow(dabsG2),]; print, "In the D-study results below, the levels of Facet 1 appear in"; print "the first column, and the levels of Facet 2 appear in the first row."; print,"D-Study Absolute Error Variances", ((0 || dfacet2) // dabsv) [format=10.3]; print,"D-Study Relative Error Variances", ((0 || dfacet2) // drelv) [format=10.3]; print,"D-Study G Coefficients", ((0 || dfacet2) // drelG) [format=10.3]; print,"D-Study Phi Coefficients", ((0 || dfacet2) // dabsG) [format=10.3]; /* saving data for plots */ if graphdat = 1 then do; create plotdata from drelv2[colname={"Facet2" "Facet1" "Gstat"}]; append from drelv2;end; if graphdat = 2 then do; create plotdata from dabsv2[colname={"Facet2" "Facet1" "Gstat"}]; append from dabsv2;end; if graphdat = 3 then do; create plotdata from drelG2[colname={"Facet2" "Facet1" "Gstat"}]; append from drelG2;end; if graphdat = 4 then do; create plotdata from dabsG2[colname={"Facet2" "Facet1" "Gstat"}]; append from dabsG2;end; quit; proc plot data = plotdata vpct=100 hpct = 100; title 'Gstat is the coefficient you specified on the GRAPHDAT= statement'; plot gstat*facet1=facet2 / vaxis=0 to 1 by .1 ;run;