Small paired samples comparison: which approach should I prefer?

Question

I have the following small dataset, that consists of scores before and after a certain treatment for 15 individuals:

df <- structure(list(before = c(4.1, 4.4, 3.7, 2.6, 4.1, 2, 3.1, 5.9, 
2.4, 6.3, 2.3, 2.3, 4.1, 3.3, 4.7), after = c(3.3, 2.3, 3.2, 
3.9, 2.4, 3.9, 0.2, 5, 3.7, 2.8, -2.6, 4.2, 2.5, 3, -1.4)), row.names = c(NA, 
15L), class = "data.frame")

Classical tests ($t$-test or Wilcoxon test), along with their associated confidence interval, do not bring sufficient evidence against the null hypothesis of no difference:

t.test(df$before, df$after, paired = TRUE)

    Paired t-test

data:  df$before and df$after
t = 2.0374, df = 14, p-value = 0.06097
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.06640299  2.58640299

wilcox.test(df$before, df$after, paired = TRUE, conf.int = TRUE)

    Wilcoxon signed rank exact test

data:  df$before and df$after
V = 90, p-value = 0.0946
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -0.2  2.6

However, a bootstrap approach seems to bring a stronger evidence (if we assume that I have correctly implemented this approach in R), with narrower 95% CI that not include zero:

my_diff <- function(data, indices) {
    data <- data[indices, ]
    return(mean(data$after - data$before))
}
set.seed(2020)
library(boot)
res <- boot(data = df, statistic = my_diff, R = 999)
boot.ci(res, type = c("bca", "perc"))

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates

CALL : 
boot.ci(boot.out = res, type = c("bca", "perc"))

Intervals : 
Level     Percentile            BCa          
95%   (-2.380, -0.153 )   (-2.381, -0.155 )  
Calculations and Intervals on Original Scale

In the case of my very small sample, which approach should be regarded as the most reliable, or at least well suited to this dataset?

Thanks!

Answer 1

There are a few issues you have going on here that area leading to these differences:

In the $t$-test you have run, you do not have the order of the before and after specified the same way as you have it specified for the bootstrap procedure (i.e. you have mixed up which term is subtracted from which). In your $t$-test, you are computing the differences as $before-after$ as opposed to $after-before$, which flips the signs and upper and lower bounds of your confidence interval. Note the differences below:

#You have computed before-after
> t.test(df$before, df$after, paired = TRUE)
    Paired t-test

data:  df$before and df$after
t = 2.0374, df = 14, p-value = 0.06097
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.06640299  2.58640299

#Rather than computing after-before
> t.test(df$after, df$before, data=df, paired = TRUE)

    Paired t-test

data:  df$after and df$before
t = -2.0374, df = 14, p-value = 0.06097
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.58640299  0.06640299
sample estimates:
mean of the differences 
                  -1.26 

To be sure, just compute differences and use a regular $t$-test

> df$diff=df$after-df$before
> t.test(df$diff)

    One Sample t-test

data:  df$diff
t = -2.0374, df = 14, p-value = 0.06097
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -2.58640299  0.06640299
sample estimates:
mean of x 
    -1.26 

Your bootstrap is using $after-before$. In addition, you've asked for Biased corrected confidence intervals and have only computed 999 bootstrap samples. If you increase the number of samples to 10,000 and ask for all bootstrap confidence intervals, you'll see your confidence intervals match better, but these confidence intervals still do not contain zero:

set.seed(2020)
library(boot)
res <- boot(data = df, statistic = my_diff, R = 10000, paired=TRUE)
boot.ci(res, type = c("all"))

> boot.ci(res, type = c("all"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = res, type = c("all"))

Intervals : 
Level      Normal              Basic         
95%   (-2.440, -0.084 )   (-2.407, -0.047 )  

Level     Percentile            BCa          
95%   (-2.473, -0.113 )   (-2.540, -0.167 )  
Calculations and Intervals on Original Scale

So the question becomes, which set of confidence intervals should you use? Personally, I'd report all confidence intervals. If you are using the confidence intervals to perform a hypothesis test, then I'd error on the side of caution and say you have mixed evidence (or perhaps no strong evidence) of a statistically significant effect depending on the test used. A thorough explanation of why you see differences in the computed confidence intervals is explained in this post. I'd advise you to read the post and the answer provided by Greg Snow. The advice he provides is quite solid, in my opinion.