The goal of bstrap is to obtain and plot a bootstrap sampling distribution of the sample mean. The idea is that if this is close to normally distributed, normal-theory tests such as t-tests should be appropriate (and if not, not).
install.packages("devtools") # if not already installed
library(devtools)
install_git("http://worktree.ca/nxskok/bstrap.git")
library(bstrap)
library(ggplot2)
library(tibble)
In the mtcars dataset, the variable disp has a right-skewed
distribution:
ggplot(mtcars, aes(x = disp)) + geom_histogram(bins = 5)
Is this skewed enough to be worth worrying about (with 32 observations)?
The function make_dist obtains a bootstrap sampling distribution of
the sample mean, based on the input number n of simulations:
make_dist(mtcars$disp, nsim = 20)
#> # A tibble: 20 × 1
#> # Rowwise:
#> x_mean
#> <dbl>
#> 1 186.
#> 2 202.
#> 3 282.
#> 4 244.
#> 5 216.
#> 6 215.
#> 7 222.
#> 8 225.
#> 9 228.
#> 10 229.
#> 11 247.
#> 12 226.
#> 13 220.
#> 14 243.
#> 15 208.
#> 16 238.
#> 17 221.
#> 18 230.
#> 19 240.
#> 20 217
The function plot_dist obtains and plots a normal quantile plot of
the bootstrap sampling distribution of the sample mean, based on ten
thousand simulations (to facilitate investigation of the tails):
plot_dist(mtcars$disp)
This distribution is very close to normal, and so we have no need to worry about using a t-test for the mean. (The sample size is large enough to overcome the non-normality, as far as the Central Limit Theorem is concerned.)
There is nothing magic about a sample size of 30, beyond which the sample size is always large enough, despite what some people seem to think. It depends on how non-normal the original data are. Here is a sample of size 30 from a distribution with long tails:
d <- tibble(x = rt(30, 3))
ggplot(d, aes(x = x)) + geom_histogram(bins = 10)
The bootstrap sampling distribution of the sample mean looks like this:
plot_dist(d$x)
This is clearly skewed to the left, because the data we got happened to have a very low outlier. Hence, sometimes a sample size of 30 is big enough and sometimes it is not, depending on how non-normal your data are.
With two or more independent samples (where one might use a two-sample t-test or analysis of variance), the same issues of normality apply: if each group’s individual sampling distribution of the sample mean is normal, the t-test or ANOVA is trustworthy. (The situation tends to be a little better than that, because data from the groups are combined into an overall test statistic, but if all the distributions are close to normal, then the test is certainly safe.)
The package includes functions to obtain bootstrap sampling distributions of the sample means for each group, and to make normal quantile plots of them.
make_several obtains bootstrap sampling distributions of the sample
mean within each group. For example, in the mtcars dataframe, each car
has a value of cyl that could be 4, 6, or 8, and we could ask whether
it is reasonable to use an ANOVA to compare the disp values among the
three groups. If, as here, the input values come from a dataframe, it is
convenient to use with as shown. Here, we obtain only 5 bootstrap
sample means from each group:
with(mtcars, make_several(disp, cyl, nsim = 5))
#> # A tibble: 15 × 2
#> group x_mean
#> <dbl> <dbl>
#> 1 4 99.3
#> 2 4 110.
#> 3 4 99.2
#> 4 4 106.
#> 5 4 94.7
#> 6 6 163.
#> 7 6 171.
#> 8 6 178.
#> 9 6 192.
#> 10 6 181.
#> 11 8 357.
#> 12 8 348.
#> 13 8 361.
#> 14 8 352.
#> 15 8 380.
It is more useful to make separate normal quantile plots of the disp
values within each cyl group. This is done by taking 10,000 bootstrap
samples within each group, and making a facetted normal quantile plot:
with(mtcars, plot_several(disp, cyl))
The only distribution with any problems is the 6-cylinder cars, whose
distribution of disp has a short lower tail. If there were a long
upper tail to go with this, the distribution would be right-skewed, and
we would be right to be cautious about running an ANOVA. As it is,
though, it looks as if we may be all right.
Each facet is given its own scale, to allow for the possibility that the
mean disp for a different number of cylinders varies; our focus here
is on getting the best possible look at the normality within each group.
A second example is based on the home_prices dataset included in the
package. These are selling prices of a total of 37 3-bedroom and
4-bedroom homes in West Lafayette, Indiana. A realtor might wish to make
an inference about the value of a fourth bedroom (all else equal) as it
effects the selling price of a home. It seems likely that distributions
of selling prices will be too right-skewed for inference based on a
two-sample $t$ procedure to be acceptable:
with(home_prices, plot_several(price, bdrms))
The bootstrap sampling distribution of the sample mean for the three-bedroom houses is skewed to the right (as expected), but for the four-bedroom houses, the distribution is close to normal.
Often, percentage changes in prices are of interest, which would imply working with logs of prices:
with(home_prices, plot_several(log(price), bdrms))
and these distributions are very close to normal. Hence, an appropriate inference is
t.test(log(price) ~ bdrms, data = home_prices)
#>
#> Welch Two Sample t-test
#>
#> data: log(price) by bdrms
#> t = -5.1887, df = 30.59, p-value = 1.296e-05
#> alternative hypothesis: true difference in means between group 3 and group 4 is not equal to 0
#> 95 percent confidence interval:
#> -0.8568304 -0.3731165
#> sample estimates:
#> mean in group 3 mean in group 4
#> 11.82912 12.44410
and thus, with 95% confidence:
exp(c(-0.857, -0.373))
#> [1] 0.4244335 0.6886652
a 3-bedroom house sells for between 42% and 69% of the price of a 4-bedroom house.