Orthogonal or Raw Polynomial regression?

Question

I am carrying out an analysis where I want to test for a relationship between two biological variables. The response is proportional (so I'm using a binomial GLM with weights) and the explanatory continuous. There are theoretical reasons to think the relationship between them should be non-linear, and specifically quadratic-like, with y increasing to a maximum and then declining, and I would like to:

a) Test whether my data supports a dome shape. b) Try and say something about the shape of this dome if so, as I have several datasets with similar domes and it would be nice to compare them (where their respective peaks are for example).

Because I have this specific idea about the shape, I am using a model with a second order polynomial, rather than a more flexible GAM.

My question is whether raw or orthogonal polynomials better suit this purpose. This question has been asked before, here, and the highest-rated response frames this as a trade-off between understanding the contribution of quadratic and linear components to the data (by using orthogonal polynomials) vs having coefficients you can actually understand and use (by using raw polynomials). I would like to know if there's a way I can have my cake and eat it - to clearly understand the contribution of linear and non-linear components and if the latter is present, its shape.

The options I have come up with are:

1. Using raw polynomials, do a model comparison between a model with just a linear x, and then with the added x^2, and take a significant improvement in model fit as an indication of my dome shape. (I don't see anyone suggesting that I should take the significance of the x^2 term itself as indicating its importance - perhaps because inference is tricky given that x and x^2 are confounded, is that right?). Then, plot the curve and calculate the maximum using the coefficients given by the model summary to tell me where the peak is.

2. As above, but replace the model comparison with a model with orthogonal polynomials, taking the significance of the quadratic term as an indication of shape, and calculating squared semipartial correlations for each term. Then, use this as a justification to run the model with raw polynomials and extract the coefficients (as hinted is okay in an answer to this question), calculate the maximum again.

3. Run the model with orthogonal polynomials, taking the significance of the quadratic term as an indication of shape. Then plot the model to show the shape, and talk qualitatively about it.

My feeling is that 1) does not necessarily provide a strong indication of the non-linear nature of the data, 2) is a bit mix and match, and feels like it complicates things, but clearly serves my purposes 3) is appealingly simple but it does not give me anything to describe the shape of the curve. In spite of there being plenty of questions on stack exchange about it, I am still really unclear on the meaning of the co-efficients generated by orthogonal polynomial regression.

What do people suggest, and why, and if so how would you report it? Is there an option that I have overlooked?

For reference, below is the code for the components of the steps described above, using some of my data.

Data

x <- c(0.000000, 3.162278, 4.472136,  6.324555,  8.366600, 10.488088, 13.038405, 15.165751, 17.029386, 18.708287)
trials <- c(22, 22, 22, 21, 23, 21, 20, 22, 21, 22)
success <- c(0.8181818, 0.8636364, 0.8636364, 0.9047619, 0.9565217, 0.9047619, 0.9000000, 0.6818182, 0.6190476, 0.7272727)

data.frame(x, success, trials) -> question_data

Models

Perform model selection to show that addition of second order polynomial improves model fit.

glm(success ~ poly(x, 1, raw = TRUE), family = binomial(link = "logit"), weights = trials, data = question_data) -> raw_poly_1
glm(success ~ poly(x, 2, raw = TRUE), family = binomial(link = "logit"), weights = trials, data = question_data) -> raw_poly_2

stats:::anova.glm(raw_poly_1, raw_poly_2, test = "Chisq")

Output:

Analysis of Deviance Table

Model 1: success ~ poly(x, 1, raw = TRUE)
Model 2: success ~ poly(x, 2, raw = TRUE)
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
1         8    10.3202                       
2         7     5.5498  1   4.7704  0.02895 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Plot the model with the second order polynomial to show structure

library(ggeffects)

ggpredict(raw_poly_2, terms = c("x")) %>% 
  plot(add.data = TRUE, jitter = FALSE) + labs(x = "\n x", y = "Proportion of Successes\n", title = NULL) +
  theme(axis.text.x = element_text(angle = 90, vjust= 0.5))

And calculate the maximum of the curve, based on the fact that it can be calculated via a quadratic equation formed from the model co-efficients.

a <- (-0.012839)
b <- 0.184377
c <- 1.513152

max <- (c - (b/(2*a)))

print(max)

Output:

[1] 8.693501

Run an orthogonal model to see whether the quadratic term is significant.

glm(success ~ poly(x, 2), family = binomial(link = "logit"), weights = trials, data = question_data) -> orth_poly

summary(orth_poly)

Output:

Call:
glm(formula = success ~ poly(x, 2), family = binomial(link = "logit"), 
    data = question_data, weights = trials)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.24589  -0.32418   0.05005   0.57034   1.10745  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.6408     0.1928   8.511   <2e-16 ***
poly(x, 2)1  -1.1888     0.5394  -2.204   0.0275 *  
poly(x, 2)2  -1.2529     0.5673  -2.209   0.0272 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 16.1744  on 9  degrees of freedom
Residual deviance:  5.5498  on 7  degrees of freedom
AIC: 40.037

Number of Fisher Scoring iterations: 4

Answer 1

The collinearity is mainly a concern when the peak (or minimum) does not occur within the range of your data. Since your example has a fair amount of data to either side of the maximum (and presumably your other data sets will be similar) this is not much of a concern and your final conclusions will be similar (possibly small round off errors) whether you use orthogonal or raw polynomials.

Answer 2

One further option is an inverse polynomial as described by John A. Nelder in https://www.jstor.org/stable/2528220

Consider the inverse quadratic $$x / y = \beta_0 + \beta_1 x + \beta_2 x^2$$ which means that for small $x$, $y \approx \beta_0^{-1} x$, and for large $x$, $y \approx (\beta_2 x)^{-1}$, which can be a hump with biologically plausible limits at either end of the range.

I used a generalized linear model in Stata with this syntax -- which didn't look too bad --

 glm success x xsq [fw=trials], link(power -1) f(binomial) vce(robust)

where xsq is x squared.

That said, the fit of the same model, except with log link, was very similar, and that with logit (which respects the bounds of the outcome in principle) was a bit different.

In short, I am pushing models that will in principle never predict negative outcomes.