Chapter 13: Comparing Two Samples

Section A.12 Chapter 13: Comparing Two Samples

Subsection A.12.1 Exercise 13.1

The score of pain before the application of the device is measured in the variable “score1". This variable is used as the response. We apply the function “t.test" in order to test the equality of the expectation of the response in the two groups. First we read in the data from the file into a data frame and then we apply the test:
🔗
```
magnets <- read.csv("_data/magnets.csv")
t.test(magnets$score1 ~ magnets$active)
```
```
## 
##  Welch Two Sample t-test
## 
## data:  magnets$score1 by magnets$active
## t = 0.41483, df = 38.273, p-value = 0.6806
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.3757896  0.5695498
## sample estimates:
## mean in group "1" mean in group "2" 
##           9.62069           9.52381
```
The computed $p$-value is 0.6806, which is above 0.05. Consequently, we do not reject the null hypothesis that the expectations in the two groups are equal. This should not come as a surprise, since patients were assigned to the groups randomly and without knowledge to which group they belong. Prior to the application of the device, the two groups should look alike.
🔗

🔗
Again, we use the variable “score1" as the response. Now apply the function “var.test" in order to test the equality of the variances of the response in the two groups:
🔗
```
var.test(magnets$score1 ~ magnets$active)
```
```
## 
##  F test to compare two variances
## 
## data:  magnets$score1 by magnets$active
## F = 0.69504, num df = 28, denom df = 20, p-value = 0.3687
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2938038 1.5516218
## sample estimates:
## ratio of variances 
##          0.6950431
```
The computed $p$-value is 0.3687, which is once more above 0.05. Consequently, we do not reject the null hypothesis that the variances in the two groups are equal. This fact is reassuring. Indeed, prior to the application of the device, the two groups have the same characteristics. Therefore, any subsequent difference between the two groups can be attributed to the difference in the treatment.
🔗

🔗
The difference in score between the treatment and the control groups is measured in the variable “change". This variable is used as the response for the current analysis. We apply the function “t.test" in order to test the equality of the expectation of the response in the two groups:
🔗
```
t.test(magnets$change ~ magnets$active)
```
```
## 
##  Welch Two Sample t-test
## 
## data:  magnets$change by magnets$active
## t = 5.9856, df = 42.926, p-value = 3.86e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  2.749137 5.543145
## sample estimates:
## mean in group "1" mean in group "2" 
##          5.241379          1.095238
```
The computed $p$-value is $3.86\times 10^{-7}\text{,}$ which is much below 0.05. Consequently, we reject the null hypothesis that the expectations in the two groups are equal. The conclusion is that, according to this trial, magnets do have an effect on the expectation of the response[^9].
🔗

🔗
Once more we consider the variable “change" as the response. We apply the function “var.test" in order to test the equality of the variances of the response in the two groups:
🔗
```
var.test(magnets$change ~ magnets$active)
```
```
## 
##  F test to compare two variances
## 
## data:  magnets$change by magnets$active
## F = 4.2062, num df = 28, denom df = 20, p-value = 0.001535
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  1.778003 9.389902
## sample estimates:
## ratio of variances 
##           4.206171
```
The computed $p$-value is 0.001535, which is much below 0.05. Consequently, we reject the null hypothesis that the variances in the two groups are equal. Hence, magnets also affect the variance of the response.
🔗

🔗

Subsection A.12.2 Exercise 13.2

We simulate the sampling distribution of the sample standard deviation for two samples, one sample of size $n_a = 29$ and the other of size $n_b=21\text{.}$ Both samples are simulated from the given Normal distribution:
🔗

🔗

mu <- 4
sig <- 4
n.a <- 29
n.b <- 21
S.a <- rep(0,10^5)
S.b <- rep(0,10^5)
for(i in 1:10^5) {
X.a <- rnorm(n.a,mu,sig)
X.b <- rnorm(n.b,mu,sig)
S.a[i] <- sd(X.a)
S.b[i] <- sd(X.b)
}
F <- S.a^2/S.b^2
mean((F < qf(0.025,n.a-1,n.b-1))|(F > qf(0.975,n.a-1,n.b-1)))

## [1] 0.04949

We compute the test statistic “F" as the ratio of the two sample standard deviations “S.a" and “S.b". The last expression computes the probability that the test statistic is either less than “qt(0.025,n.a-1,n.b-1)", or it is larger than “qt(0.975,n.a-1,n.b-1)". The term “qt(0.025,n.a-1,n.b-1)" is the 0.025-percentile of the $F$-distribution on 28 and 20 degrees of freedom and the term “qt(0.975,n.a-1,n.b-1)" is the 0.975-percentile of the same $F$-distribution. The result of the last expression is the actual significance level of the test.

🔗

We obtain that the actual significance level of the test when the measurements are Normally distributed is 0.05074, which in agreement with the nominal significance level of 5%. Indeed, the nominal significance level is computed under the assumption that the distribution of the measurement is Normal.

🔗

We repeat essentially the same simulations as before. We only change the distribution of the samples from the Normal to the Exponential distribution:
🔗
```
lam <- 1/4
n.a <- 29
n.b <- 21
S.a <- rep(0,10^5)
S.b <- rep(0,10^5)
for(i in 1:10^5) {
X.a <- rexp(n.a,lam)
X.b <- rexp(n.b,lam)
S.a[i] <- sd(X.a)
S.b[i] <- sd(X.b)
}
F <- S.a^2/S.b^2
mean((F < qf(0.025,n.a-1,n.b-1))|(F > qf(0.975,n.a-1,n.b-1)))
```
```
## [1] 0.27685
```
The actual significance level of the test is $0.27596\text{,}$ which is much higher than the nominal significance level of 5%.
🔗

Through this experiment we may see that the $F$-test is not robust to the divergence from the assumed Normal distribution of the measurement. If the distribution of the measurement is skewed (the Exponential distribution is an example of such skewed distribution) then the application of the test to the data may produce unreliable conclusions.
🔗

🔗

🔗

Subsection A.12.3 Exercise 13.3

We input the data to R and then compute the estimate:
🔗
```
s.a <- 13.4
s.b <- 16.7
s.a^2/s.b^2
```
```
## [1] 0.6438381
```
The estimate is equal to the ratio of the sample variances $s_a^2/s_b^2\text{.}$ It obtains the value $0.6438381\text{.}$ Notice that the information regarding the sample averages and the sizes of the sub-samples is not relevant for the point estimation of the parameter.
🔗

🔗
We use the formula:
🔗

\begin{equation*} \big[(s_a^2/s_b^2)/\mbox{\texttt{qf(0.975,14,14)}} , (s_a^2/s_b^2)/\mbox{\texttt{qf(0.025,14,14)}}\big] \end{equation*}

in order to compute the confidence interval:

🔗
```
n.a <- 15
n.b <- 15
(s.a^2/s.b^2)/qf(0.975,n.a-1,n.b-1)
(s.a^2/s.b^2)/qf(0.025,n.a-1,n.b-1)
```
```
## [1] 0.2161555
## [1] 1.917728
```
The confidence interval we obtain is $[0.2161555, 1.917728]\text{.}$
🔗

🔗
The estimate of the parameter is not affected by the change in the sample sizes and it is still equal to $0.6438381\text{.}$ For the confidence interval we use now the formula:
🔗

\begin{equation*} \big[(s_a^2/s_b^2)/\mbox{\texttt{qf(0.975,149,149)}} , (s_a^2/s_b^2)/\mbox{\texttt{qf(0.025,149,149)}}\big]\,: \end{equation*}

🔗
```
n.a <- 150
n.b <- 150
(s.a^2/s.b^2)/qf(0.975,n.a-1,n.b-1)
(s.a^2/s.b^2)/qf(0.025,n.a-1,n.b-1)
```
```
## [1] 0.466418
## [1] 0.8887467
```
The corrected confidence interval is $[0.466418, 0.8887467]\text{.}$
🔗

🔗

🔗

Prev Top Next