R/rhierMnlDPParallel.R

In scalablebayesm: Distributed Markov Chain Monte Carlo for Bayesian Inference in Marketing

#' @title MCMC Algorithm for Hierarchical Multinomial Logit with Dirichlet Process Prior Heterogeneity
#' @description
#' \code{rhierMnlDPParallel} is an MCMC algorithm for a hierarchical multinomial logit with a Dirichlet Process prior describing the distribution of heteorogeneity. A base normal model is used so that the DP can be interpreted as allowing for a mixture of normals with as many components as panel units. This is a hybrid Gibbs Sampler with a RW Metropolis step for the MNL coefficients for each panel unit. This procedure can be interpreted as a Bayesian semi-parametric method in the sense that the DP prior can accommodate heterogeneity of an unknown form.
#' 
#' @param Data A list of: `p`- number of choice alternatives, `lgtdata` - An \eqn{nlgt} size list of multinomial logistic data, and optional `Z`- matrix of unit characteristics.
#' @param Prior A list with one required parameter: `ncomp`, and optional parameters: `mubar`, `Amu`, `nu`, `V`, `Ad`, `deltaBar`, and `a`.
#' @param Mcmc A list with one required parameter: 'R' - number of iterations, and optional parameters: `s`, `w`, `keep`, `nprint`, and `drawcomp`. 
#' @param verbose If \code{TRUE}, enumerates model parameters and timing information.
#' 
#' @return A list containing the following elements:
#' \itemize{
#'   \item \strong{nmix}: A list with the following components:
#'   \itemize{
#'     \item \strong{probdraw}: A matrix of size \code{(R/keep) x (ncomp}, containing the probability draws at each Gibbs iteration.
#'     \item \strong{compdraw}: A list containing the drawn mixture components at each Gibbs iteration.}
#'   \item \strong{Deltadraw} (optional): A matrix of size \code{(R/keep) x (nz * nvar)}, containing the delta draws, if \code{Z} is not \code{NULL}. If \code{Z} is \code{NULL}, this element is not included.
#' }
#' 
#' @details{
#'  \subsection{Model and Priors}{
#'    \eqn{y_i} \eqn{\sim}{~} \eqn{MNL(X_i, \beta_i)} with \eqn{i = 1, \ldots, length(lgtdata)} and where \eqn{\theta_i} is \eqn{nvar x 1}
#'    
#'    \eqn{\beta_i = Z\Delta}[i,] + \eqn{u_i} \cr
#'    Note:  Z\eqn{\Delta} is the matrix \eqn{Z * \Delta}; [i,] refers to \eqn{i}th row of this product \cr
#'    Delta is an \eqn{nz x nvar} matrix 
#'    
#'    \eqn{\beta_i} \eqn{\sim}{~} \eqn{N(\mu_i, \Sigma_i)}
#'    
#'    \eqn{\theta_i = (\mu_i, \Sigma_i)} \eqn{\sim}{~} \eqn{DP(G_0(\lambda), alpha)}\cr
#'    
#'    \eqn{G_0(\lambda):}\cr
#'    \eqn{\mu_i | \Sigma_i} \eqn{\sim}{~} \eqn{N(0, \Sigma_i (x) a^{-1})}\cr
#'    \eqn{\Sigma_i} \eqn{\sim}{~} \eqn{IW(nu, nu*v*I)}\cr
#'    \eqn{delta = vec(\Delta)} \eqn{\sim}{~} \eqn{N(deltabar, A_d^{-1})}\cr
#'    
#'    \eqn{\lambda(a, nu, v):}\cr
#'    \eqn{a} \eqn{\sim}{~} uniform[alim[1], alimb[2]]\cr
#'    \eqn{nu} \eqn{\sim}{~}  dim(data)-1 + exp(z) \cr
#'    \eqn{z} \eqn{\sim}{~}  uniform[dim(data)-1+nulim[1], nulim[2]]\cr
#'    \eqn{v} \eqn{\sim}{~} uniform[vlim[1], vlim[2]]
#'    
#'    \eqn{alpha} \eqn{\sim}{~} \eqn{(1-(alpha-alphamin) / (alphamax-alphamin))^{power}} \cr
#'    alpha = alphamin then expected number of components = \code{Istarmin} \cr
#'    alpha = alphamax then expected number of components = \code{Istarmax}
#'        
#'    \eqn{Z} should NOT include an intercept and is centered for ease of interpretation. The mean of each of the \code{nlgt} \eqn{\beta}s is the mean of the normal mixture.  Use \code{summary()} to compute this mean from the \code{compdraw} output.
#'    
#'    We parameterize the prior on \eqn{\Sigma_i} such that \eqn{mode(\Sigma) = nu/(nu+2) vI}. The support of nu enforces a non-degenerate IW density; \eqn{nulim[1] > 0}.
#'    
#'    The default choices of alim, nulim, and vlim determine the location and approximate size of candidate "atoms" or possible normal components. The defaults are sensible given a reasonable scaling of the X variables. You want to insure that alim is set for a wide enough range of values (remember a is a precision parameter) and the v is big enough to propose Sigma matrices wide enough to cover the data range.  
#'    
#'    A careful analyst should look at the posterior distribution of a, nu, v to make sure that the support is set correctly in alim, nulim, vlim.  In other words, if we see the posterior bunched up at one end of these support ranges, we should widen the range and rerun.  
#'    
#'    If you want to force the procedure to use many small atoms, then set nulim to consider only large values and set vlim to consider only small scaling constants.  Set alphamax to a large number.  This will create a very "lumpy" density estimate somewhat like the classical Kernel density estimates. Of course, this is not advised if you have a prior belief that densities are relatively smooth.
#'  }
#' }
#' 
#' 
#' @author Federico Bumbaca, Leeds School of Business, University of Colorado Boulder, \email{federico.bumbaca@colorado.edu}
#' @references Bumbaca, F. (Rico), Misra, S., & Rossi, P. E. (2020). Scalable Target Marketing: Distributed Markov Chain Monte Carlo for Bayesian Hierarchical Models. Journal of Marketing Research, 57(6), 999-1018.
#' 
#' @examples
#' 
#'\donttest{
#' set.seed(66)
#' R = 500
#' 
#' p = 3                                # num of choice alterns
#' ncoef = 3  
#' nlgt = 1000                           # num of cross sectional units
#' nz = 2
#' nvar = 3
#' Z = matrix(runif(nz*nlgt), ncol=nz)
#' Z = t(t(Z)-apply(Z,2,mean))          # demean Z
#' ncomp = 3                            # no of mixture components
#' Delta=matrix(c(1,0,1,0,1,2), ncol=2)
#' 
#' comps = NULL
#' comps[[1]] = list(mu=c(0,-1,-2),   rooti=diag(rep(2,3)))
#' comps[[2]] = list(mu=c(0,-1,-2)*2, rooti=diag(rep(2,3)))
#' comps[[3]] = list(mu=c(0,-1,-2)*4, rooti=diag(rep(2,3)))
#' pvec=c(0.4, 0.2, 0.4)
#' 
#' ##  simulate from MNL model conditional on X matrix
#' simmnlwX = function(n,X,beta) {
#'   k = length(beta)
#'   Xbeta = X%*%beta
#'   j = nrow(Xbeta) / n
#'   Xbeta = matrix(Xbeta, byrow=TRUE, ncol=j)
#'   Prob = exp(Xbeta)
#'   iota = c(rep(1,j))
#'   denom = Prob%*%iota
#'   Prob = Prob / as.vector(denom)
#'   y = vector("double", n)
#'   ind = 1:j
#'   for (i in 1:n) {
#'     yvec = rmultinom(1, 1, Prob[i,])
#'     y[i] = ind%*%yvec}
#'   return(list(y=y, X=X, beta=beta, prob=Prob))
#' }
#' 
#' ## simulate data with a mixture of 3 normals
#' simlgtdata = NULL
#' ni = rep(50,nlgt)
#' for (i in 1:nlgt) {
#'   betai = Delta%*%Z[i,] + as.vector(bayesm::rmixture(1,pvec,comps)$x)
#'   Xa = matrix(runif(ni[i]*p,min=-1.5,max=0), ncol=p)
#'   X = bayesm::createX(p, na=1, nd=NULL, Xa=Xa, Xd=NULL, base=1)
#'   outa = simmnlwX(ni[i], X, betai)
#'   simlgtdata[[i]] = list(y=outa$y, X=X, beta=betai)
#' }
#' 
#' ## plot betas
#'
#'   old.par = par(no.readonly = TRUE)
#'   bmat = matrix(0, nlgt, ncoef)
#'   for(i in 1:nlgt) { bmat[i,] = simlgtdata[[i]]$beta }
#'   par(mfrow = c(ncoef,1))
#'   for(i in 1:ncoef) { hist(bmat[,i], breaks=30, col="magenta") }
#'   par(old.par)
#' 
#' ## set Data and Mcmc lists
#' keep = 5
#' Mcmc1 = list(R=R, keep=keep)
#' Prior1=list(ncomp=1)
#' Data1 = list(p=p, lgtdata=simlgtdata, Z=Z)
#' Data2 = partition_data(Data = list(Data1), s = 1)
#' 
#' out = parallel::mclapply(Data2, FUN = rhierMnlDPParallel,
#' Prior = Prior1, Mcmc = Mcmc1, mc.cores = 1, mc.set.seed = FALSE)
#'}
#' 
#' @export
rhierMnlDPParallel=function(Data,Prior,Mcmc, verbose = FALSE){
  #
  #  created 3/08 by Rossi from rhierMnlRwMixture adding DP draw for to replace finite mixture of normals
  #
  # revision history:
  #   changed 12/17/04 by rossi to fix bug in drawdelta when there is zero/one unit
  #   in a mixture component
  #   added loglike output, changed to reflect new argument order in llmnl, mnlHess 9/05
  #   changed weighting scheme to (1-w)logl_i + w*Lbar (normalized) 12/05
  #   3/07 added classes
  #   W. Taylor 4/15 - added nprint option to MCMC argument
  #
  # purpose: run hierarchical mnl logit model with mixture of normals 
  #   using RW and cov(RW inc) = (hess_i + Vbeta^-1)^-1
  #   uses normal approximation to pooled likelihood
  #
  # Arguments:
  #   Data contains a list of (p,lgtdata, and possibly Z)
  #      p is number of choice alternatives
  #      lgtdata is a list of lists (one list per unit)
  #          lgtdata[[i]]=list(y,X)
  #             y is a vector indicating alternative chosen
  #               integers 1:p indicate alternative
  #             X is a length(y)*p x nvar matrix of values of
  #               X vars including intercepts
  #             Z is an length(lgtdata) x nz matrix of values of variables
  #               note: Z should NOT contain an intercept
  #   Prior contains a list of (deltabar,Ad,lambda_hyper,Prioralpha)
  #       alpha: starting value
  #       lambda_hyper: hyperparms of prior on lambda
  #       Prioralpha: hyperparms of alpha prior; a list of (Istarmin,Istarmax,power)
  #       if elements of the prior don't exist, defaults are assumed
  #   Mcmc contains a list of (s,c,R,keep,nprint)
  #
  # Output:  as list containing
  #   Deltadraw R/keep  x nz*nvar matrix of draws of Delta, first row is initial value
  #   betadraw is nlgt x nvar x R/keep array of draws of betas
  #   probdraw is R/keep x 1 matrix of draws of probs of mixture components
  #   compdraw is a list of list of lists (length R/keep)
  #      compdraw[[rep]] is the repth draw of components for mixtures
  #   loglike  log-likelikelhood at each kept draw
  #
  # Priors:
  #    beta_i = D %*% z[i,] + u_i
  #       vec(D)~N(deltabar)
  #       u_i ~ N(theta_i)
  #       theta_i~G
  #       G|lambda,alpha ~ DP(G|G0(lambda),alpha)
  #
  #        lambda:
  #           G0 ~ N(mubar,Sigma (x) Amu^-1)
  #           mubar=vec(mubar)
  #           Sigma ~ IW(nu,nu*v*I)  note: mode(Sigma)=nu/(nu+2)*v*I
  #           mubar=0
  #           amu is uniform on grid specified by alim
  #           nu is log uniform, nu=d-1+exp(Z) z is uniform on seq defined bvy nulim
  #           v is uniform on sequence specificd by vlim
  #
  #        Prioralpha:
  #           alpha ~ (1-(alpha-alphamin)/(alphamax-alphamin))^power
  #           alphamin=exp(digamma(Istarmin)-log(gamma+log(N)))
  #           alphamax=exp(digamma(Istarmax)-log(gamma+log(N)))
  #           gamma= .5772156649015328606
  #
  # MCMC parameters
  #   s is the scaling parameter for the RW inc covariance matrix; s^2 Var is inc cov
  #      matrix
  #   w is parameter for weighting function in fractional likelihood
  #      w is the weight on the normalized pooled likelihood 
  #   R is number of draws
  #   keep is thinning parameter, keep every keepth draw
  #   nprint - print estimated time remaining on every nprint'th draw
  #--------------------------------------------------------------------------------------------------
  
  llmnlFract=
    function(beta,y,X,betapooled,rootH,w,wgt){
      z=as.vector(rootH%*%(beta-betapooled))
      return((1-w)*bayesm::llmnl(beta,y,X)+w*wgt*(-.5*(z%*%z)))
    }
  # Check arguments
  if (missing(Data)) {stop("Requires Data argument -- list of p, lgtdata, and (possibly) z", .call = FALSE)}
  if (is.null(Data$p)) {stop("Requires Data element p (# choice alternatives)", .call = FALSE)}
  p = Data$p
  if (is.null(Data$lgtdata)) {stop("Requires Data element lgtdata (list of data for each unit)", .call = FALSE)}
  lgtdata = Data$lgtdata
  nvar=ncol(Data$lgtdata[[1]]$X)
  nlgt = length(lgtdata)
  
  drawdelta = TRUE
  if (is.null(Data$Z)) {
    if(verbose) {message("Z not specified")}
    drawdelta = FALSE
  } else {
    if (nrow(Data$Z) != nlgt) {
      stop(paste("Nrow(Z) ", nrow(Z), "ne number logits ", nlgt), .call = FALSE)
    } else {
      Z = Data$Z
    }
  }
  
  if (drawdelta) {
    nz = ncol(Z)
  }
  
  # Check lgtdata for validity
  ny = 0
  nX = 0
  for (i in 1:nlgt) {
    ny = ny + length(lgtdata[[i]]$y)
    nX = nX + nrow(lgtdata[[i]]$X)
  }
  
  indy = 1
  indX = 1
  ypooled = vector(length = ny)
  Xpooled = matrix(0, nrow = nX, ncol = nvar)
  
  if (!is.null(lgtdata[[1]]$X)) {
    oldncol = ncol(lgtdata[[1]]$X)
  }
  
  for (i in 1:nlgt) {
    if (is.null(lgtdata[[i]]$y)) {stop(paste("Requires element y of lgtdata[[", i, "]]"), .call = FALSE)}
    if (is.null(lgtdata[[i]]$X)) {stop(paste("Requires element X of lgtdata[[", i, "]]"), .call = FALSE)}
    ypooled[indy:(indy + length(lgtdata[[i]]$y) - 1)] = lgtdata[[i]]$y
    indy = indy + length(lgtdata[[i]]$y)
    nrowX = nrow(lgtdata[[i]]$X)
    if ((nrowX / p) != length(lgtdata[[i]]$y)) {
      stop(paste("nrow(X) ne p*length(yi); exception at unit", i), .call = FALSE)
    }
    newncol = ncol(lgtdata[[i]]$X)
    if (newncol != oldncol) {stop(paste("All X elements must have same # of cols; exception at unit", i), .call = FALSE)}
    Xpooled[indX:(indX + nrow(lgtdata[[i]]$X) - 1), ] = lgtdata[[i]]$X
    indX = indX + nrow(lgtdata[[i]]$X)
    oldncol = newncol
  }
  
  nvar = ncol(Xpooled)
  levely = as.numeric(levels(as.factor(ypooled)))
  if (length(levely) != p) {stop(paste("y takes on ", length(levely), " values -- must be = p"))}
  bady = FALSE
  for (i in 1:p) {
    if (levely[i] != i) bady = TRUE
  }
  
  if(verbose){
  message("Table of Y values pooled over all units")
  print(table(ypooled))
  if (bady) {stop("Invalid Y")}
  }
  # check on prior
  
  alimdef=ScalableBayesmConstant.DPalimdef
  nulimdef=ScalableBayesmConstant.DPnulimdef
  vlimdef=ScalableBayesmConstant.DPvlimdef
  
  if(missing(Prior)) {Prior=NULL}
  
  if(is.null(Prior$lambda_hyper)) {lambda_hyper=list(alim=alimdef,nulim=nulimdef,vlim=vlimdef)}
  else {lambda_hyper=Prior$lambda_hyper;
  if(is.null(lambda_hyper$alim)) {lambda_hyper$alim=alimdef}
  if(is.null(lambda_hyper$nulim)) {lambda_hyper$nulim=nulimdef} 
  if(is.null(lambda_hyper$vlim)) {lambda_hyper$vlim=vlimdef}
  }
  if(is.null(Prior$Prioralpha)) {Prioralpha=list(Istarmin=ScalableBayesmConstant.DPIstarmin,Istarmax=min(50,0.1*nlgt),power=ScalableBayesmConstant.DPpower)}
  else {Prioralpha=Prior$Prioralpha;
  if(is.null(Prioralpha$Istarmin)) {Prioralpha$Istarmin=ScalableBayesmConstant.DPIstarmin} else {Prioralpha$Istarmin=Prioralpha$Istarmin}
  if(is.null(Prioralpha$Istarmax)) 
  {Prioralpha$Istarmax=min(50,0.1*nlgt)} else {Prioralpha$Istarmax=Prioralpha$Istarmax}
  if(is.null(Prioralpha$power)) {Prioralpha$power=ScalableBayesmConstant.DPpower}
  }	
  gamma= ScalableBayesmConstant.gamma
  Prioralpha$alphamin=exp(digamma(Prioralpha$Istarmin)-log(gamma+log(nlgt)))
  Prioralpha$alphamax=exp(digamma(Prioralpha$Istarmax)-log(gamma+log(nlgt)))
  Prioralpha$n=nlgt
  #
  # Check Prior arguments for validity
  if (lambda_hyper$alim[1] < 0) {stop("alim[1] must be >0", .call = FALSE)}
  if (lambda_hyper$nulim[1] < 0) {stop("nulim[1] must be >0", .call = FALSE)}
  if (lambda_hyper$vlim[1] < 0) {stop("vlim[1] must be >0", .call = FALSE)}
  if (Prioralpha$Istarmin < 1) {stop("Prioralpha$Istarmin must be >= 1", .call = FALSE)}
  if (Prioralpha$Istarmax <= Prioralpha$Istarmin) {stop("Prioralpha$Istarmin must be < Prioralpha$Istarmax", .call = FALSE)}
  
  if (is.null(Prior$Ad) & drawdelta) {Ad = ScalableBayesmConstant.A * diag(nvar * nz)} else {Ad = Prior$Ad}
  if (drawdelta) {
    if (ncol(Ad) != nvar * nz | nrow(Ad) != nvar * nz) {stop("Ad must be nvar*nz x nvar*nz", .call = FALSE)}
  }
  if (is.null(Prior$deltabar) & drawdelta) {deltabar = rep(0, nz * nvar)} else {deltabar = Prior$deltabar}
  if (drawdelta) {
    if (length(deltabar) != nz * nvar) {stop("deltabar must be of length nvar*nz", .call = FALSE)}
  }
  
  # Check on MCMC
  if (missing(Mcmc)) {
    stop("Requires Mcmc list argument", .call = FALSE)
  } else {
    if (is.null(Mcmc$s)) {s = ScalableBayesmConstant.RRScaling / sqrt(nvar)} else {s = Mcmc$s}
    if (is.null(Mcmc$w)) {w = ScalableBayesmConstant.w} else {w = Mcmc$w}
    if (is.null(Mcmc$keep)) {keep = ScalableBayesmConstant.keep} else {keep = Mcmc$keep}
    if (is.null(Mcmc$nprint)) {nprint = ScalableBayesmConstant.nprint} else {nprint = Mcmc$nprint}
    if (nprint < 0) {stop('nprint must be an integer greater than or equal to 0', .call = FALSE)}
    if (is.null(Mcmc$maxuniq)) {maxuniq = ScalableBayesmConstant.DPmaxuniq} else {maxuniq = Mcmc$maxuniq}
    if (is.null(Mcmc$gridsize)) {gridsize = ScalableBayesmConstant.DPgridsize} else {gridsize = Mcmc$gridsize}
    if (is.null(Mcmc$R)) {stop("Requires R argument in Mcmc list", .call = FALSE)} else {R = Mcmc$R}
  }
  
  if(verbose){
  # Print out problem
  cat("Starting MCMC Inference for Hierarchical Logit:")
  cat("Dirichlet Process Prior")
  cat(paste(" ", p, " alternatives; ", nvar, " variables in X"))
  cat(paste("For ", nlgt, " cross-sectional units"))
  cat(" ")
  
  cat("Prior Parameters:")
  cat("G0 ~ N(mubar, Sigma (x) Amu^-1)")
  cat("mubar = 0")
  cat("Sigma ~ IW(nu, nu*v*I)")
  cat(paste("Amu ~ uniform[", lambda_hyper$alim[1], ",", lambda_hyper$alim[2], "]"))
  cat(paste("nu ~ uniform on log grid [", nvar - 1 + exp(lambda_hyper$nulim[1]), ",", nvar - 1 + exp(lambda_hyper$nulim[2]), "]"))
  cat(paste("v ~ uniform[", lambda_hyper$vlim[1], ",", lambda_hyper$vlim[2], "]"))
  cat(" ")
  
  cat("alpha ~ (1 - (alpha - alphamin) / (alphamax - alphamin))^power")
  cat(paste("Istarmin = ", Prioralpha$Istarmin))
  cat(paste("Istarmax = ", Prioralpha$Istarmax))
  cat(paste("alphamin = ", Prioralpha$alphamin))
  cat(paste("alphamax = ", Prioralpha$alphamax))
  cat(paste("power = ", Prioralpha$power))
  cat(" ")
  
  if (drawdelta) {
    cat("deltabar")
    print(deltabar)
    cat("Ad")
    print(Ad)
  }
  
  cat(" ")
  cat("MCMC Parameters:")
  cat(paste("s =", round(s, 3), " w =", w, " R =", R, " keep =", keep, " nprint =", nprint, " maxuniq =", maxuniq, " gridsize for lambda hyperparameters =", gridsize))
  cat(" ")
  }
  
  
  # Allocate space for draws
  oldbetas = matrix(double(nlgt * nvar), ncol = nvar)
  
  # Initialize compute quantities for Metropolis
  if(verbose){
  print(paste("Initializing Metropolis candidate densities for", nlgt, "units ..."))
  }
  betainit=c(rep(0,nvar))
  #
  #  compute pooled optimum
  #
  out=optim(betainit,bayesm::llmnl,method="BFGS",control=list( fnscale=-1,trace=0,reltol=1e-6), 
            X=Xpooled,y=ypooled)
  betapooled=out$par
  H=bayesm::mnlHess(betapooled,ypooled,Xpooled)
  rm(Xpooled)
  invisible(gc())
  rootH=chol(H)
  #Revise
  hess = vector(mode="list", length=nlgt)
  #
  # initialize betas for all units
  #
  # Rico: bug fix for when hessian is not positive definite
  #
  if(is.null(Prior$hessdelta)) {hessdelta=0} else {hessdelta=Prior$hessdelta}
  for (i in 1:nlgt) 
  {
    wgt=length(lgtdata[[i]]$y)/length(ypooled)
    out=optim(betapooled,llmnlFract,method="BFGS",control=list( fnscale=-1,trace=0,reltol=1e-4), 
              X=lgtdata[[i]]$X,y=lgtdata[[i]]$y,betapooled=betapooled,rootH=rootH,w=w,wgt=wgt)
    if(out$convergence == 0) 
    { #hess=bayesm::mnlHess(out$par,lgtdata[[i]]$y,lgtdata[[i]]$X)+hessdelta*diag(nvar)
      #Revise Check:
      hess=bayesm::mnlHess(out$par,lgtdata[[i]]$y,lgtdata[[i]]$X)
      lgtdata[[i]]=c(lgtdata[[i]],list(converge=1,betafmle=out$par,hess=hess)) }
    else
    { lgtdata[[i]]=c(lgtdata[[i]],list(converge=0,betafmle=c(rep(0,nvar)),
                                       hess=diag(nvar))) }
    oldbetas[i,]=lgtdata[[i]]$betafmle
    if((i%%50==0) & verbose) cat("  completed unit #",i,fill=TRUE)
    #bayesm::fsh()
  }
  rm(ypooled)
  
  #Initialize placeholders when drawdelta == FALSE
  if(drawdelta){
    olddelta = rep(0,nz*nvar)
  } else{
    Z = matrix(0)
    deltabar = 0
    Ad = matrix(0)
  }
  
  #  return(out_s)
  draws = rhierMnlDPParallel_rcpp_loop(R, keep, nprint, lgtdata, Z, deltabar, Ad, Prioralpha,
                                       lambda_hyper, drawdelta, nvar, oldbetas, s, maxuniq,
                                       gridsize, ScalableBayesmConstant.A, ScalableBayesmConstant.nuInc,
                                       ScalableBayesmConstant.DPalpha, verbose)  
  return(draws)
}