/* For the trial analyses, the program was set to read a raw data
   file named and then process this data file within
   PROC IML -- see the "USE" command near the start of the program.
   Alternatively, the READ or EDIT statements (PROC IML commands)
   can also be used to read a data file.  */

/* data from Griffin & Gonzalez (1995) -- 1s & 2s were inserted in col 2 */
data datgg1;
  input dyad 2-3 person 6 var1 10-13 var2 19-20 var3 26-27;
  cards;
  4  1   52.5     20      5
  4  2   43.5     25      3
  6  1   48.5     64     17
  6  2   41.0     44      4
  7  1   33.5     50      2
  7  2   34.5     19      2
 12  1   45.5     38      8
 12  2   64.5     21     10
 14  1   25.0     26      0
 14  2   26.0     20      2
 17  1   10.5      4      1
 17  2   11.5      6      1
 19  1   51.0     35     12
 19  2   51.5     22      6
 21  1   55.0     27     20
 21  2   63.0     21      3
 22  1   36.5     10      3
 22  2   38.5     19      5
 23  1   59.0     27      4
 23  2   38.5     37      2
 27  1   48.0     29      0
 27  2   53.5     11     22
 29  1   45.5     25      6
 29  2   55.0     25     13
 30  1   40.0      6      3
 30  2   36.0      4      0
 31  1   27.0     76      4
 31  2   32.5     37     10
 33  1   44.0     29      5
 33  2   34.0     14      2
 34  1   45.5     64     45
 34  2   67.5     33     13
 35  1   30.0      9      0
 35  2   29.5     14      0
 37  1    0.0      0      0
 37  2    0.0      1      1
 40  1   31.0     13      0
 40  2   50.5     21      3
 42  1   41.5     26      2
 42  2   42.0     27     11
 43  1   33.5     20      9
 43  2   29.0     11      3
 45  1   62.0     25     22
 45  2   57.0     21      9
 47  1   51.0     23      1
 47  2   58.0     22      5
 48  1    1.5      2      0
 48  2    1.5      5      0
;
run;
proc means data=datgg1; var dyad person var1 var2 var3;
run;


options nocenter nodate nonumber  linesize=100 pagesize=500;  title;
proc iml;
reset noname;

/* The Correlational Analysis of Dyad-Level Data
   in the Distinguishable-Case
   Gonzalez & Griffin, 1999, Personal Relationships, 6, 449-469 */

/* Create a data matrix where rows = cases, & columns = variables;
   the first collumn of data must contain the dyad numbers
   (the dyad id# given to both dyad members); the second collumn
   contains the person id #s (e.g., sex) coded as 1 or 2;
   the third and subsequent collumns contain the variables to be analyzed.

   the data files must have been sorted by dyad id# and by
   person id # (e.g., "sort cases by dyad, person");

   Cases with missing values are not permitted in the data file;

   Use the following name for the data matrix: "data" */


/* The following command causes the program to process the SAS
   data file that was created above.  */
use datgg1;
read all var _num_ into data;


/* Data can also be entered into the program directly, as in the
   following example using the "data = {" command. */

/* data from Griffin & Gonzalez (1995) -- 1s & 2s were inserted in col 2 */
data = {
  4  1   52.5     20      5,
  4  2   43.5     25      3,
  6  1   48.5     64     17,
  6  2   41.0     44      4,
  7  1   33.5     50      2,
  7  2   34.5     19      2,
 12  1   45.5     38      8,
 12  2   64.5     21     10,
 14  1   25.0     26      0,
 14  2   26.0     20      2,
 17  1   10.5      4      1,
 17  2   11.5      6      1,
 19  1   51.0     35     12,
 19  2   51.5     22      6,
 21  1   55.0     27     20,
 21  2   63.0     21      3,
 22  1   36.5     10      3,
 22  2   38.5     19      5,
 23  1   59.0     27      4,
 23  2   38.5     37      2,
 27  1   48.0     29      0,
 27  2   53.5     11     22,
 29  1   45.5     25      6,
 29  2   55.0     25     13,
 30  1   40.0      6      3,
 30  2   36.0      4      0,
 31  1   27.0     76      4,
 31  2   32.5     37     10,
 33  1   44.0     29      5,
 33  2   34.0     14      2,
 34  1   45.5     64     45,
 34  2   67.5     33     13,
 35  1   30.0      9      0,
 35  2   29.5     14      0,
 37  1    0.0      0      0,
 37  2    0.0      1      1,
 40  1   31.0     13      0,
 40  2   50.5     21      3,
 42  1   41.5     26      2,
 42  2   42.0     27     11,
 43  1   33.5     20      9,
 43  2   29.0     11      3,
 45  1   62.0     25     22,
 45  2   57.0     21      9,
 47  1   51.0     23      1,
 47  2   58.0     22      5,
 48  1    1.5      2      0,
 48  2    1.5      5      0  };

norig = nrow(data);
tailed = 2;

/* collumn 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;

/* correlation matrix   module */
start corrcoef(matname);
cases = nrow(matname);
nm1 = cases - 1;
zscores = j(cases, ncol(matname), -999);
do loopc = 1 to ncol(matname);
mean = sum(matname[,loopc]) / cases;
sd = sqrt ( ssq((matname[,loopc]- mean)) / nm1  );
zscores[,loopc] = (matname[,loopc] - mean) / sd;
end;
rdata = (1 / nm1 ) * ( t(zscores) * zscores );
return(rdata);
finish;

/* selecting pairs that have the same dyad #s, and where sex=1 & sex=2  */
datarev=j(1,ncol(data),-9999);
do aa = 1 to nrow(data)-1;
if (data[aa,1]=data[aa+1,1] & data[aa,2]=1 & data[aa+1,2]=2) then do;
datarev = ( datarev // data[aa:(aa+1),] );
end;end;

data = datarev[2:nrow(datarev),2:ncol(data)];

n = nrow(data) / 2;

/* setting up pairwise ("double-counted") collumns of data  */
xp = j(nrow(data),(ncol(data)-1),-9999);
do aa = 1 to nrow(data)/2;
xp[aa*2-1,] = data[aa*2,2:ncol(data)];
xp[aa*2,] = data[aa*2-1,2:ncol(data)];
end;
m = ( data || xp );

cr = corrcoef(m);

/* correlations between residuals -- Velicer, 1976 */
c22 = cr[2:nrow(cr),2:ncol(cr)];
c11 = cr[1,1];
c12 = cr[1,2:ncol(cr)];
c22part = c22 - t(c12) * inv(c11) * c12 ;
d = diag(1 / (sqrt(vecdiag(c22part))));
resd = d * c22part * d ;

overall = resd[1:(nrow(resd)/2), 1:(ncol(resd)/2)];
iccross = resd[((nrow(resd)/2)+1):nrow(resd),1:(nrow(resd)/2)];

/* p-Levels of the intraclass correlations  */
intra   = vecdiag(iccross);
zintra  =  (sqrt(n) * intra) / (1 - intra##2);
pzintra = (1 - probnorm(abs(zintra))) * tailed;

/* tests of assumptions involving variances & covariances   */

/* matrix of pairwise ("double-counted") collumns of data for just one sex */
half = j(1, (ncol(m)-1), -9999);
do aa = 1 to nrow(m);
if (m[aa,1] = 1) then do;
half = ( half// m[aa,(2:ncol(m))]  );
end;end;
half = half[2:nrow(half),];

/* variance-covariance matrix & correlation matrix   */
nr2 = nrow(half);
rawsp2 = t(half) * half ;
rsums2 = t(csum(half));
corsp2 = rawsp2 - (1/nr2) * (rsums2) * t(rsums2) ;
vcv2 = corsp2 * (1/(nr2-1));
d2 = inv(diag(sqrt(vecdiag(vcv2))));
cr2 = d2 * vcv2 * d2;

first  = vcv2[1:(nrow(vcv2)/2), 1:(ncol(vcv2)/2)];
second=vcv2[(nrow(vcv2)/2+1):nrow(vcv2),(ncol(vcv2)/2+1):ncol(vcv2)];
iccross2=cr2[((nrow(cr2)/2)+1):nrow(cr2),1:(nrow(cr2)/2)];
iccrossv=vcv2[((nrow(vcv2)/2)+1):nrow(vcv2),1:(nrow(vcv2)/2)];

/* tests for equality of variances   */
vartests = j(nrow(first), 6, -9999);
do bb = 1 to nrow(first);
vartests[bb,1]=( bb );
vartests[bb,2:3]=( first[bb,bb] || second[bb,bb]  );
vartests[bb,4]=((first[bb,bb]-second[bb,bb])*sqrt(nrow(half)-2))
   / (2*sqrt(first[bb,bb]*second[bb,bb]*(1-(iccross2[bb,bb]**2))));
vartests[bb,5]=nrow(half)-2;
vartests[bb,6]=(1-probt(abs(vartests[bb,4]),vartests[bb,5]))*tailed;
end;

/* tests for equality of the within-category population covariances */
covarswi = j(1 , 6, -9999);
do cc = 1 to ncol(first)-1;
do rr = (cc+1) to nrow(first);
line     = j(1 , 6, -9999);
line[1,1:2] = ( cc || rr );
line[1,3:4] = ( first[cc,rr] || second[cc,rr] );
line[1,5]=((first[cc,rr]-second[cc,rr])*sqrt(nrow(half))) /
  sqrt(first[cc,cc]*first[rr,rr]+second[cc,cc]*second[rr,rr]+
  2*(((first[cc,rr]+second[cc,rr])/2)  **2) -
  2*(iccrossv[cc,cc]*iccrossv[rr,rr]+iccrossv[cc,rr]*iccrossv[rr,cc]));
line[1,6]=(1 - probnorm(abs(line[1,5])))*tailed;
covarswi = ( covarswi // line );
end;
end;
covarswi = covarswi[2:nrow(covarswi),];

/* tests for equality of the across-category population covariances */
covarscr = j(1 , 6, -9999);
do cc = 1 to ncol(first)-1;
do rr = (cc+1) to nrow(first);
line     = j(1 , 6, -9999);
line[1,1:2] = ( cc || rr );
line[1,3:4] = ( iccrossv[cc,rr] || iccrossv[rr,cc] );
line[1,5]=((iccrossv[cc,rr]-iccrossv[rr,cc])*sqrt(nrow(half)))/
  sqrt(first[cc,cc]*second[rr,rr]+second[cc,cc]*first[rr,rr]+
  2*(((iccrossv[cc,rr]+iccrossv[rr,cc])/2)  **2) -
  2*(first[cc,rr]*second[cc,rr]+iccrossv[cc,cc]*iccrossv[rr,rr]));
line[1,6]=(1 - probnorm(abs(line[1,5])))*tailed;
covarscr = ( covarscr // line );
end;
end;
covarscr = covarscr[2:nrow(covarscr),];


/* initializing */
nstarovr   = j(nrow(overall),ncol(overall),-9999);
zoverwi    = nstarovr;
pzoverwi   = nstarovr;
zovercr    = nstarovr;
pzovercr   = nstarovr;
indcorrs   = nstarovr;
tind       = nstarovr;
ptind      = nstarovr;
dyadcors   = nstarovr;
dstar      = nstarovr;
zdyad      = nstarovr;
pzdyad     = nstarovr;

do ii = 1 to nrow(overall);
do jj = 1 to ncol(overall);
if (ii ^= jj) then do;

/* p-Levels of the overall within-partner correlations  */
nstarovr[ii,jj] = (2*n) /
       (1+iccross[ii,ii]*iccross[jj,jj]+iccross[ii,jj]**2);
zoverwi[ii,jj] = overall[ii,jj]*sqrt(nstarovr[ii,jj]);
pzoverwi[ii,jj] = (1-probnorm(abs(zoverwi[ii,jj])))*tailed;

/* p-Levels of the overall cross-partner correlations  */
zovercr[ii,jj] = iccross[ii,jj]*sqrt(nstarovr[ii,jj]);
pzovercr[ii,jj] = (1-probnorm(abs(zovercr[ii,jj])))*tailed;

/* individual-level correlations  */
indcorrs[ii,jj]=(overall[ii,jj]-iccross[ii,jj])/
        (sqrt(1-iccross[ii,ii])*sqrt(1-iccross[jj,jj]));
if abs(indcorrs[ii,jj]) < 1 then do;
tind[ii,jj] = (indcorrs[ii,jj]*sqrt(n-2)) /
        sqrt(1-(indcorrs[ii,jj]**2));
ptind[ii,jj] = (1 - probt(abs(tind[ii,jj]),n-2))*tailed;
end;

/* dyad-level correlations  */
if (iccross[ii,ii] > 0 & iccross[jj,jj] > 0) then do;
dyadcors[ii,jj]=iccross[ii,jj]/
        (sqrt(iccross[ii,ii])*sqrt(iccross[jj,jj]));
dstar[ii,jj]  = (2*n / (1+iccross[ii,ii]*iccross[jj,jj]
        +overall[ii,jj]**2))*iccross[ii,ii]*iccross[jj,jj];
zdyad[ii,jj]  = dyadcors[ii,jj]*sqrt(dstar[ii,jj]);
pzdyad[ii,jj] = (1 - probnorm(abs(zdyad[ii,jj])))*tailed;
end;

end;
end;
end;

/* preparing the output matrices */
owpart  = j(1,5,-9999);
ocpart  = owpart;
dyadrzp = owpart;
withrzp = j(1,6,-9999);
do cc = 1 to ncol(overall)-1;
do rr = (cc+1) to nrow(overall);
owpart=( owpart // (cc || rr || overall[cc,rr] ||
        zoverwi[cc,rr] || pzoverwi[cc,rr] ) );
ocpart=( ocpart // (cc || rr || iccross[cc,rr] ||
        zovercr[cc,rr] || pzovercr[cc,rr] ));
dyadrzp=( dyadrzp // (cc || rr || dyadcors[cc,rr] ||
        zdyad[cc,rr] || pzdyad[cc,rr] ));
withrzp=( withrzp // (cc || rr || indcorrs[cc,rr] ||
        tind[cc,rr] || (n-2) || ptind[cc,rr] ));
end;
end;
owpart  = owpart[2:nrow(owpart),];
ocpart  = ocpart[2:nrow(ocpart),];
dyadrzp = dyadrzp[2:nrow(dyadrzp),];
withrzp = withrzp[2:nrow(withrzp),];


print "The Correlational Analysis of Dyad-Level Data";
print "in the Distinguishable-Case";
print "Gonzalez & Griffin, 1999, Personal Relationships, 6, 449-469.";
print "All Significance Tests Are Two-Tailed";
print "Number of Individual Cases in the Data File:", norig;
print "Number of Individual Cases That Will be Processed:", (nrow(data));

print "Variance-Covariance Matrix for the Primary Variables:";
print vcv2 [format=7.4];

print "Tests for Equality of Variances:";
labels = ( "Var.#" ||"Variance" ||"Variance" ||"t" ||"df" ||"p");
print vartests [colname=labels format=9.4];

print "Tests for Equality of the Within-Person Covariances:";
labels = ( "Var.#" || "Var.#" ||  "Covar." || "Covar." || "z" || "p");
print covarswi [colname=labels format=9.4];

labels = ("Var.#"|| "Var.#" ||  " r"|| "z"|| "p");
print "Overall Within-Partner Correlations, z-values & p-Levels:";
print owpart [colname=labels format=9.4];

labels = ("Var.#"|| "Var.#" ||  "Covar."|| "Covar."|| "z"|| "p");
print "Tests for Equality of the Across-Person Covariances:";
print covarscr  [colname=labels format=9.4];

labels = ("Var.#"|| "Var.#" ||  " r"|| "z"|| "p");
print "Overall Cross-Partner Correlations, z-values & p-Levels:";
print ocpart [colname=labels format=9.4];

labels = ("Var.#"|| "r"|| "z"|| "p");
print "Intraclass Correlations, z-values & p-Levels:";
print ( t(1:nrow(intra))||intra||zintra||pzintra) [colname=labels format=9.4];

labels = ("Var.#"|| "Var.#" ||  " r"|| "z"|| "p");
print "Dyad-Level Correlations, z-values & p-Levels:";
print dyadrzp [colname=labels format=9.4];

labels = ("Var.#"|| "Var.#" ||  " r"|| "t"|| "df"|| "p");
print "Individual-Level Correlations, t-values & p-Levels:";
print withrzp [colname=labels format=9.4];

print,"-9999 indicates that a value could not be computed.";
print "The computations sometimes produce dyad-level";
print "correlations > 1.00 or < -1.00.  This is most likely to";
print "occur when an intraclass correlations is weak.";
print "See Griffin & Gonzalez, 1995, p. 434-437 for details.";

quit;