Documentation Example 22 for PROC MCMC

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: MCMCEX22                                            */
/*   TITLE: Documentation Example 22 for PROC MCMC              */
/*          One-Compartment Model with PK Data                  */
/* PRODUCT: STAT                                                */
/*  SYSTEM: ALL                                                 */
/*    KEYS:                                                     */
/*   PROCS: MCMC                                                */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: Fang Chen                                           */
/*     REF: PROC MCMC, EXAMPLE 22                               */
/*    MISC:                                                     */
/****************************************************************/

data theoph;
   input subject time conc dose;
   datalines;
 1  0.00  0.74 4.02
 1  0.25  2.84 4.02
 1  0.57  6.57 4.02
 1  1.12 10.50 4.02
 1  2.02  9.66 4.02
 1  3.82  8.58 4.02
 1  5.10  8.36 4.02
 1  7.03  7.47 4.02
 1  9.05  6.89 4.02
 1 12.12  5.94 4.02
 1 24.37  3.28 4.02
 2  0.00  0.00 4.40
 2  0.27  1.72 4.40
 2  0.52  7.91 4.40
 2  1.00  8.31 4.40
 2  1.92  8.33 4.40
 2  3.50  6.85 4.40
 2  5.02  6.08 4.40
 2  7.03  5.40 4.40
 2  9.00  4.55 4.40
 2 12.00  3.01 4.40
 2 24.30  0.90 4.40
 3  0.00  0.00 4.53
 3  0.27  4.40 4.53
 3  0.58  6.90 4.53
 3  1.02  8.20 4.53
 3  2.02  7.80 4.53
 3  3.62  7.50 4.53
 3  5.08  6.20 4.53
 3  7.07  5.30 4.53
 3  9.00  4.90 4.53
 3 12.15  3.70 4.53
 3 24.17  1.05 4.53
 4  0.00  0.00 4.40
 4  0.35  1.89 4.40
 4  0.60  4.60 4.40
 4  1.07  8.60 4.40
 4  2.13  8.38 4.40
 4  3.50  7.54 4.40
 4  5.02  6.88 4.40
 4  7.02  5.78 4.40
 4  9.02  5.33 4.40
 4 11.98  4.19 4.40
 4 24.65  1.15 4.40
 5  0.00  0.00 5.86
 5  0.30  2.02 5.86
 5  0.52  5.63 5.86
 5  1.00 11.40 5.86
 5  2.02  9.33 5.86
 5  3.50  8.74 5.86
 5  5.02  7.56 5.86
 5  7.02  7.09 5.86
 5  9.10  5.90 5.86
 5 12.00  4.37 5.86
 5 24.35  1.57 5.86
 6  0.00  0.00 4.00
 6  0.27  1.29 4.00
 6  0.58  3.08 4.00
 6  1.15  6.44 4.00
 6  2.03  6.32 4.00
 6  3.57  5.53 4.00
 6  5.00  4.94 4.00
 6  7.00  4.02 4.00
 6  9.22  3.46 4.00
 6 12.10  2.78 4.00
 6 23.85  0.92 4.00
 7  0.00  0.15 4.95
 7  0.25  0.85 4.95
 7  0.50  2.35 4.95
 7  1.02  5.02 4.95
 7  2.02  6.58 4.95
 7  3.48  7.09 4.95
 7  5.00  6.66 4.95
 7  6.98  5.25 4.95
 7  9.00  4.39 4.95
 7 12.05  3.53 4.95
 7 24.22  1.15 4.95
 8  0.00  0.00 4.53
 8  0.25  3.05 4.53
 8  0.52  3.05 4.53
 8  0.98  7.31 4.53
 8  2.02  7.56 4.53
 8  3.53  6.59 4.53
 8  5.05  5.88 4.53
 8  7.15  4.73 4.53
 8  9.07  4.57 4.53
 8 12.10  3.00 4.53
 8 24.12  1.25 4.53
 9  0.00  0.00 3.10
 9  0.30  7.37 3.10
 9  0.63  9.03 3.10
 9  1.05  7.14 3.10
 9  2.02  6.33 3.10
 9  3.53  5.66 3.10
 9  5.02  5.67 3.10
 9  7.17  4.24 3.10
 9  8.80  4.11 3.10
 9 11.60  3.16 3.10
 9 24.43  1.12 3.10
10  0.00  0.24 5.50
10  0.37  2.89 5.50
10  0.77  5.22 5.50
10  1.02  6.41 5.50
10  2.05  7.83 5.50
10  3.55 10.21 5.50
10  5.05  9.18 5.50
10  7.08  8.02 5.50
10  9.38  7.14 5.50
10 12.10  5.68 5.50
10 23.70  2.42 5.50
11  0.00  0.00 4.92
11  0.25  4.86 4.92
11  0.50  7.24 4.92
11  0.98  8.00 4.92
11  1.98  6.81 4.92
11  3.60  5.87 4.92
11  5.02  5.22 4.92
11  7.03  4.45 4.92
11  9.03  3.62 4.92
11 12.12  2.69 4.92
11 24.08  0.86 4.92
12  0.00  0.00 5.30
12  0.25  1.25 5.30
12  0.50  3.96 5.30
12  1.00  7.82 5.30
12  2.00  9.72 5.30
12  3.52  9.75 5.30
12  5.07  8.57 5.30
12  7.07  6.59 5.30
12  9.03  6.11 5.30
12 12.05  4.57 5.30
12 24.15  1.17 5.30
;

proc fcmp outlib=sasuser.funcs.PK;
   subroutine OneComp(t,y[*],dy[*],ka,ke);
   outargs dy;
   dy[1] = -ka*y[1];
   dy[2] = ka*y[1]-ke*y[2];
   endsub;
run;


options cmplib=sasuser.funcs;

proc mcmc data=theoph nmc=10000 seed=27 outpost=theophO diag=none
   nthreads=8;
   ods select PostSumInt;
   array b[2];
   array muB[2] (0 0);
   array cov[2,2];
   array S[2,2] (1 0 0 1);

   array init[2] dose 0;
   array sol[2];

   parms beta1 -3.22 beta2 0.47 beta3 -2.45 ;
   parms cov {0.03 0 0 0.4};
   parms s2y;

   prior beta: ~ normal(0, sd=100);
   prior cov ~ iwish(2, S);
   prior s2y  ~ igamma(shape=3, scale=2);

   random b ~ mvn(muB, cov) subject=subject;
   cl   = exp(beta1 + b1);
   ka   = exp(beta2 + b2);
   ke   = exp(beta3);
   v    = cl/ke;
   call ode('OneComp',sol,init,0,time,ka,ke);
   mu   = (sol[2]/v);
   model conc ~ normal(mu,var=s2y);
run;

proc mcmc data=theoph nmc=10000 seed=22 outpost=theophC;
   array b[2];

   array mu[2] (0 0);
   array cov[2,2];
   array S[2,2] (1 0 0 1);

   parms beta1 -3.22 beta2 0.47 beta3 -2.45 ;
   parms cov {0.03 0 0 0.4};
   parms s2y;

   prior beta: ~ normal(0, sd=100);
   prior cov ~ iwish(2, S);
   prior s2y  ~ igamma(shape=3, scale=2);

   random b ~ mvn(mu, cov) subject=subject;
   cl   = exp(beta1 + b1);
   ka   = exp(beta2 + b2);
   ke   = exp(beta3);
   pred = dose*ke*ka*(exp(-ke*time)-exp(-ka*time))/cl/(ka-ke);
   model conc ~ normal(pred,var=s2y);
run;