#
# Polychoric r for RxC table
#
polychoric.r <- function(x) {
  nr<-nrow(x)
  nc<-ncol(x)
  tot<-sum(x)
  row.thresh<-qnorm(cumsum(apply(x,1,sum))[-nr]/tot)
  col.thresh<-qnorm(cumsum(apply(x,2,sum))[-nc]/tot)
  poly.r<-0
  opt<-optim(par=c(row.thresh,col.thresh,poly.r), 
             fn=table.lik, 
             gr=NULL, 
             method="Nelder-Mead",
             lower=-Inf, upper=Inf,
             control=list(fnscale=-1),
             hessian=FALSE,
             nr, nc, x)
  res<-list()
  res$Observed.table<-x
  res$Expected.table<-tot* expected.prop(nr, nc, 
                                         opt$par[1:(nr-1)], 
                                         opt$par[nr:(nr+nc-2)],
                                         opt$par[nr+nc-1])
  res$r<-opt$par[nr+nc-1]
  res$LRTS<-poisson.lik(x,x/tot)-opt$value
  res$Df<-nr*nc-nr-nc
  res$Pr<-NA ; if (res$Df>0) res$Pr<-pchisq(res$LRTS,res$Df,lower=F)
  res
}
expected.prop <- function(nr, nc, row.thresh, col.thresh, poly.r) {
  require(mvtnorm)
  res<-matrix(nr=nr,nc=nc)
  big.row<-c(-Inf, row.thresh, Inf)
  big.col<-c(-Inf, col.thresh, Inf)
  r<-diag(2)
  r[1,2]<-r[2,1]<-poly.r
  for(i in 1:nr) for(j in 1:nc) {
    res[i,j]<-pmvnorm(lower=c(big.row[i],big.col[j]),
                      upper=c(big.row[i+1],big.col[j+1]),
                      cor=r)
  }
  res
}
table.lik <- function(par, nr, nc, x) {
  phat<-expected.prop(nr, nc, par[1:(nr-1)], 
                              par[nr:(nr+nc-2)],
                              par[nr+nc-1])
  poisson.lik(x,phat)
}
poisson.lik <- function(x,p) {
  ln <- function(x) { res<-log(x); res[x==0]<-0; res }
  sum(x*ln(p))
}