The MCMC Procedure

Example 54.9 Multivariate Normal Random-Effects Model

Gelfand et al. (1990) use a multivariate normal hierarchical model to estimate growth regression coefficients for the growth of 30 young rats in a control group over a period of 5 weeks. The following statements create a SAS data set with measurements of Weight, Age (in days), and Subject.

title 'Multivariate Normal Random-Effects Model';
data rats;
   array days[5] (8 15 22 29 36);
   input weight @@;
   subject = ceil(_n_/5);
   index = mod(_n_-1, 5) + 1;
   age = days[index];
   drop index days:;
datalines;     
151 199 246 283 320 145 199 249 293 354 
147 214 263 312 328 155 200 237 272 297 
135 188 230 280 323 159 210 252 298 331 
141 189 231 275 305 159 201 248 297 338 
177 236 285 350 376 134 182 220 260 296 
160 208 261 313 352 143 188 220 273 314 
154 200 244 289 325 171 221 270 326 358 
163 216 242 281 312 160 207 248 288 324 
142 187 234 280 316 156 203 243 283 317 
157 212 259 307 336 152 203 246 286 321 
154 205 253 298 334 139 190 225 267 302 
146 191 229 272 302 157 211 250 285 323 
132 185 237 286 331 160 207 257 303 345 
169 216 261 295 333 157 205 248 289 316 
137 180 219 258 291 153 200 244 286 324 
;

The model assumes normal measurement errors,

where indexes rat (Subject variable), indexes the time period, Weight and Age denote the weight and age of the th rat in week , and is the variance in the normal likelihood. The individual intercept and slope coefficients are modeled as the following:

You can use the following independent prior distributions on , , and :

The following statements fit this multivariate normal random-effects model:

proc mcmc data=rats nmc=10000 outpost=postout
   seed=17 init=random;
   ods select Parameters REParameters PostSummaries;
   array theta[2] alpha beta;
   array theta_c[2];
   array Sig_c[2,2];
   array mu0[2] (0 0);
   array Sig0[2,2] (1000 0 0 1000);
   array S[2,2] (0.02 0 0 20);

   parms theta_c Sig_c {121 0 0 0.26} var_y;
   prior theta_c ~ mvn(mu0, Sig0);
   prior Sig_c ~ iwish(2, S);
   prior var_y ~ igamma(0.01, scale=0.01);

   random theta ~ mvn(theta_c, Sig_c) subject=subject 
      monitor=(alpha_9 beta_9 alpha_25 beta_25);
   mu = alpha + beta * age;
   model weight ~ normal(mu, var=var_y);
run;

The ODS SELECT statement displays information about model parameters, random-effects parameters, and the posterior summary statistics. The ARRAY statements allocate memory space for the multivariate parameters and hyperparameters in the model. The parameters are (theta where the variable name of each element is alpha or beta), (theta_c), and (Sig_c). The hyperparameters are (mu0), (Sig0), and (S). The multivariate hyperparameters are assigned with constant values using parentheses .

The PARMS statement declares model parameters and assigns initial values to Sig_c using braces . The PRIOR statements specify the prior distributions. The RANDOM statement defines an array random effect theta and specifies a multivariate normal prior distribution. The SUBJECT= option indicates cluster membership for each of the random-effects parameter. The MONITOR= option monitors the individual intercept and slope coefficients of subjects 9 and 25.

You can use the following syntax in the RANDOM statement to monitor all parameters in an array random effect:

monitor=(theta)

This would produce posterior summary statistics on and .

The following syntax monitors all parameters:

monitor=(alpha)

If you did not name elements of theta to be alpha and beta, the SAS System creates variable names automatically in a consecutive fashion, as in theta1 and theta2.

A seemingly confusing message might occur if a symbol name matches an internally generated variable name for elements of an array. For example, if, instead of using the symbol var_y in the SAS program for the model variance , you used s2, the SAS System produces the following error message:

ERROR: The initial value 0 for the parameter S2 is outside of the prior 
       distribution support set.

This is confusing because the program does not assign an initial value for the parameter s2 in the PARMS statement, and you might expect that PROC MCMC would not generate an invalid initial value. The confusion is caused by the ARRAY statement that defines the array variable S:

   array S[2,2] (0.02 0 0 20);

Elements of S are automatically given names s1–s4. PROC MCMC interprets s2 as an element in S that was given a value of 0, hence producing this error message.