/*------------------------------------ Problem 1. Cefamandole Data ------------------------------------*/ /* PROBLEM: Based on the following cefamandole data, the problem is to run a nonlinear mixed-effects */ /* model so as to obtain estimates of the population PK parameters assuming all random */ /* effects are uncorrelated. We also wish to compare the first-order exansion around 0 to */ /* the Laplace MLE (LMLE) estimate using METHOD=GAUSS and QPOINTS=1 options as well as to */ /* the "exact" MLE estimate using the default adaptive Gaussian quadrature method (i.e., */ /* numerical integration). We will also examine the impact of different starting values on */ /* these estimates so BE CAREFULL in your code (HINT, HINT,...) */ /*-----------------------------------------------------------------------------------------------------*/ /* Data set taken from Davidian and Giltinan, 1993 */ %include 'c:\edward\macro\nlinmix\covparms.sas' / nosource2;run; LIBNAME SC 'c:\edward\rm\short course examples\output datasets\cefamandole\'; data Cefamandole; input Subject time Conc; log_conc=log(Conc); hour=time/60; datalines; 1 10 127.00 1 15 87.00 1 20 47.40 1 30 39.90 1 45 24.80 1 60 17.90 1 75 11.70 1 90 10.90 1 120 5.70 1 150 2.55 1 180 1.84 1 240 1.50 1 300 0.70 1 360 0.34 2 10 120.00 2 15 90.10 2 20 70.00 2 30 40.10 2 45 24.00 2 60 16.10 2 75 11.60 2 90 9.20 2 120 5.20 2 150 3.00 2 180 1.54 2 240 0.73 2 300 0.37 2 360 0.19 3 10 154.00 3 15 94.00 3 20 84.00 3 30 56.00 3 45 37.10 3 60 28.90 3 75 25.20 3 90 20.00 3 120 12.40 3 150 8.30 3 180 4.50 3 240 3.40 3 300 1.70 3 360 1.19 4 10 181.00 4 15 119.00 4 20 84.30 4 30 56.10 4 45 39.80 4 60 23.30 4 75 22.70 4 90 13.00 4 120 8.00 4 150 2.40 4 180 1.60 4 240 1.10 4 300 0.48 4 360 0.29 5 10 253.00 5 15 176.00 5 20 150.00 5 30 90.30 5 45 69.60 5 60 42.50 5 75 30.60 5 90 19.60 5 120 13.80 5 150 11.40 5 180 6.30 5 240 3.80 5 300 1.55 5 360 1.22 6 10 140.00 6 15 120.00 6 20 106.00 6 30 60.40 6 45 60.90 6 60 42.20 6 75 26.80 6 90 22.00 6 120 14.50 6 150 8.80 6 180 6.00 6 240 3.00 6 300 1.30 6 360 1.03 ; /* Two-stage OLS regression approach */ proc sort data= cefamandole ; by Subject time; run; goptions reset=all hsize = 6.25 in vsize = 5.25 in htext = 3.0 pct htitle = 3.5 pct vorigin = 0 in horigin = 0 in ftext = swiss lfactor = 1; axis1 order = (0 to 300 by 30) value = (font=swissl height=1) label = (font=swissl h=1) major = (h=0.5) minor = none; axis2 order = (0 to 6 by .5) value = (font=swissl height=1) label = (font=swissl h=1) major = (h=1) minor = none; axis3 order = (0 to 300 by 30) value = none label = none major = (h=0.5) minor = none; symbol1 color=black value=none i=join l=2 r=6 w=1; symbol2 color=black value=none i=join l=1 w=5; proc gplot data=cefamandole ; plot conc*hour=Subject / nolegend vaxis=axis1 haxis=axis2; run; quit; /* STEP 1: Compare COVPARMS starting values for FIRO and LMLE to STS starting values for FIRO and LMLE */ %macro STEP1; ods output parameterestimates=sc.OLSparms_new; proc nlmixed data=cefamandole qpoints=1; parms beta1=5 beta2=1.5 beta3=4 beta4=-0.12 s2=1 ; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); resid=conc-func; var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[ 0 , 0 , 0, 0 , 0, 0, 0 , 0, 0, 0 ]) subject=Subject; predict func out=predout der; id resid; run; %covparms(parms=OLSparms_new, predout=predout, resid=resid, method=ml, random=der_b1i der_b2i der_b3i der_b4i, subject=subject, type=vc, covname=psi, output=MLEparms); proc nlmixed data=cefamandole method=firo; parms /data=MLEparms; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; T_half = log(2)/exp(beta4); func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); func_avg = exp(beta1)*exp(-exp(beta2)*hour) + exp(beta3)*exp(-exp(beta4)*hour); var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; estimate 't(1/2) - terminal' T_half; predict func out=predplot;id func func_avg; title1 'Model 1-b: First-order PA approximation (ELS) with power function of mean'; title2 'Starting values based on COVPARMS starting values'; run; proc nlmixed data=cefamandole qpoints=1; parms /data=MLEparms; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; T_half = log(2)/exp(beta4); func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); func_avg = exp(beta1)*exp(-exp(beta2)*hour) + exp(beta3)*exp(-exp(beta4)*hour); var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; estimate 't(1/2) - terminal' T_half; predict func out=predplot;id func func_avg; title1 'Model 2-b: LMLE approximation (Laplace) with power function of mean'; title2 'Starting values based on COVPARMS starting values'; run; %mend step1; /* STEP 2: Compare eyeball guesstimate starting values for FIRO and LMLE to STS starting values for FIRO and LMLE */ %macro STEP2; ods output parameterestimates=sc.FIROparms1; proc nlmixed data=cefamandole method=firo; parms beta1=5 beta2=1.5 beta3=4 beta4=-0.12 s2=1 ; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; T_half = log(2)/exp(beta4); func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); func_avg = exp(beta1)*exp(-exp(beta2)*hour) + exp(beta3)*exp(-exp(beta4)*hour); var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; estimate 't(1/2) - terminal' T_half; predict func out=predplot;id func func_avg; title1 'Model 1: First-order PA approximation (ELS) with power function of mean'; title2 'Starting values based on eyeball guesstimates'; run; ods output parameterestimates=sc.LMLEparms1; proc nlmixed data=cefamandole qpoints=1; parms beta1=5 beta2=1.5 beta3=4 beta4=-0.12 s2=1 ; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; T_half = log(2)/exp(beta4); func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); func_avg = exp(beta1)*exp(-exp(beta2)*hour) + exp(beta3)*exp(-exp(beta4)*hour); var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; estimate 't(1/2) - terminal' T_half; predict func out=predplot;id func func_avg; title1 'Model 2: LMLE approximation (Laplace) with power function of mean'; title2 'Starting values based on eyeball guesstimates'; run; %mend STEP2; /* STEP 3: Compare FIRO starting values for LMLE */ %macro step3; proc print data=sc.FIROparms1; run; proc nlmixed data=cefamandole qpoints=1; parms /data=sc.FIROparms1; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; T_half = log(2)/exp(beta4); func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); func_avg = exp(beta1)*exp(-exp(beta2)*hour) + exp(beta3)*exp(-exp(beta4)*hour); var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; estimate 't(1/2) - terminal' T_half; predict func out=predplot;id func func_avg; title1 'Model 2: LMLE approximation (Laplace) with power function of mean'; title2 'Starting values based on FIRO starting values'; run; %mend step3; /*------------------------------------------------------------------------------*/ /* STEP 4: Compare exact MLE using adaptive and qpoints=10 Guassisan quadrature */ /* In this step DO NOT run the versions that take forever to run. */ /*------------------------------------------------------------------------------*/ %macro step4; /* The following took a long time to run */ /* ods output parameterestimates=sc.FullMLEparms1; proc nlmixed data=cefamandole method=gauss qpoints=10; parms /data=sc.LMLEparms1; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; T_half = log(2)/exp(beta4); func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); func_avg = exp(beta1)*exp(-exp(beta2)*hour) + exp(beta3)*exp(-exp(beta4)*hour); var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; estimate 't(1/2) - terminal' T_half; predict func out=predplot;id func func_avg; title1 'Model 3: Full MLE with 10 point Gaussian quadrature (qpoints=10) and power function of mean'; title2 'Starting values based on LMLE estimates'; run; */ ods output parameterestimates=sc.FullMLEparms2; proc nlmixed data=cefamandole method=gauss; parms /data=sc.LMLEparms1; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; T_half = log(2)/exp(beta4); func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); func_avg = exp(beta1)*exp(-exp(beta2)*hour) + exp(beta3)*exp(-exp(beta4)*hour); var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; estimate 't(1/2) - terminal' T_half; predict func out=predplot;id func func_avg; title1 'Model 4: Full MLE with adaptive Gaussian quadrature (qpoints=adaptive) and power function of mean'; title2 'Starting values based on LMLE estimates'; run; /* The following took a long time to run */ /* ods output parameterestimates=sc.FullMLEparms3; proc nlmixed data=cefamandole method=gauss; parms beta1=5 beta2=1.5 beta3=4 beta4=-0.12 s2=1 ; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; T_half = log(2)/exp(beta4); func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); func_avg = exp(beta1)*exp(-exp(beta2)*hour) + exp(beta3)*exp(-exp(beta4)*hour); var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; estimate 't(1/2) - terminal' T_half; predict func out=predplot;id func func_avg; title1 'Model 5: Full MLE with adaptive Gaussian quadrature (qpoints=adaptive) and power function of mean'; title2 'Starting values based on eyeball guesstimates'; run; */ %mend step4; %macro MISC; /* Alternative determiniation of starting values */ ods output parameterestimates=OLSparms_new; proc nlmixed data=cefamandole qpoints=1; parms beta1=5 beta2=1.5 beta3=4 beta4=-0.12 s2=1 ; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); resid=conc-func; var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[ 0 , 0 , 0, 0 , 0, 0, 0 , 0, 0, 0 ]) subject=Subject; predict func out=predout der; id resid; run; %covparms(parms=OLSparms_new, predout=predout, resid=resid, method=ml, random=der_b1i der_b2i der_b3i der_b4i, subject=subject, type=vc, covname=psi, output=MLEparms); proc nlmixed data=cefamandole method=firo; parms /data=MLEparms; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; T_half = log(2)/exp(beta4); func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); func_avg = exp(beta1)*exp(-exp(beta2)*hour) + exp(beta3)*exp(-exp(beta4)*hour); var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; estimate 't(1/2) - terminal' T_half; predict func out=predplot;id func func_avg; title1 'Model 2: First-order PA approximation (ELS) with power function of mean'; title2; run; proc nlmixed data=cefamandole qpoints=1; parms /data=MLEparms; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; T_half = log(2)/exp(beta4); func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); func_avg = exp(beta1)*exp(-exp(beta2)*hour) + exp(beta3)*exp(-exp(beta4)*hour); var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; estimate 't(1/2) - terminal' T_half; predict func out=predplot;id func func_avg; title1 'Model 2: Second-order SS approximation (Laplace) with power function of mean'; title2; run; data MLEparms1; set MLEparms; if Parameter in ('psib4') then delete; run; proc nlmixed data=cefamandole method=firo; parms /data=MLEparms1; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4; func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); resid=conc-func; var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i ~ normal([0,0,0],[psib1, 0 , psib2, 0 , 0, psib3 ]) subject=Subject; predict func out=predout_new der; id resid; run; proc nlmixed data=cefamandole qpoints=1; parms /data=MLEparms1; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4; func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); resid=conc-func; var_conc = s2*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i ~ normal([0,0,0],[psib1, 0 , psib2, 0 , 0, psib3 ]) subject=Subject; predict func out=predout_new der; id resid; run; /* Model 3: First-order Approximation Method - Power of the mean variance model */ data MLEparms2; set MLEparms; if Parameter in ('theta') then delete; run; proc nlmixed data=cefamandole method=firo;*qpoints=1; parms /data=MLEparms2; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); var_conc = s2*(func**2); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; run; /* Model 4 Laplace Approximation Method - Constant variance model */ proc nlmixed data=cefamandole method=firo qpoints=1; parms /data=MLEparms2; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); var_conc = s2; model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; run; /* Model 5: Laplace Approximation Method - Power of the mean variance model */ proc nlmixed data=cefamandole qpoints=1; parms /data=MLEparms2; beta1i=beta1 + b1i; beta2i=beta2 + b2i; beta3i=beta3 + b3i; beta4i=beta4 + b4i; func = exp(beta1i)*exp(-exp(beta2i)*hour) + exp(beta3i)*exp(-exp(beta4i)*hour); var_conc = sigsq*(func**(2*theta)); model conc ~ normal(func, var_conc); random b1i b2i b3i b4i ~ normal([0,0,0,0],[psib1, 0 , psib2, 0 , 0 , psib3, 0 , 0 , 0 , psib4]) subject=Subject; run; %mend MISC;