How to automatically test for the best parameters for transformed independent variable in linear model

Question

Let's assume that I have a linear model with $k$ variables:

$y = \beta_0 + \beta_1\cdot x_1 + \dots + \beta_k \cdot x_k$.

Now, I want to add variable $x_{k+1}$, but, according to domain knowledge, the dependency of $y$ onto $x_{k+1}$ is not linear, but rather "S-shaped". To capture this dependency, a well-established method is to use parametrised arcus tangens function: $D \cdot \arctan{\frac{x_{k+1} - A}{B}} + C$. To include it in the linear model we can safely ignore $D$ and $C$ parameters (as they will contribute to $\beta_{k+1}$ and $\beta_0$ respectively), so eventually I'd end up with the model:

$y = \beta_0 + \beta_1\cdot x_1 + \dots + \beta_k \cdot x_k + \beta_{k+1}\cdot \arctan{\frac{x_{k+1} - A}{B}}$.

What I would like to do is to find a way to test different transformations (with different $A$ and $B$ parameters automatically. Firstly, I went for searching the grid of different $A$ and $B$ values, i.e. fitting the model to different transformations of $x_{k+1}$ and returning list of models, sorted by $R^2$, $RMSE$ or $AIC$. However, this usually favours very sharp curves:

This doesn't make much sense in the eye test though. I know that eye-test may be biased but to me, more propable dependency would be captured by a shape as below:

One could already noticed that there are plenty of observations where $x_{k+1}$ is or is near to $0$, which is actually a characteristic of those variables that will be approximated by $arctan$. That may probably affect the fitting, but still something clearly doesn't work here.

Another approach that I considered is to fit the model using gradient descent search rather than using standard software (like lm from R). The problem is that with this $arctan$ transformation, we don't have a convex objective function anymore.

How can I then automatically test for the best parameters of the $arctan$ transformation? Am I on the right track with either idea or can I approach this problem differently?r

Answer 1

One approach would be to use an algorithm designed for non-convex problems like Bayesian optimization. However, if you have already evaluated a fine grid of parameters this is unlikely to offer significant improvement. Here is an example of how you could implement Bayesian optimization for this problem.

First, we need some data. Just for fun let’s extract the data from the images you posted (in brief since this is off topic).

In Mathematica:

img = Import[NotebookDirectory[]<>"LHrXQ.png"]
img2 = ImageResize[ImageTake[img, {40, 450}, {50, 1000}], 200]
pixels = {(#[[1]]-10)*3.8,#[[2]]*9+100}&/@PixelValuePositions[img2,Black, 0.4];

ListPlot[pixels, Frame->True, ImageSize->500]

Now in order to use Bayesian optimization we need to define the objective as a function of parameters A and B. Here we will maximize the R^2 value for the fitted model.

In python:

import pandas as pd
import numpy as np
import random
from sklearn.linear_model import LinearRegression

data = pd.read_csv('mathematica_data.csv')

def objective(params):
    """Whatever you want to do for your regression."""

    # So it works with GpyOpt
    A = params[0][0]
    B = params[0][1]

    temp = data.copy()

    # Transform variable
    xt = [np.arctan((x - A)/B) for x in data['x1'].tolist()]
    temp['x1'] = xt

    # Fit a linear model
    reg = LinearRegression().fit(temp.drop('y', axis=1), temp['y'])

    # Compute scores of interest
    r2 = reg.score(temp.drop('y', axis=1), temp['y'])

    # GPyOpt will minimize so we want -f
    return - r2

Now use GPyOpt to optimize the objective.

import GPyOpt

domain = [{'name': 'A', 'type': 'continuous', 'domain': (-300.0,300.0)},
          {'name': 'B', 'type': 'continuous', 'domain': (0.1,300.0)}]


bo = GPyOpt.methods.BayesianOptimization(f=objective,
                                         domain=domain,
                                         model_type='GP',
                                         acquisition_type='EI',
                                         initial_design_numdata=10,
                                         initial_design_type='random',
                                         acquisition_jitter=0.01,
                                         num_cores=-1,
                                         de_duplication=True,
                                         exact_feval=True)

# Run optimization
bo.run_optimization(max_iter=100)

We can plot the optimizers convergence:

bo.plot_convergence()

And information about how it sampled the parameters:

bo.plot_acquisition() 

Of course using only the independent variable extracted from your plot even the best R^2 values suggests there is little to no relationship between arctan((x-A)/B).

Answer 2

You can build an optimization problem.

Your features will be the parameters, then you create a loss function and try to minimize it with gradient descent or brute force (depends on the search space of your problem)

Answer 3

I'm not sure if you are bound to the type of model presented in your question. However, an alternative would be to use generalised additive models (GAM), e.g. with regression splines or locale regression. These methods usually give a very good fit with non-linear patterns in $X$ and there is no need to provide parameterization of $X$ so that it is easy to find good models.

Here is an example based on simulated data (R code):

# Generate data
x <- -50:100
y <- 0.001*x^3
plot(x,y)
df = data.frame(y,x)

# Linear regression
reg_ols=lm(y~.,data=df)
pred_ols = predict(reg_ols, newdata=df)

# GAM with regression splined (df=3)
library(gam)
reg_gam = gam(y~s(x,3), data=df)
pred_gam = predict(reg_gam, newdata=df)

# Find opt. number of splines
library(Metrics)
for (sp in seq(1:50)){
  gamx=gam(y~s(x,sp), data=df)
  print(mse(y, predict(gamx, newdata=df)))
}

# Plot prediction and actual data
require(ggplot2)
df2 = data.frame(x,y,pred_ols, pred_gam)
ggplot(df2, aes(x)) +                    
  geom_line(aes(y=y),size=1, colour="red") +  
  geom_line(aes(y=pred_ols),size=1, colour="blue") +
  geom_line(aes(y=pred_gam),size=1, colour="black", linetype = "dashed")

As you can see, the model gives a very good fit to my non-linear function while there is no need to provide parameterization (see figure). See ISL Chapter 7.7 for an applied introduction.

Answer 4

You could use Nonlinear Least Squares, in which one of the regressors is your arctan function with two more parameters to be estimated.

In R, for example:

library(minpack.lm)

df <- datasets::airquality

my_atan <- function(x, A, B){atan((x-A)/B)}

nlsLM(Ozone ~ a + b * Temp + c * my_atan(Temp, A, B),
      data = df,
      start = list(a = 0, b = 0, c = 0, A = 0, B = 1))