Johnson-Neyman 3 Covariate Case Ancova; acvjn3c[]
 Steve Hunka, University of Alberta. 403-492-7276;steve.hunka@ualberta.ca
 Last Revision: July 19/96

This program carries out the Johnson-Neyman one-way Analysis of Covariance for the case when the regression parameters are not homogeneous among the groups. The number of covariates is restricted to 3; the number of groups must be two or more. Using contour plots, regions of statistical significance are identified as a function of covariate values to which the adjustment in group effect is made and the alpha level employed. The user may also determine the significance of group differences for a selected point set {p,q,s} of the covariates. The general linear model using a mu-model (cell-means model) design matrix is employed in the OLS estimate of the parameters ordered as   (m1,a11,a12,a13,m2,a21,...) in which the terms in m are adjusted group intercepts and the terms in a the regression coefficients. 
   
  A polynomial relating the sum of squares for a specific contrast matrix is formed by
  
            (K'B)'Inverse[K'Inverse[X'X]K]K'B-(dfk MSe F)=0
            
where K is the contrast matrix; X the design matrix; B the parameter estimates, dfk the degrees of freedom for the contrast matrix, MSe the mean square error term and F the value of the F distribution for a given level of alpha. The regions of significance are obtained by contour plotting the polynomial. Any covariate point set {p,q,s} which falls on the surface of the representation, or within its enclosing volume, will be significant at the alpha level selected or less.

  Input to the program is a data matrix of order (N,4) where N represents the number of observations and 4 the number of input variables.  The observations are ordered by group. Each row is an observation and contains in  order, the dependent variable followed by the 3 covariates.
  
  For the 1 and 2 covariate case see MathSource document 0207-953.
  
 
References
1. Johnson, P.O. & Hoyt, C. On determining three dimensional regions of significance. Journal of Experimental Education, Vol. 15, No. 4, 1947, 342-353.
2. Hunka, S. Using Mathematica to Solve Johnson-Neyman Problems, Mathematica in Education, Vol. 3, No. 3, 1994.
3. Hunka S., Identifying Regions of Signficance in ANCOVA Problems Having Nonhomogeneous Regressions, Br. Jr. Math. & Stat. Psych., 1995, 48, 161-188.
4. Hunka, S. & Leighton, J. Defining Johnson-Neyman Regions of Significance in the Three Covariate ANCOVA Using Mathematica, Jr. of Educational and Behavioural Statistics, In press 1996.

Executing the program
The program is executed by acvjn3c[d_,ng_], where d_ is the data matrix of order observations by variables. Observations are ordered by groups, and the first column is taken as the dependent variable.  ng_ is a vector of the number of observations in each group. The polynomial expression representing the relationship between the covariates and SSmcfm is  returned in case the researcher wishes to further explore regions (volumes) of significance.

Mathematica Packages required:
Statistics`ContinuousDistributions`
Graphics`Graphics3D`
Graphics`ContourPlot3D`

<<Statistics`ContinuousDistributions`
<<Graphics`ContourPlot3D`
<<Graphics`Graphics3D`

User interactive input queries
1.  Display X'X? (y/n)
2.  Plot Regression Planes? (y/n); if "y" vector input of contour 
     values for each covariate is requested. It is suggested that the means of the
     dependent variable (Y which is output) be used initially. Option to replot.
3.  Before plotting regions of significance for SSmcfm, the plotting range
     for each covariate and an alpha level is requested. 
4.  Display plots of row contrasts of SSmcfm? (y/n); if "y" the plotting range
     for each covariate and alpha level is requested; option to replot.
5.  User input of weights for symbolic contrasts? (y/n). If "y", input of 
     the number of rows in the contrast, and a matrix of weights requested. This
     option allows the user to define, for example, the Helmert contrast in
     symbolic form.  The plotting range for each covariate and an alpha level is
     requested.
6.  Test significance of specific covariate values? (y/n); if "y" the user is
     requested to input the weights to define a contrast matrix as in (5) and a 3
     element numerical vector containing the covariate values to which an
     adjustment is to be made. A F test of SS is then made. The test can be
     repeated for other values of the covariates and other contrast matrices.

Output of the program
 1) the OLS parameter estimates
 2) F test of SSmcfm
 3) F test of homogeneity of regression for each covariate; contrast matrix used
 4) Minimum and maximum of data
 5) variable means
 6) plots of regression planes (optional)
 7) symbolic contrast matrix for SSmcfm
 8) region of signficance plot for SSmcfm and optionally for each row of the
     symbolic contrast matrix used.
 9) point of significance, i.e, contour values for plots dependent on alpha level 
     which is the value of SS at which significance is just reached.
 10) plot of region of signficance for user defined contrast
 11) F test of SS for user specified contrast and covariate values

Error Messages
1. "Error: Input data must have 4 cols. Dimensions={...}". Indicates that the
     data matrix is incorrectly structured. 
2.  "Error: Sum group sizes not equal to nn". Indicates that the number of 
     observations in each group does not sum to the number of rows in the data
     matrix.

Functions Called:
kbin3[]    to input a numerical value, vector, or matrix; checks appropriatness
kbinstrg[] to input a string and check against acceptable list
getMp[]    to get a part of a matrix
pmat[]      to display a matrix using a set number of columns
ssk[]        to obtain sum of squares for a contrast; if the contrast is in
               symbolic form, a polynomial is returned.
knTrst[]    creates a contrast for SSmcfm
cradd[]     matrix circle + operator
regcnt3d[]  plots regression hyperplanes; The function y=f[x1,x2,x3] is
               evaluated using the regression eqns for each group, then a contour
               plot is done for a specific value of Y for each group given by the
               user.
in43dcp[]  plots regions of significance; 3D contour and shadow plot of the
               the polynomial produced by ssk; plot shows region of significance
chkpts[]    tests a covariate point set for significance
getCon[]    gets the number of rows and the weights to define a contrast
                matrix in symbolic form 
probF        calculates probability of F given ss,mse, and dfs

kbin3[q_,p_] function to accept numeric input from keyboard posing question in q. If p is any integer then an integer is to be input;  p is any real value then an real value is to be input; if p={n} then a n element vector is expected; if p={r,c} then a matrix with r rows and c columns is expected. If the input is inappropriate the question is posed again.  If the input is appropriate it is returned.

kbin3[q_,p_]:=
 Module[{c,ans,vrp,vrp2},
  c=0;
  vrp=" Input Error. Type a vector with # elements=";
  vrp2=" Input Error. Type a matrix of order ";
  While[c==0,
   ans=Input[q];
    Which[
      Head[ans]===Symbol,(*check for character*)
       Print[ans," Input Error; Characters not allowed"];,
       (*integer or real required*)
      Head[p]===Integer || Head[p]===Real, 
      (*chk for vec when scalar wanted*)
      If[Head[p]=!=Head[ans],
       Print[ans," Input Error. Input must be ",Head[p]];,
       c=1];,
       
      Head[p]===List && Length[p]==1, (*vector rqd*)
       If[VectorQ[ans]==False, 
       Print[ans,StringJoin[vrp,ToString[p]]];,
         If[Length[ans]!=First[p],
       Print[ans," Error: Vector requires # of elements=",
             First[p]];,
       c=1]];,
       
      Head[p]===List && Length[p]==2, (*matrix rqd*)
      If[MatrixQ[ans]==False,
      Print[ans,StringJoin[vrp2,ToString[p]]];,
          If[Dimensions[ans]!=p,
      Print[ans," Error; Matrix dimensions must be ",p];,
       c=1]]; 
     ]; (*end Which *)    
   ]; (*end While *)
 Return[ans];
 ]; (* end Module*)

kbinstrg[q,p] used to input a string in response to question in q. p contains acceptable answers in the form of a list,eg., p={{ans1},{ans2},....}. If the input is appropriate it is returned. If the input is inappropriate the question is posed again.

kbinstrg[q_,p_]:=
 Module[{c,ans},
  c=0;
  While[c==0,
   ans=InputString[q];
    Do[If[{ans}==p[[j]],
         c=1; Break[]];
      ,{j,1,Length[p]}]; (* end Do *)  
   If[c!=1,Print[ans," Input Error. Type Again"]];  
   ]; (*end While*)
 Return[ans];
 ]; (* end Module*)

getMp[p_,x_,n_]  gets the pth part of matrix x; x has {n1,n2,...np} rows
          in each part. Operates on vectors also.

getMp[p_,x_,n_]:=   
 Module[{inx,mp},
  inx=Drop[FoldList[Plus,0,n],-1];
  mp=Take[Drop[x,inx[[p]]],n[[p]]];
  Return[mp];
       ];

getCon[ng_] gets the number of rows and the weights to define a symbolic contrast matrix, then using these weights generates the symbolic contrast matrix required. Weights must be entered as a matrix not a vector even if a one row contrast is used. For example, if {{1,1,1,1,-1,-1,-1,-1,0,0,0,0}} is entered, the contrast returned is {{1,p,q,s,-1,-p,-q,-s,0,0,0,0}}.
Calls kbin3

getCon[ng_]:=
  Module[{rws,pqs,wts,cont},
  rws=kbin3["How many rows in the contrast?",0];
  wts=kbin3["Type contrast matrix weights only",
            {rws,4 ng}];
  Print["Contrast Matrix Weights Input",wts];
  pqs={1,p,q,s};
  cont=Partition[Flatten[Table[Flatten[Table[pqs,{ng}]],
       {Length[wts]}]] Flatten[wts],ng 4];
  Print["Contrast Matrix is:"];
  Print[TableForm[cont,TableSpacing->{3,0},
        TableDirections->{Row},TableHeadings->Automatic]];
  Return[cont];
  ];

pmat[m,pc] print a matrix m in units of pc columns

pmat[m_,pc_]:=
 Module[{m1,a,c,indx},
 m1=Transpose[m];
 c=Part[Dimensions[m1],1];
 a=Floor[c/pc];
 indx=Join[Table[pc,{a}],{Mod[c,pc]}];
 If[0==Part[indx,-1],
    indx=Drop[indx,-1]];
 Do[Print[""];
  Print["Table Part=",j,"/",Length[indx]];
  Print[TableForm[Transpose[Take[m1,indx[[j]]]],
       TableHeadings->Automatic,
       TableSpacing->{0,4}]];
  m1=Drop[m1,indx[[j]]];
  ,{j,1,Length[indx]}];
  
 ];

ssk[k_,xx_,pb_]
Function ssk calculates sum of squares for contrast matrix k, X'X, and parameters bp.

ssk[k_,xx_,pb_]:=
(k.pb).(Inverse[k.Inverse[xx].Transpose[k]]).
      (k.pb);

knTrst[k_,ng_]
   creates a contrast matrix for SSmcfm; k_ holds pattern e.g., for 3 groups {1,p,q,s} generates {{1,p,q,s,-1,-p,-q,-s,0,0,0,0}, {0,0,0,0,1,p,q,s,-1,-p,-q,-s}}. ng is the number of groups.

knTrst[k_,ng_]:=
 Module[{t},
  t=Table[
   Join[k,-k,Table[0,{Length[k] (ng-2)}]],
       {ng-1}];(*end Table*)
  Do[
    t[[j]]=RotateRight[t[[j]],{(j-1) (Length[k])}],
      {j,1,ng-1}];
      Return[t];
 ];

cradd[a_,b_]
 Carries out circle + operation using matrix a and b. Used to form a linear model design matrix where a and b appear as block diagonals in the result.

cradd[a_,b_]:=
Join[
Transpose[
  Join[
       Transpose[a],
       Transpose[Table[i j 0,
                     {i,Dimensions[a][[1]]},
                     {Dimensions[b][[2]]}]]]],
Transpose[
  Join[
       Table[i j 0,
             {i,Dimensions[a][[2]]},
             {Dimensions[b][[1]]}],
             Transpose[b] ]]]

probF[ss,dfs,mse,dfe] calculates probability of F. ss is the sum of squares for the hypothesis being tested, dfs is the degree of freedom for ss, mse is the mean square error term having dfe degrees of freedom.
Requires: <<Statistics`ContinuousDistributions`
Mod. Date: April 26/96

probF[ss_,dfs_,mse_,dfe_]:=
 Module[{ms,f,fd,p},
  ms=ss/dfs;
  Print["SS=",ss," df=",dfs," MS=",ms];
  Print["MSe=",mse," dfe=",dfe];
  f=ms/mse;
  fd=FRatioDistribution[dfs,dfe];
  p=1-CDF[fd,f];
  Print["F=",f," Prob=",p];
  ];

regcntr3d[eq_,rng_,ybar_] plots contours in 3d space. The function y=f[x1,x2,x3] is evaluated using the regression eqns for each group, then a contour plot is done for a specific value of y given by the user.The minimum and maximum predicted value of the dependent variable is calculated for each regression equation for the user. Using this information or the mean values of Y the user is asked to select contour values.  If the contour values are not adequate, no planes may appear, or some planes may be missing indicating that the contour value selected does not exist within the range of the plotted data for one or more groups. The plot range is determined by the range of the covariate data. The user is given the option to replot using different contour values.  
Parameters: eq is a matrix holding in each row the equations in the form of {constant, coef 1, coef 2, coef 3} with 1 row per group; rng is a matrix holding by row the minimum and maximum values of the dependent variable followed by the covariates.  
The package Graphics`ContourPlot3D` is required. Calls kbin3, kbinstrg
Mod. Date: Apr 26/96

regcntr3d[eq_,rng_]:=
Module[{ans,pobj,exprs,pr1,pr2,pr3,pr4,pr5,pr6,
        eqn,ycntr,title},
 ans="y";
 exprs=eq.{1,p,q,s};
 pr1=rng[[2]][[1]];pr2=rng[[2]][[2]];
 pr3=rng[[3]][[1]];pr4=rng[[3]][[2]];
 pr5=rng[[4]][[1]];pr6=rng[[4]][[2]];
 
 While[ans=="y",
 pobj=List[];
 ycntr=kbin3["Type vector of 3 contour values",
    {Length[eq]}];
 Print["User selected contour values:",ycntr];
 Do[ (*over number of groups*)
  (*eqn=exprs[[j]];*)
  AppendTo[pobj,
   ContourPlot3D[exprs[[j]],{p,pr1,pr2},{q,pr3,pr4},{s,pr5,pr6},
      Contours->{ycntr[[j]]},MaxRecursion->1,
      Axes->True,
      DisplayFunction->Identity]],
   {j,1,Length[eq]}]; (*end Do *) 
 title=StringInsert["Regression Planes at ",
    ToString[ycntr],-1]; 
 Show[pobj,DisplayFunction->$DisplayFunction,
      Axes->True,AxesLabel->{"p","q","s"},
      PlotLabel->title];
 ans=kbinstrg["Replot at other contour values? y/n",
       {{"y"},{"n"}}];    
  ]; (*end While *)
  ]; (*end Module *)

in43dcp[exp_,ms_,d1_,d2_]; 3D Contour and Shadow Plotting function.
The function carries out a 3D contour and shadow plot of a polynomial of the form z=f(p,q,s) for a specific value of z. It is used in the Johnson-Neyman Ancova analysis where p,q, and s are the three covariates. The argument exp is the polynomial for sum of squares expressed in terms of p,q, and s; ms is the mean square error; d1 is df1; d2 is dfe. The user is asked to input the plot range for each covariate in the form {start,end}, and the alpha level for significance. The MaxRecursion parameter is set to 1 and generally will result in a plot of sufficient resolution within a reasonable time and memory requirements. For a considerable increase in resolution, but with a considerable increase in time and memory requirements, the user can change the option MaxRecursion to 2. Calls kbinstrg, kbin3.
Requires: Statistics`ContinuousDistributions`
               Graphics`ContourPlot3D`
               Graphics`Graphics3D`
Mod. Date: Apr 30/96

in43dcp[exp_,ms_,d1_,d2_]:=
 Module[{rpq,par,sv,quest,q1,siglev,alf,gobj,
         labelalf},
  rpq="y";
  q1="For covariate type plotting range values eg.{1,10}";
  sv={"p ","q ","s "};
  While[rpq=="y", par=List[];
   Do[
     quest=StringInsert[q1,sv[[j]],15];
     AppendTo[par,
     kbin3[quest,{2}]],{j,1,3}];
  par=Flatten[par];
  Print["Current plot ranges=",par];
  alf=kbin3["Type alpha level.",0.0];
  siglev=ms d1 Quantile[FRatioDistribution[d1,d2],1-alf];
  Print["Mse=",ms," df1=",d1," df2=",d2," alpha=",alf];
  Print["SS Point of Significance=",siglev];
  labelalf=StringJoin["Region of Signficance for alpha=",
      ToString[alf]];
  gobj=ContourPlot3D[exp,
      {p,par[[1]],par[[2]]},{q,par[[3]],par[[4]]},
      {s,par[[5]],par[[6]]},MaxRecursion->1,
      Contours->{siglev},DisplayFunction->Identity,
      Axes->True];     
  Shadow[gobj,Axes->True,AxesLabel->{"p","q","s"},
      DisplayFunction->$DisplayFunction,
      PlotLabel->labelalf,
      AxesEdge->{{-1,-1},{1,-1},{-1,-1}},
      FaceGrids->{{-1,0,0},{0,1,0},{0,0,-1}}];
  rpq=kbinstrg["Do you want to replot? y/n",
         {{"y"},{"n"}}];
   ]; (*end While *)
   ]; (*end Module*)

chkPts[x_,b_,mse_,dfe_,ng_] calculates the probability in testing SS of a contrast associated with a particular set of covariate values input by user. x is X'X where X is the design matrix; b is a vector of parameter estimates; mse is the mean square error; dfe is the degree of freedom for mse; ng is the number of groups. First the symbolic contrast is constructed, then user values of the covariate given as {p,q,s} are substituted into the symbolic contrast. ssk is called to get the sum of squares, and probF to carry out the F test. 
Calls kbin3, kbinstrg, getCon, probF, ssk.
Mod. Date: April 27/96

chkPts[x_,b_,mse_,dfe_,ng_]:=
Module[{rpq,con,ss,dfc,pqs,rpcon,consub},
 rpcon="y";
 Print["\nTest of adjustment to user (p,q,s)"];
  While[rpcon=="y",
   rpq="y";
  con=getCon[ng]; (*generates symbolic contrast*)
  While[rpq=="y",
  pqs=kbin3["Type specific covariate values as vector {p,q,s}",
             {3}];
  Print["User input {p,q,s}:",pqs];
  (*now substitute specific values of p,q,s *)
  consub=con /. {p->pqs[[1]],q->pqs[[2]],s->pqs[[3]]};
   Print["Numerical Contrast matrix:"];
   Print[TableForm[consub,TableSpacing->{3,0},
        TableDirections->{Row},TableHeadings->Automatic]];
   dfc=Length[consub];
   ss=ssk[consub,x,b]; (*get SS *)
   probF[ss,dfc,mse,dfe];
   rpq=kbinstrg["For current contrast do you have more covariate values? y/n",
         {{"y"},{"n"}}];
   ];(*end While rpq==y*)
   rpcon=kbinstrg["Do you have another contrast? y/n",
         {{"y"},{"n"}}];
   ]; (*end While rpcon==y*)
   ];(*end Module*)

acvjn3c is the main Mathematica program to do the Johnson-Neyman 1-way ancova analysis. d_ is a data matrix of order (N,4) where N is the number of observations ordered by group membership; ng is a vector of the group sample sizes.
Calls: cradd, getMp, kbinstrg, kbin3, pmat, knTrst, ssk, regcntr3d, in43dcp, getCon, chkPts, probF. 
Requires: Statistics`ContinuousDistributions`, Graphics`ContourPlot3D`,
               Graphics`Graphics3D`
Mod. Date: May 1/96

acvjn3c[d_,ng_]:= (*Johnson-Neyman Ancova*)
 Module[{ans,nt,ngrp,x,y,ncv,xtx,dm,b,dfe,mse,con,ssmcfm,
      concv,ss1,f1,fdist,prob,ss2,ss3,mnx,dt,d1,dfm,
      bm,meanxy,ymns},
  If[4!=Part[Dimensions[d],2],
   Print["Error:Input data must have 4 cols. Dimensions=",
         Dimensions[d]];
   Return[]];
  nt=Length[d];(*N total sample size*) 
  If[Apply[Plus,ng]!=nt,
   Print["Error:Sum group sizes not equal to ",nt];
   Return[]];
  
  y=Flatten[Take[Transpose[d],1]];
  x=Drop[Transpose[d],1];
  ngrp=Length[ng];
  x=Transpose[Join[{Table[1,{nt}]},x]];
  
  (*now set up the design matrix*)
  dm=cradd[getMp[1,x,ng],getMp[2,x,ng]];
  Do[dm=cradd[dm,getMp[j,x,ng]],{j,3,ngrp}];
  xtx=Transpose[dm].dm;
  ans=kbinstrg["Display X'X? y/n",{{"y"},{"n"}}];
  If[ans=="y",pmat[xtx,6]];
  b=N[Inverse[xtx].Transpose[dm].y];
  Print["Mu-Model Parameter Estimates"];
  Print[TableForm[b,TableSpacing->{0,4},
        TableHeadings->Automatic]];
  ssmcfm=(b.xtx.b)-(Apply[Plus,y]^2)/nt;
  dfm=Length[xtx]-1;
  Print["\nSSmodelcfm=",ssmcfm," dfm=",dfm];
  dfe=Dimensions[dm].{1,-1};
  mse=((y.y)-b.Transpose[dm].y)/dfe;
  probF[ssmcfm,dfm,mse,dfe];
   
  (*now test each covariate regression coeff*)
  Print["\nTest of Regression Homogeneity"];
  con={1,0,0,0};
  Do[Print["\nCovariate ",j];
    con=RotateRight[con,1];
    concv=knTrst[con,ngrp];
    Print[TableForm[concv,TableSpacing->{3,0},
      TableDirections->{Row},TableHeadings->Automatic]];
    ss1=ssk[concv,xtx,b];
    d1=Length[concv];
    probF[ss1,d1,mse,dfe];
  ,{j,1,3}];(*end do*)
  
  bm=Partition[b,4];
  mnx=Table[Join[{Min[Transpose[d][[j]]]},
            {Max[Transpose[d][[j]]]}],{j,1,4}];
  Print["\nVariable by Min/Max of data:"];
  Print[TableForm[mnx,TableSpacing->{0,4},
        TableHeadings->Automatic]];
  meanxy=List[];
  Do[ (*do over number of groups*)
  	AppendTo[meanxy,
  	  N[Apply[Plus,getMp[j,d,ng]]/ng[[j]]]];
    ,{j,1,Length[ng]}]; (*end do *)
  Print["\nVariable means Y, P, Q, S:"];
  Print[TableForm[meanxy,TableSpacing->{0,4},
        TableHeadings->Automatic]];
  ymns=Flatten[Take[Transpose[meanxy],1]];
  
  ans=kbinstrg["Plot Regression planes? y/n",{{"y"},{"n"}}];
  Print["\nRegression plots: ",ans];
  If[ans=="y",
    regcntr3d[bm,mnx]];
  
  (*Now do the symbolic processing*)
  Print["\nProcessing Plot for SSmcfm"]; 
  concv=knTrst[{1,p,q,s},ngrp]; (*for SSmcfm*)
  Print["SSmcfm Symbolic Contrast Matrix"];
  Print[TableForm[concv,TableSpacing->{3,0},
        TableDirections->{Row},TableHeadings->Automatic]];
  
  ss2=Simplify[ssk[concv,xtx,b]]; 
  (*Now do the 3D contour plot of region of significance*)
  in43dcp[ss2,mse,d1,dfe];
  
  (*using concv do row contrasts separately*)
  Print["\nProcessing SSmcfm row contrasts"];
  If[ngrp>=3,
  ans=kbinstrg["Display plot of row contrasts of SSmcfm? y/n",
       {{"y"},{"n"}}];
  Print["Display plots of row contrasts of SSmcfm: ",ans];
  If[ans=="y",
   Do[Print["Row Contrast Matrix is:"];
    Print[TableForm[concv[[j]],TableSpacing->{0,3},
        TableDirections->{Column},TableHeadings->Automatic]]; 
    ss3=Simplify[ssk[{concv[[j]]},xtx,b]]; 
   in43dcp[ss3,mse,1,dfe]; 
   ,{j,1,ngrp-1}]; (*end Do*)
    ]]; (*end if ans==y;end if ngrp>=3 *)
    
  (*User contrast wt matrix input next*)
  Print["\nProcessing user weights for contrasts"];
  ans="y";
  While[ans=="y",
  ans=kbinstrg["User input of wts for symbolic contrast? y/n",
       {{"y"},{"n"}}];
  Print["User wts for symbolic contrast: ",ans]; 
  If[ans=="n",Break[]]; 
  concv=getCon[ngrp];
   ss3=Simplify[ssk[concv,xtx,b]];
   in43dcp[ss3,mse,Length[concv],dfe];
    ]; (*end while ans==y*)
    
  (*now allow checking sign.of a specific {p,q,s}*)  
  If["y"==kbinstrg["Test sign.of specific covariate values? y/n",
  {{"y"},{"n"}}],
     chkPts[xtx,b,mse,dfe,ngrp]];
  Return[ss2];
    ]; (*End Module*)

Test data as used in reference 4; the data is taken from: Huitema, B.E. (1980), The Analysis of Covariance and Alternatives, New York; Wiley.

huitdta={{10,1,9,3},{10,2,7,5},{11,2,9,4},{10,3,9,3},
         {11,4,9,3},{11,5,8,4},{10,5,7,4},{11,6,9,4},
         {11.5,6,6,3},{12,7,8,5},{12,8,7,4},{11,8,7,5},
         {11,9,6,6},{12.5,10,6,6},{12,11,6,6},
         {5,1,8,2},{6,1.5,9,4},{6,2.5,9,4},{7,3.5,8,5},
         {8,4.5,8,4},{9,4.5,8,4},{9,5,7,3},{9,6,8,3},
         {10.5,6,6,4},{11,7,7,3},{12.5,7,5,5},
         {12.5,7.5,8,5},{14,8,6,6},{14.5,9,7,5},
         {16,10,6,6}}

User input to program queries for the example given below:
1.  Display X'X?; y
2. Plot Regression planes?; y
    a) Type vector of 3 contour values; {10.75,9.05,11.8}
    b) Replot at other contour values; n
3. For plotting of the region of signficance for SSmcfm;
    a) For covariate p type plotting range values; {0,10}
    b) For covariate q type plotting range values; {0,10}
    c) For covariate s type plotting range values; {0,10}
    d) Type alpha level; .01
    e) Do you want to replot?; y
        (replotted with same ranges but alpha at .001)
4. Display plot of row contrasts of SSmcfm? y
    For each row of the contrast, the following queries:
    a) For covariate p type plotting range values; {0,20}
    b) For covariate q type plotting range values; {0,10}
    c) For covariate s type plotting range values; {0,10}
    d) Type alpha level;  .01
    e) Do you want to replot?; n
        (same input values for next row of contrast)
5. User input of weights for symbolic contrast?; y
    a) How many rows in the contrast?; 1
    b) Type contrast matrix weights only; 
         {{1,1,1,1,-.5,-.5,-.5,-.5,-.5,-.5,-.5,-.5}}
    c) For covariate p type plotting range values; {0,15}
    d) For covariate q type plotting range values; {0,10}
    e) For covariate s type plotting range values; {0,10}
    f) Type alpha level; .01
    g) Do you want to replot?; n
    h) User input of weights for symbolic contrast?; n
6.  Test sign. of specific covariate values?; y
    a) How many rows in the contrast?; 2
    b) Type contrast matrix weights only;
       {{1,1,1,1,-1,-1,-1,-1,0,0,0,0},{0,0,0,0,1,1,1,1,-1,-1,-1,-1}}
    c) Type specific covariate values as a vector {p,q,s}; {3,4,5}
    d) For current contrast do you have more covariate values?; n
    e) Do you have another contrast?; n

Processing Time: Approximately 10 minutes including user input on a Power Macintosh 6100/60.

acvjn3c[huitdta,{10,10,10}]