25.1.1 Essential points

Recall the basic steps of a Bayesian analysis from Section 2.3 (p. 25): Identify the data, define a descriptive model, specify a prior, compute the posterior distribution, interpret the posterior distribution, and, check that the model is a reasonable description of the data. Those steps are in logical order, with each step building on the previous step. That logical order should be preserved in the report of the analysis. The essential points listed below mirror the basic mechanical steps of Bayesian analysis, preceded by an important preliminary step that might recede in importance as Bayesian methods become standard procedure.

• Motivate the use of Bayesian (non-NHST) analysis: Many audiences, including journal editors and reviewers, are comfortable with NHST and are not familiar with Bayesian methods. I have found that most editors and reviewers are eager to publish research that uses the best possible methods, including Bayesian analyses, but the editors and reviewers appreciate an explanation of why you have used Bayesian analysis instead of NHST. You may motivate your use of Bayesian data analysis on several grounds, depending in part on the particular application and audience. For example, Bayesian models are designed to be appropriate to the data structure, without having to make approximation assumptions typical in NHST (e.g., homogeneity of variances across groups, and normally distributed noise). The inferences from a Bayesian analysis are richer and more informative than NHST because the posterior distribution reveals joint probabilities of combinations of parameter values. And, of course, there is no reliance on sampling distributions and p values to interpret the parameter estimates.

• Clearly describe the data structure, the model, and the model's parameters: Ultimately you want to report the meaningful parameter values, but you can't do that until you explain the model, and you can't do that until you explain the data being modeled. Therefore, recapitulate the data structure, reminding your reader of the predicted and predictor variables. Then describe the model, emphasizing the meaning of the parameters. This task of describing the model can be arduous for complex hierarchical models, but it is necessary and crucial if your analysis is to mean anything to your audience.

• Clearly describe and justify the prior: It is important to convince your audience that your prior is appropriate and does not predetermine the outcome. The prior should be amenable to a skeptical audience. The prior should be at least mildly informed to match the scale of the data being modeled. If there is copious previous research using very similar methods, it should not be ignored. Optionally, as mentioned again below, it may be helpful to try different priors and report the robustness of the posterior. When the goal is continuous parameter estimation, the posterior distribution is usually robust against reasonable changes in vague priors, but when the goal is model comparison (such as with inclusion coefficients) then the posterior probabilities of the models can be strongly affected by the choice priors.

• Report the MCMC details, especially evidence that the chains were converged and of sufficient length. Section 7.5 (p. 178) explained in detail how to examine the chains for convergence, and every program in this book produces diagnostic graphics of chains for at least some of the parameters. The programs all produce summaries of the parameters that include the ESS. Your report should indicate that the chains were checked for convergence, and your report should indicate the ESS of the relevant parameters. Usually this section of the report can be brief.

• Interpret the posterior: Many models have dozens or even hundreds of parameters, and therefore it is impossible to summarize all of them. The choice of which parameters or contrasts to report is driven by domain-specific theory and by the results themselves. You want to report the parameters and contrasts that are theoretically meaningful. You can report the posterior central tendency of a parameter and its HDI in text alone; histograms of posteriors are useful for the analyst to understand the posterior and for explanation in a textbook, but may be unnecessary in a concise report. Describe effects of shrinkage if appropriate. If your model includes interactions of predictors, be careful how you interpret lower-order effects. Finally, ifyou are using a ROPE for decisions, justify its limits.

25.1.2 Optional points

The following points are not necessarily crucial to address in every report, but should be considered. Whether or not to include these points depends on the particulars of the application domain, the points you want to make, and the audience to which the report is being addressed.

• Robustness of the posterior for different priors: When there is contention about the prior, it can be most convincing simply to conduct the analysis with different priors and demonstrate that the essential conclusions from the posterior do not change. This may be especially important when doing model comparisons, such as when using inclusion coefficients (or using Bayes factors to assess null values). Which priors should be used? Those that are meaningful and amenable to your audience, such as reviewers and editors of the journal to which the report is submitted.

• Posterior predictive check: By generating simulated data from the credible parameter values of the model, and examining the qualities of the simulated data, the veracity of the model can be further bolstered, if the simulated data do resemble the actual data. On the other hand, if the simulated data are discrepant from the actual data in systematic and interpretable ways, then the posterior predictive check can inspire new research and new models. For a perspective on posterior predictive checks, see the article by Kruschke (2013b) and Section 17.5.1 (among others) of this book.

• Power analysis: If there is only a weak effect in your results, what sample size would be needed to achieve some desired precision in the estimate? If you found a credibly nonzero difference, what was the retrospective power of your experiment and what is its replication power? This sort of information can be useful not only for the researcher's own planning, but it can also be useful to the audience of the report to anticipate potential follow-up research and to assess the robustness of the currently reported results. Chapter 13 described power analysis in depth.

25.1.3 Helpful points

Finally, to help science be cumulative, make your results available on the web:

• Post the raw data: There are at least two benefits of posting the original data. One benefit is that subsequent researchers can analyze the data with different models. New insights can be gained by alternative modeling interpretations. The longevity of the original research is enhanced. A second benefit is that if an exact or near-exact replication is conducted, the original data set can be concatenated with the new data set, to enhance sensitivity of the new data. Either way, your work gets cited!

• Post the MCMC sample of the posterior: There are at least two benefits of making the posterior publicly available. One is that other researchers can explore the posterior for effects and comparisons that were not in the report. Complex models have many parameters, and no single report can cover every possible perspective on the posteriordistribution. The longevity and impact of the research is thereby enhanced. A second benefit is that if subsequent researchers do follow-up research with a similar design and model, then the posted posterior can inform the prior of the subsequent analysis. Because the full posterior automatically incorporates all the covariations between all the parameters, the full posterior can be more useful than summaries of marginal distributions in a report. Either way, your work gets cited!

25.2.1 R code for computing hdi of a grid approximation

We can imagine the grid approximation of a distribution as a landscape ofpoles sticking up from each point on the parameter grid, with the height of each pole indicating the probability mass at that discrete point. We can imagine the highest density region by visualizing a rising tide: We gradually flood the landscape, monitoring the total mass of the poles that protrude above water, stopping the flood when 95% (say) of the mass remains protruding. The waterline at that moment defines the highest density region (e.g., Hyndman, 1996).

The function, listed below, finds the approximate highest density region in a somewhat analogous way. It uses one extra trick at the beginning, however. It first sorts all the poles in order of height, from tallest to shortest. The idea is to move down the sorted queue ofpoles until the cumulative probability has just barely exceeded 95% (or whatever mass is desired). The resulting height is the “waterline” that defines all points inside the highest density. See the comments in the top of the code for details of how to use the function.

HDIofGrid = function(probMassVec, credMass=0.95) {
 # Arguments:
 # probMassVec is a vector of probability masses at each grid point.
 # credMass is the desired mass of the HDI region.
 # Return value:
 # A list with components:
 # indices is a vector of indices that are in the HDI
 # mass is the total mass of the included indices
 # height is the smallest component probability mass in the HDI
 # Example of use: For determining HDI of a beta(30,12) distribution
 # approximated on a grid:
 # > probDensityVec = dbeta(seq(0,1,length=201), 30, 12)
 # > probMassVec = probDensityVec / sum(probDensityVec)
 # > HDIinfo = HDIofGrid(probMassVec)
 # > show(HDIinfo)
 sortedProbMass = sort(probMassVec, decreasing=TRUE)
 HDIheightIdx = min(which(cumsum(sortedProbMass) >= credMass))
 HDIheight = sortedProbMass[HDIheightIdx] HDImass = sum(probMassVec[probMassVec >= HDIheight])
 return(list(indices = which(probMassVec >= HDIheight), mass = HDImass, height = HDIheight))
}

25.2.2 HDI of unimodal distribution is shortest interval

The algorithms for computing the HDI of an MCMC sample or of a mathematical function rely on a crucial property: For a unimodal probability distribution on a single variable, the HDI of mass M is the narrowest possible interval of that mass. Figure 25.1 illustrates why this is true. Consider the 90% HDI as shown. We construct another interval of 90% mass by moving the limits of the HDI to right, such that each limit ismoved to a point that covers 4%, as marked in grey in Figure 25.1. The new interval must also cover 90%, because the 4% lost on the left is replaced by the 4% gained on the right.

Consider the grey regions in Figure 25.1. Their left edges have the same height, because the left edges are defined by the HDI. Their areas are the same, because, by definition, the areas are both 4%. Notice, however, that the left grey area is narrower than the right grey area, because the left area falls at a point where the distribution is increasing, but the right area falls at a point where the distribution is decreasing. Consequently, the distance between right edges of the two grey zones must be greater than the HDI width. (The exact widths are marked in Figure 25.1.) This argument applies for any size of grey zone, going right or left of the HDI, and for any mass HDI. The argument relies on unimodality, however.

Given the argument and diagram in Figure 25.1, it is not too hard to believe the converse: For a unimodal distribution on one variable, for any mass M,the interval containing mass M that has the narrowest width is the HDI for that mass. The algorithms described below are based on this property of HDIs. The algorithms find the HDI by searching among candidate intervals of mass M. The shortest one found is declared to be the HDI. It is an approximation, of course. See Chen and Shao (1999) for more details, and Chen, He, Shao, and Xu (2003) for dealing with the unusual situation of multimodal distributions.

25.2.3 R code for computing HDI of a MCMC sample

It is important to remember that an MCMC sample is only a random and noisy representation of the underlying distribution, and therefore the HDI found from the MCMC sample is also noisy. For details, review the discussion of HDI approximation error that accompanies Figures 7.13 (p. 185) and 7.14 (p. 186). Below is the code for finding the HDI of an MCMC sample. It is brief and hopefully self-explanatory after the discussion of the previous section.

HDIofMCMC = function(sampleVec, credMass=0.95) {
 # Computes highest density interval from a sample of representative values,
 # estimated as shortest credible interval.
 # Arguments:
 # sampleVec
 # is a vector of representative values from a probability distribution.
 # credMass
 # is a scalar between 0 and 1, indicating the mass within the credible
 # interval that is to be estimated.
 # Value:
 # HDIlim is a vector containing the limits of the HDI sortedPts = sort(sampleVec)
ciIdxInc = ceiling(credMass * length(sortedPts))
nCIs = length(sortedPts) - ciIdxInc
ciWidth = rep(0, nCIs) for(iin 1:nCIs) {
 ciWidth[i]= sortedPts[i + ciIdxInc] - sortedPts[i]
}
HDImin = sortedPts[which.min(ciWidth)]
HDImax = sortedPts[which.min(ciWidth) + ciIdxInc]
HDIlim = c(HDImin, HDImax)
return(HDIlim)
}

25.2.4 R code for computing HDI of a function

The function described in this section finds the HDI of a unimodal probability density function that is specified mathematically in R. For example, the function can find HDI's of normal densities or of beta densities or of gamma densities, because those densities are specified as functions in R.

What the program accomplishes is just a search of HDI's, converging to the shortest one, but it does this by using some commands and R functions that may not have been used elsewhere in the book. One function that the program uses is the inverse cumulative density function (ICDF) for whatever probability distribution is being targeted. We have seen one case of an ICDF previously, namely the probit function, which is the inverse of the cumulative-density function for the normal distribution. In R, the ICDF of the normal is the qnorm(x) function, where the argument x is a value between zero and one. The program for finding the HDI takes, as one of its arguments, R's name for the ICDF of the function. For example, if we want to find an HDI of a normal density, we pass in ICDFname=“qnorm”.

The crucial function called by the program is R's optimize routine. The optimize routine searches and finds the minimum of a specified function over a specified domain. In the program below, we define a function called intervalWidth that returns the width of the interval that starts at lowTailPr and has 95% mass. This intervalWidth function is repeatedly called from the optimize routine until it converges to a minimum.

HDIofICDF = function(ICDFname, credMass=0.95, tol=1e-8,…){
 # Arguments:
 # ICDFname is R's name for the inverse cumulative density function
 # of the distribution.
 # credMass is the desired mass of the HDI region.
 # tol is passed to R's optimize function.
# Return value:
 # Highest density interval (HDI) limits in a vector.
# Example of use: For determining HDI of a beta(30,12) distribution, type
 # > HDIofICDF(qbeta, shape1 = 30, shape2 = 12)
 # Notice that the parameters of the ICDFname must be explicitly named;
 # e.g., HDIofICDF(qbeta, 30, 12) does not work.
 # Adapted and corrected from Greg Snow's TeachingDemos package.
incredMass = 1.0 - credMass
 intervalWidth = function(lowTailPr, ICDFname, credMass,…){ ICDFname(credMass + lowTailPr,…)- ICDFname(lowTailPr, …)
}
optInfo = optimize(intervalWidth, c(0, incredMass), ICDFname=ICDFname, credMass=credMass, tol=tol, …)
HDIlowTailPr = optInfo$minimum
return(c(ICDFname(HDIlowTailPr,…),
 ICDFname(credMass + HDIlowTailPr,…)))
}

25.3.1 Examples

An example was given in Figure 21.11 (p. 640), which had a normal distribution on parameter β₀ transformed to a distribution over parameter μ via a logistic function. The question was, What is the implied distribution on μ?Footnote 2 on p. 639 explained the application of Equation 25.1.

As another example, we can apply Equation 25.1 to the case in which θ = f(ϕ) = logistic(ϕ) = 1/[1+ exp(−ϕ)], and the distribution on ϕ is given by p(ϕ|a, b) = (f (ϕ))^a (1 — f (ϕ))^b / B (a, b). An example of this distribution, for a = 1 and b = 1, is shown in the top-left panel of Figure 25.2. Notice that the derivative of the logistic function f is f′(ϕ) = exp(−ϕ)/[1+ exp(−ϕ)]² = f (ϕ) (1 − f (ϕ)) = θ(1 − θ). Therefore, according to Equation 25.1, the equivalent probability density at θ = f(ϕ) is p(θ) = p(ϕ)/f′(ϕ) = θ^a (1 − θ)B /[θ(1 − θ)B(θ|a, b)] = θ^(a−1) (1 − θ)^(b−1) /B(a, b) = beta(θ|a, b). The upper row of Figure 25.2 shows this situation when a = b = 1. An intuitive way to think of this situation is that the probability on ϕ is dense near ϕ = 0, but that is exactly where the logistic transformation stretches the distribution. On the other hand, the probability on ϕ is sparse at large positive or large negative values, but that is exactly where the logistic transformation compresses the distribution.

As another example, consider a case in which θ = f (ψ) = 1 − exp(−ψ), with the probability density p(ψ) = exp(−ψ). Notice that the derivative of the transformation is f′(ψ) = exp(−ψ), and therefore the equivalent density at θ = f (θ) is p(f(ψ)) = p(ψ)/f′(ψ) = 1. In other words, the equivalent density on θ is the uniform density, as shown in the lower row of Figure 25.2.

As a final example, Figure 25.3 shows a uniform distribution on standard deviation transformed to the corresponding distribution on precision. By definition, precision is the reciprocal of squared standard deviation.

25.3.2 Reparameterization of two parameters

When there is more than one parameter being transformed, the calculus becomes a little more involved. Suppose we have a probability density on a two-parameter space, p(α₁,α₂). Let β₁ = f₁(α₁,α₂) and β₂ = f₂(α₁, α₂). Our goal is to find the probability density p(β₁, β₂) that corresponds to p(α₁, α₂). We do this by considering infinitesimal corresponding regions. Consider a point α₁, α₂. The probability mass of a small region near that point is the density at that point times the area of the small region: p(α₁, α₂) dα₁dα₂. The corresponding region in the transformed parameters should have the same mass. That mass is the density at the transformed point times the area of the region mapped to from the originating region. In vector calculus textbooks, you can find discussions demonstrating that the area of the mapped-to region is$\left|\det \left(J\right)\right|\text{?}d{\alpha }_1\text{?}d{\alpha }_2$ where $J$ is the Jacobian matrix: $J_{r,c}=d{f}_r\left({\alpha }_1,\text{?}{\alpha }_2\right)/d{\alpha }_c$ and $\det \left(J\right)$ is the determinant of the Jacobian matrix. Setting the two masses equal and rearranging yields $p\left({\beta }_1,\text{?}{\beta }_2\right)=p\left({\alpha }_1,\text{?}{\alpha }_2\right)/\left|\det \left(J\right)\right|$. You can see that Equation 25.1 is a special case. As we have had no occasions to apply this transformation, no examples will be provided here.

But the method has been mentioned for those intrepid souls who may wish to venture into the wilderness ofmultivariate probability distributions with nothing more than pen and paper.