I am using Multiple Linear Regression to assess the impact of two predictors on Y and, especially, whether during a certain time there is an impact on Y that cannot be explained by X1 and X2.

  • Y, X1 and X2 are physical measurements.
  • Season is a variable for the Seasonality of Y, also derived from physical measurements.
  • Xbinary is a binary variable that is 1 for the specific time I am interested in and 0 for all other times.

The timing of Xbinary is derived from a theoretical assumption. My assumption is that during Xbinary=1 something is happening to Y that we are not measuring.

The relationship of my predictors and Y look like this:

enter image description here

The Multiple Linear Model looks fine:

lm(formula = "Y ~ I(Season^3) + X1 + I(X1^2) + X2 + I(X2^3) + Xbinary", 
    data = df)

    Min      1Q  Median      3Q     Max 
-8.5271 -1.2125 -0.0326  1.1987  5.7848 

              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -8.239e+00  4.730e-01 -17.418  < 2e-16 ***
I(Season^3)  4.992e+00  1.792e-01  27.858  < 2e-16 ***
X1           1.935e-02  1.435e-03  13.485  < 2e-16 ***
I(X1^2)     -1.554e-05  1.947e-06  -7.983 3.15e-15 ***
X2           6.238e-01  3.830e-02  16.290  < 2e-16 ***
I(X2^3)     -4.925e-04  3.448e-05 -14.284  < 2e-16 ***
Xbinary     -1.704e+00  1.415e-01 -12.037  < 2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.871 on 1275 degrees of freedom
Multiple R-squared:  0.7645,    Adjusted R-squared:  0.7634 
F-statistic: 689.8 on 6 and 1275 DF,  p-value: < 2.2e-16

But this looks troubleful:

enter image description here enter image description here enter image description here

As far as I understand, the assumption of Heteroscedasticity and Linearity are violated. As far as I understand, the first causes lower p values and imprecise coefficients, the second causes wrong coefficients?

As the data are real-world physical measurements, the true relationship is very likely to be non-linear. Still, I chose the linear regression because the coefficients are easy to interpret. I know my model is wrong, but does it mean I cannot draw any conclusion from it? Can I interpret it as a rough estimate of the real relationship? I know the assumptions are violated, but it does not look completely wrong to me.

Here is an interesting side note: As I know that the linear model is not well describing the true relationship, I put the same data into a random forest and extreme gradient boosting. I used Xbinary=0 as training data and Xbinary=1 as test data. The average predicition error of both machine learning models is roughly the coefficient of Xbinary from the linear model. That cannot be complete coincidence, can it?

  • 5
    $\begingroup$ Your model assumes very specific, restricted functional forms for the relationships between your outcome and the continuous predictors. Have you tried letting the data tell you the correct functional form by fitting with regression splines? That's still a regression model linear in coefficients, but it allows for flexible fitting of nonlinear associations between outcomes and predictors. That's nicely implemented with the rcs() function in the R rms package, which also provides useful tools for validation and display of results. $\endgroup$
    – EdM
    Commented Jan 12, 2022 at 18:48
  • $\begingroup$ @EdM you are absolutely right. When I did my research for the topic I saw the rms package before but at a first glance it looked complicated to me and I wanted to have something easy to interpret before I go into more sophisticated approaches. Now I found out how to use rms and .. I am completely flashed by how easy and informative the effect plots are. Interestingly, the coefficients from the linear model are quite close to the ols results using rcs. So, the model wasn't completely bad, but looking at the effect plots from rms tells me the real story. $\endgroup$
    – Felix Phl
    Commented Jan 13, 2022 at 10:05
  • $\begingroup$ Do you have any idea why the variation is so small as function of 'season' except when the value of season is equal to 1? Do you have any theoretical ideas about the particular shape of that curve for season? What does the data mean and how can they be connected, isn't there any better model or at least more reasonable model than a straight line? $\endgroup$ Commented Jan 13, 2022 at 11:14
  • $\begingroup$ @SextusEmpiricus image the phenological cycle of a plant, at the beginning of the vegetation season leafs will grow and at the end leafs will fall off. During that time Y is mostly dominated by this effect, during the rest of the season not. Season is meant to describe this relationship. $\endgroup$
    – Felix Phl
    Commented Jan 13, 2022 at 11:51
  • $\begingroup$ @FelixPhl is 'season' a time variable? What does the value mean? $\endgroup$ Commented Jan 13, 2022 at 12:34

4 Answers 4


A few things to start with. I agree that nature is rarely linear, but I disagree that the approach you've taken is a reasonable attempt to account for this. Polynomial terms do add non-linearity, but they are high bias models (because we can only estimate functions in the set of polynomials). A better approach to adjusting for non-linearity would be use a Spline which can a) account for non-linear effects, and b) can be estimated through a linear model like linear regression.

As to your titular question, we my personal response would be "It Depends" on what you want to do with the model. The residuals vs predicted plot clearly shows there is residual confounding. This means your estimates of the conditional mean will be biased in some places and not in others.

  • $\begingroup$ I do not understand b). Polynomial models can also be estimated through a linear model, just like splines. In fact, for splines this is only true as long as the nodes are pre-specified (or chosen from the data before the fitting). $\endgroup$
    – jarauh
    Commented Jan 13, 2022 at 16:25
  • $\begingroup$ @jarauh Polynomials are high bias. If I use a quadratic term effect I can only estimate quadratic functions. If the data generating process is not even approximately quadratic, then bias remains. This assumes we can know a priori that a quadratic effect is appropriate. Splines are superior because they have the ability to estimate a wide range of functions and need not specify anything (except the knot locations, for which there is some research to help) prior to fitting the model. $\endgroup$ Commented Jan 13, 2022 at 16:29
  • $\begingroup$ That was not my question. You claim that splines have the advantage that they "can be estimated through a linear model [...]". Just wanted to say that the same is true for polynomials. $\endgroup$
    – jarauh
    Commented Jan 13, 2022 at 16:43
  • $\begingroup$ @Demetri Pananos That is not true about splines. A spline requires that you provide an order to fit the nodes. The reason is usually because in a determinant system you want to take derivatives and require smooth function. If the data points are too far apart and you use a high order spline you also bias the resultant curve. You MUST always have a model of the system that's why social "sciences" are not science. $\endgroup$
    – wbg
    Commented Jan 13, 2022 at 16:47
  • $\begingroup$ Concerning your other argument: Splines are local polynomials. If the true function is not a local polynomial (e.g. a step function), then "model misspecification bias" remains also for splines. Also, splines have many more degrees of freedom than simple quadratic polynomial. Using a polynomial with the same number of degrees of freedom, the bias may be reduced. There's a huge literature about "fractional polynomials", which is a quite versatile class of functions for modelling. Of course, splines are more versatile, but that's mainly because the number of degrees is (almost) unbounded. $\endgroup$
    – jarauh
    Commented Jan 13, 2022 at 16:50

One troubling situation is when the regressor variables are not independent from each other. This can make it look like there is seemingly (causal) relationship that is not truly present.


Say we have a model where $y$ is a function of $x$ described by an exponential term $E[y|x] = e^{0.6x}$ and some Gaussian noise with deviation $\sigma = 0.1$

$$y \sim N(\mu = e^{0.6x}, \sigma = 0.1)$$

and we sample 10 sets of size 50 where each set is denoted by $z = 1,2,\dots, 10$.

Now, let's consider adding something tricky and let the $x$ data not be the same for each $z$ and let it following a continuous uniform distribution like

$$x \sim U(\text{lower} = 0, \text{upper} = 1+z/10)$$

Below we have plotted a simulation of this data in a plot and we have defined a binary variable based on the value of $z$. This variable might relate to your binary variable

$$z_{bin} = \begin{cases}1 &\quad \text{if z=10} \\ 0 &\quad \text{else} \end{cases}$$

The lines in the figure represent a linear fit with $x$ and $z_{bin}$ as variables. Note that in such a linear model the variable $z_{bin}$ relates to a significant effect ($\text{effect-size} = 0.09793$ and $\text{p-value} = 5.23 \cdot 10^{-7}$).

The reason for this significant effect is not because $z_{bin}$ occurs in the relationship between $y$ and $x$, but it is because we have different values of $x$ for $z_{bin} = 10$ and this will create some bias.


I believe this often happens in epidemiological research or other types of indirect observational research. A linear model is applied to the data and there is an indirect "controlling" for a variable by including it as a co-factor in the model. But, the only perfect way to control for variables is to make sure that the "independent" variables are also truly independent.

As a statistical model, it would still be not so wrong to infer certain relationships. However, one needs to be very critical/skeptical towards conclusions regarding causal relationships.

R-code to produce image

### create data
n = 50
yr = as.numeric(replicate(n,1:10))
x = runif(n*10,0,1+yr*0.1)
y = exp(0.6*x)+rnorm(n*10,0,0.1)
bin = (yr == 10)

### plotting 
col = rgb(1-bin,0.5,1-bin,0.6)
plot(x,y, pch = 21, col = rgb(0,0,0,0.6), bg = col, cex = 0.7)
title("example of observed effect for a variable \n that is not in the data generation model")

### linear model
mod = lm(y ~ x + bin)

xs = seq(0,2,0.01)
lines(xs,mod$coefficients[1] + mod$coefficients[2] * xs, col = rgb(1,0,1), lwd = 2)
lines(xs,mod$coefficients[1] + mod$coefficients[2] * xs + mod$coefficients[3], col = rgb(0,0.5,0), lwd = 2)
lines(xs,exp(0.6*xs), col = 1, lty = 2, lwd = 2)

legend(0,3, c("data generation model", "binary = 0" , "binary = 1"),
       lty = c(2,1,1), col = c(rgb(0,0,0),
                               rgb(0,0.5,0)), lwd = 2)

text(1.5,1.25, expression(y %~% N(mu == e^(0.6*x), sigma == 0.1)), pos = 3)
  • $\begingroup$ You have $y$ on both sides of an equation; I think that's a simple typo and that you intend $z$ on the right hand side instead. $\endgroup$
    – Glen_b
    Commented Jan 12, 2022 at 22:55
  • 1
    $\begingroup$ Thanks for the answer, which is certainly informative, but it's not really clear how it relates to my question. I checked the variable inflation factor and it was ~1 for all predictors, so it shouldn't be a problem here, right? I'm switching to rcs now anyway, but not sure how that handles multicollinearity. $\endgroup$
    – Felix Phl
    Commented Jan 13, 2022 at 10:08
  • $\begingroup$ @FelixPhl This is what came to my mind when I read your question, but I must admit that I wasn't very sure what your question exactly is about. Your question is very much loaded with details. Could you possibly trim down the question, or add short summary, such that it becomes more clear what is supposed to relate to it. The question title indicates something like 'drawing conclusion from a linear fit when the assumption of linearity is violated', but in the question body text it is more about consideration of non-linear modeling like using splines or quadratic functions. $\endgroup$ Commented Jan 13, 2022 at 11:01
  • $\begingroup$ I was mostly interested in whether I can "trust" the coefficients or not, in particular xbinary. But that has been covered in other answers now. I know that there is an effect during xbinary=1 which we cannot measure, but I wanted to get a value for the effect that is isolated from the effect of X1 and X2. $\endgroup$
    – Felix Phl
    Commented Jan 13, 2022 at 11:43
  • $\begingroup$ @FelixPhl I am not sure what you mean by "trust the coefficients", but this trust might be wrong if it means whether you want to use the significant linear model to suggest that there is some true effect or relationship, as in a causal relationship or whether the 'xbinary' variable would still play a role in the statistical relationship if the best fitting curve would have been used. The reason for it to be wrong is in my answer. The variable could be useful in the linear model only indirectly because it correlates with the error between the linear model and the true model. $\endgroup$ Commented Jan 13, 2022 at 12:22

As far as I understand, the assumption of Heteroscedasticity and Linearity are violated. As far as I understand, the first causes lower p values and imprecise coefficients, the second causes wrong coefficients?

If the "real" relationship isn't linear, then the linear coefficients aren't merely the "wrong coefficients", there aren't any right coefficients, because there's no line that correctly describes the data. Linear regression is about modelling the data as being linear, it isn't really assuming that the "real" relationship is actually linear. It's asking "What linear relationship minimizes squared error?" I.e. "What linear relationship is closest to the actual relationship?" It's basically projecting the space of all relationships onto the space of linear ones.

So the fact that the relationship isn't linear doesn't invalidate a linear model. What I find more concerning than the fact that the relationship is nonlinear is the fact that I'm seeing non-monotonic relationships. I would focus on addressing that issue. A linear model is one in which only first order terms are considered, but a non-monotonic relationship means that higher order terms are dominating. Using linear regression on categorical data is also a bit iffy.

Also, you might want to make the Y vs. Xbinary chart more visually informative, for instance by applying dithering to the Xbinary value.

  • 2
    $\begingroup$ And what if the 'model' is missing variables inducing apparent non-linear behavior? Understanding possible driving elements of the data should be the 1st goal as this aids in decision making, future forecasting,..., and in developing a better model, in my opinion. $\endgroup$
    – AJKOER
    Commented Jan 13, 2022 at 11:33
  • $\begingroup$ @AJKOER Unless you're talking about feature selection, I consider data collection and model building to be closely related, but separate issues. If OP can collect data on an additional variable, great, but otherwise they have to work with the data they have. $\endgroup$ Commented Jan 14, 2022 at 1:17

On the question: "Linearity assumption violated - can I still draw conclusions from my model?", my answer is a limited yes!

You may chose to verify my claim by constructing a simple non-linear model and assuming normal error generating model, project out. Further, fit a linear version of the model. I claim the Least-Square theoretical derived confidence intervals on the linear model's parameters, for which you have expressed an intuitive interpretation, should be generally well contained in the associated CIs. Why? Because you have basically encased model specification error into the error term.

As to how, and determining the best selected level for the CIs to explore, that is a simulation exercise. Hopefully, the exercise will render some meaningful guidance.

Even more work, break the curve into linear segments, and model per above, which may render further guidance as confirmed per simulation runs, as gross model mis-specification error is reduced. How do you select the point to go from one linear model to another, cut data points equally into half or...or select a random sub sample of the data which is either just inspected, or modeled for a so-called regime change. This data, however, for independence reasons, is not employed further.

This is my proposed answer to the 'any' part of your question.

Note: Not making an even bad guess on the non-linear nature of the model results in more work, but does preserve the intuitive content of the parameters, whose possible values are now expressed, at least, in a distributional sense.

  • 4
    $\begingroup$ To the original question, the phrase "my model" is loaded. If your model is typically a default linear model, expect your analyses to disagree with the way the world works. A better thought for "my model" is y ~ f(x) with our job of estimating f as others have said. $\endgroup$ Commented Jan 12, 2022 at 21:00
  • 1
    $\begingroup$ @FrankHarrell I agree. As I said, my choice to start with the linear model was based on the idea that it is the only regression form that it is easy to understand what happens (which is quite imporant in exploratory data analysis). Now that I learned about rcs I know that you can have both, a sufficiently good model and informative coefficients (or effect plots). $\endgroup$
    – Felix Phl
    Commented Jan 13, 2022 at 10:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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