pdf("g12.pdf");options(width=64)
#setwd("C:\\Users\\kolassa\\Class553")
setwd("~/Taught1/960-553/Data")
<pre>
#<a name=block1> Block 1</a>
#*******************************************************/
#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')
summary(glm(cbind(y,n-y)~li+sex+aop,data=goorin,
   family="binomial"))
library(PHInfiniteEstimates)#For network and inference; Calculate test results for exact logistic regression.
(gnet<-network(goorin[,1:3],n=goorin[,4],conditionon=1:3,
   resp=goorin[,5]))
cbind(gnet$possible,gnet$count/sum(gnet$count))
inference(gnet)
cat('\n Exact Logistic Regression  \n')
#I can't find an R routine that does exact logistic regression
#without Monte Carlo, other than my own package.
library(elrm)#For elrm; Calculate test results for exact logistic regression via Monte Carlo.
elrmout<-elrm(y/n~li+sex+aop,~aop,data=goorin,iter=1000000)
names(elrmout)
elrmout$distribution
#<a name=block2> Block 2</a>
#****************************************************/
# Housing data from Cox and Snell (1980), Example W,*/
# on satisfaction with housing in Copenhagen.  Data */
# was extracted from R package MASS.  Data may also */
# be found at                                       */
# www.stat.ucla.edu/data/cox-and-snell/exampleW.data*/
#****************************************************/
housing<-read.table("housing.dat",header=FALSE,
  col.names=c("nn","sat","infl","type","contact","freq"))
housing$nsat<-(housing$sat=="Medium")+2*(housing$sat=="High")
housing$ninf<-(housing$infl=="Medium")+2*(housing$infl=="High")
housing$ncont<-1*(housing$contact=="High")
# Some of these numeric values are zero.  Our LCA software 
# needs them to be positive.
housing$nsat<-housing$nsat+1
housing$ninf<-housing$ninf+1
housing$ncont<-housing$ncont+1
# Previously we assigned numerical codes to the above
# variables, for use in linear by linear interactions.
# We have not previously given numerical codes to housing 
# type, since it is not ordered.  Assign positive 
# integer codes now.
housing$ntype<-as.numeric(factor(housing$type))
library(poLCA)#For poLCA; Perform latent class analysis.
fullhouse<-housing[rep(1:72,housing$freq),]
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 positive 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 
# Default number of latent classes is 2.
(lcaout<-poLCA(a,fullhouse))
# Try a larger number of classes
(lcaout3<-poLCA(a,fullhouse,nclass=3))
</pre>