Ancova Program for solving the Johnson-Neyman Problem for 1 and 2 Covariates

Steve Hunka, University of Alberta, Edmonton, AB, Canada, T6G 2G5,
e-mail: steve.hunka@ualberta.ca

Program Description: The program handles the 1 and 2 covariate case of analysis of covariance when homogeneity of regression coefficients does not exist. Arguments to the program are the data matrix and a vector of the number of observations per group. Parameters are estimated using the general linear model approach. The significance of each covariate is tested, and a plot is made of the regression lines or planes. In the case of one covariate the unknown is solved directly and plotted to show the boundaries for the region of significance. In the two covariate case, the polynomial surface is plotted including a plot of the surface defining the region of significance as well as a contour plot identifying various levels of the significance regions as a function of the two covariates. 
Revised: June 1995

References:
S. Hunka. Using Mathematica to solve Johnson-Neyman Problems, Mathematica in Education, Vol. 3, No. 3, 32-36, Summer 1994.
S. Hunka. Identifying regions of significance in ANCOVA  problems having non-homogeneous regressions. Br. Jr. Math. Stat. Psych., 48, 161-188, 1995.
R.F. Potthoff. Johnson-Neyman Technique, in Encyclopedia of Statistical Sciences, Eds. S. Kotz and N.L. Johnson, 4, 299-303, John Wiley & Sons, New York, 1983.

Program Name: ancvjn[d_,ng_]
Calling arguments: d_ is the data matrix in which each row is an observation holding the dependent variable followed by one or two covariates. The rows are ordered by groups. ng_ is a vector in which each element is the number of observations in each group.
Functions called: getMp, cradd, ssk, knTrst, plotLines, plotRp, int3dplot,chkPts,  as described below.
Packages: Statistics`ContinuousDistributions`

Interactive Input: In the 1 covariate case the user inputs an alpha level for the desired region of significance and if a particular value of the covariate is to be tested for significance, the value to be used.  In the 2 covariate case the following are requested of the user:
1) an alpha value for the surface plot
2) the minimum and maximum values of each covariate for plotting purposes
3) the number of plot points
4) a vector of alpha values for the contour plot and whether shading is desired
5) to replot at a different of alpha using sequence starting at (1)
6) calculate significance of a covariate pair of values

Program Output:
1) the parameter estimates
2) a test of SSmcfm (at covariate values of zero) including contrast used
3) test of regression homogeneity for each covariate
4) symbolic contrast matrix to test SSmcfm for covariate adjustment to p and q
5) points of significance for user selected alpha values
6) plot of regression lines or planes
7) a plots of the polynomial or polynomial surfaces
8) contour plot of user selected regions of significance (2 covariate case only)
9) the polynomial representing SSmcfm for defining region of significance

Package calls, Functions, and Program follow:

<<Statistics`ContinuousDistributions`

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];
       ];

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] ]]]

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;k_ holds pattern e.g., {0,p,q}, and 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];
 ];

plotLines[p_,rng_,ttl_] plots regression lines in one plot. The arguments are defined as follows: p is a matrix in which each row represents the parameters of each line, i.e., {{c1,x1},{c2,x2}}, rng is a vector of 2 elements giving the range (min,max) for x over which each line is to be plotted, and ttl is a title for the plot. ttl can be defined as simply title="yourtitle", and then title is used in the argument list. Each line is colored differently.

plotLines[p_,rng_,ttl_]:=
 Module[{i,pobj,np,xv,ev},
   pobj=List[];
   np=Length[p]; (*# param/grp*)
   xv={1,x1}; (*symbolic const & predictor*)
   ex=p.xv;
   Do[
    AppendTo[pobj,
     Plot[ex[[i]],{x1,rng[[1]],rng[[2]]},
      DisplayFunction->Identity,
      PlotStyle->{RGBColor[Random[Real,{0,1}],
                           Random[Real,{0,1}],
                           Random[Real,{0,1}]]}]],
    {i,1,np}];
   Show[pobj,DisplayFunction->$DisplayFunction,
       AxesLabel->{"x","y"},
       PlotLabel->ttl
       ];
   ]; (*end Module*)

plotRp[p_,rng_,ttl_] plots regression planes in one 3d cube. The arguments are defined as follows: p is a matrix in which each row represents the parameters of each plane, i.e., {{c,x1,x2},{c,x1,x2}..{}}, and rng is a matrix of order (2,2) in which each row gives the range (min,max) for x1 and x2 over which each plane is to be plotted.ttl is a title,e.g, title="yourtitle".

plotRp[p_,rng_,ttl_]:=
 Module[{i,pobj,np,xv,ev},
   pobj=List[];
   np=Length[p]; (*# param/grp*)
   xv={1,x1,x2}; (*symbolic const & covars*)
   ex=p.xv;
   Do[
    AppendTo[pobj,
     Plot3D[ex[[i]],{x1,rng[[1]][[1]],rng[[1]][[2]]},
         {x2,rng[[2]][[1]],rng[[2]][[2]]},
      DisplayFunction->Identity]],
    {i,1,np}];
   Show[pobj,DisplayFunction->$DisplayFunction,
       AxesLabel->{"x1","x2","y"},
       PlotLabel->ttl
       ];
   ]; (*end Module*)

int3dplot[pmm_,ss_,mse_,dfc_,dfe_] is used to plot the expression ss having unknowns p and q. The argument pmm is a vector containing the minimum and maximum values of p and q respectively; ss is the expression; mse is means square error; dfc is degrees of freedom for contrast; dfe is error degrees of freedom. For each variable p and q (x1 and x2) the user is presented with the currently set minimum and maximum values and allowed to change them. The user could save the expression output by ancvjn2, and after defining the matrix for pmm, use the function to further explore the polynomial surface and regions of significance of ss-siglev independently of ancvjn.

int3dplot[pmm_,ss_,mse_,dfc_,dfe_]:=
 Module[{rpq,pp,vminmax,q1,mw,vw,par,q3,q2,q4,pobj,
 siglev,alf},
 rpq="y";
 pp=10;
 q4="Do you want to replot with new values? y/n";
 vminmax=pmm;
  q1="Type a value of to use.It currently is ";
  mw=Flatten[Table[{"minimum ","maximum "},{2}]];
  vw={"x1 ","x1 ","x2 ","x2 "};
  
 While[rpq=="y",
  alf=Input["Type alpha for surface plot, e.g.,.05"];
  Print["Alpha=",alf];
  siglev=mse dfc 
    Quantile[FRatioDistribution[dfc,dfe],1-alf];
  Print["Point of Signficance (c): MSe x dfc x F=",siglev];  
 par=List[];
   Do[
  q2=StringInsert[q1,vw[[i]],17];(*insert x1 or x2*)
  q2=StringInsert[q2,mw[[i]],8];(*insert min or max*)
  q2=StringJoin[q2,ToString[vminmax[[i]]]];
  par=Append[par,Input[q2]];
      ,{i,1,4}]; (*end Do *)
  vminmax=par;
  Print["Plot ranges are:",par];
  q3="How many PlotPoints do you want? Currently it is ";
  q3=StringJoin[q3,ToString[pp]];
  pp=Input[q3];
  Print["PlotPoints=",pp];
  pobj=Plot3D[ss-siglev,
    {p,par[[1]],par[[2]]},{q,par[[3]],par[[4]]},
        PlotPoints->pp,
        AxesLabel->{"x1","x2","ss-c"},
        PlotLabel->"Polynomial Surface of ss-c"];
        Show[pobj,PlotRange->{-1,0},
           PlotLabel->"Polynomial Surface at ss-c=0"];
    
   (* now do contour plot *)
   alf=Input["Type alpha for contours, e.g.,{.05,.01}"];
  Print["Alpha=",alf];
  siglev=mse dfc 
    Quantile[FRatioDistribution[dfc,dfe],1-alf];
  Print["Points of Signficance: MSe x dfc x F are"];
  Print[siglev];
   
   cs="y"==InputString["Contour shading? y/n"];
   ContourPlot[ss,
        {p,par[[1]],par[[2]]},{q,par[[3]],par[[4]]},
        Contours->siglev,
        ContourShading->cs,
        FrameLabel->{"x1","x2"},ColorFunction->Hue,
        PlotLabel->"Regions of Signficance"];
   
   rpq=InputString[q4]; (*quest to repeat *)
     ]; (*end While*)
     ]; (* end Module *)

chkPts[x,b,mse,dfe,ng] calculates the probability in testing SSmcfm 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. chkPts calls knTrst and ssk.

chkPts[x_,b_,mse_,dfe_,ng_]:=
Module[{cpat,rpq,con,ms,dfc,f,fdist,prob},
 rpq="y";
 Print[""];Print["Probability of Covariate set points"];
 While[rpq=="y",
   cpat=Join[{1},
   	Input["Type vector of covariate values,eg.{2,3}"]];
   con=knTrst[cpat,ng];
   Print["Contrast Matrix"];Print[con];
   dfc=Length[con];
   ms=ssk[con,x,b]/dfc; (*get Mean square *)
   f=ms/mse;
   Print["MS=",ms,",df=",dfc,",F=",f];
   fdist=FRatioDistribution[dfc,dfe];
   prob=1-CDF[fdist,f];
   Print["Probability=",prob];
   rpq=InputString["Do you have more values? y/n"];
   ];(*end While*)
   ];(*end Module*)

The Mathematica program.

ancvjn[d_,ng_]:= (*Johnson-Neyman Ancova*)
 Module[{nt,ngrp,x,y,ncv,xtx,dm,b,dfc,dfe,j,mse,con,
      concv,ss1,f1,fdist,prob,ss2,siglev,b1,mnx,dt,
      soln,pobj,title,alf},
  y=Flatten[Take[Transpose[d],1]];
  x=Drop[Transpose[d],1];
  nt=Length[d];(*N total*)
  ngrp=Length[ng];
  ncv=-1+Last[Dimensions[d]];(*# covariates *)
  x=Transpose[
   Partition[Join[Table[1,{nt}],Flatten[x]],nt]];
  
  (*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;
  b=N[Inverse[xtx].Transpose[dm].y];
  Print["Mu-Model Parameter Estimates"];
  Print[TableForm[b,TableSpacing->{0},
    TableHeadings->Automatic]];
  ssmcfm=(b.xtx.b)-(Apply[Plus,y]^2)/nt;
  dfm=Length[xtx]-1;Print[""];
  Print["SSmcfm=",ssmcfm," dfm=",dfm];
  dfe=Dimensions[dm].{1,-1};
  mse=((y.y)-b.Transpose[dm].y)/dfe;
  Print["Mean Square Error=",mse," dfe=",dfe];
  f1=(ssmcfm/dfm)/mse;Print["Fmodel=",f1];
  fdist=FRatioDistribution[dfm,dfe];
  prob=1-CDF[fdist,f1];
  Print["Probability=",prob];
  (*now plot the lines or planes*)
  dt=Transpose[d]; (*for convenience*)
  mnx=List[];
  Do[AppendTo[mnx,Min[dt[[i]]]];
     AppendTo[mnx,Max[dt[[i]]]],{i,2,Length[dt]}];
  Print["Covariate Range"];
  Print[mnx];  
  (*Plot regression lines or planes *)
  Which[ncv==1, title="Regression Lines";
         b1=Partition[Flatten[b],2];      
         plotLines[b1,mnx,title],
        ncv==2, title="Regression Planes";
         b1=Partition[Flatten[b],3];
         plotRp[b1,Table[mnx,{2}],title];,
        ncv>2, Print["#cov>2,so no plots"];
       ]; (*end Which*)
  
  (*now test each cov regr coeff*)
  con=Join[{1,0},Table[0,{ncv-1}]];
  Do[con=RotateRight[con,1];
    concv=knTrst[con,ngrp];
    Print["Contrast matrix covariate ",j];
    Print[TableForm[concv,TableSpacing->{0,4}]];
    ss1=ssk[concv,xtx,b];
    dfc=Length[concv];
    Print["SS covariate=",ss1," dfc=",dfc];
    f1=(ss1/dfc)/mse;
    Print["F1=",f1];
    Needs["Statistics`ContinuousDistributions`"];
    fdist=FRatioDistribution[dfc,dfe];
    prob=1-CDF[fdist,f1];
    Print["Probability=",prob];
  ,{j,1,ncv}];(*end do*)
 
  (*Now do the symbolic processing*)
  concv=knTrst[Join[{1},Take[{p,q},ncv]],ngrp];
  Print["Symbolic Contrast Matrix"];
  Print[TableForm[concv,TableSpacing->{0,4}]];
  ss2=Simplify[ssk[concv,xtx,b]]; (*symbolic sum of squares*)
 
  (*Do 2D plot & Solve, or 3D plots*)
  Which[ncv==1,
     alf=Input["Type the level of alpha, e.g.,.05"];
     Print["Alpha=",alf];
     siglev=mse dfc 
     Quantile[FRatioDistribution[dfc,dfe],1-alf];
     Print["Point of Signficance: MSe x dfc x F=",siglev];  
         Plot[ss2-siglev,{p,mnx[[1]],mnx[[2]]},
         PlotLabel->"Quadratic Plot",
         AxesLabel->{"x","ss-c"}];
         soln=Solve[ss2-siglev==0,p];
         Print["Solution to Quadratic is",soln];,
        ncv==2,
         int3dplot[mnx,ss2,mse,dfc,dfe];
        ]; (*end Which*)
        
   If["y"==InputString["Check covariate point sets? y/n"],
      chkPts[xtx,b,mse,dfe,ngrp]];
        
   Return[ss2];
       ]; (*end module *)

Data matrix: {y,cov1,cov2} one observation per row as shown below for 2 covariate case.

data={{7,4,5},{10,6,5},{9,8,9},{11,9,8},{14,12,11},
      {6,4,3},{7,6,7},{8,8,7},{8,10,11},{9,12,10},
      {6,4,3},{6,6,7},{5,7,6},{5,10,11},{4,13,12}}

ssjn=ancvjn[data,{5,5,5}] (*2 covariate case*)

Data matrix {y,cov1}; 1 covariate case

data2={{5,3},{4,3},{4,4},{2,4},{2,5},
       {3,2},{3,4},{3,5},{4,5},{4,7},
       {6,5},{4,3},{6,4},{4,5},{7,6}}

ancvjn[data2,{5,5,5}] (*one covariate case*)