summary.source_attribution: Summary function for saBayes output.

Description

This function generates a median and Chen Shao highest posterior density (HPD) region for each parameter in a source attribution model (source_attribution object). This algorithm requires unimodiality of the marginal posterior distribution for each parameter. If the marginal posterior distribution is not unimodel, the HPD region calculated may be invalid.

Usage

## S3 method for class 'source_attribution':
summary(object, alpha = 0.05, burn_in = 0, thin = 1, interval_type = "chen-shao", ...)

Arguments

object

an object that inherits from class source_attribution, such as a source_attribution object produced using saBayes.

alpha

a 100(1-alpha)% interval has a probability content of 1-alpha.

burn_in

the number iterations to remove from the beginning of the MCMC chain prior to calculating the HPD intervals. A value of 0 results in no burn in being removed.

thin

the amount to thin the MCMC chain by prior to calculating the HPD intervals. A vlaue of 1 results in no thinning.

interval_type

The type of HPD to calculate. This defaults to Chen Shao, but can be set to SPIn

...

Does nothing at the moment.

Value

Returns a list with possible names a, q, alpha, r, d, li, lj and lj_proportion. The list contains the same nesting structure as posterior.
lll{ Parameter Format and nesting Returns Source effects (a) A nested list with a matrix for each time and location. eg posterior$a[[t]][[l]][,j] gives the median, lower and upper values of the Chen Shao HPD region for a j for time t, location l. Always returned. Type effects (q) A matrix. eg posterior$q[,i] gives the the median, lower and upper values of the Chen Shao HPD region for q i for time t. Returned if names(posterior) contains q. r A nested list containing an array for each time. eg posterior$r[[t]][i,j,] gives the median, lower and upper values of the Chen Shao HPD region for r (type i, source j) at time t. Returned if names(posterior) contains r. Dispersion parameter (d) An vector. Returned if names(posterior) contains d. Rate of human infection from type i (Lambda i: li) A nested list with a matrix for each time and location. eg posterior$li[[t]][[l]][,i] gives the median, lower and upper values of the Chen Shao HPD region for li i for time t, location l. Returned if names(posterior) contains li. Rate of human infection from source j (Lambda j: lj) A nested list with a matrix for each time and location. eg posterior$lj[[t]][[l]][,j] gives the median, lower and upper values of the Chen Shao HPD region for lj j for time t, location l. Returned if names(posterior) contains lj. Proportion of human infection from source j (lj_proportion) A nested list with a matrix for each time and location. eg posterior$lj_proportion[[t]][[l]][,j] gives the median, lower and upper values of the Chen Shao HPD region for lj_proportion j for time t, location l. Returned if names(posterior) contains lj. }

References

Chen, M.-H. and Shao, Q.-M. (1998). Monte Carlo estimation of Bayesian credible and HPD intervals, Journal of Computational and Graphical Statistics, 7. Liu Y, Gelman A, Zheng T (2015). "Simulation-efficient shortest probability intervals." Statistics and Computing.

Examples

Run this code

##########################################################################
## Access simulated data set #############################################
##########################################################################
data(sim_SA)

##########################################################################
## Set priors ############################################################
##########################################################################

priors <- list(a = 1, r = 1, theta = c(0.01, 0.00001))

##########################################################################
## Run model #############################################################
##########################################################################

res <- saBayes(formula = Human~Source1+Source2+Source3+Source4+Source5, 
               time=~Time, location=~Location, type=~Type,
               data=sim_SA$data, priors = priors,
               alpha_conc = 1, prev = sim_SA$prev,
               likelihood_dist = "pois", n_iter = 20)
               
##########################################################################
#### Summary #############################################################
##########################################################################

s <- summary(res, alpha = 0.05, burn_in = 1, thin = 1)
# print the median and HPD for the type effects
s$q            
# print the median and HPD for the lambda's for each source
# for time 2, location A
s$lj_proportion$time2$locationA 
## or
s$lj_proportion[[2]][[1]]