How to compare Odds Ratios from models with different outcomes, same predictors

[Minhtrang]

Dear all,

I ran logistic procedure to get the Odds Ratios (CI) from these models

AAC=age

CAC=age

AAC: Abdominal aorta calcification (outcome 1)

CAC: Coronary artery calcification (outcome 2)

Age: predictor

I received comment from a reviewer, saying that there is a statistical method to compare Odds Ratios of Age from these 2 models.

I have no experience with this kind of comparison. It's more common to compare nested models than models with different outcomes.(Although the outcomes in this example may be similar in pathophysiology)

I would like to hear from your experience in solving this issue.

Thank you very much in advance.

[StatDave]

Since the odds ratios are just a function of the parameter estimates, testing the equality of the age odds ratios is equivalent to testing the equality of the age parameters. Since these parameters are estimated in two separate models on two different response variables, this requires a multivariate model and is one of the specialized cases that PROC CATMOD can help with. The following fits the multivariate model and uses a CONTRAST statement to compare the age parameters on the responses. The DIRECT statement is needed to treat age as continuous (CATMOD treats all predictors as categorical by default - the opposite of what LOGISTIC and other procedures do). Note that the the parameters of the multivariate model will not match those of the separate models on the response since the multivariate model takes the correlation among the responses into account. 

proc catmod;
response logits;
direct age;
model aac*cac=age / noprofile;
contrast 'age' @1 age 1 @2 age -1;
run; quit;

[Minhtrang]

Dear StatDave_sas,

Thank you very much for your reply.

I tried your suggestion and found an ERROR as

"The response functions are linearly dependent since the number of
functions per population, 2 is greater than or equal to the number of
response levels, 1, in population 1."

Would you mind showing what should be done to solve this issue?

[StatDave]

CATMOD was created many years ago and was intended as a modeling procedure for purely categorical data. Consequently, its method often has trouble with truly continuous predictors which have little or no replication of levels. This message is the typical result with such predictors. This is why PROC LOGISTIC should always be preferred to CATMOD for any sort of logistic modeling. One way around this problem in CATMOD is to categorize the continuous age variable so that most or all response combinations appear in all of the categorized age levels. The following statements round age to 5 year groups and reruns the analysis. You might have to adjust the categorization (and you don't have to do it with simple rounding) to avoid the error.

data new; set original; 
 agernd=round(age,5); 
 run;
proc catmod;
 response logits;
 direct agernd;
 model aac*cac=agernd / noprofile;
 contrast 'agernd' @1 agernd 1 @2 agernd -1;
 run; quit;

Another option is to use PROC QLIM in SAS/ETS followed by the NLEstimate macro to do the test. QLIM doesn't have the problem  that CATMOD has with continuous predictors, but QLIM can only fit a bivariate probit, not logit, model. These statements fit the bivariate probit model and save the model parameter estimates and covariance matrix which are then used in the NLEstimate macro to compare the two age parameters (which are the 2nd and 4th parameters shown in the QLIM results).

proc qlim covb;
 model aac cac=age / discrete;
 ods output parameterestimates=pe covb=c;
 run;
%NLEstimate(inest=pe, incovb=c, f=b_p2-b_p4, label=parm diff)

[Minhtrang]

Dear StatDave_sas,

I still have the same error using the suggested codes.

I found a paper which produced similar wanted results as following

I want to calculate the Homogeneity P value (the last column). 

The author mentioned about the statistical method  to get that P_value as

"To test whether the magnitude of association of each risk factor with calcification differed between the 2 sites, we
tested for equivalence of the odds ratios (ORs) by computing a Mantel-Haenszel  statistic based on the weighted sum of the
squared deviations of the stratum-specific log ORs from their weighted mean"

I would like to know if we can perform this analysis in SAS.

Thank you for your kind help and patience.

[JacobSimonsen · Accepted solution]

you can do it by first making a dataset with one response variable and a new variable indicating what type of response it is (that dataset will be twice as long as the original). Then testing interaction between the response type and age.

data mydata;
  do i=1 to 1000;
  age=rand('uniform',0,10);
    AAC=rand('bernoulli',1/(1+exp(5-0.2*age)));
    CAC=rand('bernoulli',1/(1+exp(4-0.25*age)));
  output;
  end;
ruN;
proc logistic data=mydata;
  model AAC(event="1")=age;
run;
proc logistic data=mydata;
  model CAC(event="1")=age;
run;
data mydata2;
  set mydata;
  type=1; respons=AAC; output;
  type=2; respons=CAC; output;
run;
 
*now test interaction (the type 3 test in this model): ;
proc logistic data=mydata2;
  class type/param=glm;
  model respons(event="1")=type age type*age;
run;

[Minhtrang]

Dear JacobSimonsen,

Thank you so much for your help! Now I understand the setting up of data for running this kind of analysis.

I have 1 more question only.

AAC and CAC correlates with each other ( for example, correlation coefficient=0.5).

Can we also take into account this correlation in the logistic model that you suggested?

[JacobSimonsen]

Yes, you can include the value of one response as predictor for the other response. Though, it is my gut feeling that something will go wrong if you do that both ways, because the predictor then will depend on the outcome.

If you want to use AAC as a predictor for CAC, and in same model you also want to test if age has same effect on the two outcomes, you can add one more variable, that are constant when the outcome takes the value AAC, but when the response is CAC it takes its value should depend on AAC. As below (which is an extension of the program I gave you before) the "type_" variable will estimate the log(oddsratio) of the AAC on CAC.

data mydata2;
  set mydata;
  type=1; type_=1; respons=AAC; output;
  type=2; type_=2*10+AAC; respons=CAC; output;
run;
 
*now test interaction (the type 3 test in this model): ;
proc logistic data=mydata2;
  class type type_ /param=glm;
  model respons(event="1")=type type_ age type*age;
run;