[R] How to compare two regression line slopes

Ben Bolker bolker at ufl.edu
Tue Jan 27 22:52:43 CET 2009


Etienne Toffin <etoffin <at> ulb.ac.be> writes:

> 
> I've made a research about how to compare two regression line slopes  
> (of y versus x for 2 groups, "group" being a factor ) using R.
> 
> I knew the method based on the following statement :
> t = (b1 - b2) / sb1,b2
> where b1 and b2 are the two slope coefficients and sb1,b2 the pooled  
> standard error of the slope (b)
> 
> However, I also found a procedure in Wonnacott & Wonnacott, that is  
> based on the use of a mute variable D that will have a binary value  
> according to the group to which a given point belongs (group : D=0;  
> group 2: D=1). Then the equation that is computed is as follow:
> y = b0 + b1.x + D.b2.x
> 
> which can be computed in R with:
>  > fit <- lm(y ~ group + x + x:group)
> where y is the response of the 2 groups.
> The p-value of x:group gives the probability for the two slopes to be  
> different, and the estimated values of parameters are these of both  
> populations.
> 
> These two methods have already been described in the mailing list but  
> not confronted and discussed.
> So, my questions are:
> - are these methods different ?
> - which one should be preferentially used ?
> 
   I think you're perfectly clear.
  These procedures are identical: the first has the virtue
of being very mechanical and transparent, but the
second is much easier (and easier to extend, e.g. to
multiple groups):

dat <- data.frame(x=rep(1:3,2),y=rep(1:3,2)+rnorm(6),
                 f=factor(rep(1:2,each=3)))


test1 <- function(dat) {
  fits <- lapply(split(dat,dat$f),lm,formula=y~x)
  sums <- lapply(fits,summary)
  coefs <- lapply(sums,coef)
  db <- coefs[[2]]["x","Estimate"]-coefs[[1]]["x","Estimate"]
  sd <- sqrt(sum(sapply(coefs,function(x) x["x","Std. Error"])^2))
  df <- sum(sapply(fits,"[[","df.residual"))
  td <- db/sd
  c(est=db,sd=sd,tstat=td,prt=2*pt(-abs(td),df))
}

test2 <- function(dat) {
  fit <- lm(y~x*f,data=dat)
  coef(summary(fit))["x:f2",]
}

 rbind(test1(dat),test2(dat))
           est       sd     tstat       prt
[1,] 0.5333555 1.382019 0.3859249 0.7367364
[2,] 0.5333555 1.382019 0.3859249 0.7367364
>




More information about the R-help mailing list