Re: [R] polychoric correlation error

From: John Fox <jfox_at_mcmaster.ca>
Date: Sat 05 Aug 2006 - 12:48:16 EST


Dear Janet,

Because you didn't set the value of the random-number generator seed, your example isn't precisely reproducible, but the problem is apparent anyway:

> set.seed(12345)
> n<-100
> test.x<-rnorm(n, mean=0, sd=1)
> test.c<-test.x + rnorm(n, mean=0, sd=.5)
> thresh.x<-c(-2.5, -1, -.5, .5, 1000)
> thresh.c<-c(-1, 1, 2, 3, 1000)
>
> discrete.x<-discrete.c<-vector(length=n)
>
> for (i in 1:n) {
+ discrete.x[i]<-which.min(thresh.x < test.x[i] ) + discrete.c[i]<-which.min(thresh.c < test.c[i] ) }
>
> table(discrete.x, discrete.c)

          discrete.c
discrete.x 1 2 3 4 5

         2 12  1  0  0  0
         3  3 12  0  0  0
         4  2 19  2  0  0
         5  0 18 21  9  1

>
> cor(test.x, test.c)
[1] 0.9184189
>
> pc <- polychor(discrete.x, discrete.c, std.err=T, ML=T)
Warning messages:
1: NaNs produced in: log(x) 
2: NaNs produced in: log(x) 
3: NaNs produced in: log(x) 

> pc

Polychoric Correlation, ML est. = 0.9077 (0.03314) Test of bivariate normality: Chisquare = 3.103, df = 11, p = 0.9893

  Row Thresholds

  Threshold Std.Err.
1  -1.12200   0.1609
2  -0.56350   0.1309
3   0.03318   0.1235


  Column Thresholds

  Threshold Std.Err.
1   -0.9389   0.1489
2    0.4397   0.1292
3    1.2790   0.1707
4    2.3200   0.3715

>

The variables that you've created are indeed bivariate normal, but they are highly correlated, and your choice of cut points makes it hard to estimate the correlation from the contingency tables, apparently producing some difficulty in the maximization of the likelihood. Nevertheless, the ML estimates of the correlation and thresholds for the set of data above are pretty good. (In your case, the optimization failed.)

BTW, a more straightforward way to create the categorical variables would be

discrete.x <- cut(test.x, c(-Inf, -2.5, -1, -.5, .5, Inf)) discrete.c <- cut(test.c, c(-Inf, -1, 1, 2, 3, Inf))

I hope this helps,
 John



John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox

> -----Original Message-----
> From: r-help-bounces@stat.math.ethz.ch
> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of
> Rosenbaum, Janet
> Sent: Friday, August 04, 2006 5:49 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] polychoric correlation error
>
>
> Dear all,
>
> I get a strange error when I find polychoric correlations
> with the ML method, which I have been able to reproduce using
> randomly-generated data.
>
> What is wrong?
> I realize that the data that I generated randomly is a bit
> strange, but it is the only way that I duplicate the error message.
>
>
> > n<-100
> > test.x<-rnorm(n, mean=0, sd=1)
> > test.c<-test.x + rnorm(n, mean=0, sd=.5) thresh.x<-c(-2.5, -1, -.5,
> > .5, 1000) thresh.c<-c(-1, 1, 2, 3, 1000)
> >
> > discrete.x<-discrete.c<-vector(length=n)
> >
> > for (i in 1:n) {
> + discrete.x[i]<-which.min(thresh.x < test.x[i] )
> + discrete.c[i]<-which.min(thresh.c < test.c[i] ) }
> > pc<-polychor(discrete.x, discrete.c, std.err=T, ML=T)
> Error in optim(c(optimise(f, interval = c(-1, 1))$minimum,
> rc, cc), f, :
> non-finite finite-difference value [1]
> In addition: There were 50 or more warnings (use warnings()
> to see the first 50)
> > print(pc)
> Error in print(pc) : object "pc" not found
> > warnings()
> Warning messages:
> 1: NaNs produced in: log(x)
> 2: NA/Inf replaced by maximum positive value
> 3: NaNs produced in: log(x)
>
>
> ---
>
> Thanks,
>
> Janet
>
> --------------------
>
> This email message is for the sole use of the intended\ > ...{{dropped}}



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sat Aug 05 14:06:10 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Sat 05 Aug 2006 - 18:17:57 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.