bimod: Bivariate Bayestian hierarchical model for combining...
In yy2725/bmsSum: Bayesian hierarchical models for combining the summay measures

Description Usage Arguments Value Examples

View source: R/bmsSum.R

bimod is modeling the measures and its uncertainty jointly for combining summaries from multiple sources. It returns the estimates and posterior distributions of parameters of interest fitted with data.

bimod(
  y = data$y,
  s = data$s,
  with_cov = FALSE,
  X = NULL,
  bi_option = TRUE,
  seed = 9999,
  iter = 3000,
  chain = 4,
  verbose = FALSE,
  warmup = 1000,
  thin = 3,
  adapt_delta = 0.99,
  max_treedepth = 3,
  ...
)

`y`	A n-by-1 vector, consisting of the measures of n observations.
`s`	A n-by-1 vector, consisting of the standard error of the measures for n observations.
`with_cov`	if TRUE, include covariates in the model.
`X`	A n-by-p matrix of p covariates including intercept, which is all 1 for the first column.
`bi_option`	if TRUE, use empirical value for standard deviation of uncertainty of measure
`seed`	The seed for random number generation. The default is generated from 1 to the maximum integer supported by R on the machine. Even if multiple chains are used, only one seed is needed, with other chains having seeds derived from that of the first chain to avoid dependent samples.The default is 9999.
`iter`	a number to specify the iteration times of Stan. The default is 3000.
`verbose`	TRUE or FALSE: flag indicating whether to print intermediate output from Stan on the console, which might be helpful for model debugging.
`warmup`	a number of warmup iterations to be excluded when computing the summaries. The default is 1000.
`thin`	A positive integer specifying the period for saving samples. The default value is 3.
`adapt_delta`	A parameters to control the sampler’s behavior. The default is 0.8.
`max_treedepth`	A parameters to control the sampler’s behavior in For algorithm NUTS.

a list of information of the fitted bivariate Bayestian hierarchical model including

pop_par: the estimates and credible intervals for the population-level parameters
ind_par: the estimates and credible intervals for the individual-level parameters
poster_dist: posterior distributions of all parameters fitted in the model

library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
library(MASS)
set.seed(99999)

### define data generation functions-no predictors
data.nonreg <-  function(n=50,mutheta=5,musigma=2,
                         tautheta=3, tausigma=1, r1=0.3, r2=0.3, sigmas=1,...){
  # mutheta,musigma,tautheta, tausigma -population level
  # r1 - individual level
  # r2 - population level
  theta = vector()
  sigma2 = vector()
  y = vector()
  s2 = vector()
  Sigma = list()
  indiv = matrix(rep(NA),ncol = 2, nrow = n)
  # generate individual true mean and variance from population
  mu <- c(mutheta,musigma)
  varmat <- matrix(c(tautheta^2,r2*tausigma*tautheta,r2*tausigma*tautheta,tausigma^2), nrow=2)
  temp <- MASS::mvrnorm(n, mu, varmat)
  theta <- temp[,1]
  sigma2 <- (exp(temp[,2]))^2
  # generate observed value from true individual
  for (i in 1:n){
    Sigma[[i]] <- c(sigma2[i],r1*sqrt(sigma2[i])*sigmas,r1*sqrt(sigma2[i])*sigmas,sigmas^2)
    indiv[i,] <- mvrnorm(1, temp[i,],matrix(Sigma[[i]], nrow = 2))
  }
  y <- indiv[,1]
  s <- exp(indiv[,2])
  dataset <- data.frame(cbind(y, s))
  return(list(data = dataset,
              para = c(n=n,mutheta=mutheta,musigma=musigma,
                       tautheta=tautheta, tausigma=tausigma, r1=r1, r2=r2, sigmas=sigmas)))
}


data <- data.nonreg(r1=0,r2=0)$data
res <- bimod(data$y, data$s, with_cov = FALSE, bi_option = TRUE,
             seed = 9999, iter = 3000, chain = 3, verbose = FALSE,
             warmup= 1000, thin=3, adapt_delta = 0.99)$pop_par