Anchoring a linear regression to a specific data point in r

Question

This is different, follow up question to someone else's question here:

I am very new to R and have no programming experience what so ever but I am holding my own for the data analysis we need in my lab. I am detecting changepoints in physiological data using cpt.mean() and until now I have been doing two separate linear regressions with lm(). One for the first half of my data and another for my second half of my data. I determine where to divide up the data after cpt.mean() tells me my changepoint. The first line would normally have a slope around 0 and the other could be anywhere from 0.2 to 2. I would then add those lines separately to my graph and clip them so that it appeared that I only had only line graphed but that that line had a "corner." At first it worked really well but the point where the two lines joined together was often further from the changepoint found with cpt.mean because what I graphed was determined by where the two equations mathematically intersected.

So I thought it would be better to find linear regression lines that anchor to the changepoint. That way when I graph two seperate lines, they connect at the changepoint that cpt.mean. I have read and understand the basic solution to this question linked above but I cant seem to figure out how to limit the regression so that it only fits half of the data.

In a made up scenario my changepoint could be my 26th data point (10, 20)

Normally I do something like this for the first half of the data: (PP is my X and MAP is my Y variable):

lm(MAP[1:26]~PP[1:26])

The simple solution to anchoring the best fit line to a specific point,(X0, Y0) was this

lm(I(y-y0)~I(x-x0) + 0)

with my variables and changepoint:

lm(I(MAP-20)~I(PP-10) + 0)

I figured I could just combine the two as follows:

lm(I(MAP[1:26]-20)~I(PP[1:26]-10) +0)

No error came up but nothing showed up on my scatter plot. Not sure whats happening.

Here is an example of a what my graphs normally look like:

Best, Ian

Answer 1

Your model implicitly is this. You have data $(x_1,y_1), \ldots, (x_m,y_m)$ and other data $(x_{m+1},y_{m+1}), \ldots, (x_n,y_n)$ for which $x_j \ge x_i$ whenever $m \lt j \le n$ and $1 \le i \le m$: the first data set is the "left portion" and the second is the "right portion". You wish to estimate coefficients $\alpha, \alpha^\prime, \beta, \beta^\prime$ such that the "best fit" (presumably least squares) to the left half is

$$y = \alpha + \beta x,$$

the best fit to the right half is

$$y = \alpha^\prime + \beta^\prime x,$$

and the graphs of these two lines meet at some point $x_0$ between $\max_{1\ldots m} x_i$ and $\min_{m+1\ldots n} x_j$.

The two-regression procedure you have been doing appears to be an attempt to fit the related model

$$Y_i = \beta_0 + \beta_1 x_i + \beta_2(x_i-x_0)^{+} + \epsilon_i$$

where the $Y_i$ are random variables representing how $y_i$ were obtained; $\epsilon_i$ are independent, zero-mean random variables, they all have the same variance $\sigma_l^2$ in the left portion, and they all have the same variance $\sigma_r^2$ in the right portion; and for any number $z$, $$z^{+}=\max(0, z)$$ is the "positive part" of $z$.

This model assumes a slope of $\beta_1$ underlies the left portion (so $\beta_1$ plays the role of $\beta$) and a slope of $\beta_1 + \beta_2$ underlies the right portion (so $\beta_1+\beta_2$ plays the role of $\beta^\prime$). By using a single intercept, it assumes there is no break at $x_0$: the curve is continuous there. This could be called a "linear spline with one knot."

The easiest way to fit a model like this is multiple regression of the $y_i$ against the $x_i$ and the $(x_i-x_0)^{+}$ (which you would compute for this purpose). That fit assumes $\sigma_r^2 = \sigma_l^2$ (homoscedasticity), which is slightly different from what you have been doing but might be what you intended. If you want to fit the full (heteroscedastic) model, it looks necessary to work out the Maximum Likelihood solution.

Neither procedure will be the same as intersecting two separate least squares solutions--but if that intersection happens to be close to $x_0$, they will all give very similar fits. One huge advantage of using least squares or Maximum Likelihood is that they afford estimates of uncertainty, such as confidence intervals. These will be appropriate when the change point is found independently of the data; when determined from the data, confidence intervals have to be increased to account for the loss of that "degree of freedom."

To illustrate these ideas and demonstrate that the various procedures differ, the figure shows a scatterplot of 40 data points with a changepoint at $x_0=15$ (marked with the vertical line). Before the changepoint, errors drawn from a distribution with standard deviation $\sigma_l=3$ were added; after the changepoint, the distribution had SD $\sigma_r=15$. The true underlying lines are shown in solid gray; the separately fitted least squares lines in dashed blue; the homoscedastic solution proposed here is the unbroken red curve, and the heteroscedastic solution is the black curve. Note how--by construction--the bend in the fitted (red and black) lines occurs exactly at the changepoint. Because these data are strongly heteroscedastic, we would expect the heteroscedastic model to be better--and indeed the black curve is closer to the true gray curve than the red curve is (and the dotted blue curve, for the separate regressions, turns out worst in this example).

It might also be of interest to compare the summaries of three models: the homoscedastic fit, the fit to the left portion, and the fit to the right portion. Here they are, without further comment.

Call:
lm(formula = y ~ x + x.plus, data = X)
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.3357     5.4842   0.244    0.809    
x            -1.2498     0.4950  -2.525    0.016 *  
x.plus        3.3036     0.6729   4.910 1.86e-05 ***
---
Residual standard error: 10.92 on 37 degrees of freedom
Multiple R-squared:  0.6695,    Adjusted R-squared:  0.6516 
F-statistic: 37.47 on 2 and 37 DF,  p-value: 1.275e-09

Call:
lm(formula = y ~ x, data = X, subset = left)
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.3424     1.0319   3.239  0.00646 ** 
x            -1.6260     0.1135 -14.327 2.43e-09 ***
---
Residual standard error: 1.899 on 13 degrees of freedom
Multiple R-squared:  0.9404,    Adjusted R-squared:  0.9359 
F-statistic: 205.3 on 1 and 13 DF,  p-value: 2.429e-09


Call:
lm(formula = y ~ x, data = X, subset = !left)
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -43.4024    10.9327  -3.970 0.000606 ***
x             1.9033     0.3781   5.034 4.29e-05 ***
---
Residual standard error: 13.63 on 23 degrees of freedom
Multiple R-squared:  0.5242,    Adjusted R-squared:  0.5035 
F-statistic: 25.34 on 1 and 23 DF,  p-value: 4.29e-05

---

Here is R code to implement the procedures evaluated here and reproduce the data.

n <- 40              # Sample size
x <- 1:n             # X coordinates
x.0 <- 15            # X-coordinate of the changepoint location
beta <- c(5, -2, 5)  # Parameters: intercept, left slope, right slope
sigma <- c(3, 15)    # Error variances: left, right
#
# Generate data.
#
set.seed(17)
left <- x <= x.0
m <- sum(left)
x.plus <- pmax(0, x-x.0)
X <- cbind(1, x, x.plus)
y <- rnorm(n, X %*% beta, ifelse(x <= x.0, sigma[1], sigma[2]))
plot(x,y, pch=16, col="#00000040")
abline(v=x.0)
abline(beta[1:2], lwd=2, col="Gray")
abline(c(beta[1] - beta[3]*x.0, sum(beta[2:3])), lwd=2, col="Gray")
#------------------------------------------------------------------------------#
# 
# As fitting occurs, lines are added to the plot of the data to document
# progress.
#
# Linear spline fit.
#
X <- data.frame(x=x, x.plus=x.plus)
fit <- lm(y ~ x + x.plus, X)

x.hat <- seq(min(x), max(x), length.out=1001)
x.hat.plus <- pmax(0, x.hat-x.0)
y.hat <- predict(fit, newdata=data.frame(x=x.hat, x.plus=x.hat.plus))
lines(x.hat, y.hat, col="Red", lwd=2)
#
# Separate least-squares fit to left and right portions.
#
fit.left <- lm(y ~ x, X, subset=left)
fit.right <- lm(y ~ x, X, subset=!left)
abline(fit.left, col="Blue", lty=3, lwd=2)
abline(fit.right, col="Blue", lty=3, lwd=2)
#
# MLE, heteroscedastic case.
#
Lambda. <- function(theta, y, X, left) {
  lambda <- exp(theta)
  w <- ifelse(left, lambda, 1)          # sigma.left/sigma.right
  fit <- lm(y ~ ., data=X, weights=w)
  n <- length(y)
  r2 <- (resid(fit) / w)^2
  s2 <- mean(r2)
  (n + n*log(s2))/2 - sum(left)*theta
}
Lambda <- Vectorize(Lambda., "theta")
theta <- seq(-2, 2, length.out=10)
log.L <- Lambda(theta, y, X, left)
# plot(theta, log.L)

ml <- optimize(Lambda, c(-10, 10), y=y, X=X, left=left)
lambda.hat <- exp(ml$minimum)
w <- ifelse(left, lambda.hat, 1)
fit.ml <- lm(y ~ ., data=X, weights=w)
y.hat <- predict(fit.ml, newdata=data.frame(x=x.hat, x.plus=x.hat.plus))
lines(x.hat, y.hat, col="Black", lwd=2)

Answer 2

I thought it might be useful to illustrate the application of @whuber's linear-spline model to the data presented as an example in the question. The main thing is that for a known location, $x_0$, of the corner—@Sal Mangiafico's answer deals with joint estimation of that & the other parameters— it's linear in the parameters. Ordinary least-squares estimation is far from being the only game in town for fitting linear models, but it's well-known what properties it has under what conditions, & often constitutes at least a useful exploratory analysis in which it can be examined which of those conditions obtain & which not, suggesting other methods to apply.

---

A least-squares fit of the linear-spline model to the data (taken from the graph) gives

Call:
lm(formula = y ~ x + x.plus, data = dd)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.5041  -1.6519   0.2423   1.9444  10.7108 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.6770     9.8149  -0.171   0.8648    
x             0.3780     0.1862   2.030   0.0461 *  
x.plus        2.6665     0.3233   8.247 5.82e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.984 on 71 degrees of freedom
Multiple R-squared:  0.8548,    Adjusted R-squared:  0.8507 
F-statistic:   209 on 2 and 71 DF,  p-value: < 2.2e-16

& it might be interesting to compare it to a model allowing a jump at $x_0$ rather than a corner:

$$Y_i = \beta_0 + \beta_1 x_i + \beta_2(x_i-x_0)^{+} + \beta_3\mathbf{1}(x_i>x_0) + \epsilon_i$$

Call:
lm(formula = y ~ x + x.plus + x.RH, data = dd)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.4672  -1.6247   0.2558   1.8407  10.7461 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -2.3868    11.9154  -0.200    0.842    
x             0.3930     0.2343   1.677    0.098 .  
x.plus        2.6652     0.3258   8.180 8.53e-12 ***
x.RH         -0.2047     1.9187  -0.107    0.915    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.012 on 70 degrees of freedom
Multiple R-squared:  0.8549,    Adjusted R-squared:  0.8486 
F-statistic: 137.4 on 3 and 70 DF,  p-value: < 2.2e-16

The regression lines on each side are very close to meeting where they're supposed to on the corner theory; &, from the standard error estimates above a 95% confidence interval for $\beta_3=0$ is $(-4.0,3.6)$:

Diagnostic plots of the residuals suggest, rather than heteroskedasticity, outliers or a fat-tailed error distribution, with especially influential points to the far left:

A smattering of robust regression methods—median regression (orange line), & M-estimation with Huber (red) & Tukey's biweight (green) objective functions—give rather similar fits to each other while reducing the slope on the left-hand side relative to the least-squares solution :

R code:

library(MASS)
library(quantreg)

# data got from graph
dd <- structure(list(
  x = c(41.5595, 42.9676, 46.4784, 46.8849, 46.8909, 49.8839, 49.8902, 50.857, 51.8846, 52.2893, 52.9127, 53.2899, 53.6952, 53.8567, 54.2075, 54.5852, 54.9895, 55.3148, 55.7208, 55.9117, 56.5042, 56.7189, 56.7204, 56.8006, 56.9859, 57.2308, 57.6912, 57.9051, 58.0147, 58.2289, 58.448, 58.5033, 58.7166, 58.7693, 58.8582, 58.9781, 58.9848, 59.4165, 59.5783, 59.6881, 59.7153, 59.793, 59.9585, 60.169, 60.1999, 60.2531, 60.3583, 60.7897, 60.9801, 61.0065, 61.2477, 61.2524, 61.275, 61.2997, 61.3005, 61.4096, 62.0922, 62.1683, 62.1958, 62.5183, 62.7057, 62.738, 62.8159, 63.1095, 63.1357, 63.1657, 63.1937, 63.4924, 63.5732, 63.7328, 63.7855, 63.8129, 63.9275, 64.38),
  y = c(18.8093,14.6087, 18.9369, 17.1767, 7.9153, 18.7761, 8.88974, 18.4948, 17.6454, 18.7261, 16.0574, 17.8767, 17.9915, 19.0715, 19.9248, 20.9486, 22.7111, 21.0643, 20.27, 17.4864, 20.7836, 23.2274, 20.8978, 22.2617, 28.285, 25.672, 24.1961, 27.833, 25.5606, 28.5726, 24.1414, 22.2665, 26.699, 28.8014, 16.7562, 40.3362, 29.9384, 31.076, 31.6446, 29.0313, 28.7473, 34.0884, 28.8616, 37.7259, 31.6464, 33.067, 37.556, 39.0913, 37.3305, 38.2396, 41.4221, 34.0926, 40.9677, 44.6609, 43.4109, 41.8772, 31.1404, 38.8679, 38.0726, 40.9144, 43.6422, 35.5173, 40.4607, 46.257, 47.4503, 42.9049, 41.3709, 39.1558, 39.6106, 43.5883, 45.6339, 45.1226, 35.0661, 45.8061)),
  .Names = c("x", "y"),
  row.names = c(1L, 2L, 3L, 5L, 4L, 7L, 6L, 8L, 9L, 11L, 10L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 23L, 22L, 24L, 25L, 26L, 27L, 28L, 29L,34L, 30L, 31L, 33L, 35L, 32L, 52L, 36L, 37L, 38L, 39L, 40L, 43L,41L, 51L, 42L, 44L, 50L, 53L, 49L, 54L, 62L, 45L, 61L, 65L, 64L, 63L, 46L, 55L, 56L, 60L, 66L, 47L, 59L, 72L, 73L, 67L, 68L, 57L, 58L, 69L, 74L, 70L, 48L, 71L),
  class = "data.frame"
  )

x.c <- 55.6 # the corner

# add linear spline term
dd$x.plus <- pmax(dd$x-x.c,0)
# add right-hand-side flag
dd$x.RH <- as.numeric(dd$x>x.c)

# make predictor sequence for plots
dfp <- data.frame(x=c(seq(40,x.c, by=0.1), seq(x.c+1e-7, 66, by=0.1)))
dfp$x.plus <- pmax(dfp$x-x.c,0)
dfp$x.RH <- as.numeric(dfp$x>x.c)

# fit linear spline model by least-squares
lm(y ~ x + x.plus, data=dd) -> mod.ls
# fit jump model by least squares
lm(y ~ x + x.plus + x.RH, data=dd) -> mod.ls.jump

# compare results
summary(mod.ls)
summary(mod.ls.jump)


# plot fits
png("plot1.png", width=500, height=400)
par(mfrow=c(2,1))
with(dd,plot(x, y, pch=20))
lines(dfp$x, predict(mod.ls, newdata=dfp), col="dodgerblue")
lines(dfp$x, predict(mod.ls.jump, newdata=dfp), col="darkorange")
abline(v=x.c, lty=2)
dev.off()
png("plot2.png", width=500, height=400)
with(dd, plot(x, y, pch=20, xlim=c(53,57), ylim=c(18,22)))
lines(dfp$x, predict(mod.ls, newdata=dfp), col="dodgerblue")
lines(dfp$x, predict(mod.ls.jump, newdata=dfp), col="darkorange")
abline(v=x.c, lty=2)
dev.off()

# make some diagnostic plots
png("diag1.png", width=400, height=400)
plot(mod.ls, which=1)
dev.off()
png("diag2.png", width=400, height=400)
plot(mod.ls, which=2)
dev.off()
png("diag3.png", width=400, height=400)
plot(mod.ls, which=5)
dev.off()
influence.measures(mod.ls)

# fit using a few robust methods
rq(y ~ x + x.plus, data=dd) -> mod.medreg # median regression
rlm(y ~ x + x.plus, data=dd, method="M", psi=psi.huber,  scale.est="MAD") -> mod.Huber
rlm(y ~ x + x.plus, data=dd, method="M", psi=psi.bisquare,  scale.est="MAD") -> mod.TukeyBW

png("plot3.png", width=500, height=400)
with(dd,plot(x, y, pch=20))
lines(dfp$x, predict(mod.ls, newdata=dfp), col="dodgerblue")
lines(dfp$x, predict(mod.medreg, newdata=dfp), col="darkorange")
lines(dfp$x, predict(mod.Huber, newdata=dfp), col="firebrick1")
lines(dfp$x, predict(mod.TukeyBW, newdata=dfp), col="forestgreen")
abline(v=x.c, lty=2)
dev.off()

Answer 3

This situation is similar to one traditionally encountered in agriculture: using linear-plateau models to model the response of a crop to the concentration of nutrients in the soil.

The code below uses the nls function in R. It includes one version where the first segment is held at zero slope, and one version where the first segment is allowed to have a non-zero slope.

There is no need to do any change point analysis beforehand. The nls procedure estimates the clx (critical value of x) where the two segments join.

It is helpful to start with realistic estimates of the slopes, intercepts, and clx. These are entered as e.g. a.ini in the code.

From the nls model, you can determine a p-value for the model, a pseudo r-squared for the model, and confidence intervals for the parameter estimates. See the webpage linked in the code for examples.

Note that nls can have difficulty finding solutions for difficult data sets. In this case, the package nlmrt is sometimes better.

Additional notes on the nls function:

The nls function is included with R in the native stats package. The function name stands for nonlinear least squares. It uses an iterative process to determine the values of the parameters in the model using least squares. It is used to fit non-linear models (for example if you wanted to fit an S-shaped curve or bell-shaped curve to data). It will also fit segmented models, as in this example or for linear-plateau or quadratic plateau models.

It is best to give the function realistic starting values for parameters, as it is possible for the algorithm to "wander off" and converge on locally-good but globally-poor parameter estimates. As mentioned, there are other functions in R that may work better on problematic data sets or models with more than a few parameters.

### Data from Scortchi's answer
X = c(41.5595, 42.9676, 46.4784, 46.8849, 46.8909, 49.8839, 49.8902, 50.857, 51.8846, 52.2893, 52.9127, 53.2899, 53.6952, 53.8567, 54.2075, 54.5852, 54.9895, 55.3148, 55.7208, 55.9117, 56.5042, 56.7189, 56.7204, 56.8006, 56.9859, 57.2308, 57.6912, 57.9051, 58.0147, 58.2289, 58.448, 58.5033, 58.7166, 58.7693, 58.8582, 58.9781, 58.9848, 59.4165, 59.5783, 59.6881, 59.7153, 59.793, 59.9585, 60.169, 60.1999, 60.2531, 60.3583, 60.7897, 60.9801, 61.0065, 61.2477, 61.2524, 61.275, 61.2997, 61.3005, 61.4096, 62.0922, 62.1683, 62.1958, 62.5183, 62.7057, 62.738, 62.8159, 63.1095, 63.1357, 63.1657, 63.1937, 63.4924, 63.5732, 63.7328, 63.7855, 63.8129, 63.9275, 64.38)
Y = c(18.8093,14.6087, 18.9369, 17.1767, 7.9153, 18.7761, 8.88974, 18.4948, 17.6454, 18.7261, 16.0574, 17.8767, 17.9915, 19.0715, 19.9248, 20.9486, 22.7111, 21.0643, 20.27, 17.4864, 20.7836, 23.2274, 20.8978, 22.2617, 28.285, 25.672, 24.1961, 27.833, 25.5606, 28.5726, 24.1414, 22.2665, 26.699, 28.8014, 16.7562, 40.3362, 29.9384, 31.076, 31.6446, 29.0313, 28.7473, 34.0884, 28.8616, 37.7259, 31.6464, 33.067, 37.556, 39.0913, 37.3305, 38.2396, 41.4221, 34.0926, 40.9677, 44.6609, 43.4109, 41.8772, 31.1404, 38.8679, 38.0726, 40.9144, 43.6422, 35.5173, 40.4607, 46.257, 47.4503, 42.9049, 41.3709, 39.1558, 39.6106, 43.5883, 45.6339, 45.1226, 35.0661, 45.8061)

Data = data.frame(X, Y)

####################################################
### Plateau-linear model
### Adapted from: http://rcompanion.org/handbook/I_11.html
###################################################

a.ini     = -150
b.ini     =    3
clx.ini   =   55

platlin = function(x, a, b, clx)
          {ifelse(x < clx, a + b * clx,
                           a + b * x)}

model = nls(Y ~ platlin(X, a, b, clx),
            data = Data,
            start = list(a   = a.ini,
                         b   = b.ini,
                         clx = clx.ini),
             trace = FALSE,
             nls.control(maxiter = 1000))

summary(model)

plot(X, Y, pch=16)
i=(410:650)/10
D1 = data.frame(X = i)
predy = predict(model, D1)
lines(i, predy, lwd=2, col="blue")

####################################################
### Linear-linear model
### Adapted from: http://rcompanion.org/handbook/I_11.html
###################################################

a1.ini     =    0
b1.ini     =    0.5
b2.ini     =    3
clx.ini    =   55

segment = function(x, a1, b1, b2, clx)
          {ifelse(x < clx, a1 + b1 * x,
                           (a1+b1*clx-b2*clx) + b2 * x)}

model2 = nls(Y ~ segment(X, a1, b1, b2, clx),
            data = Data,
            start = list(a1  = a1.ini,
                         b1  = b1.ini,
                         b2  = b2.ini,
                         clx = clx.ini),
             trace = FALSE,
             nls.control(maxiter = 1000))

summary(model2)

A = summary(model2)$parameters
a2 = A[1,1]+A[2,1]*A[4,1]-A[3,1]*A[4,1]
names(a2) = "a2"
a2

plot(X, Y, pch=16)
i=(410:650)/10
D1 = data.frame(X = i)
predy = predict(model2, D1)
lines(i, predy, lwd=2, col="blue")

Answer 4

As per the comments above, the question is unclear. Unless there is strong evidence that the model should consist of two line segments, it is most often a mistake to model data that way. In order to prove that such a procedure is correct, one has to offer evidence to that effect. Now the counterargument to that claim is that the change-point exists by "assumption." But is that really statistics or a mathematical exercise in futility? I claim the latter, and unless someone can show a concrete example of a physical model process that results in such a horrible assumption, I will stand my ground.

Later in this post I show, as an example, how a biexponential, i.e., a single function, fits the data, and it has only four parameters, whereas two line segments have 5, because the cut point is a variable. Now to be clear, one can continue the process of cutting the data into segments and solve for exact line segments connecting each and every adjacent point. Then the residuals are zero, the correlations are all ~1, and we have made a model that is a fit of a perfect fit with no noise. What's wrong with that picture? 1) It isn't a proper data model. 2) I could have solved for each point with a Lagrange polynomial so it is not unique. 3) It is highly unusual to have a continuous model of some physical process that has a discontinuous derivative at a point. Thus, without proof amounting to an unusual and changed model at some cut point, we have to assume that the model is posited in error. As there is no evidence provided to the effect that two separate physical processes have given rise to the data, most likely such a model is incorrect. And, if there is any such indication please provide it. Without that proof of principle, the question is ambiguous and should be closed.

Let us assume we have a 2D point, which is the literal interpretation of a "change point," i.e., $(x_c,y_c)$. Now suppose we do not and we only know $x_c$, then let $y_c=\dfrac{y_L(x_c)+y_R(x_c)}{2}$, where $y_L=m_L\,x+B_L$ for the left hand line and $y_R=m_R\,x+B_R$ is the right hand line. Suppose we do not like that, then do something else, like a weighted average. Anyway, we now assume that by hook or by crook we have $(x_c,y_c)$. Now we need new lines; $y'_{L}$ and $y'_{R}$, and since we already have one point, we only need one more point and one logical extra point is to require each respective line to go through its data centroid $(\bar{x}_{L},\bar{y}_{L})$, and $(\bar{x}_{R},\bar{y}_{R})$. Now, those centroids are merely the mean $x$ and $y$ coordinates of the left hand data, and the right hand data.

Then for the left hand new line we have simultaneous equations

$$y'_{L}-y_{c}=m'_{L}(x'_{L}-x_{c}),$$ and $$y'_{L}-\bar{y}_{L}=m'_{L}(x'_{L}-\bar{x}_{L})$$ whence $$m'_{L}=\dfrac{y_{c}-\bar{y}_{L}}{x_{c}-\bar{x}_{L}}$$thus $$y'_{L}=\dfrac{y_{c}-\bar{y}_{L}}{x_{c}-\bar{x}_{L}}(x'_{L}-\bar{x}_{L})+\bar{y}_{L}$$and similarly $$y'_{R}=\dfrac{\bar{y}_{R}-y_{c}}{\bar{x}_{R}-x_{c}}(x'_{R}-\bar{x}_{R})+\bar{y}_{R}$$

Try it and note there are many many answers to your question.

@whuber There are assumptions being made that are being ignored. The fitting in the question does not pass my eyeball test, never mind more sophisticated testing. To show what my eyeball is seeing I have put in some hand drawn eyeball regression lines. My claim is that these lines are as spurious as the actual OLS regression lines, maybe less.

Compare this illustration with the OP's and see if you can see what I see. I have two sizable problems with the original regressions shown. The first is that both regression lines have visually obviously flattened slopes. The second is that the left line does not have homoscedastic residuals. Now, to proceed with calculation of slope confidence intervals, where I, admittedly personally, have no confidence in the aptness of the regression lines shown is derivative enough to give me pause for thought. In particular, there is so-called omitted variable bias, and this is more pronounced at the end points of the regressions. I would suggest using Theil regression to reduce this bias. However, as stated the problem has no exact solution, and all that can be done is to do a median (approximate) regression for this bivariate problem. Supply the actual data, and I will show this.

The results are very dependent on what the cut point is, and the regression method used. For example, if I use @Scortchi 's data and a slightly lesser x-value cut point and Passing-Bablok conversion regression (type III) for the left hand points and Passing-Bablok regression type I for the right hand points I get significant slopes and intercepts as below:

Now, the left hand points would not have either significant slope nor intercept using ordinary least squares regression, but do using Passing-Bablok regression. If I move the cut point to the right, as for Scortchi this significance for the left hand regression disappears.

Next thought, I think the data is inadequate to prefer two lines over a single function with local terms; let us say a mixture of two exponential functions, e.g.,