What is the best linear regression model to use when the shape of the data resembles a cubic distribution?

Question

I have a data set with a distribution of one variable against the other resembling a cubic one (rises to some point and then falls to a steady level without a consequent rise). I know in which cases to use log-linear, log-lin, lin-log, and reciprocal or log reciprocal linear models, but I am not sure what to do here (I have checked all of the above and they not surprisingly turned out to be a bad fit). Is there any linear model that would help me in this case?

Answer 1

Restricted cubic splines (natural splines) are an excellent choice. These are piecewise cubic polynomials that can fit any shape given enough knots. The following code in R shows how to fit such relationships and to plot the fit with confidence bands.

require(rms)
dd <- datadist(mydata); options(datadist='dd')
f <- ols(y ~ rcs(x1, 5))  # 5 knots at default locations
f   # print model stats
plot(Predict(f))  # or plot(Predict(f, x1)) # plots over 10th smallest to 10th largest observation

Answer 2

I would have thought a "cubic regression" would work well for a cubic relationship. Call $Y_{i}$ the dependent variable, and $X_{i}$ the independent variable (or regressor). You simply use a polynomial regression:

$$Y_{i}=\left(\sum_{j=0}^{p}\beta_{j}X_{i}^{j}\right)+e_{i}$$

I would use BIC to select the value of $p$. To do this is very easy - calculate the coefficient of determination $R_{p}^{2}$ from a standard OLS regression output. Then a convenient form of BIC is given by:

$$BIC_{p}=n\log(1-R_{p}^{2})+p\log(n)$$

Although this is the standard form, with the natural logarithm's, a more convenient numerical form is given by $$BIC10_{p}=-\frac{1}{2}\log_{10}(e)BIC_{p}$$

The reason I say this is that in this form above, you get BIC expressed in based 10 log units, and this leads to a very quick interpretation of the actual number of the BIC. If BIC is positive, then the current order $p$ polynomial is more supported by the data (compared to intercept only model), and the numerical value in odds form is $10^{BIC10_{p}}$. So if $BIC10_{p}=1$, then the order $p$ polynomial is 10 times more likely than the intercept only model, if $BIC10_{p}=10$ then the order $p$ polynomial is 10 billion times more likely. BIC10 tells you how many digits are in the odds ratio. So a reasonable way to proceed is to continue to increase the order of a polynomial until $BIC10_{p}$ becomes sufficiently large.

One thing to be careful of though, is that this type of procedure is not likely to work well for extrapolation outside the range of the $X_{i}$ values. This is mainly because this is a data driven procedure.