Profile-Likelihood-Based Confidence Intervals

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: IMNLPEX5                                            */
/*   TITLE: Profile-Likelihood-Based Confidence Intervals       */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS:                                                     */
/*   PROCS: IML                                                 */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: saswmh/sasghg               UPDATE: frwick 5/2015   */
/*     REF: Technical Report Examples                           */
/*    MISC:                                                     */
/*                                                              */
/****************************************************************/

proc iml;
/*------------------------------------------------*/
/*  Profile-Likelihood Based Confidence Intervals */
/*  For Weibull (Three Parameter) Log-likelihood  */
/*------------------------------------------------*/

carcin = { 143  164  188  188  190  192  206
         209  213  216  220  227  230  234
         246  265  304  216  244 };

p = ncol(carcin); m = p - 2;
minx = carcin[><];
rang = carcin[<>] - minx;
* print m minx rang thet;

/*--- 1. Two parameter estimation: for thet = 0 ---*/

start f_weib2(x) global(carcin,thet);
   /* use x[1]=sig and x[2]=c */
   p = ncol(carcin); m = p - 2;
   temp = carcin - thet;
   sum1 = sum(log(temp[1:m]));
   sum2 = sum((temp / x[1])##x[2]);
   f = m*log(x[2]) - m*x[2]*log(x[1]) + (x[2]-1)*sum1 - sum2;
   return(f);
finish f_weib2;

start g_weib2(x) global(carcin,thet);
   /* use x[1]=sig and x[2]=c */
   p = ncol(carcin); m = p - 2;
   g = j(1,2,0.);
   temp = carcin - thet;
   sum1 = sum(log(temp[1:m]));
   sum2 = sum((temp / x[1])##x[2]);
   sum3 = sum(((temp / x[1])##x[2]) # (log(temp / x[1])));
   g[1] = -m * x[2] / x[1] + sum2 * x[2] / x[1];
   g[2] = m / x[2] - m * log(x[1]) + sum1 - sum3;
   return(g);
finish g_weib2;

n = 2; thet = 0.;
x0 = j(1,n,.5);
optn = {1 2};
con = { 1.e-6 1.e-6 ,
     .     .    };
CALL NLPTR(rc,xopt,"f_weib2",x0,optn,con,,,,"g_weib2");
fopt = f_weib2(xopt);


/*--- Compute confidence interval for x at optimum ---*/
/* compute Hessian at optimum */
xopt = xopt`;
call nlpfdd(f,g,hes2,"f_weib2",xopt,,"g_weib2");
hin2 = inv(hes2);

/* quantile of normal distribution */
prob = 0.05;
stderr = sqrt(abs(vecdiag(hin2)));
z = quantile("Normal", 1 - prob/2);
xlb = xopt - z * stderr;
xub = xopt + z * stderr;
results = xlb || xopt || xub;
*print results[L="Normal Distribution Confidence Interval"
              c={"Lower Bound" "Estimate" "UpperBound"}
              r={"sigma" "c"}];

start f_plwei2(x) global(carcin,ipar,lstar);
   /* use x[1]=sig, x[2]=c, thet */
   like = f_weib2(x);
   grad = g_weib2(x);
   grad[ipar] = like - lstar;
   return(grad`);
finish f_plwei2;

/* quantile of chi**2 distribution */
chqua = cinv(1-prob,1); lstar = fopt - .5 * chqua;
/* set number of equations = 2, and print parameter = 0 */
optn = {2 0};                    /* Dont print in NLPLM */
do ipar = 1 to 2;
   /* Compute initial step: */
   /* Choose (alfa,delt) to go in right direction */
   /* Venzon & Moolgavkar (1988), p.89 */
   if ipar=1 then ind = 2; else ind = 1;
   delt = - inv(hes2[ind,ind]) * hes2[ind,ipar];
   alfa = - (hes2[ipar,ipar] - delt` * hes2[ind,ipar]);
   if alfa > 0 then do;
   alfa = .5 * sqrt(chqua / alfa);
   end; else do;
   print "Bad alpha";
   alfa = .1 * xopt[ipar];
   end;
   if ipar=1 then delt = 1 // delt;
           else delt = delt // 1;

   /* Get upper end of interval */
   x0 = xopt + alfa * delt;
   /* set lower bound to optimal value */
   con2 = con; con2[1,ipar] = xopt[ipar];
   CALL NLPLM(rc,xres,"f_plwei2",x0,optn,con2);
   f = f_plwei2(xres); s = ssq(f);
   if (s < 1.e-6) then xub[ipar] = xres[ipar];
                else xub[ipar] = .;

   /* Get lower end of interval */
   x0 = xopt - alfa * delt;
   /* reset lower bound and set upper bound to optimal value */
   con2[1,ipar] = con[1,ipar]; con2[2,ipar] = xopt[ipar];
   CALL NLPLM(rc,xres,"f_plwei2",x0,optn,con2);
   f = f_plwei2(xres); s = ssq(f);
   if (s < 1.e-6) then xlb[ipar] = xres[ipar];
                else xlb[ipar] = .;
end;
results = xlb || xopt || xub;
print results[L="Profile-Likelihood Confidence Interval"
              c={"Lower Bound" "Estimate" "UpperBound"}
              r={"sigma" "c"}];