Case Study: Comparing Two Groups

Corresponding reading

  • Chapter 4.5 of Agresti and Kateri (2021).
  • Chapter 11 of Rice (2006) (11.3 in particular for paired samples, not covered in Agresti and Kateri (2021))

So far we have only considered creating confidence intervals for the mean of a single population. Now we are going to consider what if we have two different populations.

Consider the following examples:

  1. We are interested in the difference in wages between men and women. We collect a random sample of men and women to and ask them their annual income.
  2. We are interested the effect of rewards on behavior. We randomly select laboratory mice, with a treatment group getting rewards for successfully accomplishing a task and the control group receiving nothing.
  3. We are interested in the effect of a weight loss drug. We collect the weight of subjects before and after they use the drug.
  4. We want to compare political opinions between generations. We randomly recruit families and recruit two members to answer a survey, one of age 18-30 and the other of age 45-60.

These examples all have the common property that we are trying to compare two populations, rather than estimate something about one population.

Questions/Goals for this Module

  • How to model comparisons between two groups
  • Extend CI concepts to two groups.
  • Understand the statistical differences between paired and unpaired data

Modeling two groups

We can consider that we have data from two populations \[(X_1,\ldots,X_{n_X}) \text{ and } (Y_1,\ldots, Y_{n_Y})\] Since these are two different populations, we consider that they each have a different mean, but we’ll assume that within each group their means are the same, \[E(X_i)=\mu_X, E(Y_i)=\mu_Y\]

NoteNotation

I am deviating from the book’s notation here. The book has group 1 be \(Y_{11},\ldots,Y_{1n_1}\) and the second group \(Y_{21},\ldots,Y_{2n_2}\), with means \(\bar{Y}_1\) and \(\bar{Y}_2\). The book’s notation fits into our regression notation better, but I find from past experience that students find double indices harder to follow – and they are harder to write on the board!

What are common additional modeling choices that we might make?

  1. What do we estimate? The important question with all these parameters floating around is what are we interested in estimating? We could be interested in \(\mu_X\) and \(\mu_Y\) separately. But in many two group comparisons, we are primarily interested in comparing the two groups. This usually translates into evaluating how different these means are from each other. So we might define another parameter \[\Delta=\mu_X-\mu_Y.\] This is an example of creating a parameter of interest that isn’t explicitly a parameter of our distribution, but rather a function of our parameters.

  2. Dependence Structure: We could assume that within each group the \(X_i\) are independent from each other and the \(Y_i\) are independent from each other. This would make sense if each observation in the group comes from a different individual. What about between the \(X_i\) and \(Y_i\)? This depends on the data we are considering.

    In our first two examples, each \(X_i\) and \(Y_j\) are from separate individuals/units. It seems reasonable to assume that each response \(X_i\) is independent from each \(Y_j\). However, in the last two, the \(X_i\) and \(Y_i\) are measurements that are linked in some way: in example 3 measurement \(X_j\) is the before and \(Y_j\) the after on the same individual; in example 4, \(X_j\) is the answer of the younger person and \(Y_j\) is the answer of the older person in the same family. It does not make sense to assume they are independent. Knowing the weight of an individual before treatment gives you information about what their weight might be after treatment (you have different expectations of the final weight for someone who starts out 180lbs versus 240lbs).

    This is called paired data, because the data comes in the form of pairs (and thus \(n_X=n_Y\)). In fact we could instead think of the data as pairs indexed by the same value: \((X_1,Y_1),\ldots,(X_n,Y_n)\). In which case it might make sense to think of the pairs \((X_i,Y_i)\) as independent from each other – i.e. to rewrite our data as observations \(i\) which consist of a vector \((X_i,Y_i)\). However, note that the \(ith\) pair wouldn’t necessarily correspond to the same unit. For example, in our example 4, \(X_i\) and \(Y_i\) are from two different people but they are from the same household \(i\).

  3. We could also consider what is an appropriate parametric model for our data. If our \(X\) and \(Y\) data are both continuous, we are likely to choose a normal distribution, \[X_i \overset{i.i.d}{\sim}N(\mu_X,\sigma^2_X), Y_i \overset{i.i.d}{\sim}N(\mu_Y,\sigma^2_Y).\]

    If our data are binary, then we would often use a Bernoulli model, \[X_i \overset{i.i.d}{\sim}Be(\pi_X), Y_i \overset{i.i.d}{\sim}Be(\pi_Y).\]

    Notice that this doesn’t address correlation between the \(X_i\) and \(Y_i\), which we might have to consider if we have paired data. We might, for example, assume that they are bivariate normal, with covariance, \[cov(X_i,Y_i)=\sigma_{XY}\] However, we will see for paired data, we can avoid needing to explicitly model (or estimate) \(\sigma_{XY}\).

Relationship of these models to Regression

We have already discussed the example of two groups implicitly, when we looked at regression for the case of where our predictor \(X\) was binary. We had, \[E(Y |X)= \begin{cases} \alpha_1, & X=1\\ \alpha_0, & X=0 \end{cases} \] So in this chapter, we’ve changed the notation (\(X\) means something different in the regression frame work), but it is ultimately describing the same setting.

We also previously introduced a parametric for independent \(Y_i\) and \(X_i\) categorical, \[Y |X \sim \begin{cases} N(\alpha_1,\sigma_1^2), & X=1\\ N(\alpha_0,\sigma_0^2), & X=0 \end{cases} \] This is the standard model for independent, non-paired data, and matches our model above. This model, however, didn’t handle dependency between our data, as we have with paired data.

Moreover, if instead we consider the standard regression model for linear models, \[Y_i=\beta_0+\beta_1 X_i +\epsilon_i\] with \(\epsilon_i\) i.i.d. \(N(0,\sigma^2)\), then it doesn’t even allow for different variances in the two groups. So while the general notation and framework of regression works for two groups, the specific analysis we usually do with the standard model doesn’t match what is usually done for analyzing even two independent groups.

Exercise/Question

Convert our two-group modeling notation above into the regression notation. What is \(X\) and \(Y\) in each of our examples in the regression notation rather than the 2-group notation?

Independent Groups

Estimation

In this case, if the two populations are independent from each other, we’ve actually already calculated (in section) the MLE assuming normality of the parameters as a special case of the regression model, \[ \begin{gathered} \hat\mu_X=\bar{X},\quad \hat\mu_Y=\bar{Y}\\ \hat\sigma_X^2=\frac{1}{n_X}\sum_{i=1}^{n_X} (X_i-\bar{X})^2,\quad \hat\sigma_Y^2=\frac{1}{n_Y}\sum_{i=1}^{n_Y} (Y_i-\bar{Y})^2 \end{gathered} \] However, as is standard, we will instead use the estimates \[S^2_X=\frac{1}{n_X-1}\sum_{i=1}^{n_X} (X_i-\bar{X})^2, \quad S^2_Y=\frac{1}{n_Y-1}\sum_{i=1}^{n_Y} (Y_i-\bar{Y})^2\]

We also can state the distributions of \(\hat{\mu}_X\) and \(\hat{\mu}_Y\), either exactly – if we assume normality – or approximately if we rely on the central limit theorem, \[\hat{\mu}_X\sim N(\mu_X,\sigma^2_X/n_X), \quad \hat{\mu}_Y\sim N(\mu_Y,\sigma^2_Y/n_Y)\]

What about \(\Delta\)? The obvious estimator for their difference is1: \[\hat\Delta=\hat\mu_X-\hat\mu_Y\]

What is the distribution of \(\hat\Delta\)? We’ve already seen that linear combinations of normal variables are normal. That is also true for asymptotic distributions:

Theorem 1 If \(X\) and \(Y\) both converge asymptotically to a normal distribution, then \(aX+bY\) will also converge to asymptotically to a normal distribution.

So depending on what assumptions we want to make, we have either

  1. We assume \(X\) and \(Y\) are both i.i.d. normal, then \(\hat\Delta\) is exactly normally distributed

or

  1. Our \(n_X\) and \(n_Y\) are sufficiently large so that \(\hat\Delta\) is approximately normal by the CLT theorem.

Therefore, since \(\hat\Delta\) is normal (or approximately so), we just need to figure out what should be the mean and variance of that normal. Clearly the expectation is just, \[E(\hat\Delta)=\mu_X-\mu_Y\] The variance is given by \[var(\hat\Delta)=var(\bar{X})+var(\bar{Y})-2cov(\hat\mu_X,\hat\mu_Y)\] So for independent groups, \(cov(\hat\mu_X,\hat\mu_Y)=0\) so we have that the variance of \(\hat{\Delta}\) is, \[\sigma^2_{\hat\Delta}=var(\hat\Delta)=\frac{\sigma_X^2}{n_X}+\frac{\sigma_Y^2}{n_Y}\] We can estimate that variance by substituting in \(S^2_{X}\) and \(S^2_{Y}\) for \(\sigma^2_X\) and \(\sigma^2_Y\), which gives us an estimate of the standard error, \[S_{\hat\Delta}=\sqrt{\frac{S_X^2}{n_X}+\frac{S_Y^2}{n_Y}}\]

CI for Two-Sample Comparisons

Just like our mean, we can consider the quantity \[T_{\Delta}=\frac{\hat\Delta-\Delta}{S_{\hat{\Delta}}}=\frac{(\bar{X}-\bar{Y})-(\mu_X-\mu_Y)}{\sqrt{\frac{S_X^2}{n_X}+\frac{S_Y^2}{n_Y}}}\] This is the same form as our previous T statistic: the difference of our estimator and the parameter, divided by the estimate of standard error of the estimator.

The same logic holds as in the one sample case:

  1. If \(n\) is large, then \(T_{\Delta}\) converges to a \(N(0,1)\) distribution, \[T_{\Delta} \overset{D}{\rightarrow} N(0,1).\] and we have a CI \[\hat{\Delta}\pm z_{\alpha/2} S_{\hat{\Delta}}\]
  2. If \(n\) is small, and we are willing to assume \(X\) and \(Y\) are normally distributed, then we will use a \(t\)-distribution for \(T_{\Delta}\) and get a CI \[\hat{\Delta}\pm t_{\alpha/2,\nu} S_{\hat{\Delta}}\] \(\nu\) is a complicated function of \(n_X\) and \(n_Y\) as well as the \(S_X\) and \(S_Y\), but can often be substituted with \(min(n_X-1,n_Y-1)\) if one were to actually calculate the values by hand (rather than use a built-in function).2 This is a conservative choice, meaning it is the smallest possible value that \(\nu\) could be, so it gives you the largest possible \(t_{\alpha/2,\nu}\) (and thus the widest intervals)
Note

Technically \(T\) does not follow the \(t\) distribution exactly, even if the data are exactly normaly distributed. Rather, the use of of \(t_\nu\) gives an approximation to the distribution of \(T\), called the Welch’s t-interval.

Example

A drug company has developed a new drug to reduce hypertension. A sample of 15 patients are chosen to test the drug. Their systolic BP is reduced by 28.3 mmHg with a sample SD of 12 mmHg. The control group of 20 patients takes the standard drug and the average BP is reduced by 17.1 mmHg with the sample SD of 11.84 mmHg. Find a 95% CI for the difference in mean reduction of systolic BP, and conduct a two sided hypothesis test for the treatment making no difference.

The samples are small enough that to get a confidence interval, we are going to need to assume that \(X\), \(Y\) are normal and rely on the \(t\)-distribution.

We’ll let the treatment group be the \(Y\), so we have \(n_Y = 15\), and the control group be \(X\), \(n_X = 20\).

  • Average systolic BP for the treatment group reduced by \(\bar{x} = 28.3\) mmHg, \(S_x = 12\) mmHg.
  • Average systolic BP for the control group reduced by \(\bar{y} = 17.1\) mmHg, \(S_y = 11.84\) mmHg.

\[\hat{\Delta}=28.3-17.1=11.2\] \[S_{\hat\Delta} = \sqrt{\frac{S_X^2}{n_X}+\frac{S_Y^2}{n_Y}}=\sqrt{(12)^2/20+(11.84)^2/15}\approx 4.06\] If we want to do our 95% confidence interval by hand, we get \(df=min(15-1,20-1)=14\) which gives us

qt(1-0.05/2,14)
[1] 2.144787

and our confidence interval as \[\hat{\Delta}\pm t_{\alpha/2,\nu} S_{\hat{\Delta}}=11.2\pm 2.14 \cdot 4.06=11.2\pm 8.71\] The more careful calculation of \(\nu\), our degrees of freedom gives

[1] 30.42119

which gives

[1] 2.041091

\[\hat{\Delta}\pm t_{\alpha/2,\nu} S_{\hat{\Delta}}=11.2\pm 2.04 \cdot 4.06=11.2\pm 8.28\] So, as expected, using our “by hand” version of \(min(n_X-1,n_Y-1)\) is more conservative, meaning the CI are a bit wider than the more complicated approximation would give.

Assuming Equal Variance

If we assume equal variance \(\sigma^2\) in the \(X\) and \(Y\) groups, then you showed in homework that the MLE of \(\mu_X\) and \(\mu_Y\) doesn’t change, and the MLE of \(\sigma^2\) is \[\hat\sigma^2_{pooled}=\frac{n_X}{n_X+n_Y}\hat\sigma_X^2+\frac{n_Y}{n_X+n_Y}\hat\sigma_Y^2=\frac{1}{n_X+n_Y}\left(\sum_{i=1}^{n_X}(X_i-\bar{X})^2+\sum_{i=1}^{n_Y}(Y_i-\bar{Y})^2\right)\] a weighted combination of the MLE in each group separately.

Unsurprisingly, there is an alternative form of the estimator based on \(S_X\) and \(S_Y\), \[ \begin{aligned} S^2_{pooled}&=\frac{n_X-1}{(n_X-1)+(n_Y-1)}S_X^2+\frac{n_Y-1}{(n_X-1)+(n_Y-1)}S_Y^2\\ &= \frac{1}{n_X+n_Y-2}\left(\sum_{i=1}^{n_X}(X_i-\bar{X})^2+\sum_{i=1}^{n_Y}(Y_i-\bar{Y})^2\right).\end{aligned}\]

This is called a pooled estimate of variance. It can be interpreted as estimating the variance in each group and then combining them for an estimate of their shared variance. But we can also see that it is basically combining the data together to get a single estimate of standard deviation across all the observations pooled together.

Using this estimate, gives us a \(T\) statistic with a \(t\) distribution.

\[ T_{\Delta,pooled} = \frac{\hat\Delta-\Delta}{\hat\sigma_{\hat\Delta,pooled}}=\frac{(\bar{X}-\bar{Y}) - (\mu_X - \mu_Y)}{S_{pooled}\sqrt{\tfrac{1}{n_X}+\tfrac{1}{n_Y}}} \sim t_{n_X+n_Y-2} \] It’s distribution is clear from our previous theorems in the one-sample case (and unlike the unequal variance case, is exact if we assume normality of our data).

Theorem 2 If we have two independent samples \(X_1, \ldots, X_{n_X} \overset{i.i.d.}{\sim} N(\mu_X, \sigma_X^2)\), \(Y_1, \ldots, Y_m \overset{i.i.d.}{\sim} N(\mu_Y, \sigma_Y^2)\), and \(\sigma_X^2 = \sigma_Y^2\), , then \[ T_{\Delta,pooled} = \frac{(\bar{X}-\bar{Y}) - (\mu_X - \mu_Y)}{S_{pooled}\sqrt{\tfrac{1}{n_X}+\tfrac{1}{n_Y}}} \sim t_{n_X+n_Y-2} \]

This gives the confidence interval

\[(\bar{X}-\bar{Y}) \pm t_{\alpha/2,\; n_X+n_Y-2} \cdot \hat\sigma_{\hat{\Delta},pooled}\]

When assume equal variances?

Both of the books only present the case when the variances are equal. It does make the mathematics of the \(t\) distribution better, in the sense that the distribution of \(T_{\Delta,pooled}\) is exactly \(t\) distributed with degrees of freedom that we can easily write down and work with. The distribution of our \(T\) random variable is messy when we don’t assume equal variance, and we have to approximate it’s distribution with a \(t\). Moreover, in all realistic settings, you are measuring the same type of feature on both groups, so they should be on the same scale, and so it is not completely implausible that the variances are equal.

However, it is also the case that the that the same thing that divides them into two groups (men vs women; treatment vs control; etc) also means that may have a more variable response in one group than another. So the short answer is that there’s rarely any good reason to assume the variances are equal, and we have standard strategies that we don’t have to make that assumptions. In practice, it is rare that practitioners assume equal variance.

On the other hand, if we had more than one group we were comparing, the standard analysis (ANOVA) routinely assumes the same variance, just like we do in regression. And the pooled estimate of variance is what is used in this case. So it is useful to be familiar with the concept of pooled variance, in particular.

Exercise/Question

Show that the estimate of variance in regression when \(X\) is categorical reduces to \(S^2_{pooled}\) given above.

Example Redo

Let’s go back to our example, and consider what we would have if we assumed the same variance.

We need to calculated the pooled estimate of standard deviation, \[S_{pooled}^2 = \frac{(n-1)S_x^2 + (m-1)S_y^2}{m+n-2} = \frac{(14)(12^2) + (19)(11.84)^2}{15+20-2} \approx 141.80 \;\Rightarrow\; S_p \approx 11.91\] This is compared to the individual estimate in each group, \(S_X=12\) and \(S_y=11.84\). As expected, it is between these two estimates.

This gives us the estimate of standard error of \(\hat\Delta\) as \[S_{\hat\Delta,pooled}=S_{pooled}\sqrt{\tfrac{1}{n_X}+\tfrac{1}{n_Y}}=11.91*\sqrt{1/20+1/15}=11.91*0.34=4.07\] When we allowed for different variances we got \(S_{\hat\Delta}=4.06\), so quite similar. This isn’t surprising, since the estimates of standard deviation in each group are quite close.

Our degrees of freedom is \(n_X+n_Y-2=33\), again, close to \(\nu\) above.

qt(1-0.05/2,33)
[1] 2.034515

our confidence interval as \[\hat{\Delta}\pm t_{\alpha/2,\nu} S_{\hat{\Delta},pooled}=11.2\pm 2.03 \cdot 4.07=11.2\pm 8.28\] So we basically get all of the same answers as assuming different variances, which is expected when you start with two groups with similar standard deviations. However, with larger differences in \(S_X\) and \(S_Y\) and greater differences in sample size, you can get larger differences.

Paired Data

We have considered the case of two independent groups. But what about our examples of when our data is paired? Now we have dependency between our \(X\) data and our \(Y\) data.

The natural way of thinking about our data is to look at differences \[D_i=X_i-Y_i\] Then \[ \begin{aligned} E(D_i)&=\mu_X-\mu_Y\\ \sigma^2_D=var(D_i)&=var(X_i)+var(Y_i)-2 cov(X_i,Y_i)\\ &= \sigma_X^2+\sigma_Y^2-2 \rho \sigma_X\sigma_Y \end{aligned} \] where \(\rho\) is the correlation of \(X_i\) and \(Y_i\), if we assume the correlation is the same for all pairs \(i\).

Estimation

We still want to estimate \(\Delta=\mu_X-\mu_Y\). How does our estimate of \(\Delta\) change if we have paired data?

First we note that \[\bar{X}-\bar{Y}=\frac{1}{n}\sum_{i=1}^n(X_i-Y_i)=\bar{D}\] So in fact our previous estimator \(\hat{\Delta}=\bar{X}-\bar{Y}=\bar{D}\) will also be the average of their pair-wise differences.

Similarly, \[E(\hat{\Delta})=E(\bar{D})=E(\bar{X})-E(\bar{Y})=\mu_X-\mu_Y\] So we are clearly still using the same estimator and it is still an unbiased estimator of \(\Delta\).

However, the variance changes. For paired data \[\sigma^2_{\bar{D}}=var(\bar{D})=var(\bar{X}-\bar{Y})=\frac{1}{n}(\sigma^2_X+\sigma^2_Y-2\rho \sigma_X\sigma_Y)\] If we had instead collected the data on independent groups of units (assuming the same size \(n\) in each group) we have, \[var(\hat\Delta)=var(\bar{X}-\bar{Y})=\frac{1}{n}(\sigma_X^2+\sigma^2_Y)\]

So you can see that if there is a strong correlation between the \(X_i\) and \(Y_i\) the variance of your estimator will be quite different. In fact, if there is a strong positive correlation, then the variance of \(\hat\Delta\) will be much smaller, and thus you will have a more accurate estimate of \(\mu_X-\mu_Y\).

We frequently expect positive correlation – the pairing is usually of units that are similar to each other (or the same) so when \(X_i\) is larger on average than the population \(X\), the same is true of \(Y_i\). For this reason, if you can construct a paired study, rather than independent observations in each group, this is usually a better design for how to collect your data. This makes sense if you think about how we can write our paired analysis as an analysis of the differences; we remove all of the things that stay the same for the pair and only measure what has changed. This removes variability from other sources. Think of our example of members of the same household – they have the same environment, the same standard of living, the same familial history. So looking at the difference of a younger and older member of the household can focus in on the generational effect.

CI

Again, we will need to estimate \(\sigma^2_{\bar{D}}\). Since we are just talking about the standard error of the mean of \(D_i\), we will use our standard estimate of the standard deviation of the differences, \(\sigma_{D}^2\), \[S_{D}^2=\frac{1}{n-1}\sum_{i=1}^n(D_i-\bar{D})^2=\frac{1}{n-1}\sum_{i=1}^n((X_i-Y_i)-(\bar{X}-\bar{Y}))^2\] We then define \[T_{D}=\frac{\bar{D}-\mu_D}{S_{\bar{D}}}=\frac{(\bar{X}-\bar{Y})-(\mu_X-\mu_Y)}{S_D/\sqrt{n}}\] If we assume the \(X_i\) and \(Y_i\) are each normally distributed, then so is \(D_i\), so we have \[T_{D}\sim t_{n-1},\] and a CI given by \[\bar{D}\pm t_{\alpha/2,n-1}S_{\bar{D}}.\] Alternatively, for large sample sizes \(n\), we can use the CLT and have \[\bar{D}\pm z_{\alpha/2}S_{\bar{D}}.\]

Note

Notice how everything is exactly the same as dealing with the mean of a single population, \(\bar{X}\), only we have now just replaced \(X_i\) with \(D_i=X_i-Y_i\).

References

Agresti, Alan, and Maria Kateri. 2021. Foundations of Statistics for Data Scientists. Chapman; Hall/CRC.
Rice, John A. 2006. Mathematical Statistics and Data Analysis. 3rd ed. Duxbury Press.

Footnotes

  1. Indeed, since the MLE of a function of parameter is just that function of the MLE, \(\hat\Delta\) is indeed the MLE of \(\Delta\).↩︎

  2. \(\nu = \frac{\left(\frac{s_X^2}{n_X} + \frac{s_Y^2}{n_Y}\right)^2}{\frac{(s_X^2 / n_X)^2}{n_X - 1} + \frac{(s_Y^2 / n_Y)^2}{n_Y - 1}}\) and is bounded between \(\min(n_X-1, n_Y-1)\) and \((n_X + n_Y - 2)\). ↩︎