3
$\begingroup$

Below are two examples of two x-y data pairs:

enter image description here

The top example clearly has more "spread" from any sort of linear trend between x and y. The data points are clustered together but there are some outliers. In contrast, the lower figure has a much tighter collection of data points closer to the linear trend. I want to quantify this spread using some sort of statistical test and determine whether x-y are "different" in a "statistically significant" way. I essentially want to be able to say:

"The top example has a statistical test result of ____ which suggests that x and y are different with a significance level of ____. In contrast, the bottom example has a statistical test result of ____ which suggests that x and y are the same within a significance level of _____."

I'm not a statistician and I am getting lost in all the different possible tests that are available.

Any ideas or suggestions are appreciated.

$\endgroup$
  • 3
    $\begingroup$ The scatterplots are a good tool for studying how $x$ and $y$ covary. But if ultimately your objective is to apply a single test of whether they are "different," it sounds like your interest is focused exclusively on $x-y.$ Why, then, not use any appropriate measure of spread of the numbers $x-y$? You scarcely need a test to compare the two datasets shown, but in more delicate situations you might consider comparing their variances or even just testing the means or medians of the absolute differences $|x-y|.$ I hesitate to offer detailed advice because there's a lot going on in the first plot! $\endgroup$ – whuber Nov 29 '18 at 20:40
3
$\begingroup$

Both plots look to be like they would fit well to a linear regression model, and in that context, it sounds to me like you want to test the error variance of the two models. This could be done by formulating your problem as a linear regression where you have a continuous explanatory variable and response variable, plus another explanatory variable that is an indicator of which of the groups you are looking at. To test equality of spread you are essentially just testing for homoscedasticity of the error variance with respect to the group indicator.

Implementation: Here is an example of one way you could do this in R. We first generate some example data and plot this, with a view to mimicking the kind of data you have described in your question. We then analyse the data via a regression analysis.

#Generate some example data for analysis
set.seed(1)
beta0 <- 0;
beta1 <- 1;
n     <- 1000;
x     <- rnorm(n, 0, 5);
g     <- ifelse(round(runif(n), 0) == 0, "Group 1", "Group 2");
sigma <- 2*(1+2*(g == "Group 2"));
e     <- rnorm(n, 0, sigma);
y     <- beta0 + beta1*x + e;
DATA  <- data.frame(y = y, x = x, g = g);

#Plot example data
library(ggplot2);
ggplot(DATA, aes(x = x, y = y)) +
   geom_point() + facet_wrap(~ g, ncol = 1) +
   theme(plot.title = element_text(hjust = 0.5, face = 'bold')) +
   ggtitle('Scatterplots of data for two groups') +
   xlab('Explanatory Variable') + ylab('Response Variable');

enter image description here

Now that we have generated the example data we can analyse this using regression analysis. We will start by fitting a standard regression model that assumes homoscedasticity (i.e., constant error variance across groups). We can formally test the heteroskedasticity by testing whether the studentised residuals of the two groups have equal variance. This is done using a standard F-test for equality of variance.

#Fit regression model
#Model allows heteroskedasticity in error term across group
MODEL <- lm(y ~ x*g, data = DATA)

#Extract studentised residuals
library(MASS);
RES <- studres(MODEL);

#Perform F-test for equality of variance across groups
library(stats);
VARTEST <- var.test(RES ~ g, alternative = 'two.sided');
VARTEST$p.value;
[1] 4.164927e-119

#Plot residuals by group (reporting variance test)
ggplot(data = data.frame(Residuals = RES, Group     = g),
       aes(x = Group, y = Residuals)) +
   geom_violin(fill = 'blue') + 
   theme(plot.title    = element_text(hjust = 0.5, face = 'bold'),
         plot.subtitle = element_text(hjust = 0.5)) +
   ggtitle('Residual Plot - Regression Model') +
   labs(subtitle = paste0('(Test for homoscedasticity across groups - p-value = ', 
                          sprintf("%.4f", VARTEST$p.value), ')')) +
   xlab(NULL) + ylab('Studentised Residual');

enter image description here

From this plot and the accompanying variance test we can immediately see that there is strong evidence of heteroskedasticity in the model (i.e., error variance differs by group). This is confirmed by the extremely low p-value of the test. As a result, if you want to model the data properly, you will need to switch to a linear regression model that allows the error variance to differ for the two groups (i.e., such a model would estimate two error variances).

$\endgroup$
  • $\begingroup$ +1, but how does MODEL allow for heteroskedasticity? OPs first scatterplot looks like it violates homoskedasticity pretty bad, which makes me think something like dglm (Smyth, 1989) would be good—OP could use group as a predictor of both conditional mean and conditional variance. $\endgroup$ – Mark White Nov 30 '18 at 3:24
  • 1
    $\begingroup$ @Mark White: MODEL assumes homoscedasticity and therefore constitutes the null hypothesis for the test of heteroskedasticity. In testing for the latter we generate the p-value under the null distribution (which assumes homoscedasticity) so it is appropriate to do this by testing the variance of residuals from a homoscedastic model. If you build the heteroskedasticity into the model then the studentised residuals would have equivalent variance by virtue of the fitting mechanism. $\endgroup$ – Reinstate Monica Nov 30 '18 at 3:34
  • $\begingroup$ @Mark White: You are certainly right that to fix the model you would need to run a different model (that incorporates the heteroskedasticity). The one I have used is just for the purposes of hypothesis testing. (I have added an additional line to the end of the answer to make this clearer.) $\endgroup$ – Reinstate Monica Nov 30 '18 at 3:35
  • 1
    $\begingroup$ Apologies, I read too quickly. Wish I could upvote this answer more than once—extracting that assume homoskedasticity and doing a test on them is brilliant. $\endgroup$ – Mark White Nov 30 '18 at 3:37
  • $\begingroup$ Thanks Mark! That's very flattering! Feel free to upvote some of my other answers to quench your desire to upvote more than once. ;) $\endgroup$ – Reinstate Monica Nov 30 '18 at 3:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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