# Polynomials in a regression model (Bayesian hierarchical model)

I am not a trained statistician and am looking to get some clarification of a model from literature. The study in question is A Hierarchical Framework for Correcting Under-Reporting in Count Data. The model as defined by equations 11 to 14 (with subscripts, nonrelevant terms removed for easier interpretation): \begin{align} z_{t} \mid y_{t} &\sim \operatorname{Binomial}\left(\pi, y_t \right) \\ \log \left(\frac{\pi}{1-\pi}\right)&=\beta_{0}+g\left(u\right) \\ y_{t} &\sim \operatorname{Poisson}\left(\lambda_{t}\right) \\ \log \left(\lambda_{t}\right) &=\log \left(P_{t, s}\right)+a_{0}+f_{1}\left(x_{s}^{(1)}\right)+f_{2}\left(x_{s}^{(2)}\right) \\ &+f_{3}\left(x_{s}^{(3)}\right)+f_{4}\left(x_{s}^{(4)}\right) \end{align}

where $$z_t$$ are observed counts and $$y_t$$ are real, true counts. And the functions $$g, f_1, \ldots, f_4(\cdot)$$ are (from the paper)

orthogonal polynomials of degrees 3,2,2,2, Compared to raw polynomials, these reduce multiple-collinearity between the monomial terms (Kennedy and Gentle 1980), and were set up using the “poly” function in R

From my understanding, this model first estimates the true count $$y_t$$. The true count itself depends on a logistic regression formula where the covariates are population, and social indicators such as $$x_s^{(1)} =$$ unemployment. The covariates are used as input to orthogonal functions. Once it estimates the true count, it uses that value in a Binomial model to count the number of "successes", i.e the observed count. The probability of success in this case is given by another regression formula that also has an orthogonal function for the covariate.

My questions are rather simple:

1. What is so important about using orthogonal functions in the regression models. Why can't simple coefficients be used (and these coefficients estimated in the Bayesian implementation).

2. The interpretation of the log of $$\pi$$ and $$\lambda$$. For $$\pi$$, I am guessing, the regression formula can evaluate to numbers outside of (0, 1), so the the ilogit will transform it between 0, 1. I don't understand why the log is taking for $$\lambda$$.

Let's handle 2. first.

As you guessed, the logit transformation of $$\pi$$ is designed so that the regression formula has no restriction on its values; any value will be mapped into $$(0,1)$$. The same is true for the log transformation of $$\lambda$$: $$\lambda$$ must be positive, and using log transformation allows the regression formula to take any value, positive or negative.

The log part of both transformations also means we get a multiplicative model rather than an additive, which often makes more sense for counts and proportions.

And, on top of all that, there are mathematical reasons that these transformations for these particular distributions lead to slightly tidier computation and are the defaults, though that shouldn't be very important reason.

Now for the orthogonal functions. These aren't saying $$f_1$$ is orthogonal to $$f_2$$; that's up to the data to decide. They are saying that $$f_1$$ is a quadratic polynomial in $$x^{(1)}$$, and that it's implemented as a weighted sum of orthogonal terms rather than a weighted sum of $$x$$, $$x^2$$. What the orthogonal polynomials actually are depends on the data, but let's pretend the data are evenly spaced on $$[-1,1]$$ and they're the Chebyshev polynomials $$T_0(x)=1,\, T_1(x)=x,\, T_2(x)=2x^2-1,\, T_3(x)=4x^3-3x$$.

If we were just doing maximum likelihood this wouldn't matter at all. Suppose the ML estimate based on the powers of $$x$$ was $$-0.1+2.7x-3x^2+4.5x^3$$. We can rewrite this in terms of the orthogonal polynomials: clearly the coefficient of $$T_3$$ has to be 4.5/4 to make the $$x^3$$ match, and the rest will take calculation. It turns out to be $$-1.6T_0+6.075T_1-1.5T_2+1.125T_3$$. These are the same polynomial, it's just a different way of writing the same model, and in this case (and nearly always with modern computers) the collinearity isn't anywhere near strong enough to cause numerical rounding problems.

With Bayesian inference, though, there's the question of priors. It makes more sense to put independent priors ($$\alpha_j$$ and $$\beta_k$$ in the paper) on the coefficients of orthogonal polynomials than to put independent priors on the coefficients of $$x$$, $$x^2$$, $$x^3$$. So, my assumption is that the orthogonal polynomials were chosen so that the relatively flat ($$N(0,10^2)$$) independent priors on their coefficients made sense.

• Thanks! So what would the structure of $f_1$ look like? I also don't see why its irrelevant for ML? Is it correct in saying that using these functions is basically capturing the effects of the covariates, while not worrying/caring about the coefficient parameters as we don't really need them? – masfenix Aug 17 at 22:43
• (Updated). It's a bit stronger than that, but yes. The reason it's stronger is that you can convert from one form to the other and work out the corresponding coefficients if you do care about them. – Thomas Lumley Aug 17 at 23:29
• Okay, very helpful (Will accept the answer). Last question though is how the polynomials are determined? Are you familiar with the R package poly`? It seems like in the code, they simply provide it with the vector of covariates. Is there a general formula/structure to these polynomials? Also, is there a wiki page or a pdf I can read on when/how to use the polynomials... I have a similar problem. – masfenix Aug 17 at 23:32
• The code inside R is here: svn.r-project.org/R/trunk/src/library/stats/R/contr.poly.R. It takes the matrix with columns $(x, x^2,x^3)$ and makes the columns orthogonal. A previous question stackoverflow.com/questions/39031172/… gives more detail – Thomas Lumley Aug 17 at 23:46