calculate CIs for parameter estimates in proc nlimixed using the likelihood profile method

[rshaw1]

I am trying to calculate confidence intervals for a parameter estimate generated by proc nlmixed using the likelihood profile method.  The CIs that are currently calculated do not appear to use the likelihood profile method.  Is there a way to get proc nlmixed to do this?

Here is my code:

/*create data table that includes Q*/

PROC SQL;

   CREATE TABLE WORK.SEIDMANETAL1986_Q AS

   SELECT t1.TSFE_yrs_Range,

          t1.TSFE_yrs_Avg,

          t1.duration,

          t1.conc,

          t1.PY,

          t1.'Obs'n,

          /* Q */

            (case

            when t1.TSFE_yrs_Avg=. then .

            when t1.TSFE_yrs_Avg<10 then 0

            when t1.TSFE_yrs_Avg>(10+t1.duration) then (t1.TSFE_yrs_Avg-10)**3-(t1.TSFE_yrs_Avg-10-t1.duration)**3

            else (t1.TSFE_yrs_Avg-10)**3

            end) AS Q

      FROM '\\Esc-server1\home\rshaw\EH027\seidmanetal1986.sas7bdat' t1;

run;

/*model*/

proc nlmixed data=WORK.SEIDMANETAL1986_Q;

parms KM 1e-08; /*starting value*/

pred = conc*KM*Q*PY;

ll=LogPDF("POISSON",Obs,Pred); /*LofPDF function Returns the logarithm of a probability density (mass) function. Poisson distribution is specified.*/

model Obs ~ general(ll);

predict pred out=model1; /*generates table with predicted values and CI’s titled “Model1”*/

run;

Here is what the parameter estimates table that proc nlmixed generates looks like:

Parameter EstimatesParameter Estimate StandardError DF t Value Pr > |t| 95% Confidence Limits Gradient

KM

3.946E-8

0

7

Infty

<.0001

-Infty

Infty

-10.3456

[Rick_SAS]

Check the log and your data.

• It's hard to read your output, but is that StdErr = 0? If so, your data is constant.
• Your output shows that the gradient is not zero at the solution that you post. Therefore either the algorithm did not converge or the procedure stopped for some other reason (bad data?)

It would be helpful if you could post a PROC PRINT of the first 10 observations we can see what the data look like.

[rshaw1]

Thank you.  I believe it did converge - the following note was included in the output:

NOTE: GCONV convergence criterion satisfied.

 It's a very small data set:

Obs TSFE_yrs_Range TSFE_yrs_Avg duration conc PY Obs

1	5-10	7.5	1.5	46.9	3952	0
2	10-15	12.5	1.5	48.3	3628	0
3	15-19	17.5	1.5	44.1	3198	0
4	20-25	22.5	1.5	43.2	2656	2
5	25-30	27.5	1.5	40.3	2094	5
6	30-35	32.5	1.5	33.5	1576	8
7	35-40	37.5	1.5	31.1	1086	2

[Rick_SAS]

Is that the data being sent to PROC NLMIXED? It doesn't contain Q.

Anyway, you might try rescaling these data. Intuitively, if you define x = conc*Q*PY, then x is large ... roughly 1E8. In order to model Obs (which is in the range 0-8), by the mean of a Poisson(KM*x) distribution, the product must be about 2, which means that KM is approx 1E-8. This estimate is so small that it is  the same order as the convergence criteria for the optimization.

Try dividing x by 1E6, such as follows:

data have;
length TSFE_yrs_Range $5;
input TSFE_yrs_Range TSFE_yrs_Avg duration conc PY Obs ;
if      TSFE_yrs_Avg=. then Q = .;
else if TSFE_yrs_Avg<10 then Q = 0;
else if TSFE_yrs_Avg>(10+duration) then 
   Q = (TSFE_yrs_Avg-10)**3-(TSFE_yrs_Avg-10-duration)**3;
else Q =(TSFE_yrs_Avg-10)**3;
x = conc*Q*PY / 1e6;
   datalines;
 5-10  7.5 1.5 46.9 3952 0
10-15 12.5 1.5 48.3 3628 0
15-19 17.5 1.5 44.1 3198 0
20-25 22.5 1.5 43.2 2656 2
25-30 27.5 1.5 40.3 2094 5
30-35 32.5 1.5 33.5 1576 8
35-40 37.5 1.5 31.1 1086 2
;

proc means data=Have; var x; run;

proc nlmixed data=have;
parms KM 1e-2; /*starting value*/
pred = KM*x;
ll=LogPDF("POISSON",Obs,Pred); /*LofPDF function Returns the logarithm of a probability density (mass) function. Poisson distribution is specified.*/
model Obs ~ general(ll);
predict pred out=model1; /*generates table with predicted values and CI’s titled “Model1”*/
run;

[rshaw1]

Thanks Rick. 

Q is generated before proc nlmixed, see below code.

My goal is to generate confidence bounds for the KM parameter using the likelihood profile method.  How does what you propose doing accomplish that?

Thank you,

Rebekha

/*calculate Q*/

proc print data='\\Esc-server1\home\rshaw\EH027\seidmanetal1986.sas7bdat';

run;

PROC SQL;

CREATE TABLE WORK.SEIDMANETAL1986_Q AS

SELECT t1.TSFE_yrs_Range,

t1.TSFE_yrs_Avg,

t1.duration,

t1.conc,

t1.PY,

t1.'Obs'n,

/* Q */

(case

when t1.TSFE_yrs_Avg=. then .

when t1.TSFE_yrs_Avg<10 then 0

when t1.TSFE_yrs_Avg>(10+t1.duration) then (t1.TSFE_yrs_Avg-10)**3-(t1.TSFE_yrs_Avg-10-t1.duration)**3

else (t1.TSFE_yrs_Avg-10)**3

end) AS Q

FROM '\\Esc-server1\home\rshaw\EH027\seidmanetal1986.sas7bdat' t1;

run;

　

proc nlmixed data=WORK.SEIDMANETAL1986_Q;

parms KM 1e-8;

pred = conc*KM*Q*PY;

ll=LogPDF("POISSON",Obs,Pred);

model Obs ~ general(ll);

predict pred out=model1 alpha=0.1 /*df=5*/;

Proc print data=model1;

run;

[Rick_SAS]

> How does what you propose doing accomplish that?

If you scale the data, you can get estimates and CIs in the new scale.

I have to think about it a little more, but it looks like this problem is linear in the parameter, so you can then divide the estimate and CIs by 1E6 to obtain the estimates on the original scale.  In fact, now that I look at it, I don't think you need NLMIXED at all. I think you can use PROC GENMOD to perform the regression. It all looks linear in the parameter. Maybe someone else can verify.

[rshaw1]

I don't know much about Proc Genmod - how would I enter this into that procedure?  I'm looking to 1st estimate KM and calculate the expected cases for each group (C=KM*Q*PY) by maximizing the log likelihood.  Proc nlmixed is currently providing a value for KM and giving me the predicted values for each group.  I'll try scaling the data as you suggest, but I don't think proc nlmixed calculates CIs using the likelihood profile method.

Thank you.