postscript("l12.ps");options(width=68)
#*************************************************************/
# Mark A: Exact logistic regression.                         */
#*************************************************************/
#*******************************************************/
#Goorin, A.M., Perez--Atayde, A., Gebhardt, M., and    */
#Andersen, J. (1987) Weekly High--Dose Methotrexate and*/
#Doxorubicin for Osteosarcoma: The Dana--Farber Cancer */
#Institute/The Children's Hospital -- Study III, J.    */
#Clin. Oncol. 5, 1178--1184, quoted by Hirji, Mehta,   */
#and Patel (1987).  Survival of osteogenic sarcoma pa- */
#tients,  with covariates li (lymphatic infiltration), */
#sex, and aop (any osteoid pathology), number of sub-  */
#jects, and number of 5 year survivors.                */
#*******************************************************/
goorin<-read.table("goorin.dat",head=F)
names(goorin)<-c("li","sex","aop","n","y")
#*************************************************************/
# Examine the data.  All individuals without li survive.     */
#*************************************************************/
print(goorin)
#************************************************************/
# Fit using exact logistic regression.  The estimate of the */
# li effect is -Infinity, since the best fit sets success   */
# abilities to 1 when li=0.                                 */
#************************************************************/
cat('\n Regular Logistic Regression for Sarcoma Data  \n')
glm(cbind(y,n-y)~li+sex+aop,data=goorin,family="binomial")
cat('\n Exact Logistic Regression  \n')
#I can't find an R routine that does exact logistic regression
#without Monte Carlo.
library(elrm)
elrm(y/n~li+sex+aop,~li,data=goorin)
#*************************************************************/
# Housing data from Cox and Snell (1980), Example W, on sat- */
# isfaction with housing in Copenhagen.  Data was extracted  */
# from R package MASS.                                       */
#*************************************************************/
housing<-as.data.frame(scan("housing.dat",
  what=list(nn=0,sat="",infl="",type="",contact="",freq=0)))
housing$nsat<-1+(housing$sat=="Medium")+2*(housing$sat=="High")
housing$ninf<-1+(housing$infl=="Medium")+2*(housing$infl=="High")
housing$ncont<-1+(housing$contact=="High")
housing$ntype<-1+(housing$type=="Apartment")+
  2*(housing$type=="Atrium")+3*(housing$type=="Terrace")
fullhouse<-housing[rep(1:72,housing$freq),]
library(poLCA)
a<-cbind(nsat,ninf,ncont,ntype)~1
#Note that poLCA prints output even when you save output
#To make poLCA work, you need to
#1. Code all variables to consecutive integers,
#a. if you have an empty category, you'll get something not
#   what you anticipate.
#2. Equivariant to
#a. Ordering of categories 
#b. Ordering of variables 
#3. Equivariant to
lcaout<-poLCA(a,housing[rep(1:72,housing$freq),])
library(polycor); library(psych)
fun.testpoly<-function(rho=.5,nsamp=50,nmc=10000,progress=F){
  out<-array(NA,c(7,nmc))
  for(i in seq(dim(out)[2])){
     if(progress&(i==((round(i/100)*100)))) cat("i=",i)
     x<-rnorm(nsamp)
     y<-rho*x+sqrt(1-rho^2)*rnorm(nsamp)
#To make sure the table is square, make sure all potentail
#values are classified along each direction, by adding some
#fake observations, and then subtracting them out.
     tt<-table(c(x>0,0,1),c(y>0,0,1))-diag(rep(1,2))
     out[1,i]<-polychor(tt)
     out[5,i]<-cohen.kappa(tt)$kappa
     tt<-table(c((x> -1)+(x>1),0,1,2),c((y> -1)+(y>1),0,1,2))-
        diag(rep(1,3))
     out[2,i]<-polychor(tt)
     out[6,i]<-cohen.kappa(tt)$kappa
     tt<-table(c((x> -1)+(x>1)+(x>0),1,2,3,0),
         c((y> -1)+(y>1)+(y>0),1,2,3,0))-diag(rep(1,4))
     out[3,i]<-polychor(tt)
     out[7,i]<-cohen.kappa(tt)$kappa
     out[4,i]<-cor(x,y)
  }
  return(out)
}
fun.plottestpoly<-function(out,rho,nsamp){
  par(mfrow=c(2,1))
  db<-density(out[4,])
  plot(db,main="Correlation Estimates from Bivariate Normal",
     xlab="rho",lty=1,col=1,
     sub=paste("True correlation",rho,"Sample size",nsamp))
  legend(quantile(db$x,.01),max(db$y),lty=1:4,col=1:4,
     legend=c("Sample Correlation",paste("Poly.",2:4,"way")))
  for(j in 1:3) lines(density(out[j,]),lty=j+1,col=j+1)
  dbs<-list(density(out[5,]),density(out[6,]),density(out[7,]))
  rx<-range(dbs[[1]]$x)
  ry<-range(c(dbs[[1]]$y,dbs[[2]]$y,dbs[[3]]$y))
  plot(rx,ry,type="n", main="Kappa Estimates from Bivariate Normal",
     xlab="kappa",
     sub=paste("True correlation",rho,"Sample size",nsamp))
  legend(rx[1],ry[2],lty=2:4,col=2:4,
     legend=c("2 way","3 way","4 way"))
  for(j in 1:3) lines(dbs[[j]],lty=j+1,col=j+1)
  return(out)
}
out<-fun.testpoly(nmc=1000)
out1<-fun.plottestpoly(out,.5,50)