data("BostonHousing2", package = "mlbench")
### fit non-normal Box-Cox type linear model with two
### baseline functions (for houses near and off Charles River)
BC_BH_2 <- BoxCox(cmedv | 0 + chas ~ crim + zn + indus + nox +
rm + age + dis + rad + tax + ptratio + b + lstat,
data = BostonHousing2)
logLik(BC_BH_2)
### classical likelihood inference
summary(BC_BH_2)
### coefficients of the linear predictor
coef(BC_BH_2)
### plot linear predictor (mean of _transformed_ response)
### vs. observed values
plot(predict(BC_BH_2, type = "lp"), BostonHousing2$cmedv)
### all coefficients
coef(BC_BH_2, with_baseline = TRUE)
### compute predicted median along with 10% and 90% quantile for the first
### observations
predict(BC_BH_2, newdata = BostonHousing2[1:3,], type = "quantile",
prob = c(.1, .5, .9))
### plot the predicted density for these observations
plot(BC_BH_2, newdata = BostonHousing2[1:3, -1],
which = "distribution", type = "density", K = 1000)
### evaluate the two baseline transformations, with confidence intervals
nd <- model.frame(BC_BH_2)[1:2, -1]
nd$chas <- factor(c("0", "1"))
library("colorspace")
col <- diverge_hcl(2, h = c(246, 40), c = 96, l = c(65, 90))
fill <- diverge_hcl(2, h = c(246, 40), c = 96, l = c(65, 90), alpha = .3)
plot(BC_BH_2, which = "baseline only", newdata = nd, col = col,
confidence = "interval", fill = fill, lwd = 2,
xlab = "Median Value", ylab = expression(h[Y]))
legend("bottomright", lty = 1, col = col,
title = "Near Charles River", legend = c("no", "yes"), bty = "n")
### cars data; with quantile functions
plot(dist ~ speed, data = cars)
m <- Colr(dist ~ speed, data = cars)
q <- predict(as.mlt(m), newdata = data.frame(speed = s <- 6:25),
type = "quantile", prob = c(1, 5, 9) / 10)
lines(s, q[1,])
lines(s, q[2,])
lines(s, q[3,])
nd <- data.frame(speed = s <- as.double(1:5 * 5))
# Prob(dist at speed s > dist at speed 0)
# speed 0 is reference, not a good choice here
PI(m, newdata = nd)
# Prob(dist at speed s > dist at speed 15)
lp15 <- c(predict(m, newdata = data.frame(speed = 15)))
PI(m, newdata = nd, reference = lp15)
PI(m, newdata = nd, reference = nd[3,,drop = FALSE])
# Prob(dist at speed s' > dist at speed s)
PI(m, newdata = nd, reference = nd)
# essentially:
lp <- predict(m, newdata = nd)
PI(object = dist(lp))
# same, with simultaneous confidence intervals
PI(m, newdata = nd, reference = nd, conf.level = .95)
# plot ROC curves + confidence bands
# compare speed 20 and 25 to speed 15
plot(ROC(m, newdata = nd[4:5,,drop = FALSE],
reference = nd[3,,drop = FALSE],
conf.level = 0.95))
# Overlap of conditional densities at speed s' and s
OVL(m, newdata = nd, reference = nd)
### ROC analysis (takes too long for CRAN Windows)
if (require("mlbench") && .Platform$OS.type != "windows") {
layout(matrix(1:4, nrow = 2))
data("PimaIndiansDiabetes2", package = "mlbench")
dia <- sort(unique(PimaIndiansDiabetes2$diabetes))
nd <- data.frame(diabetes = dia,
age = 29, mass = 32) ### median values
### unconditional ROC analysis: glucose tolerance test
m0 <- Colr(glucose ~ diabetes, data = PimaIndiansDiabetes2)
# ROC curve + confidence band
plot(ROC(m0, newdata = nd[2,,drop = FALSE], conf.level = .95))
# Wald interval for AUC
PI(m0, newdata = nd[2,,drop = FALSE], conf.level = .95)
# score interval for AUC
PI(-c(coef(m0), score_test(m0)$conf.int[2:1]))
### adjusted ROC analysis for age and mass
m1 <- Colr(glucose ~ diabetes + age + mass, data = PimaIndiansDiabetes2)
# ROC curve + confidence band (this is the same for all ages /
# masses)
plot(ROC(m1, newdata = nd[2,,drop = FALSE],
reference = nd[1,,drop = FALSE],
conf.level = .95))
# Wald interval for adjusted AUC
PI(m1, newdata = nd[2,,drop = FALSE], reference = nd[1,,drop = FALSE],
conf.level = .95)
# Score interval for adjusted AUC
PI(-c(coef(m1)[1], score_test(m1, names(coef(m1))[1])$conf.int[2:1]))
### conditional ROC analysis: AUC regression ~ age + mass
m2 <- Colr(glucose ~ diabetes * (age + mass), data = PimaIndiansDiabetes2)
# ROC curve for a person with age = 29 and mass = 32
plot(ROC(m2, newdata = nd[2,,drop = FALSE],
reference = nd[1,,drop = FALSE],
conf.level = .95))
# AUC for persons ages 21:81, all with mass = 32
nd1 <- data.frame(diabetes = nd[1,"diabetes"], age = 21:81, mass = 32)
nd2 <- data.frame(diabetes = nd[2,"diabetes"], age = 21:81, mass = 32)
auc <- PI(m2, newdata = nd2, reference = nd1, one2one = TRUE,
conf.level = 0.95)
plot(nd1$age, auc[, "Estimate"], xlab = "Age (in years)", ylab =
"AUC", ylim = c(0, 1), type = "l")
lines(nd1$age, auc[, "lwr"], lty = 3)
lines(nd1$age, auc[, "upr"], lty = 3)
}
Run the code above in your browser using DataLab