This example uses the MCMC procedure to fit a Bayesian multiple linear regression (MLR) model by using a multivariate prior on the regression parameters. It demonstrates how to use existing SAS multivariate density functions for specifying prior distributions.
Researchers are interested in determining the relationship of gestational length and litter size on brain weight after accounting for body weight. They study the evalutionary benefits of larger brain weight after accounting for the coinciding disadvantages of longer gestational length and smaller litter size.
The following data set contains the average brain and body weights in grams and kilograms, respectively, the average gestational length in days, and the litter size for 95 different species as given in:
data brainweight; input body gestation litter @@; log_brain = log(brain); log_body = log(body); log_gestation = log(gestation); datalines; 17.5 3.5 26 1 3.5 0.93 34 4.6 3.15 0.15 46 3 1.14 0.049 51 1.5 ... more lines ... 200 39 180 1 210 66 158 1.2 125 49 150 2.4 106 30 151 2 ;
Figure 1 displays the distribution of each variable. The skewed distributions and the differing orders of magnitude for the minimum and maximum brain weight, body weight, and gestational length suggest the need for the log transformation.
Suppose you want to fit a Bayesian MLR model for the logarithm of brain weight with density as
where is the vector of covariates listed as for species. The African elephant has been omitted from this example because it is easily recognized as an extreme outlier even on the log scale.
The likelihood function for the logarithm of the brain weight and corresponding covariates is
where denotes a conditional probability density. The normal density is evaluated at the specified value of and the corresponding mean parameter , which is defined in Equation 1. The four regression parameters in the likelihood are through .
Suppose you had expert or prior knowledge that some of the covariates were correlated. You might want to use a multivariate prior to incorporate your information. Using the multivariate normal prior, you enable covariates to be independent or correlated a priori. You can also specify the a priori correlation that you believe to be positive or negative. Suppose you thought, a priori, that body weight and gestational age were positively correlated and that body weight and litter size were negatively correlated. Calculate the prior covariance as the product of the prior correlation and standard deviations. More formally, the formula for calculating the covariance of the and covariate is where is the prior correlation and and are the sample standard deviations of the th and th covariates, respectively.
Suppose the prior correlation between the first and second covariate was , and a priori you thought the standard deviations of the regression coefficients for log of body weight and gestational age were and , respectively. The prior covariance of the log of body weight and gestational age can be calculated as . Similarly, suppose the prior correlation between the first and third covariate was and a priori the standard deviation of the litter size was . The prior covariance of the log of body weight and litter size can be calculated as .
Suppose the following prior distributions are placed on the parameters:
where indicates a prior distribution, indicates a multivariate normal prior of dimension four with mean vector and variance matrix , and igamma indicates the density function for the inverse-gamma distribution.
The following SAS statements use the likelihood function and prior distributions to fit the Bayesian MLR model with the multivariate prior. The PROC MCMC statement invokes the procedure and specifies the input data set. The NBI= option specifies the number of burn-in iterations. The NMC= option specifies the number of posterior simulation iterations. The THIN=5 option specifies that one of every five samples is saved in the posterior sample. The SEED= option specifies a seed for the random number generator (the seed guarantees the reproducibility of the random stream). The PROPCOV=QUANEW option uses the estimated inverse Hessian matrix as the initial proposal covariance matrix. The MONITOR= option specifies a list of symbols, which can be either parameters or functions of the parameters in the model, for which inference is to be done. The symbol is a shorthand for all model parameters, in this case, beta0, beta1, beta2, beta3, sig2. The symbol ExpBeta2 and ExpBeta3 are defined in the program for and , respectively.
ods graphics on; proc mcmc data=brainweight nbi=5000 nmc=25000 thin=5 seed=1181 propcov=quanew monitor=(_parms_ ExpBeta2 ExpBeta3); array mu_pr mu_pr0-mu_pr3; array beta beta0-beta3; array data 1 log_body log_gestation litter; begincnst; call zeromatrix(mu_pr); rc = logmpdfset('lp',16,0,16,0,3,2.25,0,-2,0,1); endcnst; parms beta: 0; parms sig2 1; ExpBeta2=exp(beta2); ExpBeta3=exp(beta3); lmvn = logmpdfnormal(of beta0-beta3, of mu_pr0-mu_pr3, 'lp'); prior beta: ~ general(lmvn); prior sig2 ~ igamma(shape = 2.001, scale = 1.001); call mult(beta, data, mu); model log_brain ~ n(mu, var = sig2); run; data _null_; rc = logmpdffree(); put rc=; run; ods graphics off;
Each of the two ARRAY statements associates a name with a list of variables and constants. The first ARRAY statement is used to specify the prior mean. The second ARRAY statement specifies names for the regression coefficients. The third ARRAY statement contains all of the covariates and is used to take advantage of PROC MCMC’s ability to use matrix functions that make the code succinct.
The BEGINCNST and ENDCNST statements are used in conjunction with the DATA step functions. Programming statements enclosed by the BEGINCNST and ENDCNST statements are executed once per procedure call. You use the LOGMPDFSET function within the statements to set up a symmetric positive definite matrix from its lower triangular elements listed in row-major format. The RC assignment statement is set to 0 when the numeric arguments describe a positive definite matrix; otherwise it is set to a nonzero value.
The first PARMS statement places all regression parameters in a single block and assigns them an initial value of 0. The second PARMS statement places the variance parameter in a separate block and assigns it an initial value of 1.
ExpBeta1 and ExpBeta2 are the names for the functions of parameters for and that you want to obtain posterior summary statistics and convergence diagnostics graphs.
The LMVN assignment statement uses the LOGMPDFNORMAL function to calculate the log density of the multivariate normal density given in Equation 3. Next, the PRIOR statement for the parameters uses the GENERAL function, which indicates that you are using a SAS statement to construct a nonstandard density or distribution. The argument is an expression that takes the value of the logarithm function of the prior or likelihood distribution. In this case, the nonstandard density is the multivariate normal distribution and the argument is the logarithm of its density. The second PRIOR statement assigns the inverse-gamma prior distribution to as given in Equation 4.
The CALL statement uses the MULT matrix multiplication subroutine to calculate . The MODEL statement specifies the likelihood function as given in Equation 2.
The DATA step after the call to the MCMC procedure calls the LOGMPDFFREE function. Use the LOGMPFFREE function to free the workspace previously allocated with the LOGMPDFSET function. When called without arguments, the LOGMPDFFREE frees all the symbols previously allocated by LOGMPDFSETSQ or LOGMPDFSET. Each freed symbol is reported in the SAS log.
Figure 2 displays convergence diagnostic graphs to assess whether the Markov chain has converged. The trace plot indicates that the chain appears to have reached a stationary distribution. The chain also has good mixing and is dense.
The autocorrelation plot indicates low autocorrelation and efficient sampling. Finally, the kernel density plot shows the smooth, unimodal shape of posterior marginal distribution for . In a similar fashion, the rest of the diagnostic plots should be examined to ensure convergence of all parameters.
|Parameter||Alpha||Equal-Tail Interval||HPD Interval|
Figure 3 reports summary and interval statistics for each regression parameter and the variance parameter. Note that for the slope coefficient both equal tail interval and HPD interval contain 0.
Figure 4 displays convergence diagnostic graphs for and . The brain weight increases by a factor of (approximately 54%) for each logarithmic change of one day of gestational age. Similarly, the brain weight decreases by (11%) for each addition to the litter size.