Monday, September 19, 2011

How to get orthogonal polynomial coefficient/vector/codes

Tips - R & Stata & SAS: How to get orthogonal polynomial coefficient/vector/codes

When we do "contrast {lvl #1 #2 #3}" for trend analyses using Stata or other software for unequally spaced levels/categories, we need the orthogonal polynomial coefficient (#1 #2 #3), which is hard to be find in books. We can get these coefficients using R, Stata, or SAS. My favorite software for this purpose is R. Below I show the examples using these different kinds of software. Note: Stata has an operator (p. for orthogonal polynomial in the level values) for unequally spaced levels, for example, "contrast p.lvl".

R
  mostly I use R to get these coefficients:
  >cntr<-poly(c(1,2,5,6),3)
  >cntr

Stata
Step 1: create a dataset with one variable [lvl]:
  .input lvl
  1. 1
  2. 2
  3. 5
  4. 6
  5. end
Step 2-b: use 'orthpoly'
  .orthpoly lvl, generate(cntr1 cntr2 cntr3) degree(3)

  Now, in the dataset, you can find three new variables 'cntr1' for the orthogonal polynomial coefficients of degree 1 (linear), and 'cntr2' for the orthogonal polynomial coefficients of degree 2 (quadratic), and 'cntr3' for the orthogonal polynomial coefficients of degree 3 (cubic).

SAS
  PROC IML;
   lvl = {1 3 5 6}
   cntrl=ORPOL(lvl);
   PRINT cntrl;
  QUIT;

No comments: