Technically this is not a refereed publication (yet), but I thought it would be nice to have a little space for the function at the centre of my master’s thesis.  It deals with the Bayesian estimation of the polychoric correlation coefficient both under bivariate normal and non-normal data.  It is based on this article written by James H Albert. I was a little bit surprised that such a useful article hadn’t really got a lot of citations, so I wanted to highlight his work by implementing it into R. In my code, I figured out how to deal with both right-skewed and left-skewed latent data (Albert’s article only handles right-skewed data) and it has proven to be very useful for handling difficult datasets, such as when one is dealing with low-count or zero-count cells in contingency tables.

The main function is called polycorGibbs() because it is based of a Gibbs sampler.

And the way it works is super easy:


library(MASS)

Q <- mvrnorm(100, c(0,0), matrix(c(1,.5,.5,1),2,2))

x <- Q[,1]
y <- Q[,2]

##Make sure you declare these as "integer", not "numeric" or "float"


xb <- as.integer(cut(y, breaks=c(-Inf, mean(y)-sd(y),mean(y)+sd(y), Inf), labels=c(1:3)))
yb <- as.integer(cut(y, breaks=c(-Inf, mean(y)-sd(y),mean(y)+sd(y), Inf), labels=c(1:3)))
plyk <- polycorGibbs(xb, yb, iter = 1e4, trace=FALSE, graph=FALSE)

Summary
mean median SD Numeric SD
rho 0.481 0.487 0.00109 0.0282
c1 -1.070 -1.070 0.00885 0.1530
c2 1.020 1.020 0.01190 0.1500
d1 -1.070 -1.060 0.01090 0.1550
d2 1.010 1.010 0.01090 0.1520

Correlation
rho c1 c2 d1 d2
rho NA 0.00163 -0.0384 -0.0153 -0.0323
c1 NA NA 0.2100 0.5080 0.2260
c2 NA NA NA 0.2330 0.5220
d1 NA NA NA NA 0.2220
d2 NA NA NA NA NA

I offer some extra functionalities like plotting the time series plot for the Markov Chain and other ways to check how well everything converged and whatnot.

Overall, the function is a tad bit finicky. But then again this was my first “major”  project and I pretty much had to teach myself Bayesian statistics and MCMC from scratch to get it to work, which is why I am very proud of it. But it will get better in time.