/***********************************************************************/.
/* Multi-Response Permutation Procedure (MRPP)                         */.
/***********************************************************************/.



/***********************************************************************/.
/* Multiresponse Permutation Procedure (MRPP)                         */.
/* One-way analysis of location/scale shift for a priori groups        */.
/* For discriptions of the algorithm see Mielke, Berry & Johnson (1976)*/.
/*     Commun. Statist. - Theor. Meth., A(5), 14, pp. 1409 - 1424.     */.
/* This SPSS macro written by Li Cai (Dec., 2002 - Jun., 2004)         */.
/* L. L. Thurstone Psychometric Lab, Department of Psychology          */.
/* University of North Carolina at Chapel Hill                         */.
/***********************************************************************/.

set workspace = 20480.
set mxloop = 1000000000.

define MRPP(dist = !tokens(1) !default(1)
           /comm = !tokens(1) !default(1)
           /weight = !tokens(1) !default(1)
           /iv = !tokens(1)
           /dvs = !charend('/')).
matrix.

print /title "          Multi-Response Permutation Procedure (MRPP)".
/* Read from working SPSS file */
get rawdata /variables=!iv !dvs /missing = omit.

/* Get dimensionality of response measurement */
compute r = ncol(rawdata)-1.

/* Get raw multi-dimensional measurement without grouping information */
compute y = rawdata(:,2:(r+1)).

/* Code levels of the grouping variable using the DESIGN function */
compute levelsm = design(rawdata(:,1)).

/* Count number of columns in the design matrix to get number of groups */
compute g = ncol(levelsm).

/* Determine sample size for each group */
compute sizem = t(csum(levelsm)).

/* Determing total sample size by counting rows in the design matrix */
compute L = nrow(levelsm).
do if (L < 6).
  print /title "Combined sample size must be greater than or equal to 6.".
  print /title "Output may be not interpretable.".
end if.

compute v = !dist.
compute comm = !comm.
compute omega = ident(r).
compute weight = !weight.
do if (weight EQ 1).
  compute CiM = sizem/L.
end if.
do if (weight EQ 2).
  compute CiM = (sizem - 1) / (L - g).
end if.

/* Compute PermL(k), where k = 2...6 */
compute nume = {1:L; 1:L; 1:L; 1:L; 1:L}.
compute nume = t(nume).
compute denom = nume.
loop i = 1 to 5.
  compute denom((L-i):L,i) = make((i+1),1,1).
end loop.
compute seq = nume &/ denom.
compute result = make(1,5,1).
loop i = 1 to nrow(seq).
  compute result = result &* seq(i,:).
end loop.
compute PermL = {1,result}.

/* Compute nperm(i,j), where i = 1..g, j = 2..6 */
compute nperm = make(g,6,1).
loop i = 1 to g.
  compute n = sizem(i,1).
  compute nume = {1:n; 1:n; 1:n; 1:n; 1:n}.
  compute nume = t(nume).
  compute denom = nume.
  loop j = 1 to 5.
    do if ((n-j-1) GT 0).
      compute denom((n-j):n,j) = make((j+1),1,1).
    else if ((n-j-1) EQ 0).
      compute denom(:,j) = make(n,1,1).
    end if.
  end loop.
  compute seq = nume &/ denom.
  compute result = make(1,5,1).
  loop j = 1 to nrow(seq).
    compute result = result &* seq(j,:).
  end loop.
  compute nperm(i,:) = {1, result}.
end loop.

/* Check nperm(i,j) */
loop i = 1 to g.
  loop j = 2 to 6.
    do if (nperm(i,j) EQ 1).
      compute nperm(i,j) = 0.
    end if.
  end loop.
end loop.

/* Compute the vector of Group Sizes choose 2 */
compute n2vector = make(g,1,0).
loop i = 1 to g.
  compute n = sizem(i,1).
  compute nume = {1:n}.
  compute denom = {1:(n-2),2,1}.
  compute seq = nume &/ denom.
  compute result = 1.
  loop j = 1 to ncol(seq).
    compute result = result * seq(1,j).
  end loop.
  compute n2vector(i,1) = result.
end loop.

/* Define the matrix for distances in combined sample */
compute ydevm = make(L,L,0).

/* Handle commensuration input here */
do if (r GE 2).
  do if (comm EQ 1).
    compute phi = make(r,1,0).
    loop i = 1 to r.
      compute yt = y(:,i).
        loop j = 1 to L.
          compute ydevm(:,j) = abs((yt - yt(j,1)))&**v.
        end loop.
      compute phi(i,1) = (msum(ydevm)/2)&**(1/v).
    end loop.
    compute omega = mdiag(phi)*mdiag(phi).
    print omega /title "Euclidean commensuration invoked: print the Phi matrix".
  else if (comm EQ 2).
    compute means = csum(y)/L.
    compute d = y.
    loop i = 1 to r.
      compute d(:,i) = make(L,1,means(1,i)).
    end loop.
    compute omega = sscp(y-d).
    print omega /title "Hotelling Commensuration invoked: print SSCP matrix".
  else if (comm EQ 0).
    print /title "No commensuration invoked".
  end if.
else if (r EQ 1).
  print /title "Commensuration not needed: dimensionality of response = 1".
end if.

/* Compute distance super matrix for all observations */
compute invomega = inv(omega).
loop i = 1 to L.
  loop j = 1 to L.
    compute d = y(j,:) - y(i,:).
    compute ydevm(j,i) = d*invomega*t(d).
  end loop.
end loop.
compute ydevm = ydevm&**(v/2).

/* Define the Summary matrix for group sizes and within-group distances */
compute summarym = make(g,2,0).
compute summarym(:,1) = sizem.

/* Compute model parameters for MRPP null distribution */
compute d1j = csum(ydevm).
compute d2j = csum(ydevm&**2).
compute d3j = csum(ydevm&**3).
compute d1 = msum(d1j).
compute d2 = msum(d2j).
compute d3 = msum(d3j).

compute Do1 = (1/permL(2))*d1.
compute Do2 = (1/permL(2))*d2.
compute Do2o1 = (1/permL(3))*(msum(d1j&**2)-d2).
compute Do2o2 = (1/permL(4))*(d1&**2-4*permL(3)*Do2o1-2*d2).
compute Do3 = (1/permL(2))*d3.
compute Do3o1 = (1/permL(3))*(msum(d1j&*d2j)-d3).
compute Do3o2 = (1/permL(4))*(d1*d2-4*permL(3)*Do3o1-2*d3).

compute Do3oa = 0.
loop i = 1 to (L-2).
  loop j = (i+1) to (L-1).
    loop k = (j+1) to L.
      compute Do3oa = Do3oa + ydevm(i,j)*ydevm(i,k)*ydevm(j,k).
    end loop.
  end loop.
end loop.
compute Do3oa = (6/permL(3))*Do3oa.

compute Do3oaa = 0.
loop i = 1 to (L-1).
  loop j = (i+1) to L.
    compute Do3oaa = Do3oaa + ydevm(i,j)*d1j(1,i)*d1j(1,j).
  end loop.
end loop.

compute Do3oaa = (1/permL(4))*(2*Do3oaa-2*permL(3)*Do3o1-permL(3)*Do3oa-d3).
compute Do3oaaa = (1/permL(4))*(msum(d1j&**3)-3*permL(3)*Do3o1-d3).
compute Do3o3 = (1/permL(5))*(permL(3)*d1*Do2o1-4*permL(4)*Do3oaa-
                 2*permL(4)*Do3oaaa-4*permL(3)*Do3o1-2*permL(3)*Do3oa).
compute Do3o4 = (1/permL(6))*(permL(4)*d1*Do2o2-8*permL(5)*Do3o3-
                 4*permL(4)*Do3o2-8*permL(4)*Do3oaa).

/* Compute the first three moments of the MRPP null distribution */
compute mu = Do1.
compute variance = 2*(msum((CiM&**2)&*(nperm(:,2)&**(-1))) - permL(2)&**(-1))*(Do2-2*Do2o1+Do2o2)
                   +4*(msum((CiM&**2)&*(sizem&**(-1)))-1/L)*(Do2o1-Do2o2).
compute delta3 = 4*msum((CiM&**3)&*(nperm(:,2)&**(-2)))*Do3+8*msum((CiM&**3)&*nperm(:,3)&*(nperm(:,2)&**(-3)))*(3*Do3o1+Do3oa)
                 +8*msum((CiM&**3)&*nperm(:,4)&*(nperm(:,2)&**(-3)))*(3*Do3oaa+Do3oaaa)
                 +6*msum((CiM&**2)&*(1-CiM+CiM&*nperm(:,4)&*(nperm(:,2)&**(-2)))&*(nperm(:,2)&**(-1)))*Do3o2
                 +12*msum((CiM&**2)&*((1-CiM)&*nperm(:,3)+CiM&*nperm(:,5)&*(nperm(:,2)&**(-1)))&*(nperm(:,2)&**(-2)))*Do3o3
                 +msum(CiM&*((1-CiM)&*(1-2*CiM)+3*CiM&*(1-CiM)&*nperm(:,4)&*(nperm(:,2)&**(-2))
                       +(CiM&**2)&*nperm(:,6)&*(nperm(:,2)&**(-3))))*Do3o4.

compute skew = (delta3-3*mu*variance-mu&**3)/(sqrt(variance)&**3).

/* Compute within-group distances for each group */
loop grpscntr = 1 to g.
  compute grpsize = sizem(grpscntr,1).
  compute indexm = make(grpsize,1,0).
  compute grpdata = make(grpsize,r,0).
  compute idxcntr = 1.

  /* Construct an index vector to locate the observations for each group */
  /* This is a different approach from the indicator variable approach used */
  /* by Mielke and Berry.  However, the two are indeed equivalent. */

  /* Loop through the design matrix, whenever a '1' is met, store the index */
  loop i = 1 to L.
    do if (levelsm(i,grpscntr) EQ 1).
      compute indexm(idxcntr,1) = i.
      compute idxcntr = idxcntr + 1.
    end if.
  end loop.
  
  /* Obtaining a matrix of data just for one group by index casting */
  compute grpdata = y(indexm,:).
  /* Define the within-group distance matrix */
  compute grpdevm = make(grpsize,grpsize,0).
  loop i = 1 to grpsize.
    loop j = 1 to grpsize.
      compute d = grpdata(j,:) - grpdata(i,:).
      compute grpdevm(j,i) = d*invomega*t(d).
    end loop.
  end loop.
  /* Set distance power and store it in the summary matrix */
  compute grpdevm = grpdevm&**(v/2).
  compute summarym(grpscntr,2) = msum(grpdevm)/2.
end loop.

/* Compute the realized value of delta */
compute delta = 0.
loop grpscntr = 1 to g.
  compute delta = delta+(1/n2vector(grpscntr,1))*summarym(grpscntr,2)*CiM(grpscntr,1).
end loop.
compute Tobs = (delta-mu)/sqrt(variance).

do if (v EQ 1).
  print /title "Metric Euclidean distance function invoked".
else if (v EQ 2).
  print /title "Squared Euclidean distance function invoked".
end if.
print L /title "Total sample size" /format F12.0.
compute summarym(:,2) = summarym(:,2)&/n2vector.
compute summarym = {summarym, CiM}.
print summarym /title "Group summary" /clabels "N" "Distance" "Weights" /format F14.12.
compute moments = {mu, variance, skew}.
compute df = 8/(skew*skew).
do if (skew LT 0).
  compute pval = 1-chicdf(((-Tobs)*sqrt(2*df)+df),df).
else if (skew EQ 0).
  compute pval = cdfnorm(Tobs).
else if (skew GT 0).
  compute pval = chicdf((Tobs*sqrt(2*df)+df),df).
end if.
compute tests = {delta, Tobs, pval}.
print moments /title "Moments of the MRPP null distribution" /clabels "Mean" "Variance" "Skewness" /format F14.12.
print tests /title "Test of Significance" /clabels "Delta" "T" "Sig." /format F14.12.
end matrix.
!enddefine.

/*************************************************************/.
/* MACRO LADreg - TAILORED AFTER AMOEBA IN NUMERICAL RECIPES */.
/* USING THE DOWNHILL SIMPLEX MINIMIZATION ALGORITHM TO FIND */.
/* PARAMETER ESTIMATES IN A UNIVARIATE                       */.
/* LEAST ABSOLUTE DEVIATION (LAD) REGRESSION MODEL.          */.
/* WRITTEN BY LI CAI (JUNE, 2002 - JULY, 2004)               */.
/* L.L. THURSTONE PSYCHOMETRIC LAB, DEPARTMENT OF PSYCHOLOGY */.
/* UNIVERSITY OF NORTH CAROLINA - CHAPEL HILL                */.
/*************************************************************/.

set workspace = 20480.
set mxloop = 100000000.
define LADreg(ftol = !tokens(1) !default(.000001)
             /maxit = !tokens(1) !default(200)
             /const= !tokens(1) !default(1)
             /dv = !tokens(1) 
             /ivs = !charend('/')).
matrix.
print /title "          Least Absolute Deviation (LAD) Regression Using".
print /title "               Downhill Simplex Minimization Method".
compute const = !const.
compute ftol = !ftol.
compute maxit = !maxit.
get rawdata /variables=!dv !ivs /names=names /missing = omit.
compute dvname = names(1,1).
compute ivnames = names(1,2:ncol(names)).

compute y = rawdata(:,1).
compute N = nrow(rawdata).
compute x = rawdata(:,2:ncol(rawdata)).
compute ndim = ncol(rawdata)-1.
do if (const EQ 1).
  compute intcept = make(N,1,1).
  compute x = {intcept,x}.
  compute ndim = ncol(rawdata).
  compute ivnames = {"Intercept",ivnames}.
end if.
print dvname /title = "Criterion variable is:"/format A8.
print N /title = "Sample Size" /format F8.0.
compute b = inv(t(x)*x)*t(x)*y.
compute alpha = 1.
compute beta = .5.
compute gamma = 2.
compute npts = ndim+1.
compute fctval = make(npts,1,0).
compute Po = ident(ndim).
loop i = 1 to ndim.
  compute Po(i,:) = Po(i,:)+t(b).
end loop.
compute P = {t(b);Po}.
print P /title "Starting Simplex:" /format F8.6.
loop i = 1 to npts.
  compute fctval(i) = csum(abs(y-x*t(P(i,:)))).
end loop.
compute counter = 0.
loop i = 1 to maxit.
  compute ind = grade(fctval).
  loop j = 1 to npts.
    do if (ind(j) EQ npts).
      compute ihi = j.
    end if.
    do if (ind(j) EQ (npts-1)).
      compute inhi = j.
    end if.
    do if (ind(j) EQ 1).
      compute ilo = j.
    end if.
  end loop.
  compute rtol = 2*abs(fctval(ihi)-fctval(ilo))/abs(fctval(ihi)+fctval(ilo)).
  do if (rtol LT ftol).
    break.
  end if.
  compute Pbr = t(csum(P)-P(ihi,:))/ndim.
  compute pr = (1+alpha)*Pbr-alpha*t(P(ihi,:)).
  compute fctpr = csum(abs(y-x*pr)).
  do if (fctpr LE fctval(ilo)).
    compute prr = gamma*pr+(1-gamma)*Pbr.
    compute fctprr = csum(abs(y-x*prr)).
    do if (fctprr LT fctval(ilo)).
      compute P(ihi,:)= t(prr).
      compute fctval(ihi)= fctprr.
    else.
      compute P(ihi,:)= t(pr).
      compute fctval(ihi)= fctpr.
    end if.
  else if (fctpr GE fctval(inhi)).
    do if (fctpr LT fctval(ihi)).
      compute P(ihi,:)= t(pr).
      compute fctval(ihi)= fctpr.
    end if.
    compute prr = beta*t(P(ihi,:))+(1-beta)*Pbr.
    compute fctprr = csum(abs(y-x*prr)).
    do if (fctprr LT fctval(ihi)).
      compute P(ihi,:)= t(prr).
      compute fctval(ihi)= fctprr.
    else.
      compute Po = make(npts,ndim,0).
      loop j = 1 to npts.
        compute Po(j,:) = P(ilo,:).
      end loop.
      compute P = .5*(P+Po).
      loop j = 1 to npts.
        compute fctval(j) = csum(abs(y-x*t(P(j,:)))).
      end loop.
    end if.
  else.
    compute P(ihi,:) = t(pr).
    compute fctval(ihi) = fctpr.
  end if.
  compute counter = counter + 1.
end loop.
print P /title "Final Simplex:" /format F8.6.
print counter /title "Number of iterations:" /format F12.2.
compute finalval = csum(abs(y-x*t(P(1,:)))).
print finalval /title "Sum of absolute values of the residuals:" /format F12.2.
print P(1,:) /cnames=ivnames /title "LAD Parameter Estimates" /format F8.6.

end matrix.
!enddefine.

/* Reads in data for MRPP analyses of Multivariate Observations */.
data list free/grp v1 v2 v3.
begin data.
    1.00   20.00   31.00   26.00
    1.00   17.00   29.00   25.00
    1.00   16.00   30.00   27.00
    1.00   18.00   31.00   28.00
    1.00   15.00   28.00   25.00
    2.00   13.00   24.00   20.00
    2.00   13.00   27.00   21.00
    2.00   15.00   26.00   21.00
    2.00   15.00   27.00   18.00
    3.00   23.00   32.00   29.00
    3.00   22.00   31.00   30.00
    3.00   22.00   36.00   32.00
    3.00   23.00   36.00   30.00
    3.00   22.00   33.00   33.00
    3.00   21.00   35.00   31.00
    3.00   25.00   35.00   33.00
    3.00   24.00   34.00   30.00
    4.00   16.00   28.00   24.00
    4.00   18.00   28.00   26.00
    4.00   18.00   29.00   25.00
    4.00   20.00   31.00   25.00
end data.

/* The following command produces output in Figure 1. */.
mrpp dist=1 comm=0 weight=1 iv=grp dvs=v1 v2 v3.

/* The following command produces output in Figure 2. */.
mrpp dist=1 comm=1 weight=1 iv=grp dvs=v1 v2 v3.

/* The following command produces output in Figure 3. */.
mrpp dist=2 comm=2 weight=2 iv=grp dvs=v1 v2 v3.

/* Run a MANOVA to compare the Bartlett-Nanda-Pillai test */.
glm v1 v2 v3 by grp.

/* Reads in data for MRPP regression analysis of LAD residuals */.
data list free/y x1 x2 x3 tc a b.
begin data.
    1.00    1.00    1.00     .00   1.000    1.00    1.00
    4.00    1.00    1.00     .00   1.000    1.00    1.00
     .00    1.00    1.00     .00   1.000    1.00    1.00
    7.00    1.00    1.00     .00   1.000    1.00    1.00
   15.00   -1.00    1.00     .00   2.000    2.00    1.00
    6.00   -1.00    1.00     .00   2.000    2.00    1.00
   10.00   -1.00    1.00     .00   2.000    2.00    1.00
   13.00   -1.00    1.00     .00   2.000    2.00    1.00
   13.00    1.00     .00    1.00   3.000    1.00    2.00
    5.00    1.00     .00    1.00   3.000    1.00    2.00
    7.00    1.00     .00    1.00   3.000    1.00    2.00
   15.00    1.00     .00    1.00   3.000    1.00    2.00
    6.00   -1.00     .00    1.00   4.000    2.00    2.00
   18.00   -1.00     .00    1.00   4.000    2.00    2.00
    9.00   -1.00     .00    1.00   4.000    2.00    2.00
   15.00   -1.00     .00    1.00   4.000    2.00    2.00
    9.00    1.00   -1.00   -1.00   5.000    1.00    3.00
   16.00    1.00   -1.00   -1.00   5.000    1.00    3.00
   18.00    1.00   -1.00   -1.00   5.000    1.00    3.00
   13.00    1.00   -1.00   -1.00   5.000    1.00    3.00
   14.00   -1.00   -1.00   -1.00   6.000    2.00    3.00
    7.00   -1.00   -1.00   -1.00   6.000    2.00    3.00
    6.00   -1.00   -1.00   -1.00   6.000    2.00    3.00
   13.00   -1.00   -1.00   -1.00   6.000    2.00    3.00
end data.

/* The following command produces output in Figure 4. */.
LADREG dv=y ivs=x1 x2 x3.

/* Compute residuals from the LAD regression */.
compute res = y-(10.00000-x1-3*x2+x3).
execute.

/* The following command produces output in Figure 5. */.
mrpp iv=tc dvs=res.

/* Run an ANOVA to test the interaction effect */.
glm y by a b.