Using Mathematica to Solve Johnson-Neyman Problems

by
Steve Hunka
Centre for Research and Applied Measurement and Evaluation
Faculty of Education
University of Alberta
Edmonton, Alberta
Canada
416-492-3762

Mathematica in Education
Vol.3 No.3
Summer 1994
(c) TELOS/Springer-Verlag

Function cradd does direct sum operation

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

Function ssk calculates sum of squares for contrast
matrix k

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

Function kt creates the K' contrast matrix

kt[k_]:=
Partition[
  Flatten[
    Join[{k,-k,Table[0,{Length[k]}]},
      RotateRight[{k,-k,Table[0,{Length[k]}]}]]]
       ,9]

Enter the data, g1 for group 1, g2 for group 2, etc. The leading
 1 is included in order to estimate the constant term, the covariates follow.

g1={{1,4,5},{1,6,5},{1,8,9},{1,9,8},{1,12,11}}
g2={{1,4,3},{1,6,7},{1,8,7},{1,10,11},{1,12,10}}
g3={{1,4,3},{1,6,7},{1,7,6},{1,10,11},{1,13,12}}

Use the direct sum operator to format the data for each group to form the design matrix

x=cradd[cradd[g1,g2],g3]

xtx=Transpose[x].x

Y is the dependent variable and we solve for B using ordinary
least squares approach

y={7,10,9,11,14,6,7,8,8,9,6,6,5,5,4}

b=N[Inverse[xtx].Transpose[x].y]

mse=((y.y)-b.Transpose[x].y)/Dimensions[x].{1,-1}

Test to see if there is homogeneity of regression estimates for the first covariate; one would do the same for the second covariate.

con={0,1,0}

conk1=kt[con]

ss1=ssk[conk1,xtx,b]

f1=(ss1/2)/mse

<<Statistics`ContinuousDistributions`

fdist=FRatioDistribution[2,6]

1-CDF[fdist,22.9648]

Set up the contrast matrix to test group effects using the unknowns
p and q

con={1,p,q}

conk2=kt[con]

ss2 gives us the algebraic expression for the sum of squares for the
contrast matrix containing p and q

Bug fix: 2.0 -> 2.2

Unprotect[Transpose];
Transpose[v_?VectorQ] := v
Protect[Transpose]

ss2=ssk[conk2,xtx,b]

Expand[ss2]

Below is Figure 1; this shows us what the sum of squares function looks like

Plot3D[ss2 - 4.284,
	{p, 0, 20}, {q, -5, 20},
	PlotPoints -> 30,
	AxesLabel -> 
		{FontForm["p", {"Times-BoldItalic", 10}], 
		 FontForm["q", {"Times-BoldItalic", 10}],
		 FontForm[ "ss  ", {"Times-BoldItalic", 10}]}]

Below is Figure 2; statistical significance is shown by the area bounded by the flat surface

Show[%, PlotRange -> {-1, 0}]

Below is Figure 3; this is the function from top

ContourPlot[ss2 - 4.284,
	{p,0,20}, {q,-5,20},
	FrameLabel -> 
		{FontForm["p", {"Times-BoldItalic", 10}],
		 FontForm["q", {"Times-BoldItalic", 10}]},
	PlotPoints -> 30]

Below is Figure 4; the outer boundary defines the region of statistical significance

Show[%, 
	PlotRange -> {0,3},
	Contours -> 3]