Líneas de regresión con el mismo intercepto

Question

Lucho mucho con la Regresión. Acabo de descubrir cómo obtener 2 líneas con la misma pendiente, pero no consigo obtener 2 líneas con el mismo intercepto. He leído mucho sobre ANCOVA (porque pensé que era lo que necesitaba), pero nadie utiliza los mismos interceptos; sólo la misma pendiente. ¿Puede alguien ayudarme?

Answer 1

library(ggplot2)
set.seed(1)
x <-  1:10
dd <- rbind(data.frame(x=x,fac="a", y=x+rnorm(10)),
            data.frame(x=2*x,fac="b", y=x+rnorm(10)))
coef <- lm(y~x:fac, data=dd)$coefficients
qplot(data=dd, x=x, y=y, color=fac)+
  geom_abline(slope=coef["x:faca"], intercept=coef["(Intercept)"])+
  geom_abline(slope=coef["x:facb"], intercept=coef["(Intercept)"])

Answer 2

Aunque se trata de un hilo bastante antiguo, probablemente sea digno de mención que la respuesta de @jem77bfp parece funcionar sólo cuando el término de intercepción es cero o cercano a cero. Consideremos:

set.seed(1)
x <-  1:10
dd <- rbind(data.frame(x=10+x,fac="a", y=x+rnorm(10)),
            data.frame(x=10+2*x,fac="b", y=x+rnorm(10)))
coef <- lm(y~x:fac, data=dd)$coefficients

#(Intercept)      x:faca      x:facb 
# -7.0223128   0.8243321   0.6023107

Y aún más drásticamente:

set.seed(1)
x <-  1:10
dd <- rbind(data.frame(x=100+x,fac="a", y=x+rnorm(10)),
            data.frame(x=100+4*x,fac="b", y=x+rnorm(10)))
coef <- lm(y~x:fac, data=dd)$coefficients

# (Intercept)      x:faca      x:facb 
# -32.4026986   0.3610346   0.3122729 

Sugerencia de @Ben Bolker lm(y~x+f:x) se ajusta a dos pendientes y dos interceptos, lo que puede comprobarse al predecir "correctamente" las pendientes cuando ambos interceptos son diferentes.

No sé si hay una manera de explotar lm pero sin duda puede explotar minpack.lm::nls.lm especificando su propio modelo de error.

test <- data.frame(x = 1:10, 
                   y1 = 10 + 2*1:10 + rnorm(10, sd = 0.05), 
                   y2 = 10 + 8*1:10 + rnorm(10, sd = 0.05))

my_fun <- function(a, x, b1, b2) data.frame(y1 = a + x * b1, y2 = a + x * b2)

# this is the function which will yield the residuals; note that we need to unlist the data.frame finally
my_fun.res <- function(p, obs, x) unlist(obs - do.call(my_fun, c(list(x = x), as.list(p))))

minpack.lm::nls.lm(par = list(a = 1, b1 = 1, b2 = 1), fn = my_fun.res, 
                   obs = test[, c("y1", "y2")], x = test$x) -> pred

summary(pred)

# Parameters:
#    Estimate Std. Error t value Pr(>|t|)    
# a  10.009693   0.025435   393.5   <2e-16 ***
# b1  2.001747   0.004517   443.1   <2e-16 ***
# b2  7.996578   0.004517  1770.3   <2e-16 ***
```