Identifying Polynomial Coefficients

[Michael Betancourt]

If a covariates is positive, x > 0, then the coefficients
for polynomials, say alpha1 x + alpha2 * x^{2}, become
strongly non-identified causing no end of fitting problems.
Unfortunately I’m trying to reproduce an existing analysis
that used such a polynomial so I can’t just get rid of the
higher-order terms.

Are there any tricks to decorrelate the coefficients by
transforming the covariates? The physicist in me wants
to do some kind of Gramm-Schmidt, but then that would
require a choice of kernel.

[Ben Goodrich]

I have anecdotally observed (mostly from rstanarm) that orthogonal polynomials yield better sampling for the coefficients than raw polynomials. You can just type

poly

in R to see its source code, which uses a QR decomposition in the default case.

Ben

[Michael Betancourt]

Right, that’s what Gramm-Schmidt would yield, provided you

define a specific kernel.  Not sure which kernel would be best,

though.

Or is it going to be easier to just build everything into a design

matrix and do the QR explicitly?

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ben Goodrich]

On Friday, April 8, 2016 at 11:32:59 AM UTC-4, Michael Betancourt wrote:

Right, that’s what Gramm-Schmidt would yield, provided you

define a specific kernel.  Not sure which kernel would be best,

though.

Or is it going to be easier to just build everything into a design

matrix and do the QR explicitly?

We have an option to do the QR of a design matrix in rstanarm, and it seems to work well (but we haven't studied it that systematically). Basically you end up doing a GLM on Q instead of X = QR and then pre-multiply the resulting coefficients by the inverse of R to obtain the coefficients with respect to X as a generated quantity.

Ben

[Michael Betancourt]

Right. So instead of

y = normal( X* beta + alpha, sigma);

you’d fit

Q = qr_Q(X);
R = qr_R(X);
y = normal(Q * beta_tilde + alpha, sigma);

then

beta = mdivide_right_tri_low(beta_tilde, R); ?

[Ben Goodrich]

On Friday, April 8, 2016 at 11:45:21 AM UTC-4, Michael Betancourt wrote:

Right.  So instead of

y = normal( X* beta + alpha, sigma);

you’d fit

Q = qr_Q(X);
R = qr_R(X);
y = normal(Q * beta_tilde + alpha, sigma);

then

beta = mdivide_right_tri_low(beta_tilde, R); ?

Yeah, basically. We actually multiply Q by sqrt(N - 1) and divide R by sqrt(N - 1) to make beta_tilde have a more reasonable scale in large datasets. A lot of times ADVI only works well with this transformation.

Ben

[Michael Betancourt]

The dimensions are screwing with me here.

Let’s say that we have N data points and M covariates, where

N >> M.  In this case our original problem is

matrix[N, M] X; // data

vector[M] beta; // parameters

vector[N] alpha; // parameters

y ~ normal(X * beta + alpha, sigma);

If the covariates are highly correlated then this can cause

problems, so instead we want to use a QR decomposition,

to write

X * beta = (Q * R) * beta = Q * (R * beta) = Q * beta_tilde

and then fit the model

matrix[N, M] X; // data

matrix[N, N] Q;

matrix[N, M] R;

vector[N] beta_tilde; // parameters

vector[N] alpha; // parameters

y ~ normal(Q * beta_tilde + alpha, sigma);

beta = R \ beta_tilde;

But doesn't this replace the M original parameters with N >> M

new parameters?  Is the cost of decor relating really the

introduction of so many new parameters?

In any case qr_R is seg faulting on me so it looks like there

are bigger problems...

[Ben Goodrich]

On Friday, April 15, 2016 at 5:44:04 AM UTC-4, Michael Betancourt wrote:

matrix[N, M] X; // data

matrix[N, N] Q;

matrix[N, M] R;

In the "thin" QR decomposition of the "tall" NxM matrix X = Q * R, Q is supposed to be NxM and R is supposed to be MxM.

https://en.wikipedia.org/wiki/QR_decomposition#Rectangular_matrix
http://eigen.tuxfamily.org/dox/classEigen_1_1HouseholderQR.html#a520fc76215d7722387a76e7b7a484ee0

Is Stan messing up? In rstanarm, we do the QR decomposition in R and don't pass in X.
 

vector[N] beta_tilde; // parameters

This should then be vector[M] beta_tilde.

Ben

[Michael Betancourt]

Eigen (and hence Stan’s) QR methods don’t return the thin decompositions

but rather the full Q and R matrices.  I’ll would play around with thinning

them out by hand, but Stan’s qr_R function seg faults preventing me from

using the decomposition at all.

[Bob Carpenter]

Michael --- would you mind filing a stan-dev/math issue with an example that
seg faults?

- Bob

[Michael Betancourt]

Done.