BayesRobustProbit {BayesRGMM}

Perform MCMC algorithm to generate the posterior samples

Description

This function is used to generate the posterior samples using MCMC algorithm from the probit model with either the hypersphere decomposition or ARMA models applied to model the correlation structure in the serial dependence of repeated responses.

Usage

BayesRobustProbit(
  fixed,
  data,
  random,
  Robustness = TRUE,
  subset = NULL,
  na.action = "na.exclude",
  arma.order = NULL,
  HS.model = NULL,
  hyper.params = NULL,
  num.of.iter = 20000,
  Interactive = FALSE
)

Arguments

fixed

a two-sided linear formula object to describe fixed-effects with the response on the left of a ‘⁠~⁠’ operator and the terms separated by ‘⁠+⁠’ or ‘⁠*⁠’ operators, on the right. The specification first*second indicates the cross of first and second. This is the same as first + second + first:second.

data

an optional data frame containing the variables named in ‘⁠fixed⁠’ and ‘⁠random⁠’. It requires an “integer” variable named by ‘⁠id⁠’ to denote the identifications of subjects.

random

a one-sided linear formula object to describe random-effects with the terms separated by ‘⁠+⁠’ or ‘⁠*⁠’ operators on the right of a ‘⁠~⁠’ operator.

Robustness

logical. If 'TRUE' the distribution of random effects is assumed to be
t-distribution; otherwise normal distribution.

subset

an optional expression indicating the subset of the rows of ‘⁠data⁠’ that should be used in the fit. This can be a logical vector, or a numeric vector indicating which observation numbers are to be included, or a character vector of the row names to be included. All observations are included by default.

na.action

a function that indicates what should happen when the data contain NA’s. The default action (‘⁠na.omit⁠’, inherited from the ‘⁠factory fresh⁠’ value of
‘⁠getOption("na.action")⁠’) strips any observations with any missing values in any variables.

arma.order

a specification of the order in an ARMA model: the two integer components (p, q) are the AR order and the MA order.

HS.model

a specification of the correlation structure in HSD model:

• HS.model = ~0 denotes independence, that is, R_i is an identity matrix,

• HS.model = ~IndTime+\cdots+IndTimer denotes AR(r) correlation structure,

• HS.model = ~DiffTime1+\cdots+DiffTimer denotes correlation structure related to rth order of time difference.

hyper.params

specify the values in hyperparameters in priors.

num.of.iter

an integer to specify the total number of iterations; default is 20000.

Interactive

logical. If 'TRUE' when the program is being run interactively for progress bar and 'FALSE' otherwise.

Value

a list of posterior samples, parameters estimates, AIC, BIC, CIC, DIC, MPL, RJR, predicted values, and the acceptance rates in MH are returned.

Note

Only a model either HSD (‘⁠HS.model⁠’) or ARMA (‘⁠arma.order⁠’) model should be specified in the function. We'll provide the reference for details of the model and the algorithm for performing model estimation whenever the manuscript is accepted.

Author(s)

Kuo-Jung Lee [e-mail]

References

Lee K, Chen R, Kwak M, Lee K (2021). “Determination of correlations in multivariate longitudinal data with modified Cholesky and hypersphere decomposition using Bayesian variable selection approach.” Statistics in Medicine, 40(4), 978–997.

Lee K, Cho H, Kwak M, Jang EJ (2020). “Estimation of covariance matrix of multivariate longitudinal data using modified Choleksky and hypersphere decompositions.” Biometrics, 76(1), 75–86.

Examples

## Not run: 
library(BayesRGMM)
rm(list=ls(all=TRUE))
Fixed.Effs = c(-0.2, -0.3, 0.8, -0.4) #c(-0.2,-0.8, 1.0, -1.2)
P = length(Fixed.Effs) 
q = 1 #number of random effects
T = 5 #time points
N = 100 #number of subjects
num.of.iter = 100 #number of iterations
HSD.para = c(-0.5,  -0.3) #the parameters in HSD model
a = length(HSD.para)
w = array(runif(T*T*a), c(T, T, a)) #design matrix in HSD model

for(time.diff in 1:a)
	w[, , time.diff] = 1*(as.matrix(dist(1:T, 1:T, method="manhattan")) 
 ==time.diff)

#Generate a data with HSD model
 HSD.sim.data = SimulatedDataGenerator(
 Num.of.Obs = N, Num.of.TimePoints = T, Fixed.Effs = Fixed.Effs, 
 Random.Effs = list(Sigma = 0.5*diag(1), df=3), 
Cor.in.DesignMat = 0., Missing = list(Missing.Mechanism = 2, 
 RegCoefs = c(-1.5, 1.2)), Cor.Str = "HSD", 
 HSD.DesignMat.para = list(HSD.para = HSD.para, DesignMat = w))

hyper.params = list(
        sigma2.beta = 1,
        sigma2.delta = 1,
        v.gamma = 5,
        InvWishart.df = 5,
        InvWishart.Lambda = diag(q) )

HSD.output = BayesRobustProbit(
fixed = as.formula(paste("y~-1+", paste0("x", 1:P, collapse="+"))), 
data=HSD.sim.data$sim.data, random = ~ 1, Robustness=TRUE, 
HS.model = ~IndTime1+IndTime2, subset = NULL, na.action='na.exclude', 
hyper.params = hyper.params, num.of.iter = num.of.iter, 
 Interactive=0)

## End(Not run)