Pregunta

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.

¿Fue útil?

Solución

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.

Otros consejos

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 

If the regression coefficients in the two models are not equal (which is almost certain) then the lines would intersect. The coef function is used to extract them. The rest is high school geometry.

For Brandon: M^-1 %*% intercepts -->

M <- matrix( c(coef(m1)[2], coef(m2)[2], -1,-1), nrow=2, ncol=2 )
intercepts <- as.matrix( c(coef(m1)[1], coef(m2)[1]) )  # a column matrix
 -solve(M) %*% intercepts
#        [,1]
#[1,] 12.78597
#[2,] 16.34479

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
Licenciado bajo: CC-BY-SA con atribución
No afiliado a StackOverflow
scroll top