add co-variates in a proc NLMIXED to estimate a rate difference

[KarinGun]

Hi all!

I am modelling the correlation between a rheumatic disease (systemic sclerosis = exposure) and cancer (outcome) using a poisson regression model in proc genmod. My co-variates are age and sex. Individuals with systemic sclerosis (1) are compared to comparators (0) without the disease. The estimate I get from this model is incidence rate ratio. Sex is a binary variable and also ssc. Age is a dicrete variable. Can = cancer. 

My main model using Proc GENMOD is: 

proc genmod data=incidence5;  

class sex(ref="2") ssc(ref = "0");

model can = ssc age sex/ dist=poisson link=log offset=logpy type3;

ods output parameterestimates=results;

 run;

(Then I exponentiated my parameterestimates to get estimates that are easier to interpret). 

Then I would like to estimate the rate difference (and adjusted incidence rate) but this can't be extracted from this proc genmod estimation as I understood it. Therefore I plan to use Proc NLMIXED. I have multiplicated the estimates with 1000 to get /1000 person years as an outcome estimate.The model I have used is: 

 proc nlmixed data=incidence5;

         lambda = exp(b0 + b1*(ssc=1) + logpy);

         model can ~ poisson(lambda);

         estimate "Rate Difference" (exp(b0+b1)-exp(b0))*1000 df=1e8; 

estimate "incidence rate comparators" exp(b0)*1000 df=1e8;

estimate "incidence rate SSc" exp(b0+b1)*1000 df=1e8;

run;

It worked fine to run this model and I got the expected estimates corresponding to crude estimates calculated outside the model. However I would like to incorporate my co-variates sex and age in this model to receive adjusted estimates of rate difference and incidence rate. I have tried to read about how to do this in the proc nlmixed documentation but didn't find a solution that worked out. Unfortunately I don't have much experience in running models in SAS yet. Do you have any suggestions on how to incorporate my co-variates in the proc NLMIXED model (if possible)? Or to get these estimates in another way/another proc? 

This is how the procedure of extracting rate difference from proc nlmixed is describes on the SAS support page:

Rate difference using the PROC NLMIXED

The rate difference can also be estimated by fitting the Poisson model using PROC NLMIXED as follows. The LAMBDA= assignment statement expresses the Poisson mean parameter, lambda, as a function of age, the offset (ln), and the model parameters (b0 and b1). The MODEL statement indicates that the response count, c, is to be modeled as a Poisson variable with mean lambda. The ESTIMATE statement is used to estimate the rate difference. A large degrees of freedom value (df=1e8) is used to produce large-sample (Wald) statistics like those from the NLMeans and NLEstimate macros above. Note that the estimated log rate for age 1 is b0+b1 and the estimated log rate for age 2 is just b0. For details about using PROC NLMIXED, see the NLMIXED documentation.

      proc nlmixed data=insure;
         lambda = exp(b0 + b1*(age=1) + ln);
         model c ~ poisson(lambda);
         estimate "Rate Difference" exp(b0+b1)-exp(b0) df=1e8;
         run;

Thank you on beforehand and have a good day!

Best, 

Karin

[StatDave · Accepted solution]

For those interested, the page that is referred to and shows the NLMIXED code is this SAS note. You can add your sex and age variables as follows. Since you treated age as continuous in your GENMOD code, the same is done in NLMIXED. I assume that the non-reference value of sex is 1 and that sex is a numeric variable.

lambda = exp(b0 + b1*(ssc=1) + b2*(sex=1) + b3*age + logpy);

If the ESTIMATE statement is unchanged from what is shown in the SAS note, it will estimate the rate difference at the reference level of sex (sex=2) and age=0 since that setting adds no contributions involving sex or age.

[KarinGun]

Hi Dave!

Thanks for your suggestion! 

Sex is coded 1 or 2. Ssc is coded 0 or 1. Age is a numeric value in integers. 

I tried to run it but I got the following message in the log (and the process was not run): 

218 proc nlmixed data=incidence5;
219 lambda = exp(b0 + b1*(ssc=1) + b2*(sex =1) + b3*age + logpy); /*lambda = poisson mean parameter*/
NOTE: Character value converted to numeric for argument 1 of '=' operation.
220 model can ~ poisson(lambda);
221 estimate "Rate Difference" (exp(b0+b1)-exp(b0))*1000 df=1e8; /* *1000 to get the estimates as /1000 person years*/
222 estimate "incidence rate comparators" exp(b0)*1000 df=1e8;
223 estimate "incidence rate SSc" exp(b0+b1)*1000 df=1e8;
224
225 run;

NOTE: To assign starting values to parameters, use the PARMS statement. The default starting value of 1.0 is in effect for all parameters.
ERROR: No valid parameter points were found.
NOTE: The SAS System stopped processing this step because of errors.
NOTE: PROCEDURE NLMIXED used (Total process time):
real time 0.04 seconds
cpu time 0.04 seconds

So it sounds like I need to add a PARMS statement. But how do I choose adequate initial values of the estimates? 

Thanks again!

Best, 

Karin

[Ksharp]

https://support.sas.com/kb/37/228.html

[StatDave]

I think the first Note indicates the problem. It suggests that one of your variables is character, not numeric. If, for example, your sex variable has character values '1' and '2', then in the LAMBDA= statement you need to use ... b2(sex='1') ...

[KarinGun]

Thanks all!

Sorry I missed the first note about character/numeric variables. You were right sex is a character value and I adjusted the code accordingly. 

Now the code looks like this: 

- After SSc diangosis;

  proc nlmixed data=incidence5;

         lambda = exp(b0 + b1*(ssc=1) + b2*(sex ='1') + b3*age + logpy); /*lambda = poisson mean parameter*/

         model can ~ poisson(lambda);

PARMS b0 -7 b1 .1 b2 .1 b3 .1;

         estimate "Rate Difference" (exp(b0+b1)-exp(b0))*1000 df=1e8; /* *1000 to get the estimates as /1000 person years*/

estimate "incidence rate comparators" exp(b0)*1000 df=1e8;

estimate "incidence rate SSc" exp(b0+b1)*1000 df=1e8;

         run;

And now it works to run 🙂 Not really the expected outcome but not identical to the unadjusted estimates.  My last question: do you think that I seem to have specified the initial values in the PARMS in a correct way?  I read the documentation about this but still not confident due to lack of experience. 

Thanks again. 

/Karin

[StatDave]

That's fine, though documentation suggests it would be written: parms b0=-7, b1=.1, b2=.1, b3=.1;
But often it isn't necessary to even include the PARMS statement, so you might want to start without one and let NLMIXED use default values (which are 1 for all parameters). You can add the PARMS statement if you encounter fitting problems.

[KarinGun]

Thank you Dave, 

If I don't use the PARMS statement I don't get any output estimates. If I use the PARMS statement as suggested I get the estimates but a bit off the expected range when compared to the unadjusted estimates. 

[SteveDenham]

Time for a sensitivity analysis of the NLMIXED results. Set up a fairly narrow gridsearch in the PARMS statement so that you can get an idea of how sensitive the parameter values are to the starting values. This may take care of your differences from the expected values.

SteveDenham

[KarinGun]

Thank you Steve. I will try that. 

Best, 

Karin