HPLogistic: Odds Ratio Plot with Wald Chi-Square p-values

[PharmlyDoc]

The question of calculating a Wald Chi-Square test has been asked in the past, but I didn't see a complete answer. 

I have a large dataset of over 150,000 rows and 12 columns, and I am using HPLogistic, but HPLogistic does not have an option for an oddsratio plot.  

I want to create an odds ratio plot with a table of ORs, 95% CI, and Wald Chi-Square p-values on the right of the plot. 

I used this a resource for calculating the Wald Chi-Square:

https://support.sas.com/kb/53/376.html 

My question is if the methodology I am using to calculate odds ratios, 95% CI, and Wald Chi-Square p-values from the HPLogistic parameterestimates table is correct? (exponentiating the parameter estimates)

I am using the neuralgia dataset as an example. First with proc logistic, then with HPLogistic to compare. For this example I am not using a selection technique.

Data Neuralgia;
input Treatment $ Sex $ Age Duration Pain $ @@;
datalines;
P  F  68   1  No   B  M  74  16  No  P  F  67  30  No
P  M  66  26  Yes  B  F  67  28  No  B  F  77  16  No
A  F  71  12  No   B  F  72  50  No  B  F  76   9  Yes
A  M  71  17  Yes  A  F  63  27  No  A  F  69  18  Yes
B  F  66  12  No   A  M  62  42  No  P  F  64   1  Yes
A  F  64  17  No   P  M  74   4  No  A  F  72  25  No
P  M  70   1  Yes  B  M  66  19  No  B  M  59  29  No
A  F  64  30  No   A  M  70  28  No  A  M  69   1  No
B  F  78   1  No   P  M  83   1  Yes B  F  69  42  No
B  M  75  30  Yes  P  M  77  29  Yes P  F  79  20  Yes
A  M  70  12  No   A  F  69  12  No  B  F  65  14  No
B  M  70   1  No   B  M  67  23  No  A  M  76  25  Yes
P  M  78  12  Yes  B  M  77   1  Yes B  F  69  24  No
P  M  66   4  Yes  P  F  65  29  No  P  M  60  26  Yes
A  M  78  15  Yes  B  M  75  21  Yes A  F  67  11  No
P  F  72  27  No   P  F  70  13  Yes A  M  75   6  Yes
B  F  65   7  No   P  F  68  27  Yes P  M  68  11  Yes
P  M  67  17  Yes  B  M  70  22  No  A  M  65  15  No
P  F  67   1  Yes  A  M  67  10  No  P  F  72  11  Yes
A  F  74   1  No   B  M  80  21  Yes A  F  69   3  No
;
run;

proc logistic data=Neuralgia plots(only maxpoints=none)=(
		ODDSRATIO(logbase=10 TYPE=HORIZONTALSTAT)) namelen=30  ;
         class Treatment(ref='P') Sex(ref='F') / param=ref;
         model Pain = Treatment Duration Sex Age Sex*Age
/ clodds=wald orpvalue;
		 oddsratio treatment / cl=wald diff=ref ;
		 oddsratio Duration / cl=wald diff=ref;
         oddsratio Sex / cl=wald diff=ref;	
		 oddsratio Age / cl=wald diff=ref;
         ods output parameterestimates=logistic_paramest;
run;
quit;

PROC Logistic: Odds Ratios with 95% confidence intervals

Now, with HPLogistic, I outputted the parameterestimates and exponentiated Estimate, Lower, and Upper to calculate the odds ratios (plot) and 95% CI. And used Estimate and StdErr to calculate the WaldCHISQ and p-values. 

proc hplogistic data=neuralgia  ;
class Treatment(ref='P') Sex(ref='F') / param=ref;
model Pain = Treatment Duration Sex Age Sex*Age
/ cl association ctable=roc  link=logit lackfit rsquare; 
ods output parameterestimates=HPLogistic_paramest;
run;
quit;


data HPLogistic_paramest_OR;
set HPLogistic_paramest;
OddsRatioEst = exp(Estimate);
LowerCL = exp(Lower);
UpperCL = exp(Upper);
run;

data HPLogistic_paramest_OR_Wald;
set HPLogistic_paramest_oddsratios;
waldchisq=((Estimate)/StdErr)**2;
pvalue=1-probchi(waldchisq,1);
run;

PROC SQL;
DELETE FROM HPLogistic_paramest_OR_Wald
WHERE Parameter in ('Intercept');
quit;

proc format;
value $parameterfmt 
'Treatment A' = 'Treatment A vs P'
'Treatment B' = 'Treatment B vs P'
'Age' = 'Age at Sex=F'
'Age*Sex M'= 'Age at Sex=M'; 
run;


title "HPLogistic: Odds Ratios with 95% Confidence Limits and Wald Chi-Square p-values";
proc sgplot data=HPLogistic_paramest_OR_Wald noautolegend;
   scatter y=Parameter x=OddsRatioEst / xerrorlower=LowerCL xerrorupper=UpperCL
           markerattrs=(symbol=diamondfilled);
   refline 1 / axis=x;
   xaxis grid type=log label="Odds Ratio (log scale)";  /* specify log scale */
   yaxis grid display=(nolabel) discreteorder=data reverse;
yaxistable OddsRatioEst LowerCL UpperCL pvalue / location=inside
position=right ;
format Parameter $parameterfmt.;  
run;

 ods output parameterestimates with exponentiated Estimate as OddsRatioEst, Lower as LowerCL, and Upper as UpperCL. waldchisq and pvalue are also calculated from Estimate and StdErr where by waldchisq=((Estimate)/StdErr)**2; pvalue=1-probchi(waldchisq,1);

proc sgplot: Odds Ratio plot created from HPlogistic_paramest_OR_Wald table

Notice that the  "Age at Sex=M" odds ratio and 95% CI in the HPLogistic OR plot is different than that of the proc logistic OR plot - is the interaction still calculated at the mean of the continuous variable? Also, the HPLogistic OR plot is missing the "Sex M vs F at Age=70.05" – how do I compute this in HPLogistic? 

[SAS_Rob]

Because of the interaction, you will not get the odds ratios by simply exponentiating the parameter estimates.  You would need to include the continous variable set at its mean value.  You are essentially setting it at 0 by ignoring the interaction effect.

This usage note might be helpful as well.

https://support.sas.com/kb/24/455.html 

[PharmlyDoc]

@SAS_Rob 

Setting age to its mean? 

Proc logistic has the "At option" for the oddsratio statement: "odds ratio age /at(Age=70) cl=odds;"

So how do I calculate the odds ratio for Sex M vs F at Age=70 using the estimates from HPLogistic? 

[SAS_Rob]

It isn't really possible to get the CL from HPLOGISTIC.  The issue is that while you can compute the point estimate for the OR which would simply be:

exp(LB) where LB is b6+70.05*b7 and b6 is the estimate for age and b7 is the estimate for the interaction

you would not be able to get a standard error directly since HPLOGISTIC does not produce a table with the covariance matrix for the parameter estimates.  Without the standard error you would not be able to compute the CL.

[PharmlyDoc]

@SAS_Rob 

Do I necessarily need the covariance matrix for Odds ratio confidence intervals? I thought it was just necessary for calculating the Wald Chi-Square and p-value. I'm still trying to wrap my head around writing an equation to calculate the standard error if I had the covariance.  (aside from knowing that proc logistic has the covb option)

https://stats.stackexchange.com/questions/418231/how-does-sas-compute-standard-error 

https://stats.idre.ucla.edu/sas/seminars/analyzing-and-visualizing-interactions/

And I noticed that hplogistic allows for typing COV in the procedure statemennt...."proc HPLogistic data=neuralgia COV;" 

What does COV do? 

[SteveDenham]

In the absence of any statements like LSMEANS or LSMESTIMATE, you could explore the CODE statement.  According to the documentation, it writes DATA step code to a file that can be used to calculate predicted values (probabilities).  If you have that code, then you can plug in the values for which you want predicted probabilities, which could then be used to get ORs.  I suspect it will not easily give you a standard error, so Wald chi-squares might take more coding.  But if you have a predicted probability, you can calculate the variance using well-known formulas, and from there get the denominator as a function of a sum of the variances.

SteveDenham