9.2 - Quantitative Factor Levels: Orthogonal Polynomials

Historically, when designed experiments involved quantitative factor levels, polynomial trends in the response were evaluated by using orthogonal polynomial contrast coding. These coefficients were used to partition the SS for the Factor into linear, quadratic, cubic, etc. contributions. In this section you can see how the orthogonal polynomial contrast coefficients are generated, and the Factor SS partitioned. This method will still be required to fit polynomial regression models with terms greater than the quadratic because even after centering there will still be multicollinearity between x and \(x_3\) as well as between \(x_2\) and \(x_4\).

The following example is taken from a textbook for ANOVA, Design of Experiments: Statistical Principles of Research Design and Analysis by Robert Kuehl.

Eample 9-1: Grain yeild Section

The treatment design consisted of five plant densities (10, 20, 30, 40, and 50). Each of the five treatments was assigned randomly to three field plots in a completely randomized experiment design. The resulting grain yields are shown in the table below (Grain Data):

  Plant Density (x)
  10 20 30 40 50
  12.2 16.0 18.6 17.6 18.0
  11.4 15.5 20.2 19.3 16.4
  12.4 16.5 18.2 17.1 16.6
Means (\(\bar{y}_i\)) 12.0 16.0 19.0 18.0 17.0

Display 3.5 - Transformation of Powers of x into Orthogonal Polynomials

Mean: \(P_0=1\)

Linear: \(P_1=\lambda_1 \left( \dfrac{x-\bar{x}}{d}\right)\)

Quadratic: \(P_2=\lambda_2 \left( \left(\dfrac{x-\bar{x}}{d}\right)^2-\left( \dfrac{t^2-1}{12} \right)\right)\)

Cubic: \(P_3=\lambda_3 \left( \left(\dfrac{x-\bar{x}}{d}\right)^3- \left(\dfrac{x-\bar{x}}{d}\right)\left( \dfrac{3t^2-7}{20} \right)\right)\)

Quartic: \(P_4 = \lambda_4 \left( \left( \dfrac{x - \bar{x}}{d} \right)^4 - \left( \dfrac{x - \bar{x}}{d} \right)^2 \left( \dfrac{3t^2 - 13}{14} \right)+ \dfrac{3 \left(t^2 - 1\right) \left(t^2 - 9 \right)}{560} \right)\)

where t = number of levels of the factor, x = value of the factor level, \(\bar{x}\) = mean of the factor levels, and d = distance between factor levels

Table 3.4 - Computations for orthogonal polynomial contrasts and sums of squares

    Orthogonal Polynomial Coefficients (\(P_{ci}\))
Density(x) \(\bar{y}_i\) Mean Linear Quadratic Cubic Quartic
10 12 1 -2 2 -1 1
20 16 1 -1 -1 2 -4
30 19 1 0 -2 0 6
40 18 1 1 -1 -2 -4
50 17 1 2 2 1 1
\(\lambda_c\) - 1 1 5/6 35/12
Sum = \(\sum P_{ci}\bar{y}_{i.}\) 82 12 -14 1 7
Divisor = \(\sum P_{ci}^{2}\) 5 10 14 10 70
\(SSP_c=r(\sum P_{ci}\bar{y}_{i.})^2/\sum P_{ci}^{2}\) - 43.2 42.0 0.3 2.1
\(\hat{\alpha}c=\sum P_{ci}\bar{y}_{i.}/\sum P_{ci}^{2}\) 16.4 1.2 -1.0 0.1 0.1

To generate the polynomial coefficients, first note that the five values of \(x\) are 10, 20, 30, 40, 50. Therefore, \(\bar{x} = 30\) and the spacing \(d = 10\). This means that the five values of \(\dfrac{x - \bar{x}}{d}\) are -2, -1, 0, 1 and 2.

Linear coefficients: The polynomial \(P_1\) for linear coefficients turn out to be:

\((-2)\lambda_1, (-1)\lambda_1, (0)\lambda_1, (1)\lambda_1, (2)\lambda_1\)

To obtain the final set of coefficients we choose \(\lambda_1\) so that the coefficients are integers. Therefore, we set \(\lambda_1 = 1\) and obtain the coefficient values in Table 3.4.

Quadratic coefficients: The polynomial \(P_2\) for linear coefficients turn out to be:

\(\begin{align*}
&\left( (-2)^2 - \left( \dfrac{5^2 - 1}{12}\right)\right)\lambda_2,
\left( (-1)^2 - \left( \dfrac{5^2 - 1}{12}\right)\right)\lambda_2,
\\& \left( (0)^2 - \left( \dfrac{5^2 - 1}{12}\right)\right)\lambda_2,
\left( (1)^2 - \left( \dfrac{5^2 - 1}{12}\right)\right)\lambda_2,
\\& \left( (2)^2 - \left( \dfrac{5^2 - 1}{12}\right)\right)\lambda_2
\end{align*}\)

which simplified gives:

\((2)\lambda_2, (-1)\lambda_2, (-2)\lambda_2, (-1)\lambda_2, (2)\lambda_2\)

To obtain the final set of coefficients we choose \(\lambda_2\) so that the coefficients are integers. Therefore, we set \(\lambda_2 = 1\) and obtain the coefficient values in Table 3.4.

Cubic coefficients: The polynomial \(P_3\) for cubic coefficients turn out to be:

\(\begin{align*}
&\left( (-2)^3 - (-2) \left( \dfrac{3(5^2) - 7}{20} \right)\right)\lambda_3,
\left( (-1)^3 - (-1) \left( \dfrac{3(5^2) - 7}{20} \right)\right)\lambda_3,
\\& \left( (0)^3 - (0) \left( \dfrac{3(5^2) - 7}{20} \right)\right)\lambda_3,
\left( (1)^3 - (1) \left( \dfrac{3(5^2) - 7}{20} \right)\right)\lambda_3,
\\& \left( (2)^3 - (2) \left( \dfrac{3(5^2) - 7}{20} \right)\right)\lambda_3
\end{align*}\)

which simplified gives:

\((-\dfrac{6}{5})\lambda_3, (\dfrac{12}{5})\lambda_3, (0)\lambda_3, (\dfrac{-12}{5})\lambda_3, (\dfrac{6}{5})\lambda_3\)

To obtain the final set of coefficients we choose $\lambda_3$ so that the coefficients are integers. Therefore, we set \(\lambda_3 = 5/6\) and obtain the coefficient values in Table 3.4.

Quartic coefficients: The polynomial $P_4$ can be used obtain the quartic coefficients in the same way as above.

Table 3.5 - Analysis of variance for the orthogonal polynomial model relationship between plant density and grain yield.

Source of Variation Degrees of Freedom Sum of Squares Mean Square F Pr > F
Density 4 87.60 21.90 29.28 .000
Error 10 7.48 0.75    
           
Contrast DF Contrast SS Mean Square F Pr > F
Linear 1 43.20 43.20 57.75 .000
Quadratic 1 42.00 42.00 56.15 .000
Cubic 1 .30 .30 .40 .541
Quartic 1 2.10 2.10 2.81 .125

The polynomial relationship expressed as a function of y and x in actual units of the observed variables is more informative than when expressed in units of the orthogional polynomial,. A direct transformation to an equation in x requires the information in Display 3.5 and Table 3.4. The necessary quantities are \(\lambda_1=1, d=10, \bar{x}=30\) and t = 5.Then, substituting for the \(P_i\)

\(\begin{align}
\hat{y}&=16.4+1.2P_i-1.0P_2\\
&=16.4+1.2(1)\left( \dfrac{x-30}{10} \right)-1.0(1)\left( \dfrac{x-30}{10}^2-\dfrac{5^2-1}{12} \right)\\
\end{align}\)

and simplifying to

\(\hat{y}=5.8+0.72x-0.01x^2\)

Generating Orthogonal Polynomials in SAS Section

Below is the code for generating polynomials from the IML procedure in SAS:

/* Generating Ortho_Polynomials from IML */
proc iml;
x={10 20 30 40 50}; 
xpoly=orpol(x,4); /* the '4' is the df for the quantitative factor */
density=x`; new=density || xpoly;
create out1 from new[colname={"density" "xp0" "xp1" "xp2" "xp3" "xp4"}];
append from new; close out1; 
quit;
proc print data=out1; run;


/* Here data is sorted and then merged with the original dataset */
proc sort data=grain; by density; run;
data ortho_poly; merge out1 grain; by density; run;
proc print data=ortho_poly; run;


/* The following code will then generate the results shown in the
	Online Lesson Notes for the Kuehl example data */
proc mixed data=ortho_poly method=type3;
class;
model yield=xp1 xp2 xp3 xp4;
title 'Using Orthog polynomials from IML';
run;
Analysis of Variance
Source DF Sum of Squares Mean Square ExpectedMean Square Error Term Error DF F Value Pr > F
xp1 1 43.200000 43.200000 Var(Residual) + Q(xp1) MS(Residual) 10 57.75 < .0001
xp2 1 42.000000 42.000000 Var(Residual) + Q(xp2) MS(Residual) 10 56.15 < .0001
xp3 1 0.300000 0.300000 Var(Residual) + Q(xp2) MS(Residual) 10 0.40 0.5407
xp4 1 2.100000 2.100000 Var(Residual) + Q(xp4) MS(Residual) 10 2.81 0.1248
Residual 10 7.480000 7.480000 Var(Residual)        

Fitting a Quadradic Model with Proc Mixed

Often we can see that only a quadratic curvature is of interest in a set of data. In this case, we can plan to simply run an order 2 (quadratic) polynomial and can easily use proc mixed (the general linear model). This method just requires centering the quantitative variable levels by subtracting the mean of the levels (30), and then creating the quadratic polynomial terms.

data grain; set grain;
x=density-30;
x2=x**2;
run;

proc mixed data=grain method=type3;
class;
model yield = x x2;
run;

The output is:

Type 3 Analysis of Variance
Source DF Sum of Squares Mean Square Expedted Mean Square Error Term Error DF F Value Pr > F
x 1 43.200000 43.200000 Var(residual) + Q(x) MS(residual) 12 52.47 <.0001
x2 1 42.000000 42.000000 Var(residual) + Q(x2) MS(residual) 12 51.01 <.0001
Residual 12 9.880000 0.823333 Var(Residual)        

We can also generate the solutions (coefficients) for the model with:

proc mixed data=grain method=type3;
class;
model yeild = x x2 / solution;
run;

which gives the following output for the regression coefficients:

Solution for Fixed Effects
Effect Estimate Standard Error DF t Value Pr > |t|
Intercept 18.4000 0.351 12 50.40 <.0001
x 0.1200 0.01657 12 7.24 <.0001
x2 -0.01000 0.001400 12 -7.14 <.0001

Here we need to keep in mind that the regression was based on centered values for the predictor, so we have to back-transform to put get the coefficients in terms of the original variables. This back-transform process (from Kutner et.al) is:

Regression Function in Terms of X. After a polynomial regression model has been developed, we often wish to express the final model in terms of the original variables rather than keeping it in terms of the centered variables. This can be dne readily. For example, the fitted second-order model for one predictor variable that is is expressed in terms of centered values \(x = X-\bar{X}\):

\(\hat{Y}=b_0+b_1x+b_{11}x^2\)

becomes in terms of the original X variable:

\(\hat{Y}=b'_0+b'_1X+b'_{11}X^2\)

where:

\(\begin{align}
b'_0&=b_0-b_1\bar{X}+b_{11}\bar{X}^2\\
b'_1&=b_1-2b_{11}\bar{X}\\
b'_{11}&=b_{11}\\
\end{align}\)

In the example above, this back-transformation uses the estimates from the Solutions for Fixed Effects table above.

data backtransform;
bprime0=18.4-(.12*30)+(-.01*(30**2));
bprime1=.12-(2*-.01*30);
bprime2=-.01;
title 'bprime0=b0-(b1*meanX)+(b2*(meanX)2)';
tite2 'bprime1=b1=2*b2*meanX';
title3 'bprime2=b2';
run;

proc print data=backtransform; var bprime0 bprime1 bprime2; run;

The output is then:

Obs bprime0 bprime1 bprime2
1 5.8 0.72 -0.01
Note! The ANOVA results and the final quadratic regression equation here are identical to the results from the orthogonal polynomial coding approach.