### R code from vignette source 'mob.Rnw' ################################################### ### code chunk number 1: setup ################################################### suppressWarnings(RNGversion("3.5.2")) library("partykit") options(prompt = "R> ", continue = "+ ", digits = 4, useFancyQuotes = FALSE) ################################################### ### code chunk number 2: PimaIndiansDiabetes ################################################### data("PimaIndiansDiabetes", package = "mlbench") ################################################### ### code chunk number 3: PimIndiansDiabetes-formula ################################################### pid_formula <- diabetes ~ glucose | pregnant + pressure + triceps + insulin + mass + pedigree + age ################################################### ### code chunk number 4: logit ################################################### logit <- function(y, x, start = NULL, weights = NULL, offset = NULL, ...) { glm(y ~ 0 + x, family = binomial, start = start, ...) } ################################################### ### code chunk number 5: PimaIndiansDiabetes-mob ################################################### pid_tree <- mob(pid_formula, data = PimaIndiansDiabetes, fit = logit) ################################################### ### code chunk number 6: PimaIndiansDiabetes-print ################################################### pid_tree ################################################### ### code chunk number 7: PimaIndiansDiabetes-glmtree ################################################### pid_tree2 <- glmtree(diabetes ~ glucose | pregnant + pressure + triceps + insulin + mass + pedigree + age, data = PimaIndiansDiabetes, family = binomial) ################################################### ### code chunk number 8: PimaIndiansDiabetes-plot ################################################### plot(pid_tree) ################################################### ### code chunk number 9: PimaIndiansDiabetes-plot2 ################################################### plot(pid_tree2, tp_args = list(ylines = 1, margins = c(1.5, 1.5, 1.5, 2.5))) ################################################### ### code chunk number 10: PimaIndiansDiabetes-sctest1 ################################################### library("strucchange") sctest(pid_tree, node = 1) ################################################### ### code chunk number 11: PimaIndiansDiabetes-sctest2 ################################################### sctest(pid_tree, node = 2) ################################################### ### code chunk number 12: PimaIndiansDiabetes-sctest3 ################################################### sctest(pid_tree, node = 3) ################################################### ### code chunk number 13: PimaIndiansDiabetes-prune (eval = FALSE) ################################################### ## pid_tree3 <- mob(pid_formula, data = PimaIndiansDiabetes, ## fit = logit, control = mob_control(verbose = TRUE, ## minsize = 50, maxdepth = 4, alpha = 0.9, prune = "BIC")) ################################################### ### code chunk number 14: PimaIndiansDiabetes-info ################################################### names(pid_tree$info) ################################################### ### code chunk number 15: PimaIndiansDiabetes-info ################################################### names(pid_tree$node$info) ################################################### ### code chunk number 16: PimaIndiansDiabetes-print3 ################################################### print(pid_tree, node = 3) ################################################### ### code chunk number 17: PimaIndiansDiabetes-coef ################################################### coef(pid_tree) coef(pid_tree, node = 1) ## IGNORE_RDIFF_BEGIN summary(pid_tree, node = 1) ## IGNORE_RDIFF_END ################################################### ### code chunk number 18: mob.Rnw:785-786 ################################################### exp(coef(pid_tree)[,2]) ################################################### ### code chunk number 19: mob.Rnw:789-790 ################################################### risk <- round(100 * (exp(coef(pid_tree)[,2])-1), digits = 1) ################################################### ### code chunk number 20: PimaIndiansDiabetes-logLik ################################################### logLik(pid_tree) AIC(pid_tree) BIC(pid_tree) ################################################### ### code chunk number 21: PimaIndiansDiabetes-deviance ################################################### mean(residuals(pid_tree)^2) deviance(pid_tree)/sum(weights(pid_tree)) deviance(pid_tree)/nobs(pid_tree) ################################################### ### code chunk number 22: PimaIndiansDiabetes-predict ################################################### pid <- head(PimaIndiansDiabetes) predict(pid_tree, newdata = pid, type = "node") ################################################### ### code chunk number 23: PimaIndiansDiabetes-width ################################################### width(pid_tree) depth(pid_tree) ################################################### ### code chunk number 24: PimaIndiansDiabetes-subset ################################################### pid_tree[3] ################################################### ### code chunk number 25: mob.Rnw:877-880 ################################################### predict(pid_tree2, newdata = pid, type = "node") predict(pid_tree2, newdata = pid, type = "response") predict(pid_tree2, newdata = pid, type = "link") ################################################### ### code chunk number 26: Journals-data ################################################### data("Journals", package = "AER") Journals <- transform(Journals, age = 2000 - foundingyear, chars = charpp * pages) ################################################### ### code chunk number 27: Journals-tree ################################################### j_tree <- lmtree(log(subs) ~ log(price/citations) | price + citations + age + chars + society, data = Journals, minsize = 10, verbose = TRUE) ################################################### ### code chunk number 28: Journals-plot ################################################### plot(j_tree) ################################################### ### code chunk number 29: Journals-print ################################################### j_tree ################################################### ### code chunk number 30: Journals-methods ################################################### coef(j_tree, node = 1:3) summary(j_tree, node = 1:3) sctest(j_tree, node = 1:3) ################################################### ### code chunk number 31: BostonHousing-data ################################################### data("BostonHousing", package = "mlbench") BostonHousing <- transform(BostonHousing, chas = factor(chas, levels = 0:1, labels = c("no", "yes")), rad = factor(rad, ordered = TRUE)) ################################################### ### code chunk number 32: BostonHousing-tree ################################################### bh_tree <- lmtree(medv ~ log(lstat) + I(rm^2) | zn + indus + chas + nox + age + dis + rad + tax + crim + b + ptratio, data = BostonHousing) bh_tree ################################################### ### code chunk number 33: BostonHousing-plot ################################################### plot(bh_tree) ################################################### ### code chunk number 34: BostonHousing-AIC ################################################### mean(residuals(bh_tree)^2) logLik(bh_tree) AIC(bh_tree) ################################################### ### code chunk number 35: TeachingRatings-data ################################################### data("TeachingRatings", package = "AER") tr <- subset(TeachingRatings, credits == "more") ################################################### ### code chunk number 36: TeachingRatings-lm ################################################### tr_null <- lm(eval ~ 1, data = tr, weights = students) tr_lm <- lm(eval ~ beauty + gender + minority + native + tenure + division, data = tr, weights = students) ################################################### ### code chunk number 37: TeachingRatings-tree ################################################### (tr_tree <- lmtree(eval ~ beauty | minority + age + gender + division + native + tenure, data = tr, weights = students, caseweights = FALSE)) ################################################### ### code chunk number 38: TeachingRatings-plot ################################################### plot(tr_tree) ################################################### ### code chunk number 39: TeachingRatings-coef ################################################### coef(tr_lm)[2] coef(tr_tree)[, 2] ################################################### ### code chunk number 40: TeachingRatings-rsquared ################################################### 1 - c(deviance(tr_lm), deviance(tr_tree))/deviance(tr_null) ################################################### ### code chunk number 41: Titanic-data ################################################### data("Titanic", package = "datasets") ttnc <- as.data.frame(Titanic) ttnc <- ttnc[rep(1:nrow(ttnc), ttnc$Freq), 1:4] names(ttnc)[2] <- "Gender" ttnc <- transform(ttnc, Treatment = factor( Gender == "Female" | Age == "Child", levels = c(FALSE, TRUE), labels = c("Male&Adult", "Female|Child"))) ################################################### ### code chunk number 42: Titanic-tree ################################################### ttnc_tree <- glmtree(Survived ~ Treatment | Class + Gender + Age, data = ttnc, family = binomial, alpha = 0.01) ttnc_tree ################################################### ### code chunk number 43: Titanic-plot ################################################### plot(ttnc_tree, tp_args = list(ylines = 1, margins = c(1.5, 1.5, 1.5, 2.5))) ################################################### ### code chunk number 44: GBSG2 ################################################### data("GBSG2", package = "TH.data") GBSG2$time <- GBSG2$time/365 ################################################### ### code chunk number 45: wbreg ################################################### library("survival") wbreg <- function(y, x, start = NULL, weights = NULL, offset = NULL, ...) { survreg(y ~ 0 + x, weights = weights, dist = "weibull", ...) } ################################################### ### code chunk number 46: logLik.survreg ################################################### logLik.survreg <- function(object, ...) structure(object$loglik[2], df = sum(object$df), class = "logLik") ################################################### ### code chunk number 47: gbsg2_tree ################################################### gbsg2_tree <- mob(Surv(time, cens) ~ horTh + pnodes | age + tsize + tgrade + progrec + estrec + menostat, data = GBSG2, fit = wbreg, control = mob_control(minsize = 80)) ################################################### ### code chunk number 48: GBSG2-plot ################################################### plot(gbsg2_tree) ################################################### ### code chunk number 49: GBSG2-scatter ################################################### gbsg2node <- function(mobobj, col = "black", linecol = "red", cex = 0.5, pch = NULL, jitter = FALSE, xscale = NULL, yscale = NULL, ylines = 1.5, id = TRUE, xlab = FALSE, ylab = FALSE) { ## obtain dependent variable mf <- model.frame(mobobj) y <- Formula::model.part(mobobj$info$Formula, mf, lhs = 1L, rhs = 0L) if(isTRUE(ylab)) ylab <- names(y)[1L] if(identical(ylab, FALSE)) ylab <- "" if(is.null(ylines)) ylines <- ifelse(identical(ylab, ""), 0, 2) y <- y[[1L]] ## plotting character and response if(is.null(pch)) pch <- y[,2] * 18 + 1 y <- y[,1] y <- as.numeric(y) pch <- rep(pch, length.out = length(y)) if(jitter) y <- jitter(y) ## obtain explanatory variables x <- Formula::model.part(mobobj$info$Formula, mf, lhs = 0L, rhs = 1L) xnam <- colnames(x) z <- seq(from = min(x[,2]), to = max(x[,2]), length = 51) z <- data.frame(a = rep(sort(x[,1])[c(1, NROW(x))], c(51, 51)), b = z) names(z) <- names(x) z$x <- model.matrix(~ ., data = z) ## fitted node ids fitted <- mobobj$fitted[["(fitted)"]] if(is.null(xscale)) xscale <- range(x[,2]) + c(-0.1, 0.1) * diff(range(x[,2])) if(is.null(yscale)) yscale <- range(y) + c(-0.1, 0.1) * diff(range(y)) ## panel function for scatter plots in nodes rval <- function(node) { ## node index nid <- id_node(node) ix <- fitted %in% nodeids(mobobj, from = nid, terminal = TRUE) ## dependent variable y <- y[ix] ## predictions yhat <- if(is.null(node$info$object)) { refit.modelparty(mobobj, node = nid) } else { node$info$object } yhat <- predict(yhat, newdata = z, type = "quantile", p = 0.5) pch <- pch[ix] ## viewport setup top_vp <- viewport(layout = grid.layout(nrow = 2, ncol = 3, widths = unit(c(ylines, 1, 1), c("lines", "null", "lines")), heights = unit(c(1, 1), c("lines", "null"))), width = unit(1, "npc"), height = unit(1, "npc") - unit(2, "lines"), name = paste("node_scatterplot", nid, sep = "")) pushViewport(top_vp) grid.rect(gp = gpar(fill = "white", col = 0)) ## main title top <- viewport(layout.pos.col = 2, layout.pos.row = 1) pushViewport(top) mainlab <- paste(ifelse(id, paste("Node", nid, "(n = "), ""), info_node(node)$nobs, ifelse(id, ")", ""), sep = "") grid.text(mainlab) popViewport() plot_vp <- viewport(layout.pos.col = 2, layout.pos.row = 2, xscale = xscale, yscale = yscale, name = paste("node_scatterplot", nid, "plot", sep = "")) pushViewport(plot_vp) ## scatterplot grid.points(x[ix,2], y, gp = gpar(col = col, cex = cex), pch = pch) grid.lines(z[1:51,2], yhat[1:51], default.units = "native", gp = gpar(col = linecol)) grid.lines(z[52:102,2], yhat[52:102], default.units = "native", gp = gpar(col = linecol, lty = 2)) grid.xaxis(at = c(ceiling(xscale[1]*10), floor(xscale[2]*10))/10) grid.yaxis(at = c(ceiling(yscale[1]), floor(yscale[2]))) if(isTRUE(xlab)) xlab <- xnam[2] if(!identical(xlab, FALSE)) grid.text(xlab, x = unit(0.5, "npc"), y = unit(-2, "lines")) if(!identical(ylab, FALSE)) grid.text(ylab, y = unit(0.5, "npc"), x = unit(-2, "lines"), rot = 90) grid.rect(gp = gpar(fill = "transparent")) upViewport() upViewport() } return(rval) } class(gbsg2node) <- "grapcon_generator" plot(gbsg2_tree, terminal_panel = gbsg2node, tnex = 2, tp_args = list(xscale = c(0, 52), yscale = c(-0.5, 8.7))) ################################################### ### code chunk number 50: gbsg2_tree-methods ################################################### gbsg2_tree coef(gbsg2_tree) logLik(gbsg2_tree) ################################################### ### code chunk number 51: bt-packages ################################################### data("Topmodel2007", package = "psychotree") library("psychotools") ################################################### ### code chunk number 52: btfit1 ################################################### btfit1 <- function(y, x = NULL, start = NULL, weights = NULL, offset = NULL, ...) btmodel(y, ...) ################################################### ### code chunk number 53: bt1 ################################################### bt1 <- mob(preference ~ 1 | gender + age + q1 + q2 + q3, data = Topmodel2007, fit = btfit1) ################################################### ### code chunk number 54: btfit2 ################################################### btfit2 <- function(y, x = NULL, start = NULL, weights = NULL, offset = NULL, ..., estfun = FALSE, object = FALSE) { rval <- btmodel(y, ..., estfun = estfun, vcov = object) list( coefficients = rval$coefficients, objfun = -rval$loglik, estfun = if(estfun) rval$estfun else NULL, object = if(object) rval else NULL ) } ################################################### ### code chunk number 55: bt2 ################################################### bt2 <- mob(preference ~ 1 | gender + age + q1 + q2 + q3, data = Topmodel2007, fit = btfit2) ################################################### ### code chunk number 56: bt2-print ################################################### bt2 coef(bt2) ################################################### ### code chunk number 57: bt2-worthf (eval = FALSE) ################################################### ## worthf <- function(info) paste(info$object$labels, ## format(round(worth(info$object), digits = 3)), sep = ": ") ## plot(bt2, FUN = worthf) ################################################### ### code chunk number 58: bt2-plot ################################################### plot(bt2) ################################################### ### code chunk number 59: bt2-plot2 ################################################### worthf <- function(info) paste(info$object$labels, format(round(worth(info$object), digits = 3)), sep = ": ") plot(bt2, FUN = worthf) ################################################### ### code chunk number 60: bt2-nodeapply (eval = FALSE) ################################################### ## par(mfrow = c(2, 2)) ## nodeapply(bt2, ids = c(3, 5, 6, 7), FUN = function(n) ## plot(n$info$object, main = n$id, ylim = c(0, 0.4))) ################################################### ### code chunk number 61: bt2-nodeapply-plot ################################################### par(mfrow = c(2, 2)) nodeapply(bt2, ids = c(3, 5, 6, 7), FUN = function(n) plot(n$info$object, main = n$id, ylim = c(0, 0.4))) ################################################### ### code chunk number 62: bt2-plot3 (eval = FALSE) ################################################### ## plot(bt2, drop = TRUE, tnex = 2, ## terminal_panel = node_btplot(bt2, abbreviate = 1, yscale = c(0, 0.5))) ################################################### ### code chunk number 63: node_btplot ################################################### ## visualization function node_btplot <- function(mobobj, id = TRUE, worth = TRUE, names = TRUE, abbreviate = TRUE, index = TRUE, ref = TRUE, col = "black", linecol = "lightgray", cex = 0.5, pch = 19, xscale = NULL, yscale = NULL, ylines = 1.5) { ## node ids node <- nodeids(mobobj, terminal = FALSE) ## get all coefficients cf <- partykit:::apply_to_models(mobobj, node, FUN = function(z) if(worth) worth(z) else coef(z, all = FALSE, ref = TRUE)) cf <- do.call("rbind", cf) rownames(cf) <- node ## get one full model mod <- partykit:::apply_to_models(mobobj, node = 1L, FUN = NULL) if(!worth) { if(is.character(ref) | is.numeric(ref)) { reflab <- ref ref <- TRUE } else { reflab <- mod$ref } if(is.character(reflab)) reflab <- match(reflab, mod$labels) cf <- cf - cf[,reflab] } ## reference if(worth) { cf_ref <- 1/ncol(cf) } else { cf_ref <- 0 } ## labeling if(is.character(names)) { colnames(cf) <- names names <- TRUE } ## abbreviation if(is.logical(abbreviate)) { nlab <- max(nchar(colnames(cf))) abbreviate <- if(abbreviate) as.numeric(cut(nlab, c(-Inf, 1.5, 4.5, 7.5, Inf))) else nlab } colnames(cf) <- abbreviate(colnames(cf), abbreviate) if(index) { x <- 1:NCOL(cf) if(is.null(xscale)) xscale <- range(x) + c(-0.1, 0.1) * diff(range(x)) } else { x <- rep(0, length(cf)) if(is.null(xscale)) xscale <- c(-1, 1) } if(is.null(yscale)) yscale <- range(cf) + c(-0.1, 0.1) * diff(range(cf)) ## panel function for bt plots in nodes rval <- function(node) { ## node index id <- id_node(node) ## dependent variable setup cfi <- cf[id,] ## viewport setup top_vp <- viewport(layout = grid.layout(nrow = 2, ncol = 3, widths = unit(c(ylines, 1, 1), c("lines", "null", "lines")), heights = unit(c(1, 1), c("lines", "null"))), width = unit(1, "npc"), height = unit(1, "npc") - unit(2, "lines"), name = paste("node_btplot", id, sep = "")) pushViewport(top_vp) grid.rect(gp = gpar(fill = "white", col = 0)) ## main title top <- viewport(layout.pos.col = 2, layout.pos.row = 1) pushViewport(top) mainlab <- paste(ifelse(id, paste("Node", id, "(n = "), ""), info_node(node)$nobs, ifelse(id, ")", ""), sep = "") grid.text(mainlab) popViewport() ## actual plot plot_vpi <- viewport(layout.pos.col = 2, layout.pos.row = 2, xscale = xscale, yscale = yscale, name = paste("node_btplot", id, "plot", sep = "")) pushViewport(plot_vpi) grid.lines(xscale, c(cf_ref, cf_ref), gp = gpar(col = linecol), default.units = "native") if(index) { grid.lines(x, cfi, gp = gpar(col = col, lty = 2), default.units = "native") grid.points(x, cfi, gp = gpar(col = col, cex = cex), pch = pch, default.units = "native") grid.xaxis(at = x, label = if(names) names(cfi) else x) } else { if(names) grid.text(names(cfi), x = x, y = cfi, default.units = "native") else grid.points(x, cfi, gp = gpar(col = col, cex = cex), pch = pch, default.units = "native") } grid.yaxis(at = c(ceiling(yscale[1] * 100)/100, floor(yscale[2] * 100)/100)) grid.rect(gp = gpar(fill = "transparent")) upViewport(2) } return(rval) } class(node_btplot) <- "grapcon_generator" plot(bt2, drop = TRUE, tnex = 2, terminal_panel = node_btplot(bt2, abbreviate = 1, yscale = c(0, 0.5))) ################################################### ### code chunk number 64: tm ################################################### tm <- data.frame(age = c(60, 25, 35), gender = c("male", "female", "female"), q1 = "no", q2 = c("no", "no", "yes"), q3 = "no") tm ################################################### ### code chunk number 65: tm-predict ################################################### tm predict(bt2, tm, type = "node") predict(bt2, tm, type = function(object) t(worth(object))) predict(bt2, tm, type = function(object) t(rank(-worth(object))))