5
$\begingroup$

Suppose I have a (linear) relationship between X and Y that might look like in my poor drawing below (simplified), i.e. I see that the slope of my Y~X regression differs depending on the value/range of my independent variable X.

What would be the appropriate way/test to determine if the three "sub-slopes" are significantly different?

My first idea was to simply split the data into three distinct groups and conduct an ANOVA of Y~X. However, this only tells me if the group means are different (which is a good start), but I need to know if the slopes are different.

My second idea was to simply add the group indicator as an interaction term to the regression, but it somehow feels wrong, but maybe it's as simple as that?

Bonus: thinking about it I probably not only want to know if the slopes are different, but the intercept-slope combination. Looking at my example below, the red slope might not be significantly different from the blue slope, however, the intercept is certainly higher. I guess, I could go back to my first idea and look at the group means instead of the intercept, but I was wondering if there's a test comparing both (intercept and slope) at the same time.

enter image description here

$\endgroup$
5
  • $\begingroup$ How are you determining the three groups? Are you going to be choosing them to make the three slopes look as different as possible or do you have some other information (e.g. the three groups could represent distinct groups, or if $X$ represents time there could be two known interventions)? $\endgroup$
    – Alex
    Jul 7, 2020 at 17:41
  • $\begingroup$ (It matters because your model for the data will have more parameters if you choose the three groups yourself) $\endgroup$
    – Alex
    Jul 7, 2020 at 18:02
  • $\begingroup$ I guess it would be more appropriate to choose the cut points a priori based on some theory and hypotheses, but to be honest, in my case we are setting the cut points post-hoc, i.e. after creating a scatterplot I saw the two "change points" and am now interested if their respective slopes are different. After some further digging I found sth. called segmented regression which I guess goes into the direction I'm aiming at: en.wikipedia.org/wiki/Segmented_regression $\endgroup$
    – deschen
    Jul 7, 2020 at 18:42
  • $\begingroup$ And what software are you using to analyse the data? $\endgroup$
    – Alex
    Jul 7, 2020 at 19:01
  • $\begingroup$ I'm using R. For the segmented regression I tried the segmented package. $\endgroup$
    – deschen
    Jul 8, 2020 at 8:11

1 Answer 1

1
$\begingroup$

The simplest way is to use a selection criterion such as the Bayesian Information Criterion (Schwarz, 1978). Other methods are available, such as the supF-tests described in section 5.1 of (Bai and Perron, 2003), but using an information criterion is far simpler.

The BIC for a model $M$ is given by

$$\mathrm{BIC} = k\ln(n)-2\ln(\widehat {L}),$$ where

  • $\hat {L}$ is the maximized value of the likelihood function of the model $M$, i.e. $\hat {L}=p(x\mid {\widehat {\theta }},M)$, where $\widehat{\theta}$ are the parameter values that maximize the likelihood function for model $M$;
  • $x$ is the observed dataset;
  • $n$ is the number of data points in $x$; and
  • $k$ is the number of parameters estimated by the model.

Models with low BIC values are preferred because good models:

  1. fit the data well (so have low $-2\ln(\hat{L})$ values), and
  2. don't have many parameters (so have lower $k \ln(n)$ values).

I've simulated some data and done an analysis below:

library(segmented)
set.seed(1)

n = 300
x1 = runif(n/3, min = 0, max = 5)
y1 = 1 + 2*x1 + rnorm(n/3)
x2 = runif(n/3, min = 5, max = 15)
y2 = 6 + 1*x2 + rnorm(n/3)
x3 = runif(n/3, min = 15, max = 20)
y3 = -1.5 + 1.5*x3 + rnorm(n/3)

x = c(x1, x2, x3)
y = c(y1, y2, y3)
par(mar = c(4.1, 4.1, 0.1, 0.1))
plot(x, y, xlab = "x", ylab = "y", cex = 0.5, pch = 16)

Data with two changes

Now let's compare the model with just one regression line, versus the model with three groups:

lm1 = lm(y ~ x)
s1 = segmented(lm1, seg.Z = ~x, npsi = 2)
par(mar = c(4.1, 4.1, 0.1, 0.1))
plot(x, y, xlab = "x", ylab = "y", cex = 0.5, pch = 16)
abline(lm1, lwd = 2)
sss = seq(from = 0, to = 20, length.out = 1000)
lines(sss, predict.segmented(s1, newdata = data.frame(x = sss)), lty = 2, lwd = 2)

With two lines

The model estimated by segmented is pretty close to the truth:

summary(s1)
intercept(s1)

gives

***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm1, seg.Z = ~x, npsi = 2)

Estimated Break-Point(s):
          Est. St.Err
psi1.x  4.959  0.268
psi2.x 15.860  0.422

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.84810    0.22134   3.832 0.000156 ***
x            2.04758    0.07677  26.673  < 2e-16 ***
U1.x        -1.04908    0.08184 -12.818       NA    
U2.x         0.67839    0.09865   6.876       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.005 on 294 degrees of freedom
Multiple R-Squared: 0.985,  Adjusted R-squared: 0.9848 

Convergence attained in 5 iter. (rel. change 3.3705e-06)

$x
              Est.
intercept1  0.8481
intercept2  6.0505
intercept3 -4.7085

so the change points are estimated well, with guesses of $4.959$ and $15.860$ when the true values are $5$ and $15$. The intercepts and slopes are also done well, with $(0.85, 2.0)$ estimated for $(1, 2)$; $(6.1, 1.0)$ estimated for $(6, 1)$; and $(-4.7, 2.7)$ estimated for $(-1.5, 1.5)$.

For a linear regression with normal errors, the BIC can be calculated from the residual sum of squares (RSS) and is given by: $$\mathrm{BIC} = n\ln(\textrm{RSS}/n) + k\ln(n) + C(n),$$ where $C(n)$ does not depend on the model complexity or fit, so we ignore it.

For the linear model we calculate the BIC as:

n*sum(lm1$residuals^2/n) + 3*log(n)

which is roughly $537$. For the simple model $k = 3$ because the parameters are the intercept, the slope and the variance.

For the three groups model we calculate the BIC as:

n*sum(s1$residuals^2/n) + 9*log(n)

which is roughly $349$. For this model $k=9$ because we have three sets of intercepts and slopes, two change points and the variance. The difference between the two BICs is roughly $188$, which is massive evidence in favour of the model with three groups. A difference of BICs of more than $10$ is considered to be very strong evidence in favour of one model according to (Kass and Raftey, 1995).

Suppose that the true model is that there is just one slope, and we try to fit a segmented regression to it, the BIC would hopefully show us that the best model is the simplest one

set.seed(100)
y2 = 3 + x + rnorm(100)
lm2 = lm(y2 ~ x)
s2 = segmented(lm2, seg.Z = ~x, npsi = 2)
n*sum(lm2$residuals^2/n) + 3*log(n)
n*sum(s2$residuals^2/n) + 9*log(n)

In this case the BIC for the simple model is roughly $326$, and the BIC for the complex model is roughly $356$. So there is strong evidence that the simple model is better.

References

Bai, J. and P. Perron (2003). "Computation and analysis of multiple structural change models", Journal of Applied Econometrics 18 (1), 1–22

Kass, Robert E.; Raftery, Adrian E. (1995), "Bayes Factors", Journal of the American Statistical Association, 90 (430): 773–795

Schwarz, Gideon E. (1978). "Estimating the dimension of a model", Annals of Statistics, 6 (2): 461–464

$\endgroup$
4
  • $\begingroup$ Thanks for this detailed response! Quite helpful. A few questions, though: using the BIC implies a full model comparison (i.e. is the segmented model better/worse than a simple linear model), right? Is there any way to also compare the three slopes with each other, e.g. is slope 1 significantly different from the second one? $\endgroup$
    – deschen
    Jul 9, 2020 at 12:57
  • $\begingroup$ Is there any disadvantage/problem inherent to the approach you described? I'm asking because your answer seems to be downvoted, but there's no explanation in a comment or anywhere else why that person downvoted your answer. $\endgroup$
    – deschen
    Jul 9, 2020 at 12:58
  • 1
    $\begingroup$ You could calculate the BIC for the model with just one break and compare the BICs for all three models if you wanted. Or if you just want to enforce slope 1 = slope 2 (but potentially intercept 1 $\neq$ intercept 2), you could do that (not sure if that's possible with this package, but you could do it manually if you change the code) and then compare the models' BICs. Also U1.x in summary(s1) gives an estimate of the difference between the first slope and the second. Since the t-stat is $-12$, this is strong evidence that the difference is non-zero. $\endgroup$
    – Alex
    Jul 9, 2020 at 15:06
  • $\begingroup$ There are some other different ways you could analyse this (and I mentioned one of them) but the method I gave is flexible and widely used. The downvote happened less than 10 seconds after I posted my answer, so it can't be from someone who read the whole thing. $\endgroup$
    – Alex
    Jul 9, 2020 at 15:10

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.