Multivariate probit model in NLMIXED

[mgx]

Dear all,

I want to fit a multivariate probit model for two ordinal responses, but I do not see how I can estimate the correlation between the underlying latent continuous variables (and let this correlation depend on the covariates). My code now looks like this: 

proc nlmixed data=mvp.mydata qpoints=5 maxiter=1000 
maxfunc=10000 technique=quanew gconv=1e-8;
eta_d = beta_d_gender*gender+
beta_d_statsoc*STATSOC1+beta_d_gewest*GEWEST1;
eta_s = beta_s_gender*gender+beta_s_bmi*BMI2+
beta_s_age*AGE+
beta_s_statsoc1*STATSOC1+beta_s_statsoc2*STATSOC2;
if var='X' then do;
if response =0 then do;
lik = cdf('NORMAL',(int_drink1-eta_d));
end;
if response =1 then do;
lik = cdf('NORMAL',(int_drink2-eta_d)) -
cdf('NORMAL',(int_drink1-eta_d));
end;
if response =2 then do;
lik = cdf('NORMAL',(int_drink3-eta_d))-
cdf('NORMAL',(int_drink2-eta_d));
end;
if response =3 then do;
lik = 1 -cdf('NORMAL',(int_drink3-eta_d));
end;
ll = log(lik);
end;
if var='Y' then do;
if response =1 then do;
lik = cdf('NORMAL',(int_smoke1-eta_s));
end;
if response =2 then do;
lik = cdf('NORMAL',(int_smoke2-eta_s)) -
cdf('NORMAL',(int_smoke1-eta_s));
end;
if response =3 then do;
lik = 1-cdf('NORMAL',(int_smoke2-eta_s));
end;
ll = log(lik);
end;
model response ~ general(ll);
correlatie=corr(eta_d,eta_s);
ods output parameterestimates=parms_joint;
run;

SAS gives now the error 'ERROR: Cannot find a library containing subroutine CORR.'. Does anyone know how I can estimate the correlation between the latent variables?

Suggestions for other procedures that can do this are also welcome!

Kind regards,

Marc

[SteveDenham]

Do you have estimates of the variances and the covariance? You could then try to calculate a correlation that way by using PROC FCMP or macro code to get a value for the correlation. No guarantees, though.

SteveDenham

[Rick_SAS]

It sounds like you want the sample Pearson correlation between the latent variables in the model, not a correlation between parameters. If so, use the PREDICT statement to write the predicted value to a data set and then use PROC CORR. It will look something like this (untested)

proc nlmixed ...;
...
   predict eta_d  out=Var1;  /* output variable is named Pred */
   predict eta_s out=Var2;  /* same */
run;

data All;
merge Var1(keep=Pred rename=(Pred=eta_d)) 
      Var2(keep=Pred rename=(Pred=eta_s)); 
run;

proc corr data=All;
var eta_d eta_s;
run;

[StatDave]

You should be able to do this more simply using the joint modeling approach available in PROC GLIMMIX. See the example titled "Joint Modeling of Binary and Count Data" in the GLIMMIX documentation.