How to 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?

@Edit: As requested, the data underlying the plot is presented below. However, the plots here are just for visual purposes and I don't want this question to be focused on this concrete data example.

structure(list(y = c(413.177956367559, 336.940568728554, 271.423795696105, 
241.781630437832, 389.295192607088, 485.838148155368, 427.90090251737, 
425.318347912293, 646.237861173935, 466.654590750197, 381.459607936796, 
379.663723975475, 285.493647486784, 473.661534824836, 587.231761325654, 
422.109479708714, 366.106914377068, 310.694550139197, 615.192496132627, 
460.83370904994, 494.628049500221, 397.255757535097, 566.843613053996, 
491.613547007966, 423.409959844659, 436.515303465714, 522.414221946626, 
458.65313447119, 447.456412546129, 363.0944247836, 420.433851503218, 
437.447927051515, 380.458699893226, 380.889241204271, 449.452597470113, 
785.354871000865, 649.307156490243, 701.38941811651, 620.810032925486, 
549.994406432794, 536.476816637358, 532.827017680477, 534.569029401081, 
622.778124246893, 724.021503453207, 854.31391696348, 610.064721258919, 
578.431363007429, 661.311883547075, 663.490971124295, 542.129859228126, 
170.964424841894, 628.744421037317, 943.235729569169, 709.588445205357, 
711.43679300902, 700.552248512387, 608.720943212614, 597.994348235098, 
527.075360298016, 642.884851825923, 635.695319226458, 624.120362301625, 
528.728589597031, 456.286681464807, 697.140660423864, 895.989738745979, 
560.108415214582, 563.561490631544, 477.14359103754, 722.11913919209, 
703.691239904751, 601.026518890877, 670.386746789934, 611.816744597946, 
615.423696704836, 405.124923085792, 375.215242490828, 735.018295944114, 
614.919204190096, 630.10627231055, 590.120497911762, 589.052605761616, 
623.901396730888, 561.173063948177, 650.609184533618, 771.331844237992, 
634.092427568623, 257.698979379167, 597.685952911712, 712.852775022346, 
128.047974331485), x = c(0, 0, 0, 0, 132.32, 308.47, 400.02, 
440.5, 332.52, 166.26, 83.13, 41.57, 20.78, 276.58, 389.72, 461.04, 
230.52, 115.26, 393.34, 531.86, 601.63, 300.82, 485.6, 242.8, 
121.4, 60.7, 30.35, 15.17, 7.59, 3.79, 1.9, 0, 0, 0, 0, 337.78, 
284.58, 310.78, 439.59, 548.97, 565.67, 386.99, 193.5, 96.75, 
48.37, 24.19, 12.09, 6.05, 3.02, 1.51, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 344.58, 518.59, 601.11, 639.44, 
656.42, 662.32, 331.16, 165.58, 82.79, 41.4, 20.7, 10.35, 5.17, 
2.59, 1.29, 0, 0, 242.51, 121.26, 60.63, 415.23, 256.85)), class = "data.frame", row.names = c(NA, 
-92L))

Answer 1

This is a standard least squares problem and, as such, is solved by finding parameters that minimize the sum of squared residuals.

A general approach is first to fit the model

$$E[Y] = \beta_0 + \beta_1 x_1 + \cdots + \beta_k x_k.$$

Replace the $y_i$ by their residuals

$$e_i = y_i - (\hat\beta_0 + \hat\beta_1 x_1 + \cdots + \hat\beta_k x_k).$$

The nonlinear part of the model can be written as a function

$$f(x_{k+1}, (A,B,\beta_{k+1})) = \beta_{k+1} \arctan\left(B(x_{k+1} - A)\right).$$

Now estimate how the residuals are related to $x_{k+1}$ by minimizing the expression

$$f_{\text{obj}}(A, B, \beta_{k+1}) = f_{\text{obj}}(\mathbf{\theta}) = \frac{1}{n}\sum_{i=1}^n \left(e_i - f(x_{k+1}, \theta)\right)^2.\tag{*}$$

(To avoid potential problems with division by zero during the solution, I have replaced $B$ by $1/B$ in this formulation.)

Rather than re-inventing an efficient, accurate, stable, test method to solve $(*),$ employ one that is available, such as the nlm function in R. Here is the R code to compute this fit, beginning with a data frame X containing variables x and y and an arbitrary universal starting value $\theta_0 = (0,1,0)$ for the search:

f <- function(x, theta) theta[3] * atan(theta[2] * (x - theta[1]))
obj <- function(theta, y, x) mean((residuals(lm(y ~ 1)) - f(x, theta))^2)
fit <- with(X, nlm(obj, c(a=0, b=1, beta=0), y=y, x=x))

To demonstrate this, I simulated eight datasets like yours with $\theta=(200, 1/100, 50),$ $k=0,$ $\beta_0=550,$ and iid Normal errors of variance $\sigma^2=50^2$ (the code appears at the end of this post). Here they are along with the true model (gray line) and least-squares fits (red curve):

The main lesson is that this model has a tendency to select large values of $B$ resulting in abrupt jumps in the fit. That's not a problem with the solution or even with the model; it's a manifestation of the limited amount of data and the relatively large value of $\sigma$ they exhibit. With, say, ten times as much data, here is what a typical fit would look like:

It's good: the graph of the fit is scarcely distinguishable from the true model. The parameter estimate, however, is $\hat\theta = (198, 0.0012, 253),$ far from the true value $\theta = (200, 0.01, 50).$ This reveals a near lack of identifiability, making it difficult to pin down all three parameters. The problem is that you can achieve very nearly the same model by rescaling $B$ (to some extent) and applying a compensating change in $\beta_{k+1}.$ Geometrically, you are zooming into the graph of the arctangent function near $x=0$ which is nearly linear over large intervals surrounding that point, whence the change of scale changes nothing (to first order).

One way to resolve this problem is to specify either $B$ or $\beta_{k+1}$ and fit the remaining parameters. I recommend specifying $B,$ for then you can use standard tools like ridge regression to regularize $\beta_{k+1},$ thereby avoiding the sharp jumps in the solution. (You can use regularization in any event, but it will be easier to apply it to a linear parameter like $\beta_{k+1}$ than to the nonlinear parameters $A$ and $B.$)

You can obtain standard errors of estimate (and even a full covariance matrix) in the usual ways, because this approach yields the Maximum Likelihood estimates of all the coefficients.

---

Here is the R code used for simulation and plotting the first figure.

f <- function(x, theta) theta[3] * atan(theta[2] * (x - theta[1]))
obj <- function(theta, y, x) mean((residuals(lm(y ~ 1)) - f(x, theta))^2)
#
# Create constant values for repeated simulations.
#
x <- c(seq(0, 10, length.out=39), seq(0, 700, length.out=92-39))
x. <- seq(min(X$x), max(X$x), length.out=101)
theta <- c(A=200, `1/B`=1/100, beta=50)
mu <- 550
#
# Simulate.
#
set.seed(17)
sim <- lapply(1:8, function(i) {
  X <- data.frame(Iteration=i, x = x, 
                  y = mu + f(x, theta) + rnorm(length(x), 0, 150))
  fit <- with(X, nlm(obj, c(a=0, b=1, beta=0), y=y, x=x))
  Y <- with(fit, data.frame(Iteration=i, x = x., y = mean(X$y) + f(x., estimate)))
  list(Data=X, Fit=Y)
})
#
# Combine all simulated data into a data frame `X` and all fitted data into `Y`.
#
X <- do.call(rbind, lapply(sim, `[[`, i=1))
Y <- do.call(rbind, lapply(sim, `[[`, i=2))
Z <- data.frame(x = x., y = mu + f(x., theta))
#
# Plot everything.
#
library(ggplot2)
ggplot(X, aes(x, y)) + 
  geom_line(data=Z, size=1, color="Black", linetype=1) + 
  geom_point(alpha=0.5, shape=21, fill="Gray") +
  geom_line(data=Y, size=1, color="#e01010") + 
  xlab(expression(x[1])) + ylab(expression(y)) +
  facet_wrap(~ Iteration, nrow=2)

Answer 2

Thanks for the data. Although this is just your graph with a smoother added, I suggest that -- beyond eyeballing the data -- the graph isn't supportive of any idea that there is an underlying sigmoid. There is nothing magic about the bandwidth used, and I tried lower and higher values with the same conclusion.

Naturally I can't exclude the possibility that a pattern would emerge if the other predictors were taken into consideration, but on that front the only suggestions that I can make, as in comments, are

1. nonlinear least squares

2. using $Xb$ for the other predictors and then looking at the residuals from that.

I have to be puzzled at the idea that $Xb$ should be fine for 13 predictors, and this one is different. In my field with any bundle of 13 predictors I would expect that some would need transformation and I might need to think about interaction terms too.

I don't know what these data are, or what literature you're drawing on, but quite often something a little exotic -- here the use of an arctangent -- somehow acquires the status of what people use without a very strong rationale. To me an arctangent is tied up with the idea that the response is bounded in principle, which isn't obviously the case here. I have a prejudice against functions with more than about 3 parameters any way: not only are they often freakishly hard to fit, they are suspect without an independent rationale.

1 intercept, 13 vanilla coefficients and 4 extra (or even 2) coefficients are in total 18 (or 16) parameters; many would find 92 uncomfortably few data points for that exercise.

This is not very constructive, or even at all, but I tried. I also appreciate that you would rather discussion didn't focus on these data alone, but with no other information on the other variables, only broad comments seem possible.

Answer 3

A possible mistake that you're doing is not displaying the effect of $x_{k+1}$ after the effect of all other $x_i$. For this, rather than directly plotting $y$ against $x_{k+1}$ you should Partial Regression Plots (aka added variable plots; Wiki, or here). The three steps are (from Wiki):

1. Computing the residuals of regressing the response variable against the independent variables but omitting $X_{k+1}$

2. Computing the residuals from regressing $X_{k+1} against the remaining independent variables

3. Plotting the residuals from (1) against the residuals from (2).

Maybe plotting the parameters found by your gridsearch in a added variable plot will help you understand the found parameters.