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"}];