# R/credible.region.R In bayesSurv: Bayesian Survival Regression with Flexible Error and Random Effects Distributions

#### Documented in credible.region

#########################################################
#### AUTHOR:     Arnost Komarek                      ####
####             (2005)                              ####
####                                                 ####
#### FILE:       credible.region.R                   ####
####                                                 ####
#### FUNCTIONS:  credible.region                     ####
#########################################################

### ======================================
### credible.region
### ======================================
##
## Function to compute a simultaneous credible region for a vector parameter
## from the MCMC sample
##
## Reference: Besag, J., Green, P., Higdon, D. and Mengersen, K. (1995)
##            Bayesian computation and stochastic systems (with Discussion)
##            Statistical Science, vol. 10, 3 - 66, page 30
## and        Held, L. (2004)
##            Simultaneous inference in risk assessment; a Bayesian perspective
##            In: COMPSTAT 2004, Proceedings in Computational Statistics (J. Antoch, Ed.)
##            213 - 222, page 214
##
## \item{sample}{a data frame or matrix with sampled values (one column = one parameter)}
## \item{probs}{probabilities for which the credible regions are to be computed}
##
credible.region <- function(sample, probs=c(0.90, 0.975))
{
if (!length(dim(sample))) stop("Incorrect sample parameter supplied")
if (sum(probs < 0.50)) stop("probs must be each at least 0.50")
if (sum(probs >= 1)) stop("probs must be strictly lower than 1")

nn <- dim(sample)[1]
p <- dim(sample)[2]
k <- floor(nn*probs)

## Get ranks separately for each parameter
ranks <- apply(sample, 2, rank, ties.method="first")

## Compute sets S_j = {n + 1 - min_i r_i^{(j)}, max_i r_i^{(j)}}, j=1,...,n
## and the set S
S.left <- nn + 1 - apply(ranks, 1, min)
S.right <- apply(ranks, 1, max)
SS <- apply(cbind(S.left, S.right), 1, max)
SS <- SS[order(SS)]

## Get jstar and the credible region
jstar <- SS[k]
result <- list()
for (kk in 1:length(jstar)){
up <- sample[ranks == jstar[kk]]
low <- sample[ranks == nn + 1 - jstar[kk]]
result[[kk]] <- rbind(low, up)
rownames(result[[kk]]) <- c("Lower", "Upper")
colnames(result[[kk]]) <- colnames(sample)
}
names(result) <- paste(probs)

return(result)
}


## Try the bayesSurv package in your browser

Any scripts or data that you put into this service are public.

bayesSurv documentation built on May 2, 2019, 3:26 a.m.