hanova1: Posterior sampling for a 1-way hierarchical ANOVA
In paulnorthrop/bang: Bayesian Analysis, No Gibbs
View source: R/one_way_anova.R
hanova1 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 factor indicating the group from
which the corresponding element of resp originates.
Must have the same length as resp.
...
Optional further arguments to be passed to
ru.
prior
The log-prior for the parameters of the hyperprior
distribution. If the user wishes to specify their own prior then
prior must be an object returned from a call to
set_user_prior.
Otherwise, prior is a character scalar giving the name of the
required in-built prior.
If prior is not supplied then a default prior is used.
See Details.
hpars
A numeric vector. Used to set parameters (if any) in
an in-built prior. If prior = cauchy then hpars is
a numeric vector of length 2 giving the respective scale parameters
of the half-Cauchy priors for σ_α and σ.
param
A character scalar.
If param = "trans" (the default) then the marginal posterior
of hyperparameter vector φ is reparameterized in terms of
log σ_α, log σ.
If param = "original" the original parameterization,
i.e. σ_α, σ is used.
The former tends to make the optimizations involved in the
ratio-of-uniforms algorithm more stable and to increase the probability
of acceptance, but at the expense of slower function evaluations.
init
A numeric vector. Optional initial estimates sent to
ru in the search for the mode of the posterior
density of (perhaps a subset of) the hyperparameter vector φ.
If an in-built prior is used then ru is used to sample from the
marginal posterior density of (σ_α, σ), so
init must have length 2. Otherwise, init has length
equal to the argument anova_d supplied to
set_user_prior.
mu0, sigma0
A numeric scalar. Mean and standard deviation of a
normal prior for μ. Only used if an in-built prior is used
or if anova_d = 2 is supplied in a call to
set_user_prior to set a user-defined prior.
The default, sigma0 = Inf, sets an improper uniform prior
for μ.
nrep
A numeric scalar. If nrep is not NULL then
nrep gives the number of replications of the original dataset
simulated from the posterior predictive distribution.
Each replication is based on one of the samples from the posterior
distribution. Therefore, nrep must not be greater than n.
In that event nrep is set equal to n.
Details
Consider I independent experiments in which the ni responses
yi from experiment/group i are normally
distributed with mean θ i and standard deviation σ.
The population parameters θ1, ...,
θI are modelled as random samples from a normal
distribution with mean μ and standard deviation
σ_α. Let φ = (μ, σ_α, σ).
Conditionally on θ1, ..., θI,
y1, ..., yI
are independent of each other and are independent of φ.
A hyperprior is placed on φ.
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 φ.
Then, conditional on these values, population parameters are sampled
directly from the conditional posterior density of
θ1, ..., θI given φ 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):
π(μ, σ_α, σ) = 1/σ,
that is, a uniform prior for (μ, σ_α, log σ),
for σ_α > 0 and σ > 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":
π(μ, σ_α, σ) = 1,
that is, a uniform prior for (μ, σ_α, σ),
for σ_α > 0 and σ > 0.
[See Section 11.6 of Gelman et al. (2014).]
prior = "cauchy": independent half-Cauchy priors for
σ_α and σ with respective scale parameters
A_α and A, that is,
π(σ_α, σ) =
1 / [(1 + σ_α^2 / A_α^2) (1 + σ^2 / A^2)].
[See Gelman (2006).] The scale parameters (A_α, A)
are specified using hpars = (A_α, A).
The default setting is hpars = c(10, 10).
Parameterizations for sampling:
param = "original" is (μ, σ_α, σ),
param = "trans" (the default) is
φ1 = μ, φ2 = log σ_α, φ3 = log σ.
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 φ.
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 θ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. https://doi.org/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))
paulnorthrop/bang documentation built on April 17, 2022, 2:12 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
View source: R/one_way_anova.R
| hanova1 | 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 factor indicating the group from
which the corresponding element of |
... |
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 μ. Only used if an in-built prior is used
or if |
nrep |
A numeric scalar. If |
Details
Consider I independent experiments in which the ni responses
yi from experiment/group i are normally
distributed with mean θ i and standard deviation σ.
The population parameters θ1, ...,
θI are modelled as random samples from a normal
distribution with mean μ and standard deviation
σ_α. Let φ = (μ, σ_α, σ).
Conditionally on θ1, ..., θI,
y1, ..., yI
are independent of each other and are independent of φ.
A hyperprior is placed on φ.
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 φ.
Then, conditional on these values, population parameters are sampled
directly from the conditional posterior density of
θ1, ..., θI given φ 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):
π(μ, σ_α, σ) = 1/σ,
that is, a uniform prior for (μ, σ_α, log σ),
for σ_α > 0 and σ > 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":
π(μ, σ_α, σ) = 1,
that is, a uniform prior for (μ, σ_α, σ),
for σ_α > 0 and σ > 0.
[See Section 11.6 of Gelman et al. (2014).]
prior = "cauchy": independent half-Cauchy priors for
σ_α and σ with respective scale parameters
A_α and A, that is,
π(σ_α, σ) =
1 / [(1 + σ_α^2 / A_α^2) (1 + σ^2 / A^2)].
[See Gelman (2006).] The scale parameters (A_α, A)
are specified using hpars = (A_α, A).
The default setting is hpars = c(10, 10).
Parameterizations for sampling:
param = "original" is (μ, σ_α, σ),
param = "trans" (the default) is
φ1 = μ, φ2 = log σ_α, φ3 = log σ.
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 φ.
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 θ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. https://doi.org/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))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.