ordinal_post_beta: Data augmentation Gibbs sampling for Multivariate Ordinal...

In rochita07/BayesMVProbit: Data augmentation Gibbs sampling for Multivariate Multinomial Probit model and Multivariate ordered Probit model

Description

This function provides posterior mean and credible interval along with trace plots, density plots and carter plot for parameters using data augmentation algorithm (Albert Chib method) for posterior sampling in the Multivariate Ordinal Probit

Usage

ordinal_post_beta(
  category,
  df_t = NULL,
  iter = 5000,
  burn = 1000,
  cred_level = 0.95,
  x_list,
  sig,
  y,
  prior_delta_mean = NULL,
  prior_delta_var = NULL,
  prior_beta_mean = NULL,
  prior_beta_var = NULL
)

Arguments

category	vector of categories for each variable, for any variable the no of category should be atleast 3 to consider this model
df_t	vector of degrees of freedom for proposed student's t distribution delta (as mentioned in the paper, link attached in vignette) for each variable
iter	scalar, number of iteration for posterior to be calculated, default is 5000
burn	scalar, number of burn in for posterior to be calculated, defult is 1000
cred_level	scaler, (should be in [0,1]) specifies at what level credible interval to be calculated, default value is 0.95
x_list	must be user defined as list, each level of the list represents covariates for each level. Each variable must have atleast one covariate, subject may have level specific covariates othwerwise same value of covariates to be repeated for each ordinal variable
sig	covariance matrix of error vector, must be symmetric positive definite matrix
y	must be user defined as matrix, each column should represent for each subject
prior_delta_mean	must be user defined as list, each level is vector of prior mean for delta, by default it takes value 0 as prior mean for all delta at a specific level
prior_delta_var	must be user defined as list, each level is a square positive definite matrix of prior variance for delta, default is Identity matrix of proper dimension
prior_beta_mean	vector of prior mean for beta , by default it takes value 0 as prior mean for all beta
prior_beta_var	a square positive definite matrix of prior variance for beta, default is Identity matrix of proper dimension

Value

Posterior_mean : posterior mean of beta in vector form indicating the ordinal measure and the corresponding covariate of the ordinal measure

Credible_interval : credible interval for beta vector at specified level

trace_plot : a plot of iterations vs. sampled values for each variable in the chain, with a separate plot per variable

density_plot : the distribution of variables, with a separate plot per variable

carter_plot : plots of credible intervals for parameters from an MCMC simulation

Examples

# Data Generation:

variable_dim = 2 # no of variables, i.e dimension of yi (can be obtained if y is given)
category = c(4, 3) ## no of categories(levels) for each variable
covariate_num = c(2, 3) # No of covariates for each level
nu_tot = category + 1 ## no of cut points 
nu = nu_tot - 3 ## no of cutpoints to be considered for each variable ( nu0, nu1, nuJ are fixed)
nu_cutoff1 = c(0, 2 , 2.5)  # of length nu + 1
nu_cutoff2 = c(0, 4.5)
beta_dim = sum(covariate_num) # dimension of beta vector

n = 50 # no of subjects

# Generation of beta

beta_act1 = rep(2, covariate_num[1])
beta_act2 = rep(2,covariate_num[2])
beta_act = c(beta_act1, beta_act2)

set.seed(1287)

# Generation of x

x1 = matrix(rnorm(covariate_num[1]* n), nrow = covariate_num[1])  # each column is for each subject
x2 = matrix(rnorm(covariate_num[2]* n), nrow = covariate_num[2])  # each column is for each subject
x_list = list(x1, x2)

x = lapply(1:n, function(x) matrix(rep(0, variable_dim * beta_dim), nrow = variable_dim) ) ## initialization
a = matrix(rep(0, variable_dim * beta_dim), nrow = variable_dim) ## to store the value in the loop
col_indi = c(0, cumsum(covariate_num))  # of length variable_dim + 1
for(i in 1:n)
{
  for(l in 1: variable_dim)
  {
    a[l, (col_indi[l] + 1) : col_indi[l + 1]] = t(x_list[[l]][,i])
  }
  x[[i]] = a
}


# Generation of error 

sig = diag(variable_dim) ## error covariance matrix of order variable_dim x variable_dim ( to be mentioned in the input)

e = mvtnorm::rmvnorm(n = n, mean = rep(0, variable_dim) , sigma = sig)
e = t(e) # # each column is for each subject


# Generation of  Z

z = matrix(rep(0, variable_dim * n), nrow = variable_dim)  # matrix of order variable_dim x n # # each column is for each subject
for(i in 1:n)
{
  z[,i] =  x[[i]] %*% beta_act + e[, i]
}

# Generation of Y 

cutoff1 = c(-Inf, nu_cutoff1, Inf)  ## To also consider the fixed cutoffs 
cutoff2 = c(-Inf, nu_cutoff2, Inf)
cutoff = list(cutoff1, cutoff2)

y = matrix(rep(0, variable_dim * n), nrow = variable_dim) ## initialization of matrix of order variable_dim x n # each column is for each subject

for(i in 1:n)
{
  for(l in 1:variable_dim)
  {
    for(j in 1: length(cutoff[[l]]))
    {
      if(z[l, i] > cutoff[[l]][j] && z[l, i] <= cutoff[[l]][j + 1]) # making one side inclusive
        y[l, i] = j
    }
  }
}

y ## example of input 'y' (user should provide)


# Example 1

ordinal_post_beta(category = c(4, 3), df_t = NULL, iter = 10,  burn = 1, cred_level = 0.95, x_list, sig = diag(2), y , prior_delta_mean = NULL, prior_delta_var = NULL, prior_beta_mean = NULL, prior_beta_var = NULL)


# Example 2

 prior_delta_mean = list()  # Prior on delta
 for(i in 1: variable_dim)
 {
   prior_delta_mean[[i]] = rep(0, nu[i])
 }
 
 prior_delta_var = list()   # Prior on delta
 for(i in 1: variable_dim)
   {
     prior_delta_var[[i]] = diag(nu[i])
   }

prior_beta_mean = rep(1, beta_dim)  # Prior on beta
prior_beta_var = diag(2, beta_dim)  # Prior on beta

ordinal_post_beta(category = c(4, 3), df_t = NULL, iter = 10,  burn = 1, cred_level = 0.95, x_list, sig = diag(2), y , prior_delta_mean, prior_delta_var, prior_beta_mean, prior_beta_var)


# Interpretation of indices of beta
# Indices for a speficific beta indicate the ordinal measure and the corresponding covariate of the ordinal measure


# Warnings: There may be warning message "Error in is.positive.definite(beta_post_var) : could not find function "is.positive.definite""

# This is beacuse of mvtnorm::rmvnorm prefers sigma to be symmetric matrix. But, covariance matrix of beta not necessarily always would be symmetric matrix.
# Though this warnings does not affect sampling from mvtnorm::rmvnorm

# For example we define covariance matrix of  posterior beta from the given example as follow: 
sum_t_x_sig_x = matrix(rep(0, beta_dim^2), nrow = beta_dim) ## initialization
for(i in 1: n)
{
  sum_t_x_sig_x = sum_t_x_sig_x + t(x[[i]]) %*% sig %*% x[[i]]
}
sum_t_x_sig_x 

# Considering N(0,I) prior on beta we calculate covariance matrix of posterior beta as :
beta_post_var = solve(solve(diag(beta_dim)) + sum_t_x_sig_x )  # variance of posterior of beta
beta_post_var

# For a single sample the following code runs with out any error
mvtnorm::rmvnorm(n, mean = rep(0, nrow(beta_post_var)), beta_post_var)

rochita07/BayesMVProbit documentation built on Dec. 11, 2019, 2:07 a.m.