Resources

Documentation Example 7 for PROC NLMIXED

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: nlmex7                                              */
/*   TITLE: Documentation Example 7 for PROC NLMIXED            */
/*          Overdispersion hierarchical nonlinear mixed model   */
/* PRODUCT: STAT                                                */
/*  SYSTEM: ALL                                                 */
/*    KEYS:                                                     */
/*   PROCS: NLMIXED                                             */
/*    DATA: Aklilu et al(2013)                                  */
/*                                                              */
/* SUPPORT: Raghavendra Rao Kurada (Raghu)                      */
/*     REF:                                                     */
/*    MISC:                                                     */
/****************************************************************/

data comet;
   set sashelp.comet;
   L = 0; M = 0; H = 0; PC = 0;
   if (dose = 1.25) then L  = 1;
   if (dose = 2.50) then M  = 1;
   if (dose = 5.00) then H  = 1;
   if (dose = 200)  then PC = 1;
run;

proc nlmixed data = comet;
   parms b0 = -25  b1 = -10 b2 = -10 b3 = -10 b4 = -10 rho = 10;
   bounds sd1 >=0, sd2 >=0;
   mu = b0+ b1*L+ b2*M+ b3*H+ b4*PC + rateff + sampeff;
   term = (length**rho)*exp(mu);
   llike = log(rho) + (rho-1)*log(length) + mu - term;
   model length ~ general(llike);
   random ratEff ~ normal(0,sd1) subject = rat;
   random sampEff ~ normal(0,sd2) subject = sample(rat);
   predict gamma(1+(1/rho))/(exp(mu)**(1/rho)) out = p1;
run;

ods select ParameterEstimates;
proc nlmixed data = comet;
   parms b0 = -25  b1 = -10 b2 = -10 b3 = -10 b4 = -10 alpha = 1 rho = 10;
   bounds sd1 >=0, sd2 >=0;
   mu = b0+ b1*L+ b2*M+ b3*H+ b4*PC + rateff + sampeff;
   num = rho*((length)**(rho-1))*exp(mu);
   den = (1+(((length)**rho)*exp(mu))/alpha)**(alpha+1);
   llike = log(num/den);
   model length ~ general(llike);
   random rateff ~ normal(0,sd1) subject = rat;
   random sampeff ~ normal(0,sd2) subject = sample(rat);
   predict (gamma(alpha-(1/rho))*gamma(1+(1/rho)))
           /(gamma(alpha)*((exp(mu)/alpha)**(1/rho))) out = p2;
run;

proc means data = p1 mean noprint;
   class sample;
   var pred;
   output out = p11 (keep = sample weibull_mean)mean = weibull_mean;
run;

proc means data = p2 mean noprint;
   class sample;
   var pred;
   output out = p21 (keep = sample weibull_gamma_mean)mean = weibull_gamma_mean;
run;

proc means data = comet mean noprint;
   class sample;
   var length;
   id dose rat;
   output out = p31(keep = dose rat sample observed_mean) mean = observed_mean;
run;

data average;
   merge p11 p21 p31;
   by sample;
   if sample = . then delete;
   label observed_mean = "Observed Average";
   label weibull_gamma_mean = "Weibull Gamma Average";
   label weibull_mean = "Weibull Average";
   label sample = "Sample Index";
   if dose = 0 then dosage = "0 Vehicle Control ";
   if dose = 1.25 then dosage = "1 Low";
   if dose = 2.50 then dosage = "2 Medium";
   if dose = 5.00 then dosage = "3 High";
   if dose = 200  then dosage = "4 Positive Control";
run;

proc sgpanel data = average;
   panelby dosage/onepanel layout = columnlattice novarname uniscale = row;
   rowaxis values=(10 to 80 by 5) label="Average Tail Length";
   series x = sample y = observed_mean;
   series x = sample y = weibull_mean;
   series x = sample y = weibull_gamma_mean;
run;