43
$\begingroup$

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

library(nlme) 
## 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:

m.set$coefficients
(Intercept) Petal.Width 
  4.7771775   0.9301727
m.vir$coefficients
(Intercept) Petal.Width 
  5.2694172   0.6508306 
m.ver$coefficients
(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?

Improve this question
$\endgroup$
2

3 Answers 3

Reset to default
38
$\begingroup$

To answer these questions with R code, use the following:
1. How can I test the difference between slopes?
Answer: Examine the ANOVA p-value from the interaction of Petal.Width by Species, then compare the slopes using lsmeans::lstrends, as follows.

library(lsmeans)
m.interaction <- lm(Sepal.Length ~ Petal.Width*Species, data = iris)
anova(m.interaction)
 Analysis of Variance Table

 Response: Sepal.Length
                      Df Sum Sq Mean Sq  F value Pr(>F)    
 Petal.Width           1 68.353  68.353 298.0784 <2e-16 ***
 Species               2  0.035   0.017   0.0754 0.9274    
 Petal.Width:Species   2  0.759   0.380   1.6552 0.1947    
 Residuals           144 33.021   0.229                    
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Obtain slopes
m.interaction$coefficients
m.lst <- lstrends(m.interaction, "Species", var="Petal.Width")
 Species    Petal.Width.trend        SE  df   lower.CL upper.CL
 setosa             0.9301727 0.6491360 144 -0.3528933 2.213239
 versicolor         1.4263647 0.3459350 144  0.7425981 2.110131
 virginica          0.6508306 0.2490791 144  0.1585071 1.143154

# Compare slopes
pairs(m.lst)
 contrast                 estimate        SE  df t.ratio p.value
 setosa - versicolor    -0.4961919 0.7355601 144  -0.675  0.7786
 setosa - virginica      0.2793421 0.6952826 144   0.402  0.9149
 versicolor - virginica  0.7755341 0.4262762 144   1.819  0.1669

2. How can I test the difference between residual variances?
If I understand the question, you can compare Pearson correlations with a Fisher transform, also called a "Fisher's r-to-z", as follows.

library(psych)
library(data.table)
iris <- as.data.table(iris)
# Calculate Pearson's R
m.correlations <- iris[, cor(Sepal.Length, Petal.Width), by = Species]
m.correlations
# Compare R values with Fisher's R to Z
paired.r(m.correlations[Species=="setosa", V1], m.correlations[Species=="versicolor", V1], 
         n = iris[Species %in% c("setosa", "versicolor"), .N])
paired.r(m.correlations[Species=="setosa", V1], m.correlations[Species=="virginica", V1], 
         n = iris[Species %in% c("setosa", "virginica"), .N])
paired.r(m.correlations[Species=="virginica", V1], m.correlations[Species=="versicolor", V1], 
         n = iris[Species %in% c("virginica", "versicolor"), .N])

3. What is a simple, effective way to present these comparisons?
"We used linear regression to compare the relationship of Sepal Length to Petal Width for each Species. We did not find a significant interaction in the relationships of Sepal Length to Petal Width for I. Setosa (B = 0.9), I. Versicolor (B = 1.4), nor I. Virginica (B = 0.6); F (2, 144) = 1.6, p = 0.19. A Fisher's r-to-z comparison indicated that the Pearson correlation for I. Setosa (r = 0.28) was significantly lower (p = 0.02) than I. Versicolor (r = 0.55). Similarly, the correlation for I. Virginica (r = 0.28) was significantly weaker (p = 0.02) than the one observed for I. Versicolor."

Finally, always visualize your results!

plotly_interaction <- function(data, x, y, category, colors = col2rgb(viridis(nlevels(as.factor(data[[category]])))), ...) {
  # Create Plotly scatter plot of x vs y, with separate lines for each level of the categorical variable. 
  # In other words, create an interaction scatter plot.
  # The "colors" must be supplied in a RGB triplet, as produced by col2rgb().

  require(plotly)
  require(viridis)
  require(broom)

  groups <- unique(data[[category]])

  p <- plot_ly(...)

  for (i in 1:length(groups)) {
    groupData = data[which(data[[category]]==groups[[i]]), ]
    p <- add_lines(p, data = groupData,
                   y = fitted(lm(data = groupData, groupData[[y]] ~ groupData[[x]])),
                   x = groupData[[x]],
                   line = list(color = paste('rgb', '(', paste(colors[, i], collapse = ", "), ')')),
                   name = groups[[i]],
                   showlegend = FALSE)
    p <- add_ribbons(p, data = augment(lm(data = groupData, groupData[[y]] ~ groupData[[x]])),
                     y = groupData[[y]],
                     x = groupData[[x]],
                     ymin = ~.fitted - 1.96 * .se.fit,
                     ymax = ~.fitted + 1.96 * .se.fit,
                     line = list(color = paste('rgba','(', paste(colors[, i], collapse = ", "), ', 0.05)')), 
                     fillcolor = paste('rgba', '(', paste(colors[, i], collapse = ", "), ', 0.1)'),
                     showlegend = FALSE)
    p <- add_markers(p, data = groupData, 
                     x = groupData[[x]], 
                     y = groupData[[y]],
                     symbol = groupData[[category]],
                     marker = list(color=paste('rgb','(', paste(colors[, i], collapse = ", "))))
  }
  p <- layout(p, xaxis = list(title = x), yaxis = list(title = y))
  return(p)
}

plotly_interaction(iris, "Sepal.Length", "Petal.Width", "Species")

irisPlot

Improve this answer
$\endgroup$
0
29
$\begingroup$

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$.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – Abe
    Commented Jul 25, 2012 at 20:27
  • $\begingroup$ 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. $\endgroup$
    – Macro
    Commented Jul 25, 2012 at 21:31
  • $\begingroup$ @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)? $\endgroup$
    – COOLSerdash
    Commented May 19, 2013 at 14:06
11
$\begingroup$

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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ It looks like we've posted fairly similar answers at almost the exact same time. +1 $\endgroup$
    – Macro
    Commented Jul 25, 2012 at 15:29
  • $\begingroup$ @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. $\endgroup$
    – gung - Reinstate Monica
    Commented Jul 25, 2012 at 15:31
  • $\begingroup$ 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. $\endgroup$
    – Macro
    Commented Jul 25, 2012 at 15:44
  • 3
    $\begingroup$ From my way of thinking you both should get upvotes which is what I am doing. $\endgroup$
    – Michael R. Chernick
    Commented Jul 25, 2012 at 15:58
  • 1
    $\begingroup$ The dummy variable suggestion is a good one provided the error variance does not differ appreciably among models. Otherwise you could apply a Satterthwaite-Welch t-test (which has the singular advantage of being available when only summary statistics are known, as is often the case when reading published papers) or use weighted least squares to fit the combined model. $\endgroup$
    – whuber ♦
    Commented Oct 23, 2016 at 16:00

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