hbbr (version 1.1.2)

hbbrAug.Fit: hbbrAug.Fit (Fits processed response data to the augmented hbbr model)

Description

Fits processed benefit-risk survey data from an appropriately designed discrete choice experiment to the augmented hbbr model that includes patients' baseline characteristics. For details see article by Mukhopadhyay, S., Dilley, K., Oladipo, A., & Jokinen, J. (2019). Hierarchical Bayesian Benefit<U+2013>Risk Modeling and Assessment Using Choice Based Conjoint. Statistics in Biopharmaceutical Research, 11(1), 52-60.

Usage

hbbrAug.Fit(brdta, Z, design, tune.param = list(tau = 0.01, eta = NULL,
  df.add = 2), mcmc = list(burnin = 5000, iter = 1e+05, nc = 2, thin =
  20), verbose = TRUE)

Arguments

brdta

processed and coded survey response data to be fitted to the hbbr model. It is a data frame in which 1st two columns indicate subject id and subject response (y = 0 or 1), and remaining columns contain information on design matrix (X). See Details below for more information.

Z

matrix of observed baseline characteristics of the patients. If there are N patients responded to the survey and we have included g characteristics of each patient then Z is a matrix of (g+1) x N, with all elements of the first column equal to 1. Note that when g=0, the model reduces to regular hbbr model.

design

design information of the experiment: design = list(b, r, bl, rl, blbls, rlbls) where, b is number of benefit attributes, r is number of risk attributes, bl and rl are vectors of integers of length b, and r indicating number of levels in j-th benefit attribute and k-th risk attribute, respectively. blbls, rlbls consists of labels of benefit and risk attributes. When blbls is NULL, it uses "B1", "B2", ... and similarly for rlbls.

tune.param

a list of tuning hyper-parameters to be used; default tune.param=list(tau=0.01, eta=NULL). See Details below for more information.

mcmc

a list of mcmc parameters to be used in the Gibbs sampler to obtain posterior samples of the parameters of interests; default: mcmc=list(burnin=1000, iter=10000, nc=4, thin=10). See Details below for more information.

verbose

TRUE or FALSE: flag indicating whether to print intermediate output summary which might be helpful to see convergence results.

Value

returns a list of useful output of interest and input specifications: (del.mcmc, del.means, del.sds, summary, logL, design, model, brdata, other.inputs).

Details

brdta is a processed and coded survey response data to be fitted to the hbbr model. It is a data frame in which 1st column contains ID of respondent, 2nd column contains response (y = 0 or 1) - each value corresponds to each choice-pair card evaluated by the respondent: y =1 if the 1st choice of the pair was preferred; 0 otherwise, 3rd column onwards contain information on design matrix (X). Each row of X is a vector of indicator variables taking values 0, 1, or -1; a value of 0 is used to denote absence of an attribute level; a value of 1 or -1 is used to indicate presence of an attribute level in the 1st choice, or in the 2nd choice, respectively in the choice-pair presented to the respondent. Note that column corresponding to the 1st level for each attribute would not be included as the part-worth parameter (beta) for the 1st level of each attribute is assumed to be 0 without loss of generality. So, if there are b benefit attributes and r risk attributes, and then have bl_j and rl_k levels (j=1,...,b; k=1,...,r) then total number of columns brdta is Sum_over_j(bl_j-1) + Sum_over_k(rl_k-1). If there are B respondents each responding to k choice-pairs then brdta will have B*k rows.

tune.param is a list of tuning hyper-parameters (tau, eta) for the hbbr model. Specifically, in the hbbr model beta.h ~ MVN(beta.bar, V.beta) where the hyper-prior of beta.bar is assumed to be MVN (beta0, B) with B = 1/tau*I; and hyper-prior of V.beta is assumed to follow inverse Wishart IW(nue, V) with V = 1/eta*I. When eta is NULL then eta will take the default value of m+3 which is the DF for the Wishart distribution. If we think the respondents have very similar part-worth vectors, then use eta=1.

mcmc is a list of MCMC specification parameters to be used for rjags package: (a) burnin - contains the number of burn-in values to be generated, (b) iter - is the total number of iterations of each chain beyond burn-in, (c) nc - is the number of independent chains, and (d) thin = posterior samples to be saved for every 'thin' values of the MCMC samples in each of the 'nc' chains. For more details see rjags package help files.

Examples

# NOT RUN {
## Sample calls:
# fits simulated response data included with this package to augmented hbbr model
# and then plots the estimated part-worth utilities.

# }
# NOT RUN {
 data("simAugData")
 hbA = hbbrAug.Fit(brdta= simAugData$brdtaAug, Z=simAugData$Z,
               design=simAugData$design,
               tune.param=list(tau=0.01, eta=NULL, df.add=2),
               mcmc=list(burnin=500, iter=10000, nc=2, thin=10))

 # define an appropriate function to plot the part-worth values...
 partworth.plot = function(attr.lvl, beta.mns, nb=3, new=TRUE, pnt =15, cl=clrs)
 {
 #check dimension
   k = length(attr.lvl) # no of attributes
   bk = length(unlist(attr.lvl)) # no of levels acrosss attributes
   if (bk - k != length(beta.mns)) stop("error 1")
   mns = rep(0, length(unlist(attr.lvl)))
   cntr = 0
   for (j in 1:k)
   {
     for (i in 1:length(attr.lvl[[j]])){
       cntr = cntr +1
       if (i > 1) mns[cntr]= beta.mns[cntr-1-(j-1)]
     }
   }
   indx = list()
   j0=1
   for (j in 1:k) {
     j1 = (j0+length(attr.lvl[[j]])-1)
     indx[[j]]= j0:j1
     j0=j1+1
  }
   if (new) {
     plot(c(1,bk), c(floor(min(beta.mns)*1.2),ceiling(max(beta.mns)*1.2)),
                                        type="n", axes=FALSE, xlab="",ylab="")
     axis(2, at=0:ceiling(max(beta.mns)*1.2), las=1, cex.axis=.7)
     axis(4, at=floor(min(beta.mns)*1.2):0, las=1, cex.axis=.7)
   }
   vl=c()
   for (j in 1:k)
   {
     points(indx[[j]], mns[indx[[j]]], type="b", pch =pnt, col=cl[j])
     vl=c(vl, max(indx[[j]])+.5)
   }
   abline(v=vl,col="gray", h=0)
   box()
 }

 # Plotting estimated betas (part-worth) for some selected baseline characteristics:

 augattr.lvl = list(b1=paste("B1",1:3,sep=""),b2=paste("B2",1:3,sep=""),
          r1=paste("R1",1:3,sep=""),r2=paste("R2",1:3,sep=""))
 clrs = c("blue", "green4","orange4", "red3")

 mns = hbA$del.means
 # est. part-worth values
 betmn1 = mns %*% matrix(c(1, 0, 1), ncol=1)  # at mean age with disease staus=1
 betmn2 = mns %*% matrix(c(1, 0, -1), ncol=1) # at mean age with disease staus=-1
 betmn3 = mns %*% matrix(c(1, 1, -1), ncol=1) # at age = mean+1*SD, disease staus=-1

 partworth.plot(attr.lvl = augattr.lvl, beta.mns = betmn1)
 partworth.plot(attr.lvl = augattr.lvl, beta.mns = betmn2, new=FALSE, pnt=17)
 partworth.plot(attr.lvl = augattr.lvl, beta.mns = betmn3, new=FALSE, pnt=16)

 # Plotting true betas at those baseline characteristics
 Del = simAugData$Del
 clrs = rep("darkgrey", 4)
 # true part-worth values
 bmn1 = Del %*% matrix(c(1, 0, 1), ncol=1)  # at mean age with disease staus=1
 bmn2 = Del %*% matrix(c(1, 0, -1), ncol=1) # at mean age with disease staus=-1
 bmn3 = Del %*% matrix(c(1, 1, -1), ncol=1) # at age = mean+1*SD, disease staus=-1
 partworth.plot(attr.lvl = augattr.lvl, beta.mns = bmn1)
 partworth.plot(attr.lvl = augattr.lvl, beta.mns = bmn2, new=FALSE, pnt=17)
 partworth.plot(attr.lvl = augattr.lvl, beta.mns = bmn3, new=FALSE, pnt=16)
# }
# NOT RUN {
# }