pdf("g09.pdf");options(width=64)
#setwd("C:\\Users\\kolassa\\Class553")
setwd("~/Taught1/960-553/Data")
<pre>
#<a name=block1> Block 1</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")
housing
#*****************************************************/
# 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.
library(nnet)#For multinom; Fits the multinomial regression model.
#Note below that the levels of sat are ordered,
#but the ordering doesn't matter except to determine
#which category is baseline.
multinom(factor(sat)~contact,data=housing,weight=freq)
#<a name=block2> Block 2</a>
#**********************************************************/
# Cumulative logit model.  First reuse a known technique. */
# 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)
library(MASS)#For polr; Fits ordinal categorical regression models.
#*************************************************************/
# 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.        */
#*************************************************************/
#Note below that the levels of rsat are ordered.  Default
#is alphabetical.  That's why I changed "High" to "Top".
#You can also do this as an argument to as.factor.
cat('\n Ordinal Logistic Regression with wrong order  \n')
polr(factor(sat)~1,weights=freq,data=housing)
cat('\n Ordinal Logistic Reg. with right order, empty model\n')
# polr regresses with ordered response, and Allows various 
# 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)
housing$rsat<-as.character(housing$sat) 
housing$rsat[housing$sat=="High"]<-"Top"
housing$rsat<-as.factor(housing$rsat)
(polr1<-polr(rsat~infl+type+contact,weights=freq,
   data=housing))
library(VGAM)#For vglm; Fits an ordinal baseline logit model and some variants.
#Fit a separate set of parameters for each transition.  There
#are ways to add constraints to keep parameters constant from
#one level of response to the next.
vglm(ordered(rsat)~infl+type+contact,
   data=housing,family=cumulative,weight=freq)
#<a name=block3> Block 3</a>
vglm(ordered(rsat)~infl+type+contact,
   data=housing,family=cratio,weight=freq)
vglm(ordered(rsat)~infl+type+contact,
   data=housing,family=acat,weight=freq)
#<a name=block4> Block 4</a>
#**************************************************************/
# Beetle data, from SAS documentation for complimentary log   */
# log link, after reformatting to prepare data for treatment  */
# by cloglog analysis.  Beetles were exposed to insecticides  */
# and observed for 13 days.  The day on which the beetle died */
# is recorded, with beetles lasting all 13 days being recorded*/ 
# as 14.  Sex of the beetle (m=1, f=2), and concentrations of*./
# insecticide, were also recorded.                            */
#**************************************************************/
beetledays<-read.table("beetledays.dat",header=TRUE)
beetledays$sex<-c("M","F")[beetledays$sex]
# Starting model with -1 forces there to be no intercept, and 
# hence effects for each day are estimable.  Effects are
# discrete hazards, which equal conditional probability of
# failure at that time given survival up until tht time.
# Baseline concentration is .20, measired in mg/cm^2.
summary(glm(y~-1+factor(day)+sex+factor(conc), data=beetledays, 
   weights=freq, family=binomial(link="cloglog")))
#<a name=block5> Block 5</a>
#************************************************************/
# 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.          */
#************************************************************/
polr1
cat('\n Complimentary log log  \n')
polr(factor(rsat)~infl+type+contact,weights=freq,data=housing,
  method="cloglog")
</pre>