library(rpart)
?rpart
.First <- function(...) {
	library(Hmisc,T)
	library(Design,T)
	invisible()
}
.First()
search()
?getHdata
getHdata()
getHdata(titanic3)
page(describe(titanic3), multi=T)

na.patterns <- naclus(titanic3)
plot(na.patterns)
naplot(na.patterns)

who.na <- rpart(is.na(age) ~ sex +
 pclass + survived, data=titanic3)
plot(who.na); text(who.na)
s <- summary(is.na(age) ~
 sex + pclass + survived,
 data=titanic3)
plot(s)
m <- lrm(is.na(age) ~ sex*pclass+
 survived, data=titanic3)
page(anova(m),multi=T)
plot(anova(m))

# Day 2
set.seed(1)
w <- aregImpute(~age+sex+
        pclass+survived,
        data=titanic3,
        n.impute=5)
names(w)
w
summary(w)
plot(w)
ecdf(titanic3$age, add=T, col=3)
page(w$imputed)

# Dataset ed contributed by
# Matt Calder
dim(ed)
names(ed)
# create export file
data.dump('ed','/tmp/ed.sdd',
 oldStyle=T)
# import data:
data.restore('/winnt/temp/ed.sdd')
plot(varclus(~x1+x2+x3+x4+x5, 
     data=ed))
page(describe(ed),multi=T)
datadensity(ed)

z <- transcan(~x1+x2+x3+x4+x5,
 data=ed)
?areg.boot
args(areg.boot)
z <- areg.boot(y ~ x1+x2+x3+x4+x5,
 data=ed, B=10)
z
plot(z)

s <- spearman2(y ~ x1+x2+x3+x4+x5, 
          p=2, data=ed)
plot(s)
f <- ols(y ~ rcs(x1,4) +
   rcs(x2,7) + rcs(x3,5) +
   rcs(x4,5) + rcs(x5,4),data=ed)
f$stats
page(f)
Function(f)
dd <- datadist(ed)
page(dd)
options(datadist='dd')
plot(f)
plot(f, x2=NA)
scat1d(ed$x2)
options(width=90)
page(anova(f))
plot(anova(f))
nomogram(f)

f <- ols(y ~ rcs(x1,4) +
   rcs(x2,7) + rcs(x3,5) +
   rcs(x4,5) + rcs(x5,4) +
   rcs(x1,4) %ia% rcs(x5,4) +
   rcs(x4,5) %ia% rcs(x5,4),
   data=ed)
f$stats
plot(anova(f))
i <- plot(f, x4=NA, x5=NA, 
 method='image')
Legend(i)
anova(f,x1,x2)

getHdata(support)
attach(subset(support,
 !is.na(totcst) & totcst>0))
length(age)
describe(race)
describe(pafi)
pafi <- impute(pafi, 333.3)
describe(pafi)
alb <- impute(alb, 3.5)
race <- combine.levels(
 race, minlev=.1)
table(race)
dd <- datadist(support)
dd <- datadist(dd, race)
options(datadist='dd')
s <- spearman2(totcst ~
 age+sex+dzgroup+num.co+scoma+
 race+meanbp+hrt+temp+pafi+alb,
 p=2)
plot(s)
f <- ols(log(totcst) ~ rcs(age,3) +
 sex+dzgroup+pol(num.co,2)+
 pol(scoma,2)+race+rcs(meanbp,4)+
 rcs(hrt,4)+rcs(temp,5)+
 rcs(alb,5)+rcs(pafi,3))
f$stats
page(f)
an <- anova(f)
options(digits=3)
print(an,'dots')
plot(an, what='partial R2')
qqnorm(resid(f));qqline(resid(f))
bwplot(dzgroup ~ resid(f))

# Day 3
r <- resid(f)
length(r)
xYplot(r ~ age | dzgroup,
 method='quantile',
 abline=list(h=0,lty=2))
plot(f, ref.zero=T, 
 ylim=c(-1.5,.7))
args(plot.Design)
plot(f, dzgroup=NA,
 method='dotchart', 
 sortdot='descend',
 ref.zero=T)
s <- summary(f, antilog=T)
plot(s, main='Cost Ratio',
     log=T, cex=.75)
specs(f)

s <- summary(f, age=c(50,70),
       antilog=T)
# Did not accept dzgroup='Coma'
# It should have!
plot(s, main='Cost Ratio',
     log=T, cex=.75)

nomogram(f, 
 fun=function(x)exp(x)/1000,
 funlabel='Median Cost $/1000',
 cex.var=.6, cex.axis=.5)
f <- update(f, x=T, y=T)
page(f$x[10,])
v <- validate(f, B=30)
page(v)
cal <- calibrate(f, B=50)
plot(cal)

# Go back to Matt's data
g <- ols(y ~ rcs(x1,4) +
   rcs(x2,7) + rcs(x3,5) +
   rcs(x4,5) + rcs(x5,4),
 x=T,y=T,data=ed)
v <- validate(g, B=80)
g$stats
g
page(v)
cal <- calibrate(g, B=50)
plot(cal)

search()
detach(2)
dim(titanic3)
s <- summary(survived ~
 age + sex + pclass, data=titanic3)
plot(s)
attach(titanic3)
dev.off()
plsmo(age, survived,
 datadensity=T)
plsmo(age, survived, group=sex,
 datadensity=T, col=1:2)
plsmo(age, survived, group=pclass,
 datadensity=T, col=1:3)
plsmo(age, survived,
 group=interaction(pclass,sex),
 datadensity=T, col=1:6)
bwplot(pclass ~ age | sex,
 panel=panel.bpplot)

dd <- datadist(titanic3)
options(datadist='dd')
f1 <- lrm(survived ~ sex*pclass*
 rcs(age,5))
f1$stats
page(anova(f1))
f <- lrm(survived ~
 (sex+pclass+rcs(age,5))^2)
f$stats
a <- anova(f)
page(a)
plot(a)
plot(f, age=NA, sex=NA,
 conf.int=F, fun=plogis)
par(mfrow=c(2,1))
for(sx in levels(sex)) {
 plot(f, age=NA, pclass=NA,
  sex=sx, conf.int=F, fun=plogis)
 title(sx)
}
par(mfrow=c(1,1))
f <- update(f, x=T, y=T)
set.seed(131)
v <- validate(f, B=80)
page(v)
cal <- calibrate(f,B=20)
plot(cal)