hanova1 {bang}  R Documentation 
Produces random samples from the posterior distribution of the parameters of a 1way hierarchical ANOVA model.
hanova1(
n = 1000,
resp,
fac,
...,
prior = "default",
hpars = NULL,
param = c("trans", "original"),
init = NULL,
mu0 = 0,
sigma0 = Inf,
nrep = NULL
)
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 logprior 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 inbuilt 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 
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("bangcanovavignette", 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 halfCauchy 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
.
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
.
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), 515533. doi:10.1214/06BA117A.
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 userdefined prior.
# ======= 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)). Ratioofuniforms 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)
# Ratioofuniforms 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 ratioofuniforms
# 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 halfCauchy priors for sigma_alpha and sigma
coag3 < hanova1(resp = coagulation[, 1], fac = coagulation[, 2],
param = "original", prior = "cauchy", hpars = c(10, 1e6))