* SIMPLE-3grps: For Interactions between a continuous
                IDV & a 3-category Moderator.


****************  START OF TRIAL-RUN DATA  *****************************
   The following commands generate artificial data that can be used for
   a trial-run of the program.  Just run this whole file.  
new file.
input program.
loop #a=1 to 50.
compute group=1.
compute idv = normal (10).
compute dv  = normal (10)  - 10.
end case.
end loop.
loop #a=1 to 50.
compute group=2.
compute idv = normal (10).
compute dv  = normal (10)  + 10.
end case.
end loop.
loop #a=1 to 50.
compute group=3.
compute idv = normal (10).
compute dv  = normal (10).
end case.
end loop.
end file.
end input program.
if group=1   dv=idv *  1.55 + dv * sqrt(1 -  .55**2).
if group=2   dv=idv * -1.55 + dv * sqrt(1 - -.55**2).
if group=3   dv=idv *  6.00 + dv * sqrt(1 -  .00**2).
descriptives var = all.

* moderated regression on the artificial data.
if group = 1 or group = 3    dum1 = 0.
if group = 2                 dum1 = 1.
if group = 1 or group = 2    dum2 = 0.
if group = 3                 dum2 = 1.
compute xdum1 = idv * dum1.
compute xdum2 = idv * dum2.
regression 
 /var= idv dum1 dum2 xdum1 xdum2  dv
 /descriptives = corr
 /statistics=defaults zpp bcov 
 /dependent=dv
 /enter idv dum1 dum2
 /test (xdum1 xdum2) .

****************  END OF TRIAL-RUN DATA  ************************************.




* This version of the program reads and processes a raw data file,
  containing the variable scores. (In contrast, the previous  
  version of the program read and processed saved matrix data).

* For analyses of your own data, it may be easiest to create the
  following three variables in your data file: idv, group, & dv
  e.g., compute idv = name of the independent variable;
  e.g., compute group = name of the group variable;
  e.g., compute dv  = name of the dependent variable;  
  There must be no missing values.

* The values of group variable must be either "1", "2", or "3";
  No other values are recognized by the program.

* The program automatically computes dummy codes & product terms.



set mxloops=999999  printback = off length=none width=100.
matrix.

* The GET command below reads the data file & creates a data
  matrix ("datamat"). Enter the name of the dataset (and the
  file path) on the GET command after "file=".  Using a "*"
  instead of a file name & path causes the program to use
  the currently active SPSS data set.

* Then enter the names of the variables after "variables =" on
  the GET command.  Only three variables names are permitted,
  and they must appear in the following order: idv, group, dv.

* If the name of the idv in your dataset is something other
  than "idv", then enter the correct name in the place of "idv"
  after "variables =".

* If the name of the group variable in your dataset is
  something other than "group", then enter the correct name in
  the place of "group" after "variables =". The values of the
  group variable must be either "1", "2", or "3"; 
  No other values are recognized by the program.

* If the name of the dv in your dataset is something other than
  "dv", then enter the correct name in the place of "dv" 
  after "variables =".

get datamat / file=*  / variables = idv, group, dv  / missing=omit.

* The program automatically computes the dummy codes & product terms.

* The program generates data for plots of simple regression lines;
  The IDV range was set at 2 SDs below and 2 SDs above the IDV mean;
  Alternative low and high levels may be specified now -- simply change
  the value for multiIDV from "2.0" to whatever multiple of the SD you prefer.

compute multiIDV = 2.0  .

****************** End of User Specifications *********************.

* separating the data for the 3 groups & creating the dummy codes.
compute datagrp1 = make(1, 2, -9999).
compute datagrp2 = make(1, 2, -9999).
compute datagrp3 = make(1, 2, -9999).
compute dummys   = make(nrow(datamat), 2, -9999).
loop #luper = 1 to nrow(datamat).
do if (datamat(#luper,2)=1).
compute datagrp1 = { datagrp1; datamat(#luper,1), datamat(#luper,3) }.
compute dummys(#luper,:) = { 0, 0 }.
end if.
do if (datamat(#luper,2)=2).
compute datagrp2 = { datagrp2; datamat(#luper,1), datamat(#luper,3) }.
compute dummys(#luper,:) = { 1, 0 }.
end if.
do if (datamat(#luper,2)=3).
compute datagrp3 = { datagrp3; datamat(#luper,1), datamat(#luper,3) }.
compute dummys(#luper,:) = { 0, 1 }.
end if.
end loop.
compute datagrp1 = datagrp1(2:nrow(datagrp1),:).
compute datagrp2 = datagrp2(2:nrow(datagrp2),:).
compute datagrp3 = datagrp3(2:nrow(datagrp3),:).

* error message for illegal group values.
do if (mmin(dummys) = -9999). 
print /title="There is an illegal value for the group variable.".
end if.

* creating the data matrix with dummy codes & the product term.
compute datam={ datamat(:,1), dummys, (datamat(:,1)&*dummys(:,1)), 
               (datamat(:,1)&*dummys(:,2)), datamat(:,3) }.

* n, mean, sd, & correlation matrix (Bernstein, p. 77-79).
compute n = nrow(datam).
compute rawsp = t(datam) * datam .
compute rsums = t(csum(datam)).
compute mn = t(rsums) / n.
compute corsp = rawsp - (1/n) * (rsums) * t(rsums) .
compute vcv = corsp * (1/(n-1)).
compute sd = t(sqrt(diag(vcv))).
compute d = inv(mdiag(sqrt(diag(vcv)))).
compute cr = d * vcv * d.

* n, mean, sscp, & sd for Group 1.
compute n1 = nrow(datagrp1).
compute rawsp1 = t(datagrp1) * datagrp1 .
compute rsums1 = t(csum(datagrp1)).
compute mn1 = t(rsums1) / n1.
compute corsp1 = rawsp1 - (1/n1) * (rsums1) * t(rsums1) .
compute vcv1 = corsp1 * (1/(n1-1)).
compute sd1 = t(sqrt(diag(vcv1))).
compute d1 = inv(mdiag(sqrt(diag(vcv1)))).
compute cr1 = d1 * vcv1 * d1.

* n, mean, sscp, & sd for Group 2.
compute n2 = nrow(datagrp2).
compute rawsp2 = t(datagrp2) * datagrp2 .
compute rsums2 = t(csum(datagrp2)).
compute mn2 = t(rsums2) / n2.
compute corsp2 = rawsp2 - (1/n2) * (rsums2) * t(rsums2) .
compute vcv2 = corsp2 * (1/(n2-1)).
compute sd2 = t(sqrt(diag(vcv2))).
compute d2 = inv(mdiag(sqrt(diag(vcv2)))).
compute cr2 = d2 * vcv2 * d2.

* n, mean, sscp, & sd for Group 3.
compute n3 = nrow(datagrp3).
compute rawsp3 = t(datagrp3) * datagrp3 .
compute rsums3 = t(csum(datagrp3)).
compute mn3 = t(rsums3) / n3.
compute corsp3 = rawsp3 - (1/n3) * (rsums3) * t(rsums3) .
compute vcv3 = corsp3 * (1/(n3-1)).
compute sd3 = t(sqrt(diag(vcv3))).
compute d3 = inv(mdiag(sqrt(diag(vcv3)))).
compute cr3 = d3 * vcv3 * d3.

* Overall regression coeffs.
compute beta = inv(cr(1:5,1:5)) * cr(1:5,6) .
compute b = (sd(1,6) &/ sd(1,1:5))   &* t(beta) .
compute a = mn(1,6) - ( rsum ( mn(1,1:5) &* b ) ) .
compute r2  = t(beta) * cr(1:5,6) .
compute r2main = t(inv(cr(1:3,1:3))*cr(1:3,6))*cr(1:3,6).
compute r2chXn = r2 - r2main.
compute fsquared = (r2 - r2main) / (1 - r2) .
compute F = ((r2-r2main)/2) / ((1-r2)/(n-5-1)).
compute dferror = n - 5 - 1.
compute pF = 1 - fcdf(F,2,dferror) .

print {r2chXn,F,{2},dferror,fsquared,pF} /format="f12.3" 
  /title="Coefficients for the Interaction"
  / space=2 /clabels="Rsq. ch." "F" "df num." "df denom." "fsquared" "Sig. F".
print {t(b),beta}/format="f12.3" 
  /title="Beta weights for the full equation:"
  /rlabels="idv" "dum1" "dum2" "xdum1" "xdum2" /clabels="raw b" "std.beta"  .
print a /title="The intercept is:" .


* simple slope info.
compute dum1 = { 0 ; 1 ; 0 }.
compute dum2 = { 0 ; 0 ; 1 }.
compute slopes={ b(1,1) ; (b(1,1)+b(1,4)) ; (b(1,1)+b(1,5)) }. 
compute aslopes={ a ; (b(1,2)+a) ; (b(1,3)+a) }.
compute mse = (n/(n-5))*(sd(1,6)**2)*(1-r2).
compute Sb=mse*inv((mdiag(sd(1,1:5))*cr(1:5,1:5)*mdiag(sd(1,1:5)))*(n-1)) .
compute SEslopes={ 
  (sqrt ( {1,0,0,dum1(1,1),dum2(1,1)} * Sb * t({1,0,0,dum1(1,1),dum2(1,1)}) )) ;
  (sqrt ( {1,0,0,dum1(2,1),dum2(2,1)} * Sb * t({1,0,0,dum1(2,1),dum2(2,1)}) )) ;
  (sqrt ( {1,0,0,dum1(3,1),dum2(3,1)} * Sb * t({1,0,0,dum1(3,1),dum2(3,1)}) ))   }.
compute tslopes = slopes &/ SEslopes .
compute df = { (n-5-1) ;  (n-5-1) ;  (n-5-1) }.
compute zslopes =slopes &*{sd1(1,1)/sd1(1,2);sd2(1,1)/sd2(1,2);sd3(1,1)/sd3(1,2)}.
compute zSE  =SEslopes &*{sd1(1,1)/sd1(1,2) ;sd2(1,1)/sd2(1,2);sd3(1,1)/sd3(1,2)}.  
compute dfs =  n-5-1 .
compute pslopes = (1 - tcdf(abs(tslopes),dfs)) * 2.

* df & t values -- from Darlington p 516 & Howell 87 p 586 --  p = 05 two-tailed .
compute dft={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,22,24,26,28,
 30,32,34,36,38,40,43,46,49,52,56,60,65,70,75,80,85,90,95,100,110,120,130,
 150,175,200,250,300,400,500,600,700,800,900,1000,1000000000;
 12.706,4.303,3.182,2.776,2.571,2.447,2.365,2.306,2.262,2.228,2.201,2.179,
  2.160,2.145,2.131,2.120,2.110,2.101,2.093,2.086,2.074,2.064,2.056,2.048,
  2.042,2.037,2.032,2.028,2.024,2.021,2.017,2.013,2.010,2.007,2.003,2.000,
  1.997,1.994,1.992,1.990,1.988,1.987,1.985,1.984,1.982,1.980,1.978,1.976,
  1.974,1.972,1.969,1.968,1.966,1.965,1.964,1.963,1.963,1.963,1.962,1.962 }.
compute tabledT = 0.
loop #a = 1 to 59  .
do if (dfs ge dft(1,#a) and dfs < dft(1,#a+1)) .
compute tabledT = dft(2,#a) .
end if.
end loop if (tabledT > 0).
compute confidLo = (zslopes - (tabledT &* zSE))   .
compute confidHi = (zslopes + (tabledT &* zSE))  .

print { aslopes , slopes , tslopes , df , pslopes}  /format="f12.3"  / space=2
 /title="Simple Slope Coefficients for the DV on the IDV for "
+ "Individual Groups:" 
 /rlabels="Group 1"  "Group 2"   "Group 3"
 /clabels="a" "raw b" "t-test" "df" "Sig. T".
print { zslopes , zSE , confidLO, confidHI } /format="f12.3" 
 /title="Standardized Simple Slopes & 95% Confidence Intervals: "
 /rlabels="Group 1"  "Group 2"   "Group 3"
 /clabels="std.beta" "SE" "95%  Low" "95%  Hi".
print {(aslopes(1,1)-aslopes(2,1)) / (slopes(2,1) - slopes(1,1))} 
 /format="f12.3" / space=2
 /title="The simple regression lines for Group 1 & Group 2 intersect at IDV =".
print {(aslopes(1,1)-aslopes(3,1)) / (slopes(3,1) - slopes(1,1))}
 /format="f12.3" 
 /title="The simple regression lines for Group 1 & Group 3 intersect at IDV =".
print {(aslopes(2,1)-aslopes(3,1)) / (slopes(3,1) - slopes(2,1))}
 /format="f12.3" 
 /title="The simple regression lines for Group 2 & Group 3 intersect at IDV =".


* Simultaneous Regions of Significance. 
compute sscp1 = mdiag(sd1(1,1:2)) * cr1 * mdiag(sd1(1,1:2)) * (n1-1).
compute sscp2 = mdiag(sd2(1,1:2)) * cr2 * mdiag(sd2(1,1:2)) * (n2-1).
compute sscp3 = mdiag(sd3(1,1:2)) * cr3 * mdiag(sd3(1,1:2)) * (n3-1).
compute ssreg1 = ((sscp1(1,2)**2)/sscp1(1,1)) .
compute ssreg2 = ((sscp2(1,2)**2)/sscp2(1,1)) .
compute ssreg3 = ((sscp3(1,2)**2)/sscp3(1,1)) .
compute Xhi = { -9999 ; -9999 ; -9999 }.
compute Xlo = { -9999 ; -9999 ; -9999 }.
compute N = n1 + n2 + n3.

* df & Bonferroni F; Huitema, 1980, pp. 293 & 387-388; #comparisons = 3.
compute df = N - 6.
compute dfF={3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,22,24,26,28,30,
 40,60,120,500,1000;
 21.49, 13.49, 10.36, 8.74, 7.78, 7.13, 6.68, 6.34, 6.08, 5.87, 5.70, 5.56,
 5.45, 5.35, 5.26, 5.18, 5.12, 5.06, 4.96, 4.88, 4.81, 4.76, 4.71, 4.54,
 4.39, 4.24, 4.13, 4.11}.
loop #a = 1 to 27  .
do if (df ge dff(1,#a) and df < dff(1,#a+1)).
compute tf = dff(2,#a) .
end if.
end loop.

* comparing groups 1 & 2.
compute ssresd =(sscp1(2,2)-ssreg1) + (sscp2(2,2)-ssreg2) .
compute aa = aslopes(1,1) - aslopes(2,1).
compute bb =  slopes(1,1) -  slopes(2,1).
compute A=((2*tf*-1)/(N-4)) * ssresd * ((1/ssreg1)+(1/ssreg2)) + (bb**2).
compute B=(2*tf/(N-4))*ssresd*((mn1(1,1)/ssreg1)+(mn2(1,1)/ssreg2))+(aa*bb).
compute C=((2*tf*-1)/(N-4)) * ssresd * ( (N/(n1*n2))+
          ((mn1(1,1)**2)/ssreg1)+((mn2(1,1)**2)/ssreg2) ) + (aa**2).
do if ( (B**2 - A*C) gt 0).
compute hi = (-B + (sqrt(B**2 - A*C)) ) / A.
compute lo = (-B - (sqrt(B**2 - A*C)) ) / A.
do if (hi > lo).
compute Xhi(1,1) = hi.
end if.
do if (lo < hi).
compute Xlo(1,1) = lo.
end if.
end if.

* comparing groups 1 & 3.
compute ssresd =(sscp1(2,2)-ssreg1) + (sscp3(2,2)-ssreg3) .
compute aa = aslopes(1,1) - aslopes(3,1).
compute bb =  slopes(1,1) -  slopes(3,1).
compute A=((2*tf*-1)/(N-4)) * ssresd * ((1/ssreg1)+(1/ssreg3)) + (bb**2).
compute B=(2*tf/(N-4))*ssresd*((mn1(1,1)/ssreg1)+(mn3(1,1)/ssreg3))+(aa*bb).
compute C=((2*tf*-1)/(N-4)) * ssresd * ( (N/(n1*n3))+
          ((mn1(1,1)**2)/ssreg1)+((mn3(1,1)**2)/ssreg3) ) + (aa**2).
do if ( (B**2 - A*C) gt 0).
compute hi = (-B + (sqrt(B**2 - A*C)) ) / A.
compute lo = (-B - (sqrt(B**2 - A*C)) ) / A.
do if (hi > lo).
compute Xhi(2,1) = hi.
end if.
do if (lo < hi).
compute Xlo(2,1) = lo.
end if.
end if.

* comparing groups 2 & 3.
compute ssresd =(sscp2(2,2)-ssreg2) + (sscp3(2,2)-ssreg3) .
compute aa = aslopes(2,1) - aslopes(3,1).
compute bb =  slopes(2,1) -  slopes(3,1).
compute A=((2*tf*-1)/(N-4)) * ssresd * ((1/ssreg2)+(1/ssreg3)) + (bb**2).
compute B=(2*tf/(N-4))*ssresd*((mn2(1,1)/ssreg2)+(mn3(1,1)/ssreg3))+(aa*bb).
compute C=((2*tf*-1)/(N-4)) * ssresd * ( (N/(n2*n3))+
          ((mn2(1,1)**2)/ssreg2)+((mn3(1,1)**2)/ssreg3) ) + (aa**2).
do if ( (B**2 - A*C) gt 0).
compute hi = (-B + (sqrt(B**2 - A*C)) ) / A.
compute lo = (-B - (sqrt(B**2 - A*C)) ) / A.
do if (hi > lo).
compute Xhi(3,1) = hi.
end if.
do if (lo < hi).
compute Xlo(3,1) = lo.
end if.
end if.

print /title=
 "Simultaneous Regions of Significance -- Johnson-Neyman Technique:" / space=2.
print /title=
   "Using Bonferroni F, p = .05, 2-tailed; see Huitema, 1980, p. 293.".
print /title=
 "The regression lines for group comparisons are significantly different".
print {xlo, xhi} /format="f12.3" 
  /title="at IDV scores < Lo Value & > Hi Value:"
  /rlabels="Grps 1&2"  "Grps 1&3"  "Grps 2&3"
  /clabels="Lo Value" "Hi Value".
print /title="-9999 indicates that meaningful values could not be computed.".


* data for plot.
compute idvlo1 = mn1(1,1) - (sd1(1,1) * multiIDV).
compute idvhi1 = mn1(1,1) + (sd1(1,1) * multiIDV).
compute idvlo2 = mn2(1,1) - (sd2(1,1) * multiIDV).
compute idvhi2 = mn2(1,1) + (sd2(1,1) * multiIDV).
compute idvlo3 = mn3(1,1) - (sd3(1,1) * multiIDV).
compute idvhi3 = mn3(1,1) + (sd3(1,1) * multiIDV).
compute idv    = {idvlo1; idvhi1; idvlo2; idvhi2; idvlo3; idvhi3} .
compute dv=({slopes(1,1);slopes(1,1);slopes(2,1);slopes(2,1);
             slopes(3,1);slopes(3,1)} &* idv)+{aslopes(1,1);
    aslopes(1,1);aslopes(2,1);aslopes(2,1);aslopes(3,1);aslopes(3,1)}.
compute group    = { 1 ; 1 ; 2 ; 2 ; 3 ; 3 }.
compute data = { group , idv , dv }.
print data /format="f12.3"  / space=2
  /title="Data for simple slope plots:" /clabels="Group" "IDV" "DV".
save data /outfile=*  / var=group idv dv .

end matrix.

* The following PLOT command can be used instead of the GRAPH 
   command in SPSS version 11 and earlier.
* plot vsize=15/hsize=50/ format=contour(2) / plot=dv with idv by group.

* The SPSS GRAPH command has few options for controlling the graphs;
  Lines for groups sometimes have gaps, in which case you should assume
  that separated lines of the same color should be joined together;
  Enter the "Data for simple slope plots:" into a more flexible graphing
  program, such as DeltaGraph, for better displays.
  
graph   / line = mean(dv) by idv by group .

graph   / scatter = idv with dv by group.