hanova1 {bang} | R Documentation |
Posterior sampling for a 1-way hierarchical ANOVA
Description
Produces random samples from the posterior distribution of the parameters of a 1-way hierarchical ANOVA model.
Usage
hanova1(
n = 1000,
resp,
fac,
...,
prior = "default",
hpars = NULL,
param = c("trans", "original"),
init = NULL,
mu0 = 0,
sigma0 = Inf,
nrep = NULL
)
Arguments
n |
A numeric scalar. The size of posterior sample required. |
resp |
A numeric vector. Response values. |
fac |
A vector of class |
... |
Optional further arguments to be passed to
|
prior |
The log-prior for the parameters of the hyperprior
distribution. If the user wishes to specify their own prior then
|
hpars |
A numeric vector. Used to set parameters (if any) in
an in-built prior. If |
param |
A character scalar.
If |
init |
A numeric vector. Optional initial estimates sent to
|
mu0 , sigma0 |
A numeric scalar. Mean and standard deviation of a
normal prior for |
nrep |
A numeric scalar. If |
Details
Consider I
independent experiments in which the ni
responses
y
i
from experiment/group i
are normally
distributed with mean \theta i
and standard deviation \sigma
.
The population parameters \theta
1, ...,
\theta
I
are modelled as random samples from a normal
distribution with mean \mu
and standard deviation
\sigma_\alpha
. Let \phi = (\mu, \sigma_\alpha, \sigma)
.
Conditionally on \theta
1, ..., \theta
I
,
y
1, ..., y
I
are independent of each other and are independent of \phi
.
A hyperprior is placed on \phi
.
The user can either choose parameter values of a default hyperprior or
specify their own hyperprior using set_user_prior
.
The ru
function in the rust
package is used to draw a random sample from the marginal posterior
of the hyperparameter vector \phi
.
Then, conditional on these values, population parameters are sampled
directly from the conditional posterior density of
\theta
1, ..., \theta
I
given \phi
and the data.
See the vignette("bang-c-anova-vignette", package = "bang")
for details.
The following priors are specified up to proportionality.
Priors:
prior = "bda"
(the default):
\pi(\mu, \sigma_\alpha, \sigma) = 1/\sigma,
that is, a uniform prior for (\mu, \sigma_\alpha, log \sigma)
,
for \sigma_\alpha > 0
and \sigma > 0
.
The data must contain at least 3 groups, that is, fac
must have
at least 3 levels, for a proper posterior density to be obtained.
[See Sections 5.7 and 11.6 of Gelman et al. (2014).]
prior = "unif"
:
\pi(\mu, \sigma_\alpha, \sigma) = 1,
that is, a uniform prior for (\mu, \sigma_\alpha, \sigma)
,
for \sigma_\alpha > 0
and \sigma > 0
.
[See Section 11.6 of Gelman et al. (2014).]
prior = "cauchy"
: independent half-Cauchy priors for
\sigma_\alpha
and \sigma
with respective scale parameters
A_\alpha
and A
, that is,
\pi(\sigma_\alpha, \sigma) =
1 / [(1 + \sigma_\alpha^2 / A_\alpha^2) (1 + \sigma^2 / A^2)].
[See Gelman (2006).] The scale parameters (A_\alpha
, A
)
are specified using hpars
= (A_\alpha
, A
).
The default setting is hpars = c(10, 10).
Parameterizations for sampling:
param = "original"
is (\mu, \sigma_\alpha, \sigma
),
param = "trans"
(the default) is
\phi1 = \mu, \phi2 = log \sigma_\alpha, \phi3 = log \sigma
.
Value
An object (list) of class "hef"
, which has the same
structure as an object of class "ru" returned from ru
.
In particular, the columns of the n
-row matrix sim_vals
contain the simulated values of \phi
.
In addition this list contains the arguments model
, resp
,
fac
and prior
detailed above and an n
by I
matrix theta_sim_vals
: column i
contains the simulated
values of \theta
i
. Also included are
data = cbind(resp, fac)
and summary_stats
a list
containing: the number of groups I
; the numbers of responses
each group ni
; the total number of observations; the sample mean
response in each group; the sum of squared deviations from the
group means s
; the arguments to hanova1
mu0
and
sigma0
; call
: the matched call to hanova1
.
References
Gelman, A., Carlin, J. B., Stern, H. S. Dunson, D. B., Vehtari, A. and Rubin, D. B. (2014) Bayesian Data Analysis. Chapman & Hall / CRC.
Gelman, A. (2006) Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1(3), 515-533. doi:10.1214/06-BA117A.
See Also
The ru
function in the rust
package for details of the arguments that can be passed to ru
via
hanova1
.
hef
for hierarchical exponential family models.
set_user_prior
to set a user-defined prior.
Examples
# ======= Late 21st Century Global Temperature Data =======
# Extract data for RCP2.6
RCP26_2 <- temp2[temp2$RCP == "rcp26", ]
# Sample from the posterior under the default `noninformative' flat prior
# for (mu, sigma_alpha, log(sigma)). Ratio-of-uniforms is used to sample
# from the marginal posterior for (log(sigma_alpha), log(sigma)).
temp_res <- hanova1(resp = RCP26_2[, 1], fac = RCP26_2[, 2])
# Plot of sampled values of (sigma_alpha, sigma)
plot(temp_res, params = "ru")
# Plot of sampled values of (log(sigma_alpha), log(sigma))
# (centred at (0,0))
plot(temp_res, ru_scale = TRUE)
# Plot of sampled values of (mu, sigma_alpha, sigma)
plot(temp_res)
# Estimated marginal posterior densities of the mean for each GCM
plot(temp_res, params = "pop", which_pop = "all", one_plot = TRUE)
# Posterior sample quantiles
probs <- c(2.5, 25, 50, 75, 97.5) / 100
round(t(apply(temp_res$sim_vals, 2, quantile, probs = probs)), 2)
# Ratio-of-uniforms information and posterior sample summaries
summary(temp_res)
# ======= Coagulation time data, from Table 11.2 Gelman et al (2014) =======
# With only 4 groups the posterior for sigma_alpha has a heavy right tail if
# the default `noninformative' flat prior for (mu, sigma_alpha, log(sigma))
# is used. If we try to sample from the marginal posterior for
# (sigma_alpha, sigma) using the default generalized ratio-of-uniforms
# runing parameter value r = 1/2 then the acceptance region is not bounded.
# Two remedies: reparameterize the posterior and/or increase the value of r.
# (log(sigma_alpha), log(sigma)) parameterization, ru parameter r = 1/2
coag1 <- hanova1(resp = coagulation[, 1], fac = coagulation[, 2])
# (sigma_alpha, sigma) parameterization, ru parameter r = 1
coag2 <- hanova1(resp = coagulation[, 1], fac = coagulation[, 2],
param = "original", r = 1)
# Values to compare to those in Table 11.3 of Gelman et al (2014)
all1 <- cbind(coag1$theta_sim_vals, coag1$sim_vals)
all2 <- cbind(coag2$theta_sim_vals, coag2$sim_vals)
round(t(apply(all1, 2, quantile, probs = probs)), 1)
round(t(apply(all2, 2, quantile, probs = probs)), 1)
# Pairwise plots of posterior samples from the group means
plot(coag1, which_pop = "all", plot_type = "pairs")
# Independent half-Cauchy priors for sigma_alpha and sigma
coag3 <- hanova1(resp = coagulation[, 1], fac = coagulation[, 2],
param = "original", prior = "cauchy", hpars = c(10, 1e6))