#########################################
# Day 1 of Regression Modeling Strategies
#########################################

# Assignment 3

library(Hmisc,T)
library(Design,T)
data.restore('/tmp/support.sdd')
support <- cleanup.import(support)
page(contents(support),multi=T)
page(describe(support),multi=T)
library(rpart)
datadensity(support)
hist.data.frame(support)
ecdf(support, datadensity='hist',side=3)
n <- naclus(support)
plot(n)
naplot(n)
sub <- subset(support, is.na(totcst) | totcst>0)
dim(sub)
w <- is.na(totcst) ~ age + sex + 
 dzgroup + scoma + meanbp + hospdead + charges
s <- summary(w, data=sub)
options(digits=3)
s
page(s)
plot(s)
f <- rpart(w, data=sub)
plot(f);text(f)
f$cptable
sum(sub$age > 24.53 & sub$age < 26.2)

# Estimate costs from charges (previous regression)
sub <- upData(sub, 
 totcst.i=ifelse(is.na(totcst),
                 exp(.2550+.9244*log(charges)),
                 totcst),
 labels=c(totcst.i='Total Cost Imputed from charges')
   )
page(contents(sub),multi=T)

# Assignment 4 Prob. 1
g <- function(y) c(Mean=mean(y), Median=median(y))
s <- summary(totcst.i ~ age+sex+dzgroup+num.co+scoma+
   race+meanbp, data=sub, fun=g)
page(s)
plot(s, which=1:2, pch=1:2, xlab='$')
# Use show.pch() to see pch defns

w <- spearman2(totcst.i ~ age+sex+dzgroup+num.co+scoma+
        race+meanbp,data=sub,p=2)
w
plot(w)
tapply(sub$totcst.i, sub$dzgroup, median, na.rm=T)


##########################
# Day 2
##########################
library(Hmisc,T)
library(Design,T)

sample(c(10,11,13,17), 4, T)
x <- 1:20
smean.cl.boot(x)   # in Hmisc library
smean.cl.normal(x)

# Assignment 5: RMS Chp 7 Prob 2(a)
v <- varclus(~., data=sub)
plot(v)
# (b)
s <- spearman2(totcst.i ~ age+sex+dzgroup+
 num.co+scoma+race+meanbp+hrt+temp+pafi+alb,
 data=sub)
plot(s)
# (c) done on Day 1
# (d)
library(rpart)
f <- rpart(is.na(pafi) ~ ., data=sub)
plot(f);text(f)
levels(sub$dzgroup)
print(paste(letters[1:8],
   levels(sub$dzgroup),sep=':'),quote=F)
describe(sub$race)

# (2)e-f
sub$race2 <- combine.levels(as.character(
   impute(sub$race)), minlev=.1)
sub$pafi <- oldUnclass(impute(sub$pafi, 333.33))
sub$alb <- oldUnclass(impute(sub$alb, 3.5))
describe(sub[c('race2','pafi','alb')])
class(sub$pafi)
attributes(sub$pafi)

dd <- datadist(sub)
options(datadist='dd')

f <- ols(totcst ~ rcs(age,3)+sex+dzgroup+
 pol(num.co,2)+pol(scoma,2)+race2+rcs(meanbp,4)+
 rcs(hrt,4)+rcs(temp,5)+rcs(alb,5)+rcs(pafi,3),
 data=sub)
qqnorm(resid(f)) # terrible non-normality!
f <- update(f, log(totcst) ~ .)
r <- resid(f)
qqnorm(r) # almost a straight line
qqline(r)
bwplot(dzgroup ~ r, data=sub)
xYplot(r ~ meanbp, method='quantile',
 data=sub,nx=100)
table(sub$num.co);table(sub$scoma)
a <- anova(f)
page(a)
plot(a,what='partial R2')
par(mfrow=c(3,4))
plot(f, ref.zero=T, ylim=c(-1.5,.7))

par(mfrow=c(1,1))
s <- summary(f, antilog=T)
plot(s, main='Cost Ratio', log=T)
nomogram(f, 
 fun=function(x)exp(x)/1000,
 funlabel='Median Cost ($/1000)',
 cex.var=.6, cex.axis=.5)

Dialog(fitPar('f', lp=F,   
 fun=list('Median Cost ($/1000)'=function(y)exp(y)/1000),  
  fun.round=0))
 runmenu.f()

f <- update(f, x=T, y=T)
f$x[1:5,]
v <- validate(f, B=50)
page(v)
page(f)
cal <- calibrate(f, B=50)
plot(cal)
g <- Function(f)
page(g)
exp(g(age=c(30,50,70, 90),sex='female'))


#################################
# Day 3
#################################

library(Hmisc,T)
library(Design,T)
?Design.trans
?scored
?Overview

data.restore('/tmp/titanic3.sdd')
titanic3 <- cleanup.import(titanic3)
page(describe(titanic3),multi=T)
dd <- datadist(titanic3)
options(datadist='dd')
options(digits=2)
attach(titanic3)
rm(sex)
s <- summary(survived ~ age+sex+pclass+sibsp+parch)
plot(s)
plsmo(age, survived, datadensity=T)
plsmo(age, survived, group=sex,
      col=1:2, datadensity=T)
plsmo(age, survived, group=pclass,
      col=1:3, datadensity=T)
plsmo(age, survived, 
      group=interaction(pclass,sex),,
      col=1:6, datadensity=T)
plsmo(age, survived,
      group=interaction(pclass,
           cut2(fare,g=3)),
      datadensity=T, col=1:9)
v <- varclus(~age+sex+pclass+fare)
plot(v)
bwplot(pclass ~ fare,
 panel=panel.bpplot)

w <- pclass=='1st'
table(w)
plsmo(age[w], survived[w],
  group=cut2(fare[w],g=4))
plsmo(fare[w], survived[w],
  group=cut2(age[w],g=4),datadensity=T)
ecdf(fare)
library(rpart)
r <- rpart(fare ~ sex+age+pclass)
plot(r);text(r)
r <- rpart(survived ~ sex+age+
  pclass+fare)
plot(r);text(r)

f1 <- lrm(survived ~ sex*pclass*
 rcs(age,5)+rcs(age,5)*(sibsp+parch))
page(anova(f1))
f <- lrm(survived ~ 
 (sex + pclass + rcs(age,5))^2 +
 rcs(age,5)*sibsp)
page(f)
a <- anova(f)
page(a); plot(a)
print(a, 'sub')
for(sx in levels(sex)) {
	plot(f, age=NA, pclass=NA, sex=sx)
	title(sx)
}

plot(f) # somewhat arbitrary b/c interactions
f <- update(f, x=T,y=T)
v <- validate(f, B=30)
page(v)
cal <- calibrate(f, B=10)
plot(cal)
f$stats
(554-26)/554

who.na <- rpart(is.na(age) ~
 sex+pclass+survived+sibsp+parch,
 minbucket=15)
plot(who.na);text(who.na)
m <- lrm(is.na(age) ~ sex*pclass +
 survived + sibsp + parch + 
 rcs(fare,4))
plot(anova(m))
plot(m, fare=NA)
xYplot(age ~ fare | pclass*sex)
spearman2(fare ~ age)

set.seed(17)
xtrans <- aregImpute(
 ~ age + sex + pclass + sibsp +
   parch + survived, n.impute=5)
page(xtrans)
plot(xtrans)
traceback()
?aregImpute
.R.; .SV4.
f.mi <- fit.mult.impute(formula(f),
 lrm, xtrans, data=titanic3)
page(anova(f.mi))
plot(f.mi, age=NA, pclass=NA)
s <- summary(f.mi, age=c(1,30),
 sibsp=0:1)
plot(s,log=T)
g <- Function(f.mi)
page(g)
page(sascode(g))
plogis(g(age=c(2,21,50),sex='male',
 pclass='3rd'))
drep <- dataRep(~roundN(age,10)+
 sex+pclass+roundN(sibsp,clip=0:1))
page(drep)
Dialog(fitPar('f.mi',lp=F,
 fun=list('Prob[Survival]'=plogis)),
 limits='data',
 basename='Titanic',
 vary=list(sex=levels(sex)),
 datarep=drep)
runmenu.Titanic()


data.restore('/tmp/diabetes.sdd')
diabetes <- cleanup.import(diabetes)
plot(v <- varclus(~.,data=diabetes))
search()[1:5]
page(v)
detach(2)
attach(diabetes)
r <- rpart(frame ~ gender+height+
 weight+waist+hip)
probs <- predict(r, diabetes)
mpc <- (t(apply(-probs,1,order)))[,1]
frame.pred <- levels(frame)[mpc]
table(frame,frame.pred)
fr <- ifelse(is.na(frame), 
 as.character(frame.pred),
 as.character(frame))
table(fr)
frame <- fr
na <- naclus(diabetes)
naplot(na)

bp.1s <- impute(bp.1s)
chol <- impute(chol)
weight <- impute(weight)
hip <- impute(hip)
set.seed(721)
f <- areg.boot(glyhb ~ 
 monotone(age) + bp.1s + chol +
 frame + monotone(weight) +
 monotone(hip),B=60)
page(f)
plot(f)
plot(fitted(f), resid(f))
qqnorm(resid(f))
page(f)
funs <- Function(f)
names(funs)
page(funs$glyhb)
plot(glyhb, funs$glyhb(glyhb))
plot(1/glyhb, funs$glyhb(glyhb))
s <- summary(f, values=list(
 age=c(20,30,40,50,60,70,80)))