Comparing Two Samples |
Often we wish to compare two samples, in order to make inferences about possible differences between the two sampled populations, or differences between subjects in different experimental conditions. It is important to remember that distributions may differ in many respects: location, spread, and shape. The standard t-test assumes that the two distributions are identical in every way except possibly for differences between mean values.
Before testing the difference between means, it is sensible to investigate the possibility that the two samples might differ in some other respect. The quickest way to compare all aspects of two samples is with a qqplot, though if you are comparing several samples you may also want to make simultaneous boxplots.
If the two samples do not appear to have the same shape, or the same spread, you should think hard about why you are using a hypothesis test to compare their means. For example, if one distribution has a standard deviation 10, and the other 5, then even if there is no difference between the means the first distribution will have much greater probability of getting a value more than 10 units above the common mean (15.8% vs 2.3%).
Supposing that you have decided to test the null hypothesis that there is no difference between the means of the two groups, then there are several options, depending on the details of the situation.
The data used in the following examples is from a longitudinal study of children in Berkeley, CA. Follow the link to the Berkeley dataset, and paste in the entire "Berkeley" data structure - everything between the horizontal lines.
Now, let's compare the mean height in centimeters of two year old boys and girls at age 2 using the Berkeley dataset. The following commands select out the heights at age two for males and females.
male <- Berkeley[,"sex"] == "Male"
m.ht2 <- Berkeley[male,"ht2"]
f.ht2 <- Berkeley[!male,"ht2"]
The exclamation point in front of "male" in the last line means
"not": "not male" is of course female.
As usual, there are other ways to separate the data into the two
groups. If you print the whole dataset, you will see that the
first 26 cases are males, and the remaining cases are females.
You could use that information to select the two subsets, but
it is good to have the more general tool for selecting from files
where the subsets aren't already segregated.
Look at the qqplot for the two groups. Try to predict the outcome of the t-test before you continue.
You might wish to plot the confidence intervals for the two groups. See this script. Copy it into an R document, then substitute your variables for x and y everywhere.
Now let's run the t-test:
t.test(m.ht2,f.ht2,var.equal=TRUE)
If all is well, the output should look like this:
Two Sample t-test data: m.ht2 and f.ht2 t = 1.0709, df = 56, p-value = 0.2888 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.8135485 2.6822985 sample estimates: mean of x mean of y 88.40000 87.46563You can see that there is hardly any difference between the two sample means - males were on average one cm. taller in the sample. The t-statistic is 1.08. The p-value is big: about 0.29, and the 95 percent confidence interval includes 0. In other words, the observed difference could easily occur due solely to sampling variation - random noise. We say that the difference is "not statistically significant" as a shorthand for "not big enough to rule out chance variation as an explanation".
power.t.test(n=26,sd=3.5,delta=1)
This is equivalent to
power.t.test(n=26,sd=3.5,delta=1,type="two.sample")
Where did "sd=3.5" come from? In general, before we do an experiment, we
don't know what the standard deviation of the data will be! Here
we already have a sample, so we can use the standard deviation based
on our dataset. When designing an experiment, it is often helpful to
look at the published literature in the field to get a sense of the
level of variability to expect.
How big a sample would we need to achieve power of .9?
power.t.test(sd=3.5,delta=1,power=.9,type="two.sample")
The power.t.test function takes the given list of values
(n,sd,delta,sig.level,power) and computes whichever
ones are missing without a default specification. "sd" defaults
to 1, and "sig.level" defaults to .05, so you need to explicitly
set them to "NULL" if you wish to have the function compute them.
Remember that the "n" computed is the sample size for each group,
not the total.
wilcox.test(m.ht2,f.ht2)
Wilcoxon rank sum test with continuity correction
data: m.ht2 and f.ht2
W = 503, p-value = 0.1756
alternative hypothesis: true mu is not equal to 0
Warning message:
Cannot compute exact p-value with ties in: wilcox.test.default(m.ht2, f.ht2)
The warning is not important - it simply announces that the
p-value is an approximate p-value, not an exact p-value. That is
the least of our worries!
Note that the alternative hypothesis is stated as though the test is specific for comparing the means of the two distributions: "true mu is not equal to 0". Be careful here! The null hypothesis can be rejected because of differing shape even when the two distributions have the same mean. It is always important to check the assumptions underlying the test!
randTest <- rep(0,10000)
x <- c(m.ht2,f.ht2)
for(i in 1:10000)
{
y <- sample(x)
randTest[i] <- mean(y[1:26])-mean(y[27:58])
}
### compute the p-value ###
sum(abs(randTest) > 0.934375)/10000
boot <- rep(0,10000)
for(i in 1:10000)
{
m <- sample(m.ht2,replace=T)
f <- sample(f.ht2,replace=T)
boot[i] <- mean(m)-mean(f)
}
The 95% confidence interval for the difference between the means is
quantile(boot,c(.025,.975))
Compare the bootstrap CI to the standard t based CI (-0.8135, 2.682),
and to the Welch CI (-0.7876, 2.656);
it should be similar.