This vignette introduces the basic functionality of the RoBTT package using a classic example from Ronald Fishers Design of Experiments (Fisher, 1971). Fisher illustrates the Student’s \(t\)-test on an example of self fertilized and cross fertilized plants as coded by Darwin. In other words, one group of plants was fertilized by gametes (sex-cells) from the same individual (self), while the other group of plants was fertilized by gametes from the same species (cross). Darwin himself concludes:
I may premise that if we took by chance a dozen men belonging to two nations and measured them, it would I presume be very rash to form any judgment from such small number on their average heights. But the case is somewhat different with crossed and self-fertilized plants, as they were exactly the same age, were subjected from first to last to the same conditions, and were descended from the same parents. p. 28 (Fisher, 1971).
Not knowing much about statistics he gives his data to Francis Galton, who proceeds to order the plants by height as a method of statistical analysis and concludes:
The observations as I received them […] certainly have no prima facie appearance of regularity. But as soon as we arrange them in the order of their magnitudes […] the case is materially altered. We now see, with few exceptions, that the largest plant on the crossed side in each pot exceeds the largest plant on the self-fertilised side, that the second exceeds the second, the third the third, and so on. . p. 29 (Fisher, 1971).
Galton then forwarded the data to Fisher, who later analyzed it with a frequentist Student’s \(t\)-test and states:
The observed value of \(t\), 2.148, thus just exceeds the 5 per cent point, and the experimental result may be judged significant, though barely so. .
We can load the included data with the code below
We reanalyze the data using a Bayesian equal variance \(t\)-test as follows.
## fit equal variance t-test
fit_01 <- RoBTT(
x1 = fertilization$Self,
x2 = fertilization$Crossed,
parallel = TRUE,
prior_delta = prior("cauchy", list(0, 1/sqrt(2))),
prior_rho = NULL, # this indicates no prior on the variance allocation factor -> equal variance test
prior_nu = NULL, # this indicates no prior on the degrees of freedom -> normal distribution test
seed = 0
)We can retrieve a summary of the estimates with the code below. This function gives us the posterior mean and median estimate for the effect in terms of Cohen’s \(d\) and 95% central credible interval. In addition, it shows the estimates for mean and variance in each of the groups.
# get overall summary
summary(fit_01)
#> Call:
#> RoBTT(x1 = fertilization$Self, x2 = fertilization$Crossed, prior_delta = prior("cauchy",
#> list(0, 1/sqrt(2))), prior_rho = NULL, prior_nu = NULL, parallel = TRUE,
#> seed = 0)
#>
#> Robust Bayesian t-test
#> Components summary:
#> Models Prior prob. Post. prob. Inclusion BF
#> Effect 1/2 0.500 0.744 2.902
#>
#> Model-averaged estimates:
#> Mean Median 0.025 0.975
#> delta 0.529 0.537 0.000 1.398We can also look at the estimates for the individual models and retrieve their prior distributions. The Bayesian equal variance \(t\)-test uses two models: One assuming absence of an effect (as indicated by prior \(\delta = \text{Spike}(0)\)) assuming presence of an effect (as indicated by prior \(\delta = \text{Cauchy}(0, 0.71)\)). Both models assume equal variances. The inclusion Bayes factor shows the evidence for the inclusion of the individual models. In this case the Bayes factor for model 2 (assuming the difference in mean) is equivalent to the Bayes factor for the presence of an effect.
# get individual model summaries
summary(fit_01, type = "models")
#> Call:
#> RoBTT(x1 = fertilization$Self, x2 = fertilization$Crossed, prior_delta = prior("cauchy",
#> list(0, 1/sqrt(2))), prior_rho = NULL, prior_nu = NULL, parallel = TRUE,
#> seed = 0)
#>
#> Robust Bayesian t-test
#> Models overview:
#> Model Distribution Prior delta Prior prob. log(marglik) Post. prob.
#> 1 normal Spike(0) 0.500 -76.82 0.256
#> 2 normal Cauchy(0, 0.71) 0.500 -75.76 0.744
#> Inclusion BF
#> 0.345
#> 2.902Finally, we can use the conditional = TRUE argument to
retrieve the estimates only for the models assuming an effect (i.e.,
without the spike).
summary(fit_01, conditional = TRUE)
#> Call:
#> RoBTT(x1 = fertilization$Self, x2 = fertilization$Crossed, prior_delta = prior("cauchy",
#> list(0, 1/sqrt(2))), prior_rho = NULL, prior_nu = NULL, parallel = TRUE,
#> seed = 0)
#>
#> Robust Bayesian t-test
#> Components summary:
#> Models Prior prob. Post. prob. Inclusion BF
#> Effect 1/2 0.500 0.744 2.902
#>
#> Model-averaged estimates:
#> Mean Median 0.025 0.975
#> delta 0.529 0.537 0.000 1.398
#>
#> Conditional estimates:
#> Mean Median 0.025 0.975
#> delta 0.713 0.703 0.008 1.454We can also analyze the example using an unequal variance \(t\)-test as follows.
## fit equal unvariance t-test
fit_10 <- RoBTT(
x1 = fertilization$Self,
x2 = fertilization$Crossed,
parallel = TRUE,
prior_delta = prior("cauchy", list(0, 1/sqrt(2))),
prior_rho = prior("beta", list(3, 3)), #prior on variance allocation
prior_rho_null = NULL, # remove models assuming equal variance
prior_nu = NULL,
seed = 0
)We again retrieve the summary; however, this time we also get a Bayes factor for unequal variances (heterogeneity) and an estimate for \(\rho\) (the precision proportion).
summary(fit_10)
#> Call:
#> RoBTT(x1 = fertilization$Self, x2 = fertilization$Crossed, prior_delta = prior("cauchy",
#> list(0, 1/sqrt(2))), prior_rho = prior("beta", list(3, 3)),
#> prior_nu = NULL, prior_rho_null = NULL, parallel = TRUE,
#> seed = 0)
#>
#> Robust Bayesian t-test
#> Components summary:
#> Models Prior prob. Post. prob. Inclusion BF
#> Effect 1/2 0.500 0.746 2.940
#> Heterogeneity 2/2 1.000 1.000 Inf
#>
#> Model-averaged estimates:
#> Mean Median 0.025 0.975
#> delta 0.539 0.546 0.000 1.421
#> rho 0.682 0.693 0.458 0.851We can again receive the individual model estimates. The
Prior rho column now contains a beta distribution rather
than spike at zero. This indicates that we used the models assuming
differences in variances.
summary(fit_10, type = "models")
#> Call:
#> RoBTT(x1 = fertilization$Self, x2 = fertilization$Crossed, prior_delta = prior("cauchy",
#> list(0, 1/sqrt(2))), prior_rho = prior("beta", list(3, 3)),
#> prior_nu = NULL, prior_rho_null = NULL, parallel = TRUE,
#> seed = 0)
#>
#> Robust Bayesian t-test
#> Models overview:
#> Model Distribution Prior delta Prior rho Prior prob. log(marglik)
#> 1 normal Spike(0) Beta(3, 3) 0.500 -75.96
#> 2 normal Cauchy(0, 0.71) Beta(3, 3) 0.500 -74.88
#> Post. prob. Inclusion BF
#> 0.254 0.340
#> 0.746 2.940The benefit of RoBTT in comparison to other commonly used \(t\)-tests is that it does not require us to rely on a single model but instead allows researchers to take different models into account at the same time. In other words, RoBTT uses Bayesian model-averaging to allow the data to guide the inference to be based most strongly on those models that predicted the data best. We can fit a \(t\)-test that model-averages over equal and unequal variances with the code provided below.
# fit t-test model averaging over equal and unequal variances
fit_1 <- RoBTT(
x1 = fertilization$Self,
x2 = fertilization$Crossed,
parallel = TRUE,
prior_delta = prior("cauchy", list(0, 1/sqrt(2))),
prior_rho = prior("beta", list(3, 3)),
prior_nu = NULL,
seed = 1,
control = set_control(adapt_delta = 0.95)
)We retrieve the overall estimates as before using the summary function.
summary(fit_1)
#> Call:
#> RoBTT(x1 = fertilization$Self, x2 = fertilization$Crossed, prior_delta = prior("cauchy",
#> list(0, 1/sqrt(2))), prior_rho = prior("beta", list(3, 3)),
#> prior_nu = NULL, parallel = TRUE, control = set_control(adapt_delta = 0.95),
#> seed = 1)
#>
#> Robust Bayesian t-test
#> Components summary:
#> Models Prior prob. Post. prob. Inclusion BF
#> Effect 2/4 0.500 0.746 2.936
#> Heterogeneity 2/4 0.500 0.705 2.393
#>
#> Model-averaged estimates:
#> Mean Median 0.025 0.975
#> delta 0.540 0.553 0.000 1.414
#> rho 0.630 0.632 0.486 0.845When looking at the individual models, we see that four models (instead of two) are incorporated. These are the two equal variance models (from the first test in this vignette) and the two unequal variance models (from the second test in this vignette). The inclusion Bayes factors tell us the evidence for each of these individual models. We see that the model assuming an effect and unequal variances has the highest predictive accuracy.
summary(fit_1, type = "models")
#> Call:
#> RoBTT(x1 = fertilization$Self, x2 = fertilization$Crossed, prior_delta = prior("cauchy",
#> list(0, 1/sqrt(2))), prior_rho = prior("beta", list(3, 3)),
#> prior_nu = NULL, parallel = TRUE, control = set_control(adapt_delta = 0.95),
#> seed = 1)
#>
#> Robust Bayesian t-test
#> Models overview:
#> Model Distribution Prior delta Prior rho Prior prob. log(marglik)
#> 1 normal Spike(0) Spike(0.5) 0.250 -76.82
#> 2 normal Spike(0) Beta(3, 3) 0.250 -75.96
#> 3 normal Cauchy(0, 0.71) Spike(0.5) 0.250 -75.76
#> 4 normal Cauchy(0, 0.71) Beta(3, 3) 0.250 -74.88
#> Post. prob. Inclusion BF
#> 0.076 0.245
#> 0.179 0.652
#> 0.219 0.842
#> 0.527 3.338A common problem for t-tests that assume the normality of residuals is outliers. When outliers are present, the estimates of these methods can be driven overly strongly by those outliers. A common way to accommodatee outliers is by also model-averaging over robust \(t\)-likelihoods. The t-distribution has thicker tails than the normal distribution; therefore, the estimate under the \(t\)-likelihood will be affected less due to outliers than the estimate using a normal likelihood. We can fit RoBTT that also model-averages over \(t\)-likelihoods as follows.
fit_2 <- RoBTT(
x1 = fertilization$Self,
x2 = fertilization$Crossed,
parallel = TRUE,
prior_delta = prior("cauchy", list(0, 1/sqrt(2))),
prior_rho = prior("beta", list(3, 3)),
prior_nu = prior("exp", list(1)), # prior on degrees of freedom
seed = 2
)When looking at the individual model estimates incorporating \(t\)-likelihoods, we find that the four models from the previous t-test are incorporated a second time, however, with t-likelihood in the likelihood column. In other words, we now use eight models overall.
summary(fit_2, type = "models")
#> Call:
#> RoBTT(x1 = fertilization$Self, x2 = fertilization$Crossed, prior_delta = prior("cauchy",
#> list(0, 1/sqrt(2))), prior_rho = prior("beta", list(3, 3)),
#> prior_nu = prior("exp", list(1)), parallel = TRUE, seed = 2)
#>
#> Robust Bayesian t-test
#> Models overview:
#> Model Distribution Prior delta Prior rho Prior nu Prior prob.
#> 1 normal Spike(0) Spike(0.5) 0.125
#> 2 t Spike(0) Spike(0.5) Exponential(1) 0.125
#> 3 normal Spike(0) Beta(3, 3) 0.125
#> 4 t Spike(0) Beta(3, 3) Exponential(1) 0.125
#> 5 normal Cauchy(0, 0.71) Spike(0.5) 0.125
#> 6 t Cauchy(0, 0.71) Spike(0.5) Exponential(1) 0.125
#> 7 normal Cauchy(0, 0.71) Beta(3, 3) 0.125
#> 8 t Cauchy(0, 0.71) Beta(3, 3) Exponential(1) 0.125
#> log(marglik) Post. prob. Inclusion BF
#> -76.82 0.006 0.042
#> -77.80 0.002 0.016
#> -75.96 0.014 0.101
#> -77.33 0.004 0.025
#> -75.76 0.017 0.124
#> -72.35 0.525 7.751
#> -74.88 0.042 0.305
#> -72.65 0.389 4.463
#> [0;31mModel (6): There were 2 divergent transitions.[0mIt is sometimes argued that point null hypotheses are not meaningful to test as evidence for the alternative could already be provided by deviations from the null that would be considered small or uninteresting. RoBTT addresses this concern by allowing researchers to specify perinull hypotheses (i.e., a prior distribution with small sd instead of the point null).
fit_3 <- RoBTT(
x1 = fertilization$Self,
x2 = fertilization$Crossed,
parallel = TRUE,
prior_delta = prior("cauchy", list(0, 1/sqrt(2)), list(0, Inf)),
prior_rho = prior("beta", list(3, 3)),
prior_nu = prior("exp", list(1)),
prior_delta_null = prior("normal", list(0, 0.15), list(0, Inf)), #prior distribution and truncation
seed = 3
)When looking at the individual models we see that the models that previously had Spike(0.5) are now using normal(0, 0.15) prior from 0 to infinity. This indicates the shift from point null to perinull testing.
summary(fit_3, type = "models")
#> Call:
#> RoBTT(x1 = fertilization$Self, x2 = fertilization$Crossed, prior_delta = prior("cauchy",
#> list(0, 1/sqrt(2)), list(0, Inf)), prior_rho = prior("beta",
#> list(3, 3)), prior_nu = prior("exp", list(1)), prior_delta_null = prior("normal",
#> list(0, 0.15), list(0, Inf)), parallel = TRUE, seed = 3)
#>
#> Robust Bayesian t-test
#> Models overview:
#> Model Distribution Prior delta Prior rho Prior nu
#> 1 normal Normal(0, 0.15)[0, Inf] Spike(0.5)
#> 2 t Normal(0, 0.15)[0, Inf] Spike(0.5) Exponential(1)
#> 3 normal Normal(0, 0.15)[0, Inf] Beta(3, 3)
#> 4 t Normal(0, 0.15)[0, Inf] Beta(3, 3) Exponential(1)
#> 5 normal Cauchy(0, 0.71)[0, Inf] Spike(0.5)
#> 6 t Cauchy(0, 0.71)[0, Inf] Spike(0.5) Exponential(1)
#> 7 normal Cauchy(0, 0.71)[0, Inf] Beta(3, 3)
#> 8 t Cauchy(0, 0.71)[0, Inf] Beta(3, 3) Exponential(1)
#> Prior prob. log(marglik) Post. prob. Inclusion BF
#> 0.125 -76.06 0.006 0.040
#> 0.125 -73.76 0.058 0.428
#> 0.125 -75.20 0.014 0.097
#> 0.125 -73.91 0.049 0.364
#> 0.125 -75.10 0.015 0.107
#> 0.125 -71.65 0.473 6.274
#> 0.125 -74.21 0.037 0.266
#> 0.125 -71.95 0.349 3.757This vignette showcased the applications of the RoBTT R package. For methodological background see Maier et al. (2022) (https://doi.org/10.31234/osf.io/d5zwc).