Pregibon link test for generalized linear models

[susannemd]

Dear all,

I am conducting an analysis of cost data using proc genmod. Before I reduce the model, I want to make sure that I specify the correct distribution and link-function. I have assumed the response distribution to be gamma (after modified park method and graphical assessment). However, I am having trouble assessing the link-function using the Pregibon link test.

Do any of you have any experience with the Pregibon link test for GLMs analysed with proc genmod? I would be very grateful for any kind of help!

Kind regards,
Susanne

[SteveDenham]

I guess Google is my friend. A search on "Pregibon link test" leads to this site https://statcompute.wordpress.com/2016/12/04/pregibon-test-for-goodness-of-link-in-sas/ , where there is macro code for implementing the test for PROC LOGISTIC.  A global replace on that to GENMOD and correctly specifying your model should make it possible to test which link is optimum. What link functions are you considering? The log is the canonical link. but I see some sites that include identity and inverse links as possible link functions.

SteveDenham

[susannemd]

Dear Steve,

Thank you for your reply.

I have tested the macro-language, unfortunately, without succes:

%macro pregibon;

options mprint mlogic nocenter;

%let links = log identity;

%let loop = 1;

%let data = test;

%let y = total;

%let x = sex age location;

proc sql noprint;

           select n(&data) - 3 into :df from &data;

quit;

%do %while (%scan(&links, &loop) ne %str());

           %let link = %scan(&links, &loop);

proc genmod data= &data ;

           model &y= &x / dist=gamma link = &link;

           output out = _out1 p=p_1;

run;

data _out2;

           set _out1(rename = (p_1 = p1));

           p2 = p1*p1;

run;

ods listing close;

ods output ParameterEstimates = _parm;

proc genmod data= _out2;

           model &y = p1 p2 / dist=gamma link = &link;

run;

ods listing;

%if &loop = 1 then %do;

data _parm1;

           format link $10.;

           set _parm(where = (variable = "p2"));

           link= upcase("&link");

run;

%end;

%else %do;

data _parm1;

set _parm1 _parm(where = (variable = "p2") in = new);

if new then link = upcase("&link");

run;

%end;

data _parm2(drop=variable);

set _parm1;

_t= estimate / stderr;

_df= &df;

_p = (1 - probt(abs(_t), _df)) * 2;

run;

%let loop = %eval(&loop + 1);

%end;

title;

proc report data= _last_ spacing = 1 headline nowindows split = "*";

           column("   * PREGIBON TEST FOR GOODNESS OF LIN

                                 * H0: THE LINK FUNCTION IS SPECIFIED CORRECTLY * "

                                 link _t _df _p);

           define link / "Link function" width = 15 order order =data;

           define _t / "T-value" width = 15 format=12.4;

           define _df / "DF" width = 10;

           define _p/ "P-value" width = 15 format=12.4;

run;

%mend;

%filasign

Kind regards,
Susanne

[SteveDenham · Accepted solution]

HI @susannemd 

When you say "without success," everyone wonders what might have happened.  Please rerun the macro, with the options line:

options mlogic mprint symbolgen;

and provide the log.  Maybe we can get this to work (I hope).

Oh, I just saw something. Your x variable includes categorical variables sex and location.  You should include a CLASS statement in the PROC GENMOD call declaring them as class variables.  If that solves the problem, then the above stuff probably doesn't apply.  If it doesn't solve it, then we need to look at the log to find out what might be happening.

[susannemd]

Dear Steve,

Thank you! That was what was missing.

The code is now up and running.

[susannemd]

Hi again,

I managed to make the wrong post the "Solution", and I don't know how to fix it. I am sorry!

[SteveDenham]

No problem Susan. I have plenty of solutions marked to my credit. You deserve this one because of how much I learned that really applies to a lot of what I do right now.

SteveDenham

[FreelanceReinh]

Hi @susannemd,

---

@susannemd wrote:
I managed to make the wrong post the "Solution", and I don't know how to fix it.

---

Here's how:

Select a different post as the solution after clicking "Not the Solution" in the option menu (see icon below) of the current solution.

[StatDave]

It should be noted that there is a test for the link function, as well as for the specified form of any predictor, in the ASSESS statement in PROC GENMOD. You can find examples of its use in the GENMOD documentation.

[SteveDenham]

HI @StatDave ,

The ASSESS option looks quite good for model determination by looking at residuals.  I am curious as to what output looks like for using the link= option.  I also get the feeling that it is going to take a fair amount of experience with this to interpret the graphical results, while the Pregibon test essentially just looks to see if the candidate link functions are linear.  Is that a fair evaluation of the two methods?

SteveDenham

[StatDave]

It isn't really necessary to interpret the graph. By adding the RESAMPLE option, a p-value for the test of link adequacy is provided.

assess link / resample;

[SteveDenham]

@StatDave ,

The Details on the ASSESS functionality says this with regards to the link option:

You can check the link function instead of the jth covariate by using values of the linear predictor (X'beta)

in place of values of the jth covariate.

So does the statement accept a vector input for this?  Would you have to run the model under a given link to generate the predicted values, save that output, and then copy it into the statement?  What if you have a long (say 10K entries) predicted vector?  I don't see a way to provide this input.  Is it as simple as your code:

assess link/resample ?

SteveDenham

[StatDave]

I think that sentence was intended to describe what the LINK option in the ASSESS statement does internally for you. So, yes, the statement I showed is all you need.

[SteveDenham]

that is great news. Thank you.