```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.                */
#*******************************************************/
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)
```