setwd("I:\\classes\\regression 2010\\data")

data <- read.csv("senicfull.csv")
data$logLOS <- log(data$LOS)

pairs(~INFRISK+BEDS+logLOS, data=data, pch=16)

# adjusted variable plot approach
# look at the association between INFRISK and logLOS, 
# adjusting for BEDS

reg.xy <- lm(logLOS ~ BEDS, data=data)
res.xy <- reg.xy$residuals

reg.xz <- lm(INFRISK ~ BEDS, data=data)
res.xz <- reg.xz$residuals

plot(res.xz, res.xy, pch=16)
reg.res <- lm(res.xy ~ res.xz)
abline(reg.res, lwd=2)
reg.infrisk.beds <- lm(logLOS ~ BEDS + INFRISK, data=data)

############
# why NURSE is not associated, after adjustment for BEDS?

reg.nurse <- lm(logLOS ~ NURSE, data=data)
reg.nurse.beds <- lm(logLOS ~ NURSE + BEDS, data=data)

reg.xy <- lm(logLOS ~ BEDS, data=data)
res.xy <- reg.xy$residuals

reg.xz <- lm(NURSE ~ BEDS, data=data)
res.xz <- reg.xz$residuals

plot(res.xz, res.xy, pch=16)
reg.res <- lm(res.xy ~ res.xz)
abline(reg.res, lwd=2)

#######################
# what about the other way?  what about why BEDS is assoc after adjustment
# for NURSE?

reg.xy <- lm(logLOS ~ NURSE, data=data)
res.xy <- reg.xy$residuals

reg.xz <- lm(BEDS ~ NURSE, data=data)
res.xz <- reg.xz$residuals

plot(res.xz, res.xy, pch=16)
reg.res <- lm(res.xy ~ res.xz)
abline(reg.res, lwd=2)
reg.nurse.beds <- lm(logLOS ~ NURSE + BEDS, data=data)

###############

reg.infrisk.beds <- lm(logLOS ~ BEDS + INFRISK, data=data)
summary(reg.infrisk.beds)

data$beds100 <- data$BEDS/100
reg.infrisk.beds100 <- lm(logLOS ~ beds100 + INFRISK, data=data)
summary(reg.infrisk.beds100)


data$beds100 <- data$BEDS/100
reg.infrisk.beds <- lm(LOS ~ beds100 + INFRISK, data=data)
summary(reg.infrisk.beds)


##################################
# fitted model for logLOS ~ BEDS + INFRISK transformed to LOS scale
# holds INFRISK constant at 2 and at 5

# estimate fitted values
coefs <- reg.infrisk.beds100$coefficients
beds.range <- seq(1, 5000,1)/100
logLOShat2 <- coefs[1] + coefs[2]*beds.range + coefs[3]*2
logLOShat5 <- coefs[1] + coefs[2]*beds.range + coefs[3]*5
plot(data$beds100, data$logLOS)
lines(beds.range, logLOShat2, col=2, lwd=2)
lines(beds.range, logLOShat5, col=5, lwd=2)
legend(4,3, c("INFRISK=2","INFRISK=5"), col=c(2,5), lwd=c(2,2))


plot(data$beds100, data$LOS)
lines(beds.range, exp(logLOShat2), col=2, lwd=2)
lines(beds.range, exp(logLOShat5), col=5, lwd=2)
legend(4,20, c("INFRISK=2","INFRISK=5"), col=c(2,5), lwd=c(2,2))


plot(data$beds100, data$LOS, xlim=c(0,50), ylim=c(6.5,30))
lines(beds.range, exp(logLOShat2), col=2, lwd=2)
lines(beds.range, exp(logLOShat5), col=5, lwd=2)
legend(20,10, c("INFRISK=2","INFRISK=5"), col=c(2,5), lwd=c(2,2))