How to determine the best relationship (linear, log, etc.) between input predictor variable(s) and output variable for multiple linear regression?

Question

I am trying to determine the most accurate relationship between two variables (each predictor versus the output eventually). I want to know if the relationship is linear, or log-linear, or log-log, or some other form. What is the systematic approach to determining the underlying relationship that isn't full-out obvious and sometimes seens non-existant.

In some regression models I have seen, the author has had 10 predictors, and then changed the whole function from linear to log and ended up doubling the predictability (R^2 from ~40% to 80%). I'm fully aware that R^2 or adjusted R^2 isn't the best for determining the relationship but it is just for example.

Thanks in advance!

Answer 1

One possibility (although "best" while popular in questions here is always a bit vague...) would be using Generalized additive models (GAM) which are of the form

$g(\operatorname{E}(Y))=\beta_0 + f_1(x_1) + f_2(x_2)+ \cdots + f_m(x_m)$

with the definitions just as in Generalized Linear Models and $f_i(x_i)$ being functions estimated from the data. Basically the functional relationship of each predictors and the linked response get estimated.

In R you can use e.g., the gam, the mgcv and the gamlss packages to fit GAMs and variants.

An example would be to fit a GAM for the daily ozone measurements in New York, May to September 1973 explained by solar radiation, wind and temperature. Each predictor's functional relationship is estimated with nonparametric smoothing splines:

require(gam)
data(airquality)
mod1<-gam(Ozone^(1/3) ~ s(Solar.R) + s(Wind) + s(Temp), data=airquality,na=na.gam.replace)
summary(mod1)

Call: gam(formula = Ozone^(1/3) ~ s(Solar.R) + s(Wind) + s(Temp), data = airquality, 
    na.action = na.gam.replace)
Deviance Residuals:
    Min      1Q  Median      3Q     Max 
-1.1620 -0.2788 -0.0484  0.3321  1.2043 

(Dispersion Parameter for gaussian family taken to be 0.219)

    Null Deviance: 90.72 on 115 degrees of freedom
Residual Deviance: 22.52 on 103 degrees of freedom
AIC: 167 

Number of Local Scoring Iterations: 2 

DF for Terms and F-values for Nonparametric Effects

            Df Npar Df Npar F  Pr(F)   
(Intercept)  1                         
s(Solar.R)   1       3   1.60 0.1932   
s(Wind)      1       3   4.52 0.0051 **
s(Temp)      1       3   5.65 0.0013 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

But perhaps it's best to plot the estimated functions

par(mfrow=c(1,3))
plot(mod1,se=TRUE)

As you can see the functional relationship looks pretty different for all predictors and the functions are all non-linear. They are "best" in the sense of the fit criteria laid out in detail in e.g., the original paper.