Statalist The Stata Listserver


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: trend test after logistic regression.


From   rgates@stata.com (Richard Gates)
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: trend test after logistic regression.
Date   Thu, 26 Jan 2006 07:29:50 -0600

I sent Ricardo <ovaldia@yahoo.com> some Mata code that 
computes the orthogonal polynomial codings for an integer 
variable and I thought that some people learning Mata might 
find it a useful.


-Rich
rgates@stata.com

. clear

. set obs 9
obs was 0, now 9

. gen int cat = mod(_n-1,3)+1

. mata:orthogonal_poly("cat")
  2
  3

. list

     +-------------------------------+
     | cat       _cat_1       _cat_2 |
     |-------------------------------|
  1. |   1   -.70710678    .40824829 |
  2. |   2            0   -.81649658 |
  3. |   3    .70710678    .40824829 |
  4. |   1   -.70710678    .40824829 |
  5. |   2            0   -.81649658 |
     |-------------------------------|
  6. |   3    .70710678    .40824829 |
  7. |   1   -.70710678    .40824829 |
  8. |   2            0   -.81649658 |
  9. |   3    .70710678    .40824829 |
     +-------------------------------+

. clear

. set obs 12
obs was 0, now 12

. gen int cat = mod(_n-1,4)+1

. mata:orthogonal_poly("cat")
  2
  3
  4

. list

     +----------------------------------------+
     | cat       _cat_1   _cat_2       _cat_3 |
     |----------------------------------------|
  1. |   1   -.67082039       .5    -.2236068 |
  2. |   2    -.2236068      -.5    .67082039 |
  3. |   3     .2236068      -.5   -.67082039 |
  4. |   4    .67082039       .5     .2236068 |
  5. |   1   -.67082039       .5    -.2236068 |
     |----------------------------------------|
  6. |   2    -.2236068      -.5    .67082039 |
  7. |   3     .2236068      -.5   -.67082039 |
  8. |   4    .67082039       .5     .2236068 |
  9. |   1   -.67082039       .5    -.2236068 |
 10. |   2    -.2236068      -.5    .67082039 |
     |----------------------------------------|
 11. |   3     .2236068      -.5   -.67082039 |
 12. |   4    .67082039       .5     .2236068 |
     +----------------------------------------+


mata

void function orthogonal_poly(string scalar svar) 
{
	real scalar m, ni
	string scalar names, namej
	real colvector v, lev
	real matrix X
	pointer (real colvector) colvector ind

	pragma unset v
	st_view(v,.,svar)
	lev = sort(uniqrows(v),1)
	m = length(lev)
	if (m<2) {
		errprintf("too few levels %d\n", m)
		exit(198)
	}
	if (m>100) {
		errprintf("too many levels %d\n", m)
		exit(198)
	}
	ind = J(m,1,NULL)

	X = orthogonal_poly_X(m)[.,2::m]
	for (i=1; i<=m; i++) {
		ind[i] = &orthogonal_poly_index(v,lev[i])
	}
	names = ""
	for (j=1; j<m; j++) {
		namej = "_"+svar+"_"+strofreal(j)
		stata("cap drop "+namej)
		st_addvar("double",namej)
		names = names + " " + namej
		for (i=1; i<=m; i++) {
			ni = length(*ind[i])
			if (ni > 0) st_store(*ind[i],namej,J(ni,1,X[i,j]))
		}
	}
	st_global("r(polyvars)",names)
}

/* generate orthogonal polynomials */
real matrix orthogonal_poly_X(real scalar n)
{
	real scalar i
	real colvector z
	real matrix X, Q, R

	z = (1::n):-mean((1::n),J(n,1,1))
	X = J(n,n,1)
	for (i=2; i<=n; i++) X[.,i] = X[.,i-1]:*z
	
	pragma unset Q
	pragma unset R	
	qrd(X,Q,R)
	X = Q*diag(R)
	
	X = X:/sqrt(colsum(X:^2))
	X = X:*(abs(X):>1.0e-15)
	return(X)
}

real colvector orthogonal_poly_index(real colvector v, real scalar lev)
{
	real scalar n, k, i, j
	real colvector ind, iv

	n = length(v)
	iv = (v:==lev)
	k = sum(iv)
	if (k == 0) return(J(0,1,.))

	ind = J(k,1,.)
	j = 0
	for (i=1; i<=n; i++) {
		if (iv[i]>0) ind[++j] = i
	}
	return(ind)
}

end 
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2021 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index