postscript("l09.ps");options(width=68)
#****************************************************/
# 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<-as.data.frame(scan("housing.dat",
  what=list(nn=0,sat="",infl="",type="",contact="",freq=0)))
housing$nsat<-(housing$sat=="Medium")+2*(housing$sat=="High")
housing$ninf<-(housing$infl=="Medium")+2*(housing$infl=="High")
housing$ncont<-(housing$contact=="Medium")+
  2*(housing$contact=="High")
fullhouse<-housing[rep(1:72,housing$freq),]
housing$rsat<-as.character(housing$sat) 
housing$rsat[housing$sat=="Medium"]<-"Intermediate"
housing$rsat<-as.factor(housing$rsat)
#*************************************************************/
# Mark C: Polytomous Regression.  Baseline logit model first.*/
#*************************************************************/
cat('\n Baseline logit model  \n')
#To make the numbers come close to those we get from SAS,
#rename the categories to change the ordering. Resulting 
#effect coefficients will be double those of SAS, since  
#proc logistic by default uses the -1,1 parameterization.
housing$rsat<-as.character(housing$sat) 
housing$rsat[housing$sat=="High"]<-"Top"
housing$rsat<-as.factor(housing$rsat)
library(nnet)
multinom(rsat~contact,data=housing,weight=freq)
#*************************************************************/
# Mark D: Cumulative logit model.  First reuse a known tech- */
# nique.  Simple model exhibits different parameterizations. */
#*************************************************************/
cat('\n Binary response for comparison purpose  \n')
housing$oksat<-housing$nsat>1
glm(oksat~factor(contact),family=binomial,data=housing,
        weight=freq)
#*************************************************************/
# The following data set will not give us what we want, be-  */
# cause it will order sat incorrectly.  Recode data to put   */
# medium between high and low, and add descending to make    */
# low the baseline.  First fit the simplest model.  This     */
# should give logits of raw proportions.  Check this.        */
#*************************************************************/
cat('\n Ordinal Logistic Regression with wrong order  \n')
cat('\n Ordinal Logistic Reg. with right order, empty model  \n')
library(MASS) 
# Venables and Ripley, Modern Applied Statistics with S    
# polr does regression with ordered response.   Allows var-
# ious link functions.                                     
polr(factor(nsat)~1,weights=freq,data=housing)
cat('\n Ordinal Logist. Reg. with added covariates  \n')
polr(factor(nsat)~infl+type+contact,weights=freq,data=housing)
#*************************************************************/
# Mark E: Cf. the complimentary log log.  Odds ratios are no */
# longer given, since they don't come naturally out of the   */
# model.  Intercepts have a completely different meaning and */
# won't be compared.  Linear parameters are multiplied by a  */
# constant, approx., and z statistics are similar.           */
#*************************************************************/
cat('\n Complimentary log log  \n')
polr(factor(nsat)~infl+type+contact,weights=freq,data=housing,
  method="cloglog")