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

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.

• For an interesting discussion of thorny issues surrounding this topic, see stats.stackexchange.com/questions/9334/… – rolando2 Oct 11 '12 at 23:23
• And don't forget to make sure the results are CrossValidated ... seriously :) – Bitwise Oct 12 '12 at 1:26

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.

• Is this method essentially attempting to find the functional form of each variable relative to the output (predictors vs prediction respectively?) The functional form could end up being y = a + b*log(x1) + x1^(1/2) Is the output of the model you quoted an equation of that form "g(E(Y))=β0 +f1(x1)+f2(x2)+⋯+fm(xm)" Where does it solve for f1(x1), f2(x2) etc? If the equation is y = a + f1(x1) + f2(x2) + ... fm(xm) Then is this similar to a regression model (multi variable linear regression)? – Pat Keough Oct 11 '12 at 23:50
• Question i+ii) It attempts to estimate an additive model for the (linked) mean. The additive model relates the predictors to the mean via smooth functions. They can be estimated in a parametric, semi- or non-parametric fashion (depending on specification) iii) what do you mean by "output of the model" iv) check the paper for the fitting algorithm v) it is regression model but linear only if the f_i were linear. – Momo Oct 12 '12 at 0:37
• I looked at rolando's suggestion and went to that link. I downloaded Eureqa and it seems to be the exact algorithm that is binding the best fit to my data. It seems extremely powerful in the sense that it does not bog your computer down and try to do 10000 calculations per second. The results so far are great.. it basically finds the g(E(Y)) that you suggested. Thanks Momo and Rolando. – Pat Keough Oct 12 '12 at 1:25