BayesRGMM (version 2.2)

BayesCumulativeProbitHSD: Perform MCMC algorithm to generate the posterior samples for longitudinal ordinal data

Description

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

Usage

BayesCumulativeProbitHSD(
  fixed,
  data,
  random,
  Robustness,
  subset,
  na.action,
  HS.model,
  hyper.params,
  num.of.iter,
  Interactive
)

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<U+2019>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.

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 \(r\)th 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.

References

Lee:etal:2021BayesRGMM

Lee:etal:2020BayesRGMM

Examples

# NOT RUN {
library(BayesRGMM)
rm(list=ls(all=TRUE))

Fixed.Effs = c(-0.1, 0.1, -0.1) #c(-0.8, -0.3, 1.8, -0.4) 
P = length(Fixed.Effs) 
q = 1 #number of random effects
T = 7 #time points
N = 100 #number of subjects
Num.of.Cats = 3 #in KBLEE simulation studies, please fix it. 
num.of.iter = 1000 #number of iterations

HSD.para = c(-0.9, -0.6) #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)

x = array(0, c(T, P, N))
for(i in 1:N){
    #x[,, i] = t(rmvnorm(P, rep(0, T), AR1.cor(T, Cor.in.DesignMat)))
    x[, 1, i] = 1:T
    x[, 2, i] = rbinom(1, 1, 0.5)
    x[, 3, i] = x[, 1, i]*x[, 2, i]
}

DesignMat = x

#Generate a data with HSD model


#MAR
CPREM.sim.data = SimulatedDataGenerator.CumulativeProbit(
 Num.of.Obs = N, Num.of.TimePoints = T, Num.of.Cats = Num.of.Cats, 
 Fixed.Effs = Fixed.Effs, Random.Effs = list(Sigma = 0.5*diag(1), df=3), 
 DesignMat = DesignMat, Missing = list(Missing.Mechanism = 2, 
 MissingRegCoefs=c(-0.7, -0.2, -0.1)), 
 HSD.DesignMat.para = list(HSD.para = HSD.para, DesignMat = w))

print(table(CPREM.sim.data$sim.data$y))
print(CPREM.sim.data$classes)

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

BCP.Est.output = BayesRobustProbitSummary(BCP.output)
# }