9

I have two linear fits that I've gotten from lm calls in my R script. For instance...

fit1 <- lm(y1 ~ x1)
fit2 <- lm(y2 ~ x2)

I'd like to find the (x,y) point at which these two lines (fit1 and fit2) intersect, if they intersect at all.

0

3 Answers 3

Reset to default

Trending sort

Trending sort is based off of the default sorting method — by highest score — but it boosts votes that have happened recently, helping to surface more up-to-date answers.

It falls back to sorting by highest score if no posts are trending.

13

Here's some high school geometry then ;-)

# First two models
df1 <- data.frame(x=1:50, y=1:50/2+rnorm(50)+10)
m1 <- lm(y~x, df1)

df2 <- data.frame(x=1:25, y=25:1*2+rnorm(25)-10)
m2 <- lm(y~x, df2)

# Plot them to show the intersection visually    
plot(df1)
points(df2)

# Now calculate it!    
a <- coef(m1)-coef(m2)
c(x=-a[[1]]/a[[2]], y=coef(m1)[[2]]*x + coef(m1)[[1]])

Or, to simplify the solve-based solution by @Dwin:

cm <- rbind(coef(m1),coef(m2)) # Coefficient matrix
c(-solve(cbind(cm[,2],-1)) %*% cm[,1])
# [1] 12.68034 16.57181 
5
  • 1
    Or just use solve. Although maybe that's more of a college solution? ;)
    – joran
    Aug 18, 2011 at 22:51
  • lol, I could have done the high school geometry. just wanted to see if there was a better way!!
    – CodeGuy
    Aug 18, 2011 at 22:55
  • Can you show a solution using solve? I'd be interested to see it. I don't use it very frequently. Aug 18, 2011 at 22:58
  • 1
    @DWin: I incorporated a slightly more compact variant of your solve solution in my answer.
    – Tommy
    Aug 18, 2011 at 23:50
  • 2
    -1 for "taking the OP back to school" - that's not really constructive
    – DQdlM
    Aug 20, 2011 at 11:27
9

One way to avoid the geometry is to re-parameterize the equations as:

y1 = m1 * (x1 - x0) + y0
y2 = m2 * (x2 - x0) + y0

in terms of their intersection point (x0, y0) and then perform the fit of both at once using nls so that the returned values of x0 and y0 give the result:

# test data
set.seed(123)
x1 <- 1:10
y1 <- -5 + x1 + rnorm(10)
x2 <- 1:10
y2 <- 5 - x1 + rnorm(10)
g <- rep(1:2, each = 10) # first 10 are from x1,y1 and second 10 are from x2,y2

xx <- c(x1, x2)
yy <- c(y1, y2)
nls(yy ~ ifelse(g == 1, m1 * (xx - x0) + y0, m2 * (xx - x0) + y0),
    start = c(m1 = -1, m2 = 1, y0 = 0, x0 = 0))

EDIT: Note that the lines xx<-... and yy<-... are new and the nls line has been specified in terms of those and corrected.

21
  • this is perfect. is it possible for you to show me, using this code you already have, how I can specify that the following requirement: the solution (the two lines) must have an intersection point between z1 and z2 (two values I specify)
    – CodeGuy
    Aug 19, 2011 at 14:59
  • @CodeGuy, Specify the arguments: algorithm = "port", lower = ...whatever..., upper = ...whatever... as per ?nls. Aug 19, 2011 at 15:59
  • I'm not familiar with what you mean. Can you add it to the code and show me?
    – CodeGuy
    Aug 19, 2011 at 16:27
  • @CodeGuy, Assuming you want to restrict x0 to lie between 2 and 4: nls(c(y1, y2) ~ ifelse(g == 1, b1 * (x1 - x0) + y0, b2 * (x2 - x0) + y0), start = c(b1 = -1, b2 = 1, y0 = 0, x0 = 3), algorithm = "port", lower = c(b1 = -Inf, b2 = -Inf, y0 = -Inf, x0 = 2), upper = c(b1 = Inf, b2 = Inf, y0 = Inf, x0 = 4)) Aug 19, 2011 at 16:55
  • wait a sec. the input is just x and y values. I don't have two sets of (x,y) values. Just one set.
    – CodeGuy
    Aug 19, 2011 at 16:58
3

I am a little surprised there isn't a built in function for this.

Here is a rudimentary function (for lm results), using the same general method as Tommy above. This uses the simple substitution method for two lines in the form "y=mx+b" to find the common intersection at y (y1=y2 ; m1*x + b1 = m2*x + b2) and solves for x:

Function definition

# Linear model Intercept function
lmIntx <- function(fit1, fit2, rnd=2) {
  b1<- fit1$coefficient[1]  #y-int for fit1
  m1<- fit1$coefficient[2]  #slope for fit1
  b2<- fit2$coefficient[1]  #y-int for fit2
  m2<- fit2$coefficient[2]  #slope for fit2
  if(m1==m2 & b1==b2) {print("Lines are identical")
  } else if(m1==m2 & b1 != b2) {print("Lines are parallel")
  } else {
    x <- (b2-b1)/(m1-m2)      #solved general equation for x
    y <- m1*x + b1            #plug in the result
    data.frame(x=round(x, rnd), y=round(y, rnd))
  }
}

Test:

line1 <- data.frame(x=c(0,1), y=c(0,2))
line2 <- data.frame(x=c(0,1), y=c(1,3))
line3 <- data.frame(x=c(0,1), y=c(1,5))

lmIntx(lm(line1$y~line1$x), lm(line2$y~line2$x))
[1] "Lines are parallel"


lmIntx(lm(line1$y~line1$x), lm(line1$y~line1$x))
[1] "Lines are identical"

lmIntx(lm(line1$y~line1$x), lm(line3$y~line3$x))
               x  y
(Intercept) -0.5 -1

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.