BookmarkSubscribeRSS Feed
☑ This topic is solved. Need further help from the community? Please sign in and ask a new question.
Rick_SAS
SAS Super FREQ

Please post the program you have written. That will enable us to suggest modifications.

Rick_SAS
SAS Super FREQ

>  it produced a constant value for the integrand. however, i need the integrand as a function of rho

No, the program I posted is a function of rho. You can see how the function changes by plotting the double integral versus rho:

 

rhoVals = do(-0.9, 0.6, 0.05);
IntegVals = j(1, ncol(rhoVals), .);
do i = 1 to ncol(rhoVals);
   rho = rhoVals[i];
   call quad(v, "outer", xinterval);
   IntegVals[i] = v;
end;

title "Double Integral as a Function of rho";
call series(rhoVals, IntegVals) grid={x y};

 

SGPlot45.png

 

As I have previously explained, this function is unbounded as rho -> +1. Therefore, it has no maximum on the interval (-1, 1).  

Do you have a reference for the MLE equations for weighted polychoric correlation?  Perhaps you copied an equation incorrectly?

 

hrmannan
Calcite | Level 5

Dear Dr Rick

the methodology for estimating weighted polychoric correlation taken from the doc

written by Paul Bailey, Ahmad Emad, Ting Zhang, Qingshu Xie,Source: vignettes/wCorrFormulas.Rmd.

the link for the doc is-

https://american-institutes-for-research.github.io/wCorr/articles/wCorrFormulas.html

 

 i couldnt paste the formulas as the format is not readable.  kindly see where i went wrong with the formulas . First the double integral of a bivariate normal density over finite limits (both are proportions so below 1) has to be solved and then the loglikelihood is to be maximized, the log likelihood is-

=sum(weights*log of the double integral)

you showed that this double integral doesnt converge. then where is the error.

looking forward to your kind help.

regards

Haider Mannan

 

hrmannan
Calcite | Level 5

dear Dr. Rick, I have sent you link for the paper. after taking log of the function and double-integrating it still seems to give unbounded function. the graph below shows this. how do I get a bounded function so I can maximise sum(weights*integrand) with respect to rho using my data when weights are survey weights? kindly help. The SAS codes & outputs are:

/* double integration for loglikelihood for estimating weighted polychoric correlation */
proc iml;
start inner(y) global(g_x, rho);
x = g_x;
return(-.5*x#x/(1-rho**2)-.5*y#y+rho#x#y/(1-rho**2)) ;
finish;

/* TEST THE INNER MODULE */
rho=0.5; /* global value of RHO */
g_x = 0; /* for testing: evaluate inner integral when x=g_x */
yinterval = {.5 .8};
call quad(z, "inner", yinterval);
print z;
free g_x; /* undefine g_x; it'll be set in the OUTER module */

/* g_x is the global variable that holds the value of x for the outer integral.
This is a constant when we evaluate the inner integral. */
start outer(x) global(g_x, rho);
g_x = x; /* save copy of x so the inner integral can use it */
yinterval = 1.5 || 3;
/* evaluate w = inner integral when x is held constant */
call quad(w, "inner", yinterval);
return (w);
finish;

/* TEST THE OUTER MODULE */
xinterval= {.6 .9};
call quad(v, "outer", xinterval);
print v;

rhoVals = do(-0.9, .6, 0.05);
IntegVals = j(1, ncol(rhoVals), .);
do i = 1 to ncol(rhoVals);
rho = rhoVals[i];
call quad(v, "outer", xinterval);
IntegVals[i] = v;
end;

title "Double Integral as a Function of rho";
call series(rhoVals, IntegVals) grid={x y};

 

Double Integral as a Function of rho

z
-0.0645

v
-0.846

hrmannan_1-1699951773556.png

 



 

 

hrmannan
Calcite | Level 5

Dear Dr Rick

the methodology for estimating weighted polychoric correlation taken from the doc

written by Paul Bailey, Ahmad Emad, Ting Zhang, Qingshu Xie,Source: vignettes/wCorrFormulas.Rmd.

the link for the doc is-

https://american-institutes-for-research.github.io/wCorr/articles/wCorrFormulas.html

 

 i couldnt paste the formulas as the format is not readable.  kindly see where i went wrong with the formulas . First the double integral of a bivariate normal density over finite limits (both are proportions so below 1) has to be solved and then the loglikelihood is to be maximized, the log likelihood is-

=sum(weights*log of the double integral)

you showed that this double integral doesnt converge. then where is the error.

looking forward to your kind help.

regards

Haider Mannan

Rick_SAS
SAS Super FREQ

Thanks for the link. You should remove the "Problem is solved" indicator if you want people to work on this problem. If you mark it as solved, people assume your issue is resolved.

hrmannan
Calcite | Level 5

dear Dr. Rick

maybe I didnt send the full link for the R doc on calculating weighted polychoric correlation. here it is-

file://ad.uws.edu.au/dfshare/HomesCMB$/30043787/My%20Documents/Phillipa's%20ABS%20data%20based%20study%20aug%202022/wCorr%20Formulas%20%E2%80%A2%20wCorr.html

Rick_SAS
SAS Super FREQ

Let me make sure I understand what you want. In the article

Bailey, Emad, Zhang, and Xie (2023)
"wCorr Formulas"
https://american-institutes-for-research.github.io/wCorr/articles/wCorrFormulas.html

you want to calculate the probability given by the following double integral inside the red box?

 

doubleintegral.png

 

hrmannan
Calcite | Level 5
Dear Dr Rick
Using the article formulas I want to determine the weighted polychoric correlation (rho). The first step is to find the joint probability function which is a double integral of a bivariate density of two standard normal (latent) variables over finite boundaries as given which have to be determined comparing the observed ordinal variables against sum of survey weights in my data. The next step is to use this result of double integral to my data to maximise the log likelihood function with respect to rho when loglikelihood is sum(weight*log(double integral)). In the lit search I found Simpsons one third rule being suggested to find the double integral of a bivariate normal density. And there are other ways.thanks.
Haider
hrmannan
Calcite | Level 5
Yes the formulas you mention from the article are the ones to be used. I have just explained the steps as I understood. Thanks
Rick_SAS
SAS Super FREQ

Great. The integrand in that formula is the standard bivariate normal density function (PDF). You do not need to use the QUAD function to compute the double integral. Instead, you can use the PROBBNRM function in SAS to integrate the bivariate PDF. The PROBBNRM function assumes that the lower limits of integration are minus infinity (just like the CDF function does), but you can use the rules of integration to obtain the integral over a rectangular region.

 

Let (a,b) be the limits of integration for the X variables. (These are the "theta prime" limits in the paper.) Let (c,d) be the limits of integration for the Y variable. (These are the "theta" limits in the paper.) Then you can use the following function to evaluate the double integral on the rectangle [a,b]x[c,d]:

proc iml;
/* compute the probability under the bivariate normal cdf with 
   correlation parameter rho for the rectangular region [a,b]x[c,d] */
start BivarProb(a, b, c, d, rho);
   Prob = PROBBNRM(b,d,rho) - PROBBNRM(a,d,rho) - PROBBNRM(b,c,rho) + PROBBNRM(a,c,rho);
   return Prob;
finish;

/* define region of integration to be the rectangular region ([a,b]x[c,d] */
a = -1;   b = 0;
c =  0;   d = 0.5;

/* prob when rho = 0 */
rho = 0;
Prob = BivarProb(a, b, c, d, rho);
print rho Prob;
/* prob when rho = -0.5 */
rho = -0.5;
Prob = BivarProb(a, b, c, d, rho);
print rho Prob;

/* compute the prob for a range of rho values */
rhos = do(-0.9, 0.9, 0.05);
Probs = j(1, ncol(rhos), .);
do i = 1 to ncol(rhos);
   rho = rhos[i];
   Probs[i] = BivarProb(a, b, c, d, rho);
end;
title "Pr[ (x,y) in [a,b]x[c,d] ]";
title2 "(a,b,c,d) = (-1,0,0,0.5)";
call series(rhos, Probs) grid={x y};

For the full program, call the BivarProb function from inside the objective function that you are optimizing. 

 

If this is the information you want, then please mark this response as the correct solution so that future programmers will be able to find the solution, too.

hrmannan
Calcite | Level 5

Dear Dr Rick, thanks very much. I am trying to understand the outputs. Does the graph mean that the ML estimate of rho is aroundr -0.9 where the probability is maximum? 

hrmannan_0-1700229482288.png

another point is that the regions of integration you have assigned are (-1,0) and (0,0.5). According to the paper i sent, shouldn't these inner and outer limits of integration be non-negative fractions? please let me know.

Haider

Rick_SAS
SAS Super FREQ

My response shows how you can compute the bivariate normal probability by using the PROBBNRM function in SAS. You can use this function instead of the QUAD function to compute the double integrals in your problem. My response does not do any optimization, nor does it compute the polychoric correlation. You will have to write that part of the program yourself. 

 

> Does the graph mean that the ML estimate of rho is aroundr -0.9 where the probability is maximum? 

No. This is not the graph of the LL. It is a graph of the bivariate normal probability (the double integral) in the rectangle [-1, 0] x [0, 0.5] for various choices of rho.

 

> According to the paper i sent, shouldn't these inner and outer limits of integration be non-negative fractions?

No, I don't think so. The limits of integration are the cut points of the underlying (latent) normal variate which gives rise to the (observed) ordinal variable. Yes, they are based on the observed cumulative probabilities of the marginal counts, which are nonnegative fractions, but the actual cut points are of the form quantile("Normal", cumprob[i]), where cumprob[i] are the cumulative probabilities. So, the cut points are in (-infinity, infinity).

hrmannan
Calcite | Level 5
Dear Dr Rick
Thanks. Please explain how you got the limits (-0.5,0) and (0,1). The paper stated two ways or is it just one method. Kindly explain.
Haider
hrmannan
Calcite | Level 5
Correction I meant how did you find the limits (-1,0) and (0,0.5) for the probability integral of a bivariate normal density?
Haider

sas-innovate-2024.png

Don't miss out on SAS Innovate - Register now for the FREE Livestream!

Can't make it to Vegas? No problem! Watch our general sessions LIVE or on-demand starting April 17th. Hear from SAS execs, best-selling author Adam Grant, Hot Ones host Sean Evans, top tech journalist Kara Swisher, AI expert Cassie Kozyrkov, and the mind-blowing dance crew iLuminate! Plus, get access to over 20 breakout sessions.

 

Register now!

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 42 replies
  • 2787 views
  • 20 likes
  • 3 in conversation