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

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

• Did you check generalised additive models? They are extremely flexible in finding a good parameterisation more or less endogenously... so no need for manual feature engineering. – Peter Jan 7 '20 at 12:35

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

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).

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)

• I actually mentioned that at the end, but please note that this is not a convex problem anymore due to partial derivatives of $arctan$ – jakes Jan 7 '20 at 18:13
• Do you care about optimality? Just do a brute force search across all values and if not possible, maybe genetic algorithm? – Carlos Mougan Jan 7 '20 at 18:27
• It doesn't have to be perfect optimality. But I've been doing grid search as mentioned in the post and the result were clearly wrong. I am not sure, but I thing this is probably what you call brute force search? – jakes Jan 7 '20 at 19:42

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.

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))