Search code examples
rpackagebreakpoints

Equal and opposite slopes in segmented package

Hi I am trying to use a segmented package in R to fit a piecewise linear regression model to estimate break point in my data. I have used the following code to get this graph.

library(segmented)
set.seed(5)
x <- c(1:10, 13:22)
y <- numeric(20)
## Create first segment
y[1:10] <- 20:11 + rnorm(10, 0, 1.5)
## Create second segment
y[11:20] <- seq(11, 15, len=10) + rnorm(10, 0, 1.5)

## fitting a linear model
lin.mod <- lm(y~x)
segmented.mod <- segmented(lin.mod, seg.Z = ~x, psi=15)
summary(segmented.mod)
plot(x,y, pch=".",cex=4,xlab="x",ylab="y")
plot(segmented.mod, add=T, lwd = 3,col = "red")

enter image description here

My theoretical calculations suggests that the slopes of the two lines about the breakpoint should be equal in magnitude but opposite in sign. I am beginner to lm and glms. I was hoping if there is a way to estimate breakpoints with slopes constrained by the relation, slope1=-slope2

enter image description here

Solution

  • This is not supported in the segmented package.

    nls2 with "plinear-brute" algorithm could be used. In the output .lin1 and .lin2 are the constant term and the slope respectively. This tries each value in the range of x as a possible bp fitting a linear regression to each.

    library(nls2)

st <- data.frame(bp = seq(min(x), max(x)))
nls2(y ~ cbind(1, abs(x - bp)), start = st, alg = "plinear-brute")

    giving:

    Nonlinear regression model
  model: y ~ cbind(1, abs(x - bp))
   data: parent.frame()
       bp     .lin1     .lin2 
14.000000  9.500457  0.709624 
 residual sum-of-squares: 45.84213

Number of iterations to convergence: 22 
Achieved convergence tolerance: NA

    Here is another example which may clarify this since it generates the data from the same model as fit:

    library(nls2)

set.seed(123)
n <- 100
bp <- 25
x <- 1:n
y <- rnorm(n, 10 + 2 * abs(x - bp))

st <- data.frame(bp = seq(min(x), max(x)))
fm <- nls2(y ~ cbind(1, abs(x - bp)), start = st, alg = "plinear-brute")

    giving:

    > fm
Nonlinear regression model
  model: y ~ cbind(1, abs(x - bp))
   data: parent.frame()
    bp  .lin1  .lin2 
25.000  9.935  2.005 
 residual sum-of-squares: 81.29

Number of iterations to convergence: 100 
Achieved convergence tolerance: NA

    Note: In the above we assumed that bp is an integer in the range of x but we can relax that if such condition is not desired by using the result of nls2 as the starting value of an nls optimization, i.e. nls(y ~ cbind(1, abs(x - bp)), start = coef(fm)[1], alg = "plinear") .