# # wisc_CDA.q # # PURPOSE: fit POM models to the diabetic retinopathy data # # AUTHOR: P. Heagerty # # DATE: 01/01/08 # updated 15 Jan 2007 # #--------------------------------------------------------------------------- # data <- matrix( scan("wisc.data"), ncol=23, byrow=T ) # wisc.all.data <- data.frame( r.ret = data[,1], l.ret = data[,2], c.ret = data[,3], r.edema = data[,4], l.edema = data[,5], r.refrac = data[,6], l.refrac = data[,7], r.iop = data[,8], l.iop = data[,9], age.diag = data[,10], duration = data[,11], glyc.hemo = data[,12], sys.bp = data[,13], dia.bp = data[,14], bmi = data[,15], pulse = data[,16], gender = data[,17], prot = data[,18], dose = data[,19], county = data[,20], r.level = data[,21], l.level = data[,22], c.level = data[,23] ) # wisc.data <- data.frame( y = wisc.all.data[,"r.level"], glyc.hemo = wisc.all.data[,"glyc.hemo"], log.duration = log( wisc.all.data[,"duration"] ), age.diag = wisc.all.data[,"age.diag"], dia.bp = wisc.all.data[,"dia.bp"] ) wisc.center.data <- data.frame( y = wisc.all.data[,"r.level"], glyc.hemo = wisc.all.data[,"glyc.hemo"] - 12, log.duration = log( wisc.all.data[,"duration"] ) - 2, age.diag = wisc.all.data[,"age.diag"] - 15, dia.bp = wisc.all.data[,"dia.bp"] - 80 ) # ### ### fit simple POM model (dia.bp80) ### # options( contr="contr.treatment" ) # ###### ###### load the necessary MASS fnxs ###### # ##### (a) OPTION 1 = read all of the needed files # source("polr.q") source("vcov.q") source("misc.q") # ##### (b) OPTION 2 = simply use the MASS library if installed # #library(MASS) # ##### create factor response # wisc.data$y <- factor( wisc.data$y ) wisc.center.data$y <- factor( wisc.center.data$y ) # attach( wisc.data ) # ##### dichotomized dia.bp # wisc.data$dia.bp80 <- as.integer( dia.bp >= 80 ) # ##### set minimum level of Y at 1 # fit0 <- polr( y ~ dia.bp80, data=wisc.data ) summary( fit0 ) # ##### continuous dia.bp # fit1 <- polr( y ~ dia.bp, data=wisc.center.data ) summary( fit1 ) # ## ### #### Fit with multiple predictors ### ## # # fit2 <- polr( y ~ glyc.hemo + log.duration + age.diag + dia.bp, data = wisc.center.data ) summary( fit2 ) # # # end-of-file...