hbbr (version 1.1.2)

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

Description

Fits processed benefit-risk survey data from an appropriately designed discrete choice experiment to the hbbr (Hierarchical Bayesian Benefit-Risk) model. 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

hbbr.Fit(brdta, 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.

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 for 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 paramaters of interests; default: mcmc=list(burnin=5000, iter=100000, nc=2, thin=20). 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: (bbar.mcmc, bbar.means, bbar.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 in X 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: (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 R2jags package help files.

Examples

# NOT RUN {
## Sample calls: fits pilot response data included with the package
# }
# NOT RUN {
  data(hbbrPilotResp)
  hbfit = hbbr.Fit(brdta=hbbrPilotResp$brdta, design=hbbrPilotResp$design,
                   mcmc=list(burnin=500, iter=10000, nc=2, thin=10))
  hb = hbfit$bbar.mcmc
  dgn = hbfit$design
  mns = hbfit$bbar.means
  sds = hbfit$bbar.sd # same as apply(hbfit$bbar.mcmc, 2, sd)

  ## Plots of MCMC draws ---------------------------------------
  op=par(mfrow=c(1,2), mar = c(4,2,3,1),oma=c(.1,.1,2,.1))
  matplot(hb,type="l",xlab="Iterations",ylab="",
          main=paste("Average Part-Worths (beta-bars)"),
          cex.main=.8, cex.lab=0.8, axes=FALSE)
  axis(1, at=seq(0,dim(hb)[1],length.out = 6),
          labels= paste(seq(0,5,1)*dim(hb)[1]/5 *hbfit$other.inputs$thin, sep=""),
          cex.axis=0.8)
  axis(2, cex.axis=0.8,las=1)

  plot(hbfit$logL, type="l",main="Log Likelihood", axes=FALSE,xlab="Iterations",ylab="",
      cex.main=.8,cex.lab=0.8)
  axis(1, at=seq(0,dim(hb)[1],length.out = 6),
      labels= paste(seq(0,5,1)*dim(hb)[1]/5 *hbfit$other.inputs$thin, sep=""),
       cex.axis=0.8)
  axis(2, cex.axis=0.8,las=1)
  title(outer=TRUE, main = paste("MCMC draws plotted at every ",
       hbfit$other.inputs$thin,"-th Iteration",sep=""),cex.main=.9)
  par(op)
  
  ## Plots for mean estimated part-worth utilities ------------------
  require(ggplot2)
  require(gridExtra)

  b.mns = c()
  b.sds = c()
  b.atr = c()
  b.lvl = c()
  j.now=1
  for (j in 1:dgn$b) {
    b.mns = c(b.mns,0, mns[j.now:(j.now-1+dgn$bl[j]-1)])
    b.sds = c(b.sds,0, sds[j.now:(j.now-1+dgn$bl[j]-1)])
    b.atr = c(b.atr, rep(dgn$blbls[j], dgn$bl[j]))
    b.lvl = c(b.lvl, paste("E", 1:dgn$bl[j],sep=""))
    j.now = j.now-1+dgn$bl[j]
  }

  r.mns = c()
  r.sds = c()
  r.atr = c()
  r.lvl = c()
  k.now=j.now
  for (k in 1:dgn$r) {
    r.mns = c(r.mns,0,mns[k.now:(k.now-1+dgn$rl[k]-1)])
    r.sds = c(r.sds,0, sds[k.now:(k.now-1+dgn$rl[k]-1)])
    r.atr = c(r.atr, rep(dgn$rlbls[k], dgn$rl[k]))
    r.lvl = c(r.lvl, paste("H", 1:dgn$rl[k],sep=""))
    k.now = k.now-1+dgn$rl[k]
  }

  d0.b = data.frame(Attributes =b.atr, lvl=b.lvl, util = b.mns, se = b.sds)
  d0.r = data.frame(Attributes =r.atr, lvl=r.lvl, util = r.mns, se = r.sds)
  y.max = max(abs(mns) + max(sds))
  pd <- position_dodge(0.2) # move them .2 to the left and right

  pb = ggplot(data = d0.b, aes(x=lvl, y=util, group=Attributes,color=Attributes)) +
    ylim(0, y.max) +
    geom_hline(yintercept = 0) +
    geom_line(size=1.5, position=pd) +
    geom_point(size=4, shape=22, fill="green",color="darkgreen", position=pd) +
    geom_errorbar(aes(ymin=util-se, ymax=util+se), width=0.2, position=pd) +
    xlab("Benefit-Attribute Levels") + ylab("Estimated Utility") +
    ggtitle("Estimated Partworth Utilities of Benefits") +
    scale_color_manual(values=c("deepskyblue3" , "#9999CC", "cyan3" )) +
    theme(legend.position="bottom",plot.title = element_text(size = 10))

  pr = ggplot( data = d0.r, aes(x=lvl, y=util, group=Attributes,color=Attributes)) +
    ylim(-y.max,0)+
    geom_hline(yintercept = 0) +
    geom_line(size=1.5, position=pd) +
    geom_point(size=4, shape=22, fill="pink",color="darkred", position=pd) +
    geom_errorbar(aes(ymin=util-se, ymax=util+se), width=0.2, position=pd) +
    xlab("Risk-Attribute Levels") + ylab("Estimated Utility") +
    ggtitle("Estimated Partworth Utilities of Risks") +
    scale_color_manual(values=c("orange" , "maroon" )) +
    theme(legend.position="bottom",plot.title = element_text(size = 10))

  grid.arrange(pb, pr, nrow = 1)

##------------------------------------------------------------------
# }
# NOT RUN {
# }