Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I would like to test the difference in response of two variables to one predictor. Here is a minimal reproducible example.

## gls is used in the application; lm would suffice for this example
m.set <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "setosa")
m.vir <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "virginica")
m.ver <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "versicolor")

I can see that the slope coefficients are different:

(Intercept) Petal.Width 
  4.7771775   0.9301727
(Intercept) Petal.Width 
  5.2694172   0.6508306 
(Intercept) Petal.Width 
   4.044640    1.426365 

I have three questions:

  1. How can I test the difference between slopes?
  2. How can I test the difference between residual variances?
  3. What is a simple, effective way to present these comparisons?

A related question, Method to compare variable coefficient in two regression models, suggests re-running the model with a dummy variable to differentiate the slopes, are there options that would allow the use of independent data sets?

share|improve this question
In regards to the first question see…. – rpierce Sep 24 '14 at 19:57
up vote 18 down vote accepted

How can I test the difference between slopes?

Include a dummy for species, let it interact with $P_i$, and see if this dummy is significant. Let $L_i$ be the sepal length and $P_i$ be the pedal width and $S_1, S_2, S_3$ be the dummy variables for the three species. The compare the model

$$ E(L_i) = \beta_0 + \beta_1 P_i $$

with the model that allows the effect of $P_i$ to be different for each species:

$$ E(L_i) = \alpha_0 + \alpha_1 S_2 + \alpha_2 S_3 + \alpha_4 P_i + \alpha_5 P_iS_2 + \alpha_6 P_i S_3 $$

The GLS estimators are MLEs and the first model is a submodel on the second, so you can use the likelihood ratio test here. The likelihoods can be extracted using the logLik function and the degrees of freedom for the test will be $4$ since you've deleted $4$ parameters to arrive at the submodel.

What is a simple, effective way to present the comparison?

I think the most appealing way would be to plot the regression lines for each species all on the same axes, maybe with error bars based on the standard errors. This would make the difference (or non-difference) between the species and their relationship to $P_i$ very apparent.

Edit: I noticed another question has been added to the body. So, I'm adding an answer to that:

How can I test the difference between residual variances?

For this, you'll need to stratify the data set and fit separate models since, the interaction-based model I suggested will constraint the residual variance to be the same in every group. If you fit separate models, this constraint goes away. In that case, you can still use the likelihood ratio test (the likelihood for the larger model is now calculated by summing the likelihoods from the three separate models). The "null" model depends on what you want to compare it with

  • if you only want to test the variance, while leaving the main effects in, then the "null" model should be the model with the interactions I've written above. The degrees of freedom for the test are then $2$.

  • If you want to test the variance jointly with the coefficients, then the null model should be the first model I've written above. The degrees of freedom for the test is then $6$.

share|improve this answer
Why is there no $S_1$ in the second model? Is gls(Sepal.Length ~ species:Petal.Width, data = iris) the correct implementation of the model in R? – Abe Jul 25 '12 at 20:27
Hi @Abe. $S_1$ is the "reference" species - the regression line for that species is given by $\alpha_0 + \alpha_4 P_i$. If species is a categorical variable then I think gls(Sepal.Length ~ species*Petal.Width, data=iris) would be the syntax. – Macro Jul 25 '12 at 21:31
@Macro Nice answer (+1)! I wonder whether you could fit the gls model but allowing for different residual variances for each Species with the option weights=varIdent(form=~1|Species) (regarding the second question)? – COOLSerdash May 19 '13 at 14:06

I agree with the previous suggestion. You should fit a multiple regression model with a dummy variable for each data set. This will allow you to test whether the intercepts differ. If you also want to know if the slopes differ, then you need to also include interactions between the dummies and the variable in question. There is no problem with the fact that the data are independent. Note that if they are both independent and (for example) different species, then you would not be able to tell whether the difference you find is due to the differing species or the differing data sets, as they are perfectly confounded. However, there is no test / get-out-of-jail-free card that will get you around that problem without gathering a new sample and running your study again.

share|improve this answer
It looks like we've posted fairly similar answers at almost the exact same time. +1 – Macro Jul 25 '12 at 15:29
@Macro, yes, but yours is mostly better (+1 earlier); you addressed all 3 questions which I missed on my first (not thorough) reading of the question. My contribution here is the part about confounding. – gung Jul 25 '12 at 15:31
yes that is a good point. I suppose if you were making this inquiry at all you'd have to be operating under the assumption that the data sets were measuring the same thing, etc... with the only difference being that the species were different. – Macro Jul 25 '12 at 15:44
From my way of thinking you both should get upvotes which is what I am doing. – Michael Chernick Jul 25 '12 at 15:58

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.