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:

enter image description here

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:

enter image description here

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, 
  • $\begingroup$ It's hard to see anything here. Your particular sigmoid doesn't track the data at all well, although at least some of the problem is that we can't see the other variables. But in principle nonlinear least squares seems a possible fitting method. I don't use R routinely but it is certainly supported there. $\endgroup$ – Nick Cox Jan 7 '20 at 10:55
  • $\begingroup$ I agree the presented variable case is not particularly obvious. But if it $arctan$ doesn't track the data well, we should just reject this variable and not include it in the model. Also I'm rather seeking general approach than "fitting" to this case, that's why I didn't include other variables as they are out of interest here. Yes, I considered non-linear least squares as well, but it was extremely hard to set starting parameters so that R nls function would find the solution. $\endgroup$ – jakes Jan 7 '20 at 11:15
  • $\begingroup$ Your question needs rewriting then. If you are focusing on the relationship with this predictor alone, the lead in with discussion of other variables is irrelevant. I have noted too with different software that you need good initial choices for nonlinear least squares. Four parameters is also a lot to play with. I don't understand why you start the curve at the origin. Why not post the data? $\endgroup$ – Nick Cox Jan 7 '20 at 11:25
  • $\begingroup$ It's not like the other variables are irrelevant. They consitute sort of base model, and testing for $x_{k-1}$ should be done on this base model so it's not like I am interested in univariate relationship between $y$ and $x_{k+1}$. Starting at the origin was just technical here, but I don't think it affects coefficients other than $\beta_0$, right? $\endgroup$ – jakes Jan 7 '20 at 11:32
  • 1
    $\begingroup$ Also, your argument about the intercept is wrong. If you add any other predictor to the model the intercept could, and usually will, change. The limiting case to consider is that without any predictors at all, the intercept is just the mean response. Saying that adding other predictors leaves it unchanged is not correct. If you try adding them one at a time you will find that the intercept jumps around and may even go outside the range of the response. Equally adding another predictor is likely to change the coefficients of predictors in the model. $\endgroup$ – Nick Cox Jan 7 '20 at 11:39

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:

Figure 2

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

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.

enter image description here

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.


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.