#' Simulate hyperpriors from an GMRF
#'
#'
#' @param R Desired prior odds ratio. Default to 2, i.e., a 95\% prior interval for the residual odds ratios lies in the interval (R, 1/R).
#' @param nsamp Sample to simulate for scaling factor
#' @param nsamp.check Sample to simulate for checking range
#' @param Amat Adjacency matrix of the areas in the data.
#' @param nperiod numerical value of how many time periods in the data
#' @param only.iid Indicator for whether or not only IID hyperpriors are simulated
#' @author Zehang Richard Li, Laina Mercer
#' @references Wakefield, J. Multi-level modelling, the ecologic fallacy, and hybrid study designs. \emph{International Journal of Epidemiology}, 2009, vol. 38 (pg. 330-336).
#' @examples
#' \dontrun{
#' data(DemoMap)
#' mat <- DemoMap$Amat
#' priors <- simhyper(R = 2, nsamp = 1e+05, nsamp.check = 5000, Amat = mat)
#' }
#' @export
simhyper <- function(R = 2, nsamp = 1e+05, nsamp.check = 5000, Amat = NULL, nperiod = 6, only.iid = TRUE) {
#################################################################### (R,1/R) is the range of the residual odds ratios gives a=d/2 where d = degrees of freedom of marginal Student’s t
d <- 1
a <- d/2
p <- 0.025
b <- (log(R))^2 * d/(2 * stats::qt(p, df = d)^2)
a.iid <- a
b.iid <- b
if (!only.iid) {
##################################################################### Check range using simulation #
tausamp <- stats::rgamma(nsamp, a, b)
Usamp <- stats::rnorm(nsamp, mean = 0, sd = 1/sqrt(tausamp))
# The adjacency matrix for 6 years
m2 = c(1, rep(2, nperiod - 2), 1)
# create the adjacency #
before <- 1:(nperiod - 2)
after <- 3:(nperiod)
alternate <- c(rbind(before, after))
adj2 <- c(2, alternate, nperiod - 1)
make.Q <- function(num.neighbors, neighbors, omega.sq = 1) {
n <- length(num.neighbors)
mat <- matrix(0, ncol = n, nrow = n)
diag(mat) <- num.neighbors
mat[cbind(rep(1:n, num.neighbors), neighbors)] <- -1
mat/omega.sq
}
vars.Q <- function(eigenvalues, eigenvectors) {
margsum <- 0
nloop <- length(eigenvalues) - 1
for (i in 1:nloop) {
ev <- eigenvectors[, i]
margsum <- margsum + ev %*% t(ev)/eigenvalues[i]
}
margvars <- diag(margsum)
margvars
}
#
sim.Q <- function(Q) {
eigenQ <- eigen(Q)
rankQ <- qr(Q)$rank
sim <- as.vector(eigenQ$vectors[, 1:rankQ] %*% matrix(stats::rnorm(rep(1, rankQ), rep(0, rankQ), 1/sqrt(eigenQ$values[1:rankQ])),
ncol = 1))
sim
}
#
Q <- make.Q(m2, adj2, 1)
eigentemp <- eigen(Q)
eigenvaluesQ <- eigentemp$values
eigenvectorsQ <- eigentemp$vectors
rankQ <- qr(Q)$rank # 5
margy <- mean(vars.Q(eigenvaluesQ, eigenvectorsQ))
a.rw1 = a
b.rw1 = b/margy
taustarsamp <- stats::rgamma(nsamp.check, a.rw1, b.rw1)
Ustarsamp <- matrix(nrow = nsamp.check, ncol = nperiod)
for (i in 1:nsamp.check) {
Qstar <- Q * taustarsamp[i]
Ustarsamp[i, ] <- sim.Q(Qstar)
}
check.rw1 = stats::quantile(exp(Ustarsamp), p = c(0.025, 0.5, 0.975))
# tausamp <- rgamma(nsamp,a,b) Usamp <- rnorm(nsamp,mean=0,sd=1/sqrt(tausamp)) quantile(exp(Usamp),p=c(0.025,0.5,0.975)) m2
# <- c(2, 3, rep(4, nperiod - 4), 3, 2) # create the adjacency # before2 <- 1:(nperiod - 4) before1 <- 2:(nperiod - 3) after1
# <- 4:(nperiod - 1) after2 <- 5:(nperiod) alternate <- c(rbind(before1, before2, after1, after2))
# # adj2<-c(2, 3, # 1, 3, 4, # alternate, # nperiod - 3, nperiod - 3, nperiod, # nperiod - 2, nperiod - 1) # Q <- make.Q(m2,
# adj2, 1) The adjacency matrix for 6 years
if (nperiod > 4) {
Q <- matrix(0, nrow = nperiod, ncol = nperiod)
Q[1, 1:3] <- c(1, -2, 1)
Q[2, 1:4] <- c(-2, 5, -4, 1)
for (j in 3:(nperiod - 2)) {
Q[j, (j - 2):(j + 2)] <- c(1, -4, 6, -4, 1)
}
Q[nperiod, (nperiod - 2):nperiod] <- c(1, -2, 1)
Q[nperiod - 1, (nperiod - 3):nperiod] <- c(1, -4, 5, -2)
} else {
stop("RW2 prior not specified for n < 5")
}
vars.Q2 <- function(eigenvalues, eigenvectors, rankQ) {
margsum <- 0
nloop <- rankQ
for (i in 1:nloop) {
ev <- eigenvectors[, i]
margsum <- margsum + ev %*% t(ev)/eigenvalues[i]
}
margvars <- diag(margsum)
margvars
}
sim.Q <- function(Q) {
eigenQ <- eigen(Q)
rankQ <- qr(Q)$rank
sim <- as.vector(eigenQ$vectors[, 1:rankQ] %*% matrix(stats::rnorm(rep(1, rankQ), rep(0, rankQ), 1/sqrt(eigenQ$values[1:rankQ])),
ncol = 1))
sim
}
eigentemp <- eigen(Q)
eigenvaluesQ <- eigentemp$values
eigenvectorsQ <- eigentemp$vectors
rankQ <- qr(Q)$rank # 4
margy <- mean(vars.Q2(eigenvaluesQ, eigenvectorsQ, rankQ))
a.rw2 = a
b.rw2 = b/margy
taustarsamp <- stats::rgamma(nsamp.check, a.rw2, b.rw2)
Ustarsamp <- matrix(nrow = nsamp.check, ncol = nperiod)
for (i in 1:nsamp.check) {
Qstar <- Q * taustarsamp[i]
Ustarsamp[i, ] <- sim.Q(Qstar)
}
check.rw2 = stats::quantile(exp(Ustarsamp), p = c(0.025, 0.5, 0.975))
####################################################################
tausamp <- stats::rgamma(nsamp, a, b)
Usamp <- stats::rnorm(nsamp, mean = 0, sd = 1/sqrt(tausamp))
m2 <- apply(Amat, 1, sum)
# create the adjacency list #
nums <- c(1:dim(Amat)[1])
adj2 <- NULL
for (i in 1:dim(Amat)[1]) {
adj2 <- c(adj2, nums[as.numeric(Amat[i, ]) == 1])
}
make.Q <- function(num.neighbors, neighbors, omega.sq = 1) {
n <- length(num.neighbors)
mat <- matrix(0, ncol = n, nrow = n)
diag(mat) <- num.neighbors
mat[cbind(rep(1:n, num.neighbors), neighbors)] <- -1
mat/omega.sq
}
vars.Q3 <- function(eigenvalues, eigenvectors, rankQ) {
margsum <- 0
nloop <- rankQ
for (i in 1:nloop) {
ev <- eigenvectors[, i]
margsum <- margsum + ev %*% t(ev)/eigenvalues[i]
}
margvars <- diag(margsum)
margvars
}
#
sim.Q <- function(Q) {
eigenQ <- eigen(Q)
rankQ <- qr(Q)$rank
sim <- as.vector(eigenQ$vectors[, 1:rankQ] %*% matrix(stats::rnorm(rep(1, rankQ), rep(0, rankQ), 1/sqrt(eigenQ$values[1:rankQ])),
ncol = 1))
sim
}
#
Q <- make.Q(m2, adj2, 1)
eigentemp <- eigen(Q)
eigenvaluesQ <- eigentemp$values
eigenvectorsQ <- eigentemp$vectors
rankQ <- qr(Q)$rank # 20
margy <- mean(vars.Q3(eigenvaluesQ, eigenvectorsQ, rankQ))
a.icar = a
b.icar = b/margy
taustarsamp <- stats::rgamma(nsamp.check, a.icar, b.icar)
Ustarsamp <- matrix(nrow = nsamp.check, ncol = dim(Amat)[1])
for (i in 1:nsamp.check) {
Qstar <- Q * taustarsamp[i]
Ustarsamp[i, ] <- sim.Q(Qstar)
}
check.icar = stats::quantile(exp(Ustarsamp), p = c(0.025, 0.5, 0.975))
}else{
a.rw1 <- a.rw2 <- a.icar <- a.iid
b.rw1 <- b.rw2 <- b.icar <- b.iid
check.rw1 <- NULL
check.rw2 <- NULL
check.icar <- NULL
}
return(list(a.iid = a.iid, b.iid = b.iid, a.rw1 = a.rw1, b.rw1 = b.rw1, check.rw1 = check.rw1, a.rw2 = a.rw2, b.rw2 = b.rw2,
check.rw2 = check.rw2, a.icar = a.icar, b.icar = b.icar, check.icar = check.icar))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.