R/DAMisc_functions.R

In DAMisc: Dave Armstrong's Miscellaneous Functions

binfit <-
function (mod) 
{
    mod <- update(mod, x = TRUE, y = TRUE)
    y <- mod$y
    null.mod <- update(mod, ".~1", data=model.frame(mod))
    b <- mod$coef[-1]
    var.ystar <- t(b) %*% var(model.matrix(mod)[, -1]) %*% b
    G <- -2 * (logLik(null.mod) - logLik(mod))
    res.col1 <- c(logLik(null.mod), deviance(mod), NA, 1 - (logLik(mod)/logLik(null.mod)), 
        1 - exp(2*(logLik(null.mod) - logLik(mod))/length(mod$residuals)), 
        var.ystar/(var.ystar + switch(mod$family[[2]], logit = pi^2/3, 
            probit = 1)), mean(mod$y == as.numeric(fitted(mod) > 
            0.5)), BIC(mod))
    res.col2 <- c(logLik(mod), G, pchisq(G, 8, lower.tail = F), 
        1 - ((logLik(mod) - mod$rank)/logLik(null.mod)), res.col1[5]/(1 - 
            (exp(2*logLik(null.mod)/length(mod[["residuals"]])))), 
        1 - (sum((y - fitted(mod))^2)/sum((y - mean(y))^2)), 
        (sum(mod$y == as.numeric(fitted(mod) > 0.5)) - max(table(y)))/(length(mod$residuals) - 
            max(table(y))), AIC(mod))
    res.vec <- c(res.col1, res.col2)[-3]
    res.col1 <- sprintf("%3.3f", res.col1)
    res.col1[3] <- ""
    res.df <- data.frame(Names1 = c("Log-Lik Intercept Only:", 
        paste("D(", mod$df.residual, "):", sep = ""), " ", "McFadden's R2:", 
        "ML (Cox-Snell) R2:", "McKelvey & Zavoina R2:", "Count R2:", 
        "BIC:"), vals1 = res.col1, Names2 = c("Log-Lik Full Model:", 
        paste("LR(", mod$rank - 1, "):", sep = ""), "Prob > LR:", 
        "McFadden's Adk R2:", "Cragg-Uhler (Nagelkerke) R2:", 
        "Efron's R2:", "Adj Count R2:", "AIC:"), vals2 = sprintf("%3.3f", 
        res.col2))
    return(res.df)
}
btscs <-
function (data, event, tvar, csunit, pad.ts = FALSE) 
{
    data$orig_order <- 1:nrow(data)
    data <- data[order(data[[csunit]], data[[tvar]]), ]
    spells <- function(x) {
        tmp <- rep(0, length(x))
        runcount <- 0
        for (j in 2:length(x)) {
            if (x[j] == 0 & x[(j - 1)] == 0) {
                tmp[j] <- runcount <- runcount + 1
            }
            if (x[j] != 0 & x[(j - 1)] == 0) {
                tmp[j] <- runcount + 1
                runcount <- 0
            }
            if (x[j] == 0 & x[(j - 1)] != 0) {
                tmp[j] <- runcount <- 0
            }
        }
        tmp
    }
    sp <- split(data, data[[csunit]])
    if (pad.ts) {
        sp <- lapply(sp, function(x) x[match(seq(min(x[[tvar]], 
            na.rm = T), max(x[[tvar]], na.rm = T)), x[[tvar]]), 
            ])
        for (i in 1:length(sp)) {
            if (any(is.na(sp[[i]][[event]]))) {
                sp[[i]][[event]][which(is.na(sp[[i]][[event]]))] <- 1
            }
            if (any(is.na(sp[[i]][[tvar]]))) {
                sp[[i]][[tvar]] <- seq(min(sp[[i]][[tvar]], na.rm = T), 
                  max(sp[[i]][[tvar]], na.rm = T))
            }
            if (any(is.na(sp[[i]][[csunit]]))) {
                sp[[i]][[csunit]][which(is.na(sp[[i]][[csunit]]))] <- mean(sp[[i]][[csunit]], 
                  na.rm = T)
            }
        }
    }
    sp <- lapply(1:length(sp), function(x) {
        cbind(sp[[x]], data.frame(spell = spells(sp[[x]][[event]])))
    })
    data <- do.call(rbind, sp)
    if (!pad.ts) {
        if (any(is.na(data$orig_order))) {
            data <- data[-which(is.na(data$orig_order)), ]
        }
        data <- data[data$orig_order, ]
    }
    else {
        data <- data[order(data[[csunit]], data[[tvar]]), ]
    }
    invisible(data)
}
DAintfun <-
function (obj, varnames, theta = 45, phi = 10, xlab=NULL, ylab=NULL, zlab=NULL,...) 
{
    if (length(varnames) != 2) {
        stop("varnames must be a vector of 2 variable names")
    }
    else {
        v1 <- varnames[1]
        v2 <- varnames[2]
        ind <- unique(c(grep(v1, names(obj$coef)), grep(v2, names(obj$coef))))
        ind <- ind[order(ind)]
        b <- obj$coef[ind]
        mod.x <- model.matrix(obj)
        not.ind <- c(1:ncol(mod.x))[!(c(1:ncol(mod.x)) %in% ind)]
        mod.x[, not.ind] <- 0
        dens <- sm.density(mod.x[, varnames], display = "none")
        b <- obj$coef[ind]
        v1.seq <- dens$eval.points[, 1]
        v2.seq <- dens$eval.points[, 2]
        eff.fun <- function(x1, x2) {
            b[1] * x1 + b[2] * x2 + b[3] * x1 * x2
        }
        hcols <- paste("gray", seq(from = 20, to = 80, length = 4), 
            sep = "")
        predsurf <- outer(v1.seq, v2.seq, eff.fun)
        cutoff <- quantile(dens$estimate, prob = c(0.25, 0.5, 
            0.75))
        pred1 <- predsurf
        pred1[dens$estimate < cutoff[1]] <- NA
        pred2 <- predsurf
        pred2[dens$estimate < cutoff[2]] <- NA
        pred3 <- predsurf
        pred3[dens$estimate < cutoff[3]] <- NA
        persp(v1.seq, v2.seq, predsurf, 
			xlab = ifelse(is.null(xlab), toupper(v1), xlab), 
			ylab = ifelse(is.null(ylab), toupper(v2), ylab), 
            zlab = ifelse(is.null(zlab), toupper("Predictions"), zlab), 
			col = hcols[1], theta = theta, phi = phi,...)
        par(new = TRUE)
        persp(v1.seq, v2.seq, pred1, col = hcols[2], axes = FALSE, 
            xlab = "", ylab = "", zlab = "", theta = theta, phi = phi, 
            zlim = c(min(c(predsurf)), max(c(predsurf))), ylim = c(min(v2.seq), 
                max(v2.seq)), xlim = c(min(v1.seq), max(v1.seq)))
        par(new = TRUE)
        persp(v1.seq, v2.seq, pred2, col = hcols[3], axes = FALSE, 
            xlab = "", ylab = "", zlab = "", theta = theta, phi = phi, 
            zlim = c(min(c(predsurf)), max(c(predsurf))), ylim = c(min(v2.seq), 
                max(v2.seq)), xlim = c(min(v1.seq), max(v1.seq)))
        par(new = TRUE)
        persp(v1.seq, v2.seq, pred3, col = hcols[4], axes = FALSE, 
            xlab = "", ylab = "", zlab = "", theta = theta, phi = phi, 
            zlim = c(min(c(predsurf)), max(c(predsurf))), ylim = c(min(v2.seq), 
                max(v2.seq)), xlim = c(min(v1.seq), max(v1.seq)))
	    invisible(list(x1=v1.seq, x2=v2.seq, pred=predsurf))
    }
}
DAintfun2 <-
function (obj, varnames, rug = TRUE, ticksize = -0.03, hist = FALSE, 
    hist.col = "gray75", nclass = c(10, 10), scale.hist = 0.5, 
    border = NA, name.stem = "cond_eff", 
	xlab = NULL, ylab=NULL, plot.type = "screen") 
{
    rseq <- function(x) {
        rx <- range(x, na.rm = TRUE)
        seq(rx[1], rx[2], length = 25)
    }
    if (!("model" %in% names(obj))) {
        obj <- update(obj, model = T)
    }
    v1 <- varnames[1]
    v2 <- varnames[2]
    ind1 <- grep(v1, names(obj$coef))
    ind2 <- grep(v2, names(obj$coef))
    s1 <- rseq(model.matrix(obj)[, v1])
    s2 <- rseq(model.matrix(obj)[, v2])
    a1 <- a2 <- matrix(0, nrow = 25, ncol = ncol(model.matrix(obj)))
    a1[, ind1[1]] <- 1
    a1[, ind1[2]] <- s2
    a2[, ind2[1]] <- 1
    a2[, ind2[2]] <- s1
    eff1 <- a1 %*% obj$coef
    se.eff1 <- sqrt(diag(a1 %*% vcov(obj) %*% t(a1)))
    low1 <- eff1 - qt(0.975, obj$df.residual) * se.eff1
    up1 <- eff1 + qt(0.975, obj$df.residual) * se.eff1
    eff2 <- a2 %*% obj$coef
    se.eff2 <- sqrt(diag(a2 %*% vcov(obj) %*% t(a2)))
    low2 <- eff2 - qt(0.975, obj$df.residual) * se.eff2
    up2 <- eff2 + qt(0.975, obj$df.residual) * se.eff2
    if (!plot.type %in% c("pdf", "png", "eps", "screen")) {
        print("plot type must be one of - pdf, png or eps")
    }
    else {
        if (plot.type == "pdf") {
            pdf(paste(name.stem, "_", v1, ".pdf", sep = ""), 
                height = 6, width = 6)
        }
        if (plot.type == "png") {
            png(paste(name.stem, "_", v1, ".png", sep = ""))
        }
        if (plot.type == "eps") {
            old.psopts <- ps.options()
            setEPS()
            postscript(paste(name.stem, "_", v1, ".eps", sep = ""))
        }
        if (plot.type == "screen") {
            oldpar <- par()
            par(mfrow = c(1, 2))
        }
        plot(s2, eff1, type = "n", ylim = range(c(low1, up1)), 
            xlab = ifelse(is.null(xlab), toupper(v2), xlab[1]), ylab = ifelse(is.null(ylab), paste("Conditional Effect of ", 
                toupper(v1), " | ", toupper(v2), sep = ""), ylab[1]))
        if (hist == TRUE) {
            rng <- diff(par()$usr[3:4])
            h2 <- hist(obj$model[[v2]], nclass = nclass[1], plot = FALSE)
            prop2 <- h2$counts/sum(h2$counts)
            plot.prop2 <- (prop2/max(prop2)) * rng * scale.hist + 
                par()$usr[3]
            av2 <- pretty(prop2, n = 3)
            axis(4, at = (av2/max(prop2)) * rng * scale.hist + 
                par()$usr[3], labels = av2)
            br2 <- h2$breaks
            for (i in 1:(length(br2) - 1)) {
                polygon(x = c(br2[i], br2[(i + 1)], br2[(i + 
                  1)], br2[i], br2[i]), y = c(par()$usr[3], par()$usr[3], 
                  plot.prop2[i], plot.prop2[i], par()$usr[3]), 
                  col = hist.col, border = border)
            }
        }
        if (rug == TRUE) {
            rug(obj$model[[v2]], ticksize = ticksize)
        }
        if (par()$usr[3] < 0 & par()$usr[4] > 0) {
            abline(h = 0, col = "gray50")
        }
        lines(s2, eff1)
        lines(s2, low1, lty = 2)
        lines(s2, up1, lty = 2)
        if (plot.type != "screen") {
            dev.off()
        }
        if (plot.type == "pdf") {
            pdf(paste(name.stem, "_", v2, ".pdf", sep = ""), 
                height = 6, width = 6)
        }
        if (plot.type == "png") {
            png(paste(name.stem, "_", v2, ".png", sep = ""))
        }
        if (plot.type == "eps") {
            postscript(paste(name.stem, "_", v2, ".eps", sep = ""))
        }
        plot(s1, eff2, type = "n", ylim = range(c(low2, up2)), 
            xlab = ifelse(is.null(xlab), toupper(v1), xlab[2]), 
			ylab = ifelse(is.null(ylab), paste("Conditional Effect of ", 
                toupper(v2), " | ", toupper(v1), sep = ""), ylab[2]))
        if (hist == TRUE) {
            rng <- diff(par()$usr[3:4])
            h1 <- hist(obj$model[[v1]], nclass = nclass[2], plot = FALSE)
            prop1 <- h1$counts/sum(h1$counts)
            plot.prop1 <- (prop1/max(prop1)) * rng * scale.hist + 
                par()$usr[3]
            av1 <- pretty(prop1, n = 3)
            axis(4, at = (av1/max(prop1)) * rng * scale.hist + 
                par()$usr[3], labels = av1)
            br1 <- h1$breaks
            for (i in 1:(length(br1) - 1)) {
                polygon(x = c(br1[i], br1[(i + 1)], br1[(i + 
                  1)], br1[i], br1[i]), y = c(par()$usr[3], par()$usr[3], 
                  plot.prop1[i], plot.prop1[i], par()$usr[3]), 
                  col = hist.col, border = border)
            }
        }
        if (rug == TRUE) {
            rug(obj$model[[v1]], ticksize = ticksize)
        }
        if (par()$usr[3] < 0 & par()$usr[4] > 0) {
            abline(h = 0, col = "gray50")
        }
        lines(s1, eff2)
        lines(s1, low2, lty = 2)
        lines(s1, up2, lty = 2)
        if (plot.type != "screen") {
            dev.off()
        }
        if (plot.type == "eps") {
            ps.options <- old.psopts
        }
        if (plot.type == "screen") {
            par <- oldpar
        }
    }
}
glmChange <-
function (obj, data, typical.dat = NULL, diffchange = c("range", "sd", "unit"), sim=FALSE, R=1000) 
{
    vars <- names(attr(terms(obj), "dataClasses"))[-1]
    pols <- grep("poly", vars)
    if (length(pols) > 0) {
        poly.split <- strsplit(vars[pols], split = "")
        start <- lapply(poly.split, function(x) grep("(", x, 
            fixed = T) + 1)
        stop <- lapply(poly.split, function(x) grep(",", x, fixed = T) - 
            1)
        pol.vars <- sapply(1:length(poly.split), function(x) paste(poly.split[[x]][start[[x]]:stop[[x]]], 
            collapse = ""))
        vars[pols] <- pol.vars
    }
    rn <- vars
    var.classes <- sapply(vars, function(x) class(data[[x]]))
    minmax <- lapply(vars, function(x) c(NA, NA))
    meds <- lapply(vars, function(x) NA)
    names(minmax) <- names(meds) <- vars
    levs <- obj$xlevels
    if (length(levs) > 0) {
        for (i in 1:length(levs)) {
            tmp.levs <- paste(names(levs)[i], unlist(levs[i]), 
                sep = "")
            col.inds <- match(tmp.levs, names(obj$coef))
            if (length(grep("1$", names(obj$coef)[col.inds])) > 
                0) {
                col.inds <- c(col.inds[which(is.na(col.inds))], 
                  col.inds[grep("1$", names(col.inds))])
                names(col.inds) <- gsub("1$", "", names(col.inds))
                col.inds <- col.inds[match(tmp.levs, names(col.inds))]
            }
            tmp.coefs <- obj$coef[col.inds]
            tmp.coefs <- obj$coef[match(tmp.levs, names(obj$coef))]
            tmp.coefs[which(is.na(tmp.coefs))] <- 0
            mm <- c(which.min(tmp.coefs), which.max(tmp.coefs))
            minmax[[names(levs)[i]]] <- factor(levs[[i]][mm], 
                levels = levs[[i]])
            tmp.tab <- table(data[[names(levs)[i]]])
            meds[[names(levs)[i]]] <- factor(names(tmp.tab)[which.max(tmp.tab)], 
                levels = levs[[i]])
        }
    }
    vars <- vars[sapply(minmax, function(x) is.na(x[1]))]
	mmc <- match.arg(diffchange)
    for (i in 1:length(vars)) {
		if(mmc == "range"){
        minmax[[vars[i]]] <- range(data[[vars[i]]], na.rm = T)
		}
		if(mmc == "sd"){
        minmax[[vars[i]]] <- median(data[[vars[i]]], na.rm = T) + c(-.5,.5)*sd(data[[vars[i]]], na.rm=T)
		}
		if(mmc == "unit"){
        minmax[[vars[i]]] <- median(data[[vars[i]]], na.rm = T) + c(-.5,.5)
		}
        meds[[vars[i]]] <- median(data[[vars[i]]], na.rm = T)
    }
    tmp.df <- do.call(data.frame, lapply(meds, function(x) rep(x, 
        length(meds) * 2)))
    if (!is.null(typical.dat)) {
        notin <- which(!(names(typical.dat) %in% names(tmp.df)))
        if (length(notin) > 0) {
            cat("The following variables in typical.dat were not found in the prediction data: ", 
                names(typical.dat)[notin], "\n\n", sep = "")
            typical.dat <- typical.dat[, -notin]
        }
        for (j in 1:ncol(typical.dat)) {
            tmp.df[[names(typical.dat)[j]]] <- typical.dat[1, 
                j]
            meds[names(typical.dat)[j]] <- as.numeric(typical.dat[1, 
                j])
        }
    }
    inds <- seq(1, nrow(tmp.df), by = 2)
    for (j in 1:length(minmax)) {
        tmp.df[inds[j]:(inds[j] + 1), j] <- minmax[[j]]
    }
    preds <- matrix(predict(obj, newdata = tmp.df, type = "response"), 
        ncol = 2, byrow = TRUE)
    diffs <- cbind(preds, apply(preds, 1, diff))
    colnames(diffs) <- c("min", "max", "diff")
    rownames(diffs) <- rn
    minmax.mat <- do.call(cbind, minmax)
    minmax.mat <- rbind(c(unlist(meds)), minmax.mat)
    rownames(minmax.mat) <- c("typical", "min", "max")
	if(sim){
    preds <- predict(obj, newdata = tmp.df, type = "link", se.fit=TRUE)
	res <- sapply(1:R, function(x)apply(matrix(family(obj)$linkinv(with(preds, rnorm(length(preds$fit), fit, se.fit))), ncol=2, byrow=TRUE), 1, diff))
	cis <- t(apply(res, 1, quantile, c(.025, .975)))
	colnames(cis) <- c("lower", "upper")
	diffs <- cbind(diffs, cis)
	}
    ret <- list(diffs = diffs, minmax = minmax.mat)
    class(ret) <- "change"
    return(ret)
}
intEff <-
function (obj, vars, data) 
{
    if (obj$family$family != "binomial") 
        stop("intEff only works for binomial GLMs")
    dat.cl <- attr(obj$terms, "dataClasses")[vars]
    if (dat.cl[vars[1]] == "factor" & length(unique(data[[vars[1]]])) == 
        2) {
        vars[1] <- paste(vars[1], colnames(contrasts(data[[vars[1]]])), 
            sep = "")
    }
    if (dat.cl[vars[1]] == "factor" & length(unique(data[[vars[1]]])) > 
        2) {
        stop("factor variables must only have two unique values, violated by v1")
    }
    if (dat.cl[vars[2]] == "factor" & length(unique(data[[vars[2]]])) == 
        2) {
        vars[2] <- paste(vars[2], colnames(contrasts(data[[vars[2]]])), 
            sep = "")
    }
    if (dat.cl[vars[2]] == "factor" & length(unique(data[[vars[2]]])) > 
        2) {
        stop("factor variables must only have two unique values, violated by v2")
    }
    v1 <- paste(vars, collapse = ":")
    v2 <- paste(vars[c(2, 1)], collapse = ":")
    inb <- any(c(v1, v2) %in% names(obj$coef))
    X <- model.matrix(obj, obj$model)
    if (!inb) 
        stop("Specified variables not interacted in model\n")
    if (v2 %in% names(obj$coef)) {
        vars <- vars[c(2, 1)]
    }
    int.var <- paste(vars, collapse = ":")
    b <- obj$coef
    lens <- apply(X[, vars], 2, function(x) length(unique(x)))
    if (any(dat.cl == "numeric")) 
        dat.cl[which(dat.cl == "numeric")] <- "c"
    if (any(dat.cl == "factor")) 
        dat.cl[which(dat.cl == "factor")] <- "d"
    if (any(lens == 2)) 
        dat.cl[which(lens == 2)] <- "d"
    type.int <- paste(sort(dat.cl), collapse = "")
    if (obj$family$link == "logit") 
        type.int <- paste("l", type.int, sep = "")
    if (obj$family$link == "probit") 
        type.int <- paste("p", type.int, sep = "")
    if (!(obj$family$link %in% c("probit", "logit"))) 
        stop("Link must be either logit or probit")
    out.dat <- switch(type.int, lcc = logit_cc(obj = obj, int.var = int.var, 
        vars = vars, b = b, X = X), lcd = logit_cd(obj = obj, 
        int.var = int.var, vars = vars, b = b, X = X), ldd = logit_dd(obj = obj, 
        int.var = int.var, vars = vars, b = b, X = X), pcc = probit_cc(obj = obj, 
        int.var = int.var, vars = vars, b = b, X = X), pcd = probit_cd(obj = obj, 
        int.var = int.var, vars = vars, b = b, X = X), pdd = probit_dd(obj = obj, 
        int.var = int.var, vars = vars, b = b, X = X))
    invisible(out.dat)
}
logit_cc <-
function (obj = obj, int.var = int.var, vars = vars, b = b, X = X) 
{
    phat <- predict(obj, type = "response")
    xb <- predict(obj, type = "link")
    phi <- phat * (1 - phat)
    linear <- b[int.var] * phi
    logitcc <- b[int.var] * phi + (b[vars[1]] + b[int.var] * 
        X[, vars[2]]) * (b[vars[2]] + b[int.var] * X[, vars[1]]) * 
        phi * (1 - (2 * phat))
    d2f <- phat * (1 - phat) * (1 - 2 * phat)
    d3f <- phat * (1 - phat) * (1 - 6 * phat + 6 * phat^2)
    b1b4x2 <- b[vars[1]] + b[int.var] * X[, vars[2]]
    b2b4x1 <- b[vars[2]] + b[int.var] * X[, vars[1]]
    deriv11 <- b[int.var] * d2f * X[, vars[1]] + b2b4x1 * d2f + 
        b1b4x2 * b2b4x1 * d3f * X[, vars[1]]
    deriv22 <- b[int.var] * d2f * X[, vars[2]] + b1b4x2 * d2f + 
        b1b4x2 * b2b4x1 * X[, vars[2]] * d3f
    deriv44 <- phi + b[int.var] * d2f * X[, vars[1]] * X[, vars[2]] + 
        X[, vars[2]] * b2b4x1 * d2f + X[, vars[1]] * b1b4x2 * 
        d2f + b1b4x2 * b2b4x1 * X[, vars[1]] * X[, vars[2]] * 
        d3f
    derivcc <- b[int.var] * d2f + b1b4x2 * b2b4x1 * d3f
    others <- X[, -c(1, match(c(vars, int.var), names(b)))]
    if (!("matrix" %in% class(others))) {
        others <- matrix(others, nrow = nrow(X))
    }
    colnames(others) <- colnames(X)[-c(1, match(c(vars, int.var), 
        names(b)))]
    nn <- apply(others, 2, function(x) b[int.var] * d2f * x + 
        b1b4x2 * b2b4x1 * x * d3f)
    mat123 <- cbind(deriv11, deriv22, deriv44, nn, derivcc)
    colnames(mat123) <- c(vars, int.var, colnames(nn), "(Intercept)")
    mat123 <- mat123[, match(colnames(X), colnames(mat123))]
    logit_se <- sqrt(diag(mat123 %*% vcov(obj) %*% t(mat123)))
    logit_t <- logitcc/logit_se
    out <- data.frame(int_eff = logitcc, linear = linear, phat = phat, 
        se_int_eff = logit_se, zstat = logit_t)
    invisible(out)
}
logit_cd <-
function (obj = obj, int.var = int.var, vars = vars, b = b, X = X) 
{
    phat <- predict(obj, type = "response")
    phi <- phat * (1 - phat)
    linear <- b[int.var] * phi
    dum <- vars[which(sapply(apply(X[, vars], 2, table), length) == 
        2)]
    cont <- vars[which(vars != dum)]
    X1 <- X2 <- X
    X1[, dum] <- 1
    X1[, int.var] <- X1[, cont] * X1[, dum]
    phat1 <- plogis(X1 %*% b)
    phi1 <- phat1 * (1 - phat1)
    d2f1 <- phi1 * (1 - 2 * phat1)
    ie1 <- (b[cont] + b[int.var]) * phi1
    X2[, dum] <- 0
    X2[, int.var] <- X2[, cont] * X2[, dum]
    phat2 <- plogis(X2 %*% b)
    phi2 <- phat2 * (1 - phat2)
    d2f2 <- phi2 * (1 - 2 * phat2)
    ie2 <- b[cont] * phi2
    logitcd <- ie1 - ie2
    deriv1 <- phi1 - phi2 + b[cont] * X[, cont] * (d2f1 - d2f2) + 
        b[int.var] * X[, cont] * d2f1
    deriv2 <- (b[cont] + b[int.var]) * d2f1
    deriv3 <- phi1 + (b[cont] + b[int.var]) * d2f1 * X[, cont]
    deriv0 <- (b[cont] + b[int.var]) * d2f1 - b[cont] * d2f2
    others <- X[, -c(1, match(c(vars, int.var), names(b)))]
    if (!("matrix" %in% class(others))) {
        others <- matrix(others, nrow = nrow(X))
    }
    colnames(others) <- colnames(X)[-c(1, match(c(vars, int.var), 
        names(b)))]
    nn <- apply(others, 2, function(x) ((b[cont] + b[int.var]) * 
        d2f1 - b[cont] * d2f2) * x)
    mat123 <- cbind(deriv1, deriv2, deriv3, nn, deriv0)
    colnames(mat123) <- c(vars, int.var, colnames(nn), "(Intercept)")
    mat123 <- mat123[, match(colnames(X), colnames(mat123))]
    logit_se <- sqrt(diag(mat123 %*% vcov(obj) %*% t(mat123)))
    logit_t <- logitcd/logit_se
    out <- data.frame(int_eff = logitcd, linear = linear, phat = phat, 
        se_int_eff = logit_se, zstat = logit_t)
    invisible(out)
}
logit_dd <-
function (obj = obj, int.var = int.var, vars = vars, b = b, X = X) 
{
    phat <- predict(obj, type = "response")
    phi <- phat * (1 - phat)
    linear <- b[int.var] * phi
    X11 <- X01 <- X10 <- X00 <- X
    X11[, vars[1]] <- 1
    X11[, vars[2]] <- 1
    X10[, vars[1]] <- 1
    X10[, vars[2]] <- 0
    X01[, vars[1]] <- 0
    X01[, vars[2]] <- 1
    X00[, vars[1]] <- 0
    X00[, vars[2]] <- 0
    X00[, int.var] <- X00[, vars[1]] * X00[, vars[2]]
    X11[, int.var] <- X11[, vars[1]] * X11[, vars[2]]
    X01[, int.var] <- X01[, vars[1]] * X01[, vars[2]]
    X10[, int.var] <- X10[, vars[1]] * X10[, vars[2]]
    phat11 <- plogis(X11 %*% b)
    phat00 <- plogis(X00 %*% b)
    phat10 <- plogis(X10 %*% b)
    phat01 <- plogis(X01 %*% b)
    phi11 <- phat11 * (1 - phat11)
    phi10 <- phat10 * (1 - phat10)
    phi01 <- phat01 * (1 - phat01)
    phi00 <- phat00 * (1 - phat00)
    logitdd <- (phat11 - phat10) - (phat01 - phat00)
    deriv1 <- phi11 - phi10
    deriv2 <- phi11 - phi01
    deriv3 <- phi11
    deriv0 <- (phi11 - phi01) - (phi10 - phi00)
    others <- X[, -c(1, match(c(vars, int.var), names(b)))]
    if (!("matrix" %in% class(others))) {
        others <- matrix(others, nrow = nrow(X))
    }
    colnames(others) <- colnames(X)[-c(1, match(c(vars, int.var), 
        names(b)))]
    nn <- apply(others, 2, function(x) ((phi11 - phi01) - (phi10 - 
        phi00)) * x)
    mat123 <- cbind(deriv1, deriv2, deriv3, nn, deriv0)
    colnames(mat123) <- c(vars, int.var, colnames(nn), "(Intercept)")
    mat123 <- mat123[, match(colnames(X), colnames(mat123))]
    logit_se <- sqrt(diag(mat123 %*% vcov(obj) %*% t(mat123)))
    logit_t <- logitdd/logit_se
    out <- data.frame(int_eff = logitdd, linear = linear, phat = phat, 
        se_int_eff = logit_se, zstat = logit_t)
    invisible(out)
}
mnlSig <-
function (obj, pval = 0.05, two.sided = TRUE, flag.sig = TRUE, 
    insig.blank = FALSE) 
{
    smulti <- summary(obj)
    multi.t <- smulti$coefficients/smulti$standard.errors
    multi.p <- (2^as.numeric(two.sided)) * pnorm(abs(multi.t), 
        lower.tail = F)
    b <- matrix(sprintf("%.3f", smulti$coefficients), ncol = ncol(multi.t))
    sig.vec <- c(" ", "*")
    sig.obs <- as.numeric(multi.p < pval) + 1
    if (flag.sig) {
        b <- matrix(paste(b, sig.vec[sig.obs], sep = ""), ncol = ncol(multi.t))
    }
    if (insig.blank) {
        b[which(multi.p > pval, arr.ind = T)] <- ""
    }
    rownames(b) <- rownames(multi.t)
    colnames(b) <- colnames(multi.t)
    b <- as.data.frame(b)
    return(b)
}
multiChange <-
function (obj, data, typical.dat = NULL, diffchange=c("range", "sd", "unit")) 
{
    vars <- names(attr(terms(obj), "dataClasses"))[-1]
    pols <- grep("poly", vars)
    if (length(pols) > 0) {
        poly.split <- strsplit(vars[pols], split = "")
        start <- lapply(poly.split, function(x) grep("(", x, 
            fixed = T) + 1)
        stop <- lapply(poly.split, function(x) grep(",", x, fixed = T) - 
            1)
        pol.vars <- sapply(1:length(poly.split), function(x) paste(poly.split[[x]][start[[x]]:stop[[x]]], 
            collapse = ""))
        vars[pols] <- pol.vars
    }
    rn <- vars
    var.classes <- sapply(vars, function(x) class(data[[x]]))
    minmax <- lapply(vars, function(x) c(NA, NA))
    meds <- lapply(vars, function(x) NA)
    names(minmax) <- names(meds) <- vars
    levs <- obj$xlevels
    if (length(levs) > 0) {
        for (i in 1:length(levs)) {
            tmp.levs <- paste(names(levs)[i], unlist(levs[i]), 
                sep = "")
            col.inds <- sapply(tmp.levs, function(x) grep(x, 
                colnames(coefficients(obj))))
            col.inds[which(sapply(col.inds, length) == 0)] <- NA
            col.inds <- unlist(col.inds)
            if (length(grep("1$", names(col.inds))) > 0) {
                col.inds <- c(col.inds[which(is.na(col.inds))], 
                  col.inds[grep("1$", names(col.inds))])
                names(col.inds) <- gsub("1$", "", names(col.inds))
                col.inds <- col.inds[match(tmp.levs, names(col.inds))]
            }
            tmp.coefs <- coefficients(obj)[, col.inds]
            tmp.coefs[which(is.na(tmp.coefs))] <- 0
            tmp.coefs <- apply(tmp.coefs, 2, sum)
            mm <- c(which.min(tmp.coefs), which.max(tmp.coefs))
            minmax[[names(levs)[i]]] <- factor(levs[[i]][mm], 
                levels = levs[[i]])
            tmp.tab <- table(data[[names(levs)[i]]])
            meds[[names(levs)[i]]] <- factor(names(tmp.tab)[which.max(tmp.tab)], 
                levels = levs[[i]])
        }
    }
    vars <- vars[sapply(minmax, function(x) is.na(x[1]))]
	mmc <- match.arg(diffchange)
    for (i in 1:length(vars)) {
		if(mmc == "range"){
        minmax[[vars[i]]] <- range(data[[vars[i]]], na.rm = T)
		}
		if(mmc == "sd"){
        minmax[[vars[i]]] <- median(data[[vars[i]]], na.rm = T) + c(-.5,.5)*sd(data[[vars[i]]], na.rm=T)
		}
		if(mmc == "unit"){
        minmax[[vars[i]]] <- median(data[[vars[i]]], na.rm = T) + c(-.5,.5)
		}
        meds[[vars[i]]] <- median(data[[vars[i]]], na.rm = T)
    }
    tmp.df <- do.call(data.frame, lapply(meds, function(x) rep(x, 
        length(meds) * 2)))
    if (!is.null(typical.dat)) {
        notin <- which(!(names(typical.dat) %in% names(tmp.df)))
        if (length(notin) > 0) {
            cat("The following variables in typical.dat were not found in the prediction data: ", 
                names(typical.dat)[notin], "\n\n", sep = "")
            typical.dat <- typical.dat[, -notin]
        }
        for (j in 1:ncol(typical.dat)) {
            tmp.df[[names(typical.dat)[j]]] <- typical.dat[1, 
                j]
            meds[names(typical.dat)[j]] <- as.numeric(typical.dat[1, 
                j])
        }
    }
    inds <- seq(1, nrow(tmp.df), by = 2)
    for (j in 1:length(minmax)) {
        tmp.df[inds[j]:(inds[j] + 1), j] <- minmax[[j]]
    }
    preds <- predict(obj, newdata = tmp.df, type = "probs")
    preds.min <- preds[seq(1, nrow(preds), by = 2), ]
    preds.max <- preds[seq(2, nrow(preds), by = 2), ]
    diffs <- preds.max - preds.min
    rownames(preds.min) <- rownames(preds.max) <- rownames(diffs) <- rn
    minmax.mat <- do.call(data.frame, minmax)
    minmax.mat <- rbind(do.call(data.frame, meds), minmax.mat)
    rownames(minmax.mat) <- c("typical", "min", "max")
    ret <- list(diffs = diffs, minmax = minmax.mat, minPred = preds.min, 
        maxPred = preds.max)
    class(ret) <- "change"
    return(ret)
}
ordChange <-
function (obj, data, typical.dat = NULL, diffchange=c("range", "sd", "unit")) 
{
    vars <- names(attr(terms(obj), "dataClasses"))[-1]
    pols <- grep("poly", vars)
    if (length(pols) > 0) {
        poly.split <- strsplit(vars[pols], split = "")
        start <- lapply(poly.split, function(x) grep("(", x, 
            fixed = T) + 1)
        stop <- lapply(poly.split, function(x) grep(",", x, fixed = T) - 
            1)
        pol.vars <- sapply(1:length(poly.split), function(x) paste(poly.split[[x]][start[[x]]:stop[[x]]], 
            collapse = ""))
        vars[pols] <- pol.vars
    }
    rn <- vars
    var.classes <- sapply(vars, function(x) class(data[[x]]))
    minmax <- lapply(vars, function(x) c(NA, NA))
    meds <- lapply(vars, function(x) NA)
    names(minmax) <- names(meds) <- vars
    levs <- obj$xlevels
    if (length(levs) > 0) {
        for (i in 1:length(levs)) {
            tmp.levs <- paste(names(levs)[i], unlist(levs[i]), 
                sep = "")
            col.inds <- match(tmp.levs, names(obj$coef))
            if (length(grep("1$", names(obj$coef)[col.inds])) > 
                0) {
                col.inds <- c(col.inds[which(is.na(col.inds))], 
                  col.inds[grep("1$", names(col.inds))])
                names(col.inds) <- gsub("1$", "", names(col.inds))
                col.inds <- col.inds[match(tmp.levs, names(col.inds))]
            }
            tmp.coefs <- obj$coef[col.inds]
            tmp.coefs[which(is.na(tmp.coefs))] <- 0
            mm <- c(which.min(tmp.coefs), which.max(tmp.coefs))
            minmax[[names(levs)[i]]] <- factor(levs[[i]][mm], 
                levels = levs[[i]])
            tmp.tab <- table(data[[names(levs)[i]]])
            meds[[names(levs)[i]]] <- factor(names(tmp.tab)[which.max(tmp.tab)], 
                levels = levs[[i]])
        }
    }
    vars <- vars[sapply(minmax, function(x) is.na(x[1]))]
	mmc <- match.arg(diffchange)
    for (i in 1:length(vars)) {
		if(mmc == "range"){
        minmax[[vars[i]]] <- range(data[[vars[i]]], na.rm = T)
		}
		if(mmc == "sd"){
        minmax[[vars[i]]] <- median(data[[vars[i]]], na.rm = T) + c(-.5,.5)*sd(data[[vars[i]]], na.rm=T)
		}
		if(mmc == "unit"){
        minmax[[vars[i]]] <- median(data[[vars[i]]], na.rm = T) + c(-.5,.5)
		}
        meds[[vars[i]]] <- median(data[[vars[i]]], na.rm = T)
    }
    tmp.df <- do.call(data.frame, lapply(meds, function(x) rep(x, 
        length(meds) * 2)))
    if (!is.null(typical.dat)) {
        notin <- which(!(names(typical.dat) %in% names(tmp.df)))
        if (length(notin) > 0) {
            cat("The following variables in typical.dat were not found in the prediction data: ", 
                names(typical.dat)[notin], "\n\n", sep = "")
            typical.dat <- typical.dat[, -notin]
        }
        for (j in 1:ncol(typical.dat)) {
            tmp.df[[names(typical.dat)[j]]] <- typical.dat[1, 
                j]
            meds[names(typical.dat)[j]] <- as.numeric(typical.dat[1, 
                j])
        }
    }
    inds <- seq(1, nrow(tmp.df), by = 2)
    for (j in 1:length(minmax)) {
        tmp.df[inds[j]:(inds[j] + 1), j] <- minmax[[j]]
    }
    preds <- predict(obj, newdata = tmp.df, type = "probs")
    preds.min <- preds[seq(1, nrow(preds), by = 2), ]
    preds.max <- preds[seq(2, nrow(preds), by = 2), ]
    diffs <- preds.max - preds.min
    rownames(preds.min) <- rownames(preds.max) <- rownames(diffs) <- rn
    minmax.mat <- do.call(data.frame, minmax)
    minmax.mat <- rbind(do.call(data.frame, meds), minmax.mat)
    rownames(minmax.mat) <- c("typical", "min", "max")
    ret <- list(diffs = diffs, minmax = minmax.mat, minPred = preds.min, 
        maxPred = preds.max)
    class(ret) <- "change"
    return(ret)
}
poisGOF <-
function (obj) 
{
    if (!("glm" %in% class(obj))) {
        stop("poisGOF only works on objects of class glm (with a poisson family)\n")
    }
    if (!("y" %in% names(obj))) {
        obj <- update(obj, y = TRUE)
    }
    ind.chisq <- (((obj$y - obj$fitted)^2)/obj$fitted)
    dev <- obj$deviance
    df <- obj$df.residual
    p.chisq <- pchisq(sum(ind.chisq), df = df, lower.tail = FALSE)
    p.dev <- pchisq(dev, df = df, lower.tail = FALSE)
    vec <- sprintf("%.3f", c(sum(ind.chisq), dev, p.chisq, p.dev))
    mat <- matrix(vec, ncol = 2, nrow = 2)
    rownames(mat) <- c("Chi-squared", "Deviance")
    colnames(mat) <- c("Stat", "p-value")
    mat <- as.data.frame(mat)
    return(mat)
}
pre <-
function (mod1, mod2 = NULL, sim = FALSE, R = 2500) 
{
    if (!is.null(mod2)) {
        if (mean(class(mod1) == class(mod2)) != 1) {
            stop("Model 2 must be either NULL or of the same class as Model 1\n")
        }
    }
    if (!any(class(mod1) %in% c("polr", "multinom", "glm"))) {
        stop("pre only works on models of class glm (with binomial family), polr or multinom\n")
    }
    if ("glm" %in% class(mod1)) {
        if (!(family(mod1)$link %in% c("logit", "probit", "cloglog", 
            "cauchit"))) {
            stop("PRE only calculated for models with logit, probit, cloglog or cauchit links\n")
        }
        if (is.null(mod2)) {
            y <- mod1[["y"]]
            mod2 <- update(mod1, ". ~ 1", data=model.frame(mod1))
        }
        pred.mod2 <- as.numeric(predict(mod2, type = "response") >= 
            0.5)
        pmc <- mean(mod2$y == pred.mod2)
        pred.y <- as.numeric(predict(mod1, type = "response") >= 
            0.5)
        pcp <- mean(pred.y == mod1$y)
        pre <- (pcp - pmc)/(1 - pmc)
        pred.prob1 <- predict(mod1, type = "response")
        pred.prob2 <- predict(mod2, type = "response")
        epcp <- (1/length(pred.prob1)) * (sum(pred.prob1[which(mod1$y == 
            1)]) + sum(1 - pred.prob1[which(mod1$y == 0)]))
        epmc <- (1/length(pred.prob2)) * (sum(pred.prob2[which(mod2$y == 
            1)]) + sum(1 - pred.prob2[which(mod2$y == 0)]))
        epre <- (epcp - epmc)/(1 - epmc)
        if (sim) {
            b1.sim <- mvrnorm(R, coef(mod1), vcov(mod1))
            b2.sim <- mvrnorm(R, coef(mod2), vcov(mod2))
            mod1.probs <- family(mod1)$linkinv(model.matrix(mod1) %*% 
                t(b1.sim))
            mod2.probs <- family(mod2)$linkinv(model.matrix(mod2) %*% 
                t(b2.sim))
            pmcs <- apply(mod2.probs, 2, function(x) mean(as.numeric(x > 
                0.5) == mod2$y))
            pcps <- apply(mod1.probs, 2, function(x) mean(as.numeric(x > 
                0.5) == mod1$y))
            pre.sim <- (pcps - pmcs)/(1 - pmcs)
            epmc.sim <- apply(mod2.probs, 2, function(x) (1/length(x)) * 
                (sum(x[which(mod2$y == 1)]) + sum(1 - x[which(mod2$y == 
                  0)])))
            epcp.sim <- apply(mod1.probs, 2, function(x) (1/length(x)) * 
                (sum(x[which(mod1$y == 1)]) + sum(1 - x[which(mod1$y == 
                  0)])))
            epre.sim <- (epcp.sim - epmc.sim)/(1 - epmc.sim)
        }
    }
    if ("multinom" %in% class(mod1)) {
        mod1 <- update(mod1, ".~.", model = TRUE, trace = FALSE)
        if (is.null(mod2)) {
            mod2 <- update(mod1, ". ~ 1", data = mod1$model, 
                trace = FALSE)
        }
        pred.prob1 <- predict(mod1, type = "prob")
        pred.prob2 <- predict(mod2, type = "prob")
        pred.cat1 <- apply(pred.prob1, 1, which.max)
        pred.cat2 <- apply(pred.prob2, 1, which.max)
        pcp <- mean(as.numeric(mod1$model[, 1]) == pred.cat1)
        pmc <- max(table(as.numeric(mod2$model[, 1]))/sum(table(mod2$model[, 
            1])))
        pre <- (pcp - pmc)/(1 - pmc)
        epcp <- mean(pred.prob1[cbind(1:nrow(pred.prob1), as.numeric(mod1$model[, 
            1]))])
        tab <- table(mod1$model[, 1])/sum(table(mod1$model[, 
            1]))
        epmc <- mean(tab[as.numeric(mod2$model[, 1])])
        epre <- (epcp - epmc)/(1 - epmc)
        if (sim) {
            b1.sim <- mvrnorm(R, c(t(coef(mod1))), vcov(mod1))
            b2.sim <- mvrnorm(R, c(t(coef(mod2))), vcov(mod2))
            mod.levs <- rownames(coef(mod1))
            var.levs <- rownames(contrasts(mod1$model[, 1]))
            tmp.reord <- c(match(mod.levs, var.levs), (1:length(var.levs))[-match(mod.levs, 
                var.levs)])
            reord <- match(1:length(var.levs), tmp.reord)
            tmp <- lapply(1:nrow(b1.sim), function(x) matrix(b1.sim[x, 
                ], ncol = ncol(coef(mod1)), byrow = T))
            tmp2 <- lapply(tmp, function(x) cbind(model.matrix(mod1) %*% 
                t(x), 0))
            tmp2 <- lapply(tmp2, function(x) x[, reord])
            mod1.probs <- lapply(tmp2, function(x) exp(x)/apply(exp(x), 
                1, sum))
            tmp <- lapply(1:nrow(b2.sim), function(x) matrix(b2.sim[x, 
                ], ncol = ncol(coef(mod2)), byrow = T))
            tmp2 <- lapply(tmp, function(x) cbind(model.matrix(mod2) %*% 
                t(x), 0))
            tmp2 <- lapply(tmp2, function(x) x[, reord])
            mod2.probs <- lapply(tmp2, function(x) exp(x)/apply(exp(x), 
                1, sum))
            pred.cat1 <- lapply(mod1.probs, function(x) apply(x, 
                1, which.max))
            pred.cat2 <- lapply(mod2.probs, function(x) apply(x, 
                1, which.max))
            pcp.sim <- sapply(pred.cat1, function(x) mean(x == 
                as.numeric(mod1$model[, 1])))
            pmc.sim <- sapply(pred.cat2, function(x) mean(x == 
                as.numeric(mod1$model[, 1])))
            pre.sim <- (pcp.sim - pmc.sim)/(1 - pmc.sim)
            epcp.sim <- sapply(mod1.probs, function(z) mean(z[cbind(1:nrow(mod1$model), 
                as.numeric(mod1$model[, 1]))]))
            epmc.sim <- sapply(mod2.probs, function(z) mean(z[cbind(1:nrow(mod1$model), 
                as.numeric(mod1$model[, 1]))]))
            epre.sim <- (epcp.sim - epmc.sim)/(1 - epmc.sim)
        }
    }
    if ("polr" %in% class(mod1)) {
        if (is.null(mod1$Hess)) {
            mod1 <- update(mod1, Hess = TRUE)
        }
        if (is.null(mod2)) {
            mod2 <- update(mod1, ". ~ 1", data = mod1$model, 
                model = TRUE, Hess = TRUE)
        }
        pred.prob1 <- predict(mod1, type = "prob")
        pred.prob2 <- predict(mod2, type = "prob")
        pred.cat1 <- apply(pred.prob1, 1, which.max)
        pred.cat2 <- apply(pred.prob2, 1, which.max)
        pcp <- mean(as.numeric(mod1$model[, 1]) == pred.cat1)
        pmc <- max(table(as.numeric(mod2$model[, 1]))/sum(table(mod2$model[, 
            1])))
        pre <- (pcp - pmc)/(1 - pmc)
        epcp <- mean(pred.prob1[cbind(1:nrow(pred.prob1), as.numeric(mod1$model[, 
            1]))])
        tab <- table(mod1$model[, 1])/sum(table(mod1$model[, 
            1]))
        epmc <- mean(tab[as.numeric(mod2$model[, 1])])
        epre <- (epcp - epmc)/(1 - epmc)
        if (sim) {
            b1 <- c(coef(mod1), mod1$zeta)
            b2 <- c(coef(mod2), mod2$zeta)
            v1 <- invisible(vcov(mod1))
            v2 <- invisible(vcov(mod2))
            b1.sim <- mvrnorm(R, b1, v1)
            b2.sim <- mvrnorm(R, b2, v2)
            mod1.probs <- lapply(1:nrow(b1.sim), function(x) simPredpolr(mod1, 
                b1.sim[x, ], n.coef = length(coef(mod1))))
            mod2.probs <- lapply(1:nrow(b2.sim), function(x) simPredpolr(mod2, 
                b2.sim[x, ], n.coef = length(coef(mod2))))
            pred.cat1 <- lapply(mod1.probs, function(x) apply(x, 
                1, which.max))
            pred.cat2 <- lapply(mod2.probs, function(x) apply(x, 
                1, which.max))
            pcp.sim <- sapply(pred.cat1, function(x) mean(x == 
                as.numeric(mod1$model[, 1])))
            pmc.sim <- sapply(pred.cat2, function(x) mean(x == 
                as.numeric(mod1$model[, 1])))
            pre.sim <- (pcp.sim - pmc.sim)/(1 - pmc.sim)
            epcp.sim <- sapply(mod1.probs, function(z) mean(z[cbind(1:nrow(mod1$model), 
                as.numeric(mod1$model[, 1]))]))
            epmc.sim <- sapply(mod2.probs, function(z) mean(z[cbind(1:nrow(mod1$model), 
                as.numeric(mod1$model[, 1]))]))
            epre.sim <- (epcp.sim - epmc.sim)/(1 - epmc.sim)
        }
    }
    ret <- list()
    ret$pre <- pre
    ret$epre <- epre
    form1 <- formula(mod1)
    form2 <- formula(mod2)
    ret$m1form <- paste(form1[2], form1[1], form1[3], sep = " ")
    ret$m2form <- paste(form2[2], form2[1], form2[3], sep = " ")
    ret$pcp <- pcp
    ret$pmc <- pmc
    ret$epmc <- epmc
    ret$epcp <- epcp
    if (sim) {
        ret$pre.sim <- pre.sim
        ret$epre.sim <- epre.sim
    }
    class(ret) <- "pre"
    return(ret)
}
print.pre <-
function (x, ..., sim.ci = 0.95) 
{
    cat("mod1: ", as.character(x$m1form), "\n")
    cat("mod2: ", as.character(x$m2form), "\n\n")
    cat("Analytical Results\n")
    cat(" PMC = ", sprintf("%2.3f", x$pmc), "\n")
    cat(" PCP = ", sprintf("%2.3f", x$pcp), "\n")
    cat(" PRE = ", sprintf("%2.3f", x$pre), "\n")
    cat("ePMC = ", sprintf("%2.3f", x$epmc), "\n")
    cat("ePCP = ", sprintf("%2.3f", x$epcp), "\n")
    cat("ePRE = ", sprintf("%2.3f", x$epre), "\n\n")
    low <- (1 - sim.ci)/2
    up <- 1 - low
    if (exists("pre.sim", x)) {
        pre.ci <- sprintf("%2.3f", quantile(x$pre.sim, c(0.5, 
            low, up)))
        epre.ci <- sprintf("%2.3f", quantile(x$epre.sim, c(0.5, 
            low, up)))
        tmp <- rbind(pre.ci, epre.ci)
        rownames(tmp) <- c(" PRE", "ePRE")
        colnames(tmp) <- c("median", "lower", "upper")
        cat("Simulated Results\n")
        print(tmp, quote = F)
    }
}
probit_cc <-
function (obj = obj, int.var = int.var, vars = vars, b = b, X = X) 
{
    phat <- predict(obj, type = "response")
    xb <- X %*% b
    phi <- dnorm(xb)
    linear <- b[int.var] * phi
    probitcc <- (b[int.var] - (b[vars[1]] + b[int.var] * X[, 
        vars[2]]) * (b[vars[2]] + b[int.var] * X[, vars[1]]) * 
        xb) * phi
    d2f <- -xb * phi
    d3f <- (xb^2 - 1) * phi
    b1b4x2 <- b[vars[1]] + b[int.var] * X[, vars[2]]
    b2b4x1 <- b[vars[2]] + b[int.var] * X[, vars[1]]
    deriv11 <- b[int.var] * d2f * X[, vars[1]] + b2b4x1 * d2f + 
        b1b4x2 * b2b4x1 * d3f * X[, vars[1]]
    deriv22 <- b[int.var] * d2f * X[, vars[2]] + b1b4x2 * d2f + 
        b1b4x2 * b2b4x1 * X[, vars[2]] * d3f
    deriv44 <- phi + b[int.var] * d2f * X[, vars[1]] * X[, vars[2]] + 
        X[, vars[2]] * b2b4x1 * d2f + X[, vars[1]] * b1b4x2 * 
        d2f + b1b4x2 * b2b4x1 * X[, vars[1]] * X[, vars[2]] * 
        d3f
    derivcc <- b[int.var] * d2f + b1b4x2 * b2b4x1 * d3f
    others <- X[, -c(1, match(c(vars, int.var), names(b)))]
    if (!("matrix" %in% class(others))) {
        others <- matrix(others, nrow = nrow(X))
    }
    colnames(others) <- colnames(X)[-c(1, match(c(vars, int.var), 
        names(b)))]
    nn <- apply(others, 2, function(x) b[int.var] * d2f * x + 
        b1b4x2 * b2b4x1 * x * d3f)
    mat123 <- cbind(deriv11, deriv22, deriv44, nn, derivcc)
    colnames(mat123) <- c(vars, int.var, colnames(nn), "(Intercept)")
    mat123 <- mat123[, match(colnames(X), colnames(mat123))]
    probit_se <- sqrt(diag(mat123 %*% vcov(obj) %*% t(mat123)))
    probit_t <- probitcc/probit_se
    out <- data.frame(int_eff = probitcc, linear = linear, phat = phat, 
        se_int_eff = probit_se, zstat = probit_t)
    invisible(out)
}
probit_cd <-
function (obj = obj, int.var = int.var, vars = vars, b = b, X = X) 
{
    phat <- predict(obj, type = "response")
    xb <- predict(obj, type = "link")
    phi <- dnorm(xb)
    linear <- b[int.var] * phi
    dum <- vars[which(sapply(apply(X[, vars], 2, table), length) == 
        2)]
    cont <- vars[which(vars != dum)]
    X1 <- X2 <- X
    X1[, dum] <- 1
    X1[, int.var] <- X1[, cont] * X1[, dum]
    phi1 <- dnorm(X1 %*% b)
    d2f1 <- -(X1 %*% b) * phi1
    ie1 <- (b[cont] + b[int.var]) * phi1
    X2[, dum] <- 0
    X2[, int.var] <- X2[, cont] * X2[, dum]
    phi2 <- dnorm(X2 %*% b)
    d2f2 <- -(X2 %*% b) * phi2
    ie2 <- b[cont] * phi2
    probitcd <- ie1 - ie2
    deriv1 <- phi1 - phi2 + b[cont] * X[, cont] * (d2f1 - d2f2) + 
        b[int.var] * X[, cont] * d2f1
    deriv2 <- (b[cont] + b[int.var]) * d2f1
    deriv3 <- phi1 + (b[cont] + b[int.var]) * d2f1 * X[, cont]
    deriv0 <- (b[cont] + b[int.var]) * d2f1 - b[cont] * d2f2
    others <- X[, -c(1, match(c(vars, int.var), names(b)))]
    if (!("matrix" %in% class(others))) {
        others <- matrix(others, nrow = nrow(X))
    }
    colnames(others) <- colnames(X)[-c(1, match(c(vars, int.var), 
        names(b)))]
    nn <- apply(others, 2, function(x) ((b[cont] + b[int.var]) * 
        d2f1 - b[cont] * d2f2) * x)
    mat123 <- cbind(deriv1, deriv2, deriv3, nn, deriv0)
    colnames(mat123) <- c(vars, int.var, colnames(nn), "(Intercept)")
    mat123 <- mat123[, match(colnames(X), colnames(mat123))]
    probit_se <- sqrt(diag(mat123 %*% vcov(obj) %*% t(mat123)))
    probit_t <- probitcd/probit_se
    out <- data.frame(int_eff = probitcd, linear = linear, phat = phat, 
        se_int_eff = probit_se, zstat = probit_t)
    invisible(out)
}
probit_dd <-
function (obj = obj, int.var = int.var, vars = vars, b = b, X = X) 
{
    phat <- predict(obj, type = "response")
    xb <- predict(obj, type = "link")
    phi <- dnorm(xb)
    linear <- b[int.var] * phi
    X11 <- X01 <- X10 <- X00 <- X
    X11[, vars[1]] <- 1
    X11[, vars[2]] <- 1
    X10[, vars[1]] <- 1
    X10[, vars[2]] <- 0
    X01[, vars[1]] <- 0
    X01[, vars[2]] <- 1
    X00[, vars[1]] <- 0
    X00[, vars[2]] <- 0
    X00[, int.var] <- X00[, vars[1]] * X00[, vars[2]]
    X11[, int.var] <- X11[, vars[1]] * X11[, vars[2]]
    X01[, int.var] <- X01[, vars[1]] * X01[, vars[2]]
    X10[, int.var] <- X10[, vars[1]] * X10[, vars[2]]
    Xb11 <- X11 %*% b
    Xb00 <- X00 %*% b
    Xb10 <- X10 %*% b
    Xb01 <- X01 %*% b
    phat11 <- pnorm(Xb11)
    phat00 <- pnorm(Xb00)
    phat10 <- pnorm(Xb10)
    phat01 <- pnorm(Xb01)
    phi11 <- dnorm(Xb11)
    phi00 <- dnorm(Xb00)
    phi10 <- dnorm(Xb10)
    phi01 <- dnorm(Xb01)
    probitdd <- (phat11 - phat10) - (phat01 - phat00)
    deriv1 <- phi11 - phi10
    deriv2 <- phi11 - phi01
    deriv3 <- phi11
    deriv0 <- (phi11 - phi01) - (phi10 - phi00)
    others <- X[, -c(1, match(c(vars, int.var), names(b)))]
    if (!("matrix" %in% class(others))) {
        others <- matrix(others, nrow = nrow(X))
    }
    colnames(others) <- colnames(X)[-c(1, match(c(vars, int.var), 
        names(b)))]
    nn <- apply(others, 2, function(x) ((phi11 - phi01) - (phi10 - 
        phi00)) * x)
    mat123 <- cbind(deriv1, deriv2, deriv3, nn, deriv0)
    colnames(mat123) <- c(vars, int.var, colnames(nn), "(Intercept)")
    mat123 <- mat123[, match(colnames(X), colnames(mat123))]
    probit_se <- sqrt(diag(mat123 %*% vcov(obj) %*% t(mat123)))
    probit_t <- probitdd/probit_se
    out <- data.frame(int_eff = probitdd, linear = linear, phat = phat, 
        se_int_eff = probit_se, zstat = probit_t)
    invisible(out)
}
searchVarLabels <-
function (dat, str) 
{
    if ("var.labels" %in% names(attributes(dat))) {
        vlat <- "var.labels"
    }
    if ("variable.labels" %in% names(attributes(dat))) {
        vlat <- "variable.labels"
    }
    ind <- sort(union(grep(str, attr(dat, vlat), ignore.case = T), grep(str, names(dat), ignore.case = T)))
    vldf <- data.frame(ind = ind, label = attr(dat, vlat)[ind])
    rownames(vldf) <- names(dat)[ind]
    vldf
}
simPredpolr <-
function (object, coefs, n.coef) 
{
    X <- model.matrix(object)
    xint <- match("(Intercept)", colnames(X), nomatch = 0L)
    if (xint > 0L) 
        X <- X[, -xint, drop = FALSE]
    n <- nrow(X)
    q <- length(coefs) - n.coef
    eta <- {
        if (n.coef > 0) {
            drop(X %*% coefs[1:n.coef])
        }
        else {
            rep(0, n)
        }
    }
    pfun <- switch(object$method, logistic = plogis, probit = pnorm, 
        cauchit = pcauchy)
    cumpr <- matrix(pfun(matrix(coefs[(n.coef + 1):length(coefs)], 
        n, q, byrow = TRUE) - eta), , q)
    Y <- t(apply(cumpr, 1L, function(x) diff(c(0, x, 1))))
    return(Y)
}
ziChange <-
function (obj, data, typical.dat = NULL, type = "count") 
{
    vars <- as.character(formula(obj))[3]
    vars <- gsub("*", "+", vars, fixed = T)
    vars <- strsplit(vars, split = "|", fixed = T)[[1]]
    vars <- c(unlist(strsplit(vars, "+", fixed = T)))
    vars <- unique(trim(vars))
    pols <- grep("poly", vars)
    if (length(pols) > 0) {
        poly.split <- strsplit(vars[pols], split = "")
        start <- lapply(poly.split, function(x) grep("(", x, 
            fixed = T) + 1)
        stop <- lapply(poly.split, function(x) grep(",", x, fixed = T) - 
            1)
        pol.vars <- sapply(1:length(poly.split), function(x) paste(poly.split[[x]][start[[x]]:stop[[x]]], 
            collapse = ""))
        vars[pols] <- pol.vars
    }
    rn <- vars
    vars.type <- as.character(formula(obj))[3]
    vars.type <- gsub("*", "+", vars.type, fixed = T)
    vars.type <- strsplit(vars.type, split = "|", fixed = T)[[1]][ifelse(type == 
        "count", 1, 2)]
    vars.type <- c(unlist(strsplit(vars.type, "+", fixed = T)))
    vars.type <- unique(trim(vars.type))
    pols <- grep("poly", vars.type)
    if (length(pols) > 0) {
        poly.split <- strsplit(vars.type[pols], split = "")
        start <- lapply(poly.split, function(x) grep("(", x, 
            fixed = T) + 1)
        stop <- lapply(poly.split, function(x) grep(",", x, fixed = T) - 
            1)
        pol.vars.type <- sapply(1:length(poly.split), function(x) paste(poly.split[[x]][start[[x]]:stop[[x]]], 
            collapse = ""))
        vars.type[pols] <- pol.vars.type
    }
    var.classes <- sapply(vars, function(x) class(data[[x]]))
    minmax <- lapply(vars, function(x) c(NA, NA))
    meds <- lapply(vars, function(x) NA)
    names(minmax) <- names(meds) <- vars
    levs <- obj$levels
    if (length(levs) > 0) {
        for (i in 1:length(levs)) {
            tmp.levs <- paste(names(levs)[i], unlist(levs[i]), 
                sep = "")
            tmp.coefs <- obj$coef[[type]][match(tmp.levs, names(obj$coef[[type]]))]
            tmp.coefs[which(is.na(tmp.coefs))] <- 0
            mm <- c(which.min(tmp.coefs), which.max(tmp.coefs))
            minmax[[names(levs)[i]]] <- factor(levs[[i]][mm], 
                levels = levs[[i]])
            tmp.tab <- table(data[[names(levs)[i]]])
            meds[[names(levs)[i]]] <- factor(names(tmp.tab)[which.max(tmp.tab)], 
                levels = levs[[i]])
        }
    }
    vars <- vars[sapply(minmax, function(x) is.na(x[1]))]
    for (i in 1:length(vars)) {
        minmax[[vars[i]]] <- range(data[[vars[i]]], na.rm = T)
        meds[[vars[i]]] <- median(data[[vars[i]]], na.rm = T)
    }
    tmp.df <- do.call(data.frame, lapply(meds, function(x) rep(x, 
        length(meds) * 2)))
    if (!is.null(typical.dat)) {
        notin <- which(!(names(typical.dat) %in% names(tmp.df)))
        if (length(notin) > 0) {
            cat("The following variables in typical.dat were not found in the prediction data: ", 
                names(typical.dat)[notin], "\n\n", sep = "")
            typical.dat <- typical.dat[, -notin]
        }
        for (j in 1:ncol(typical.dat)) {
            tmp.df[[names(typical.dat)[j]]] <- typical.dat[1, 
                j]
            meds[names(typical.dat)[j]] <- as.numeric(typical.dat[1, 
                j])
        }
    }
    inds <- seq(1, nrow(tmp.df), by = 2)
    for (j in 1:length(minmax)) {
        tmp.df[inds[j]:(inds[j] + 1), j] <- minmax[[j]]
    }
    preds <- matrix(predict(obj, newdata = tmp.df, type = type), 
        ncol = 2, byrow = T)
    diffs <- cbind(preds, apply(preds, 1, diff))
    colnames(diffs) <- c("min", "max", "diff")
    rownames(diffs) <- rn
    diffs <- diffs[which(rownames(diffs) %in% vars.type), ]
    minmax.mat <- do.call(cbind, minmax)
    minmax.mat <- rbind(c(unlist(meds)), minmax.mat)
    rownames(minmax.mat) <- c("typical", "min", "max")
    ret <- list(diffs = diffs, minmax = minmax.mat)
    class(ret) <- "change"
    return(ret)
}
cat2Table <- function(eff.obj, digits=2, rownames = NULL, colnames=NULL){
	print.digs <- paste("%.", digits, "f", sep="")
	cis <- paste("(", 
		sprintf(print.digs, eff.obj$lower), 
		",",
		sprintf(print.digs, eff.obj$upper), 
		")", sep="")
	cis.mat <- matrix(cis, 
		nrow = length(levels(eff.obj$x[[1]])),
		ncol = length(levels(eff.obj$x[[2]])))

	out.mat <- NULL
	est.mat <- matrix(sprintf(print.digs, eff.obj$fit),  
		nrow = length(levels(eff.obj$x[[1]])),
		ncol = length(levels(eff.obj$x[[2]])))
	for(i in 1:nrow(est.mat)){
		out.mat <- rbind(out.mat, est.mat[i,], cis.mat[i,])
	}
	if(is.null(colnames)){
		colnames(out.mat) <- levels(eff.obj$x[[2]])
	}
	if(is.null(rownames)){
		rn <- levels(eff.obj$x[[1]])
	}else{
		rn <- rownames
		}
		rn2 <- NULL
		for(i in 1:length(rn)){
			rn2 <- c(rn2, rn[i], paste(rep(" ", i), collapse=""))
		}
		rownames(out.mat) <- rn2
	out.mat
}
BGMtest <- function(obj, vars, digits = 3, level = 0.05, two.sided=T){
	cl <- attr(terms(obj), "dataClasses")
    if(any(cl[vars] != "numeric"))stop("Both variables in vars must be numeric")
    facs <- apply(attr(terms(obj), "factors"), 1, sum)
	if(any(facs[vars] != 2))stop("Each variable in vars must be involved in only 2 terms")
    r.x <- range(obj$model[, vars[1]])
	r.z <- range(obj$model[, vars[2]])
	b <- coef(obj)
	V <- vcov(obj)
	nb <- names(b)
	inds <- lapply(vars, function(x)grep(x, names(b)))
	A <- matrix(0, nrow=5, ncol=length(b))
    A[1:2, inds[[1]]] <- cbind(1, r.z)
	A[3:4, inds[[2]]] <- cbind(1, r.x)
    A[5, do.call(intersect, inds)] <- 1
	cond.b <- A %*% b        
	sp1 <- do.call(max, lapply(strsplit(gsub("-", "", as.character(cond.b)), split=".", fixed=T), function(x)nchar(x[1])))
	cond.se <- sqrt(diag(A %*% V %*% t(A)))
	sp2 <- do.call(max, lapply(strsplit(as.character(cond.se), split=".", fixed=T), function(x)nchar(x[1])))
	cond.t <- cond.b/cond.se
	sp3 <- do.call(max, lapply(strsplit(gsub("-", "", as.character(cond.t)), split=".", fixed=T), function(x)nchar(x[1])))
	cond.p <- (2^two.sided)*pt(abs(cond.t), obj$df.residual, lower.tail=F)
	pdigs1 <- paste("%", (sp1), ".", (digits), "f", sep="")     
	pdigs2 <- paste("%", (sp2), ".", (digits), "f", sep="")     
	pdigs3 <- paste("%", (sp3), ".", (digits), "f", sep="") 
	pdigs4 <- paste("%.", digits, "f", sep="") 
	rn <- c("P(X|Zmin)", "P(X|Zmax)", "P(Z|Xmin)", "P(Z|Xmax)", "P(XZ)")
	out <- cbind(sprintf(pdigs1, cond.b), sprintf(pdigs2, cond.se), 
		sprintf(pdigs3, cond.t), sprintf(pdigs4, cond.p))
	if(any(out[,1] < 0)){
		indp <- which(out[,1] > 0)
		out[indp,1] <- paste(" ", out[indp,1], sep="")
		out[indp,3] <- paste(" ", out[indp,3], sep="")
	}
	pad <- apply(out, 2, function(x)max(nchar(x)))
	nchars <- nchar(out)
	newchars <- t(apply(nchars, 1, function(x)pad-x))
	tmp <- apply(newchars, c(1,2), function(x)ifelse(x > 0, rep(" ", x), ""))
	out <- matrix(paste(tmp, out, sep=""), nrow=nrow(out), ncol=ncol(out))
	rownames(out) <- rn
	colnames(out) <- c(
		paste(paste(rep(" ", pad[1]-3), collapse=""), "est", sep=""), 
		paste(paste(rep(" ", pad[2]-2), collapse=""), "se", sep=""), 
		paste(paste(rep(" ", pad[3]-1), collapse=""), "t", sep=""), 
		ifelse(pad[4] > 6, 
			paste(paste(rep(" ", pad[4]-6), collapse=""), "p-value", collapse=""), 
			"p-value"))
	return(noquote(out))
}
intQualQuant <- function(obj, vars, level = .95 , 
	labs = NULL, n = 10 , onlySig = FALSE, plot = FALSE, 
	vals = NULL){
cl <- attr(terms(obj), "dataClasses")[vars]
if(length(cl) != 2){
	stop("vars must identify 2 and only 2 model terms")
}
if(!all(c("numeric", "factor") %in% cl)){
	stop("vars must have one numeric and one factor")
}
facvar <- names(cl)[which(cl == "factor")]
quantvar <- names(cl)[which(cl == "numeric")]
faclevs <- obj$xlevels[[facvar]]
if(is.null(labs)){
	labs <- faclevs
}
if(!is.null(vals)){
	quantseq <- vals
}
else{
	qrange <- range(obj$model[[quantvar]], na.rm=TRUE)
	quantseq <- seq(qrange[1], qrange[2], length=n)
}
b <- coef(obj)
faccoef <- paste(facvar, faclevs, sep="")
main.ind <- sapply(faccoef, function(x)
	grep(paste("^", x, "$", sep=""), names(b)))
main.ind <- sapply(main.ind, function(x)
	ifelse(length(x) == 0, 0, x))
int.ind1 <- sapply(faccoef, function(x){
	g1 <- grep(paste("[a-zA-Z0-9]*\\:", x, "$", sep=""), 
		names(b))
	ifelse(length(g1) == 0, 0, g1)
})
int.ind2 <- sapply(faccoef, function(x){
	g2 <- grep(paste("^", x, "\\:[a-zA-Z0-9]*", sep=""), 	
		names(b))
	ifelse(length(g2) == 0, 0, g2)
})
{if(sum(int.ind1) != 0){int.ind <- int.ind1}
else{int.ind <- int.ind2}}

inds <- cbind(main.ind, int.ind)
nc <- ncol(inds)
dn <- dimnames(inds)
if(!("matrix" %in% class(inds))){inds <- matrix(inds, nrow=1)}
outind <- which(main.ind == 0)
inds <- inds[-outind, ]
if(length(inds) == nc){
	inds <- matrix(inds, ncol=nc)
}
rownames(inds) <- dn[[1]][-outind]
colnames(inds) <- dn[[2]]


if(length(faclevs) < 2){stop("Factor must have at least two unique values")}
{if(length(faclevs) > 2){
combs <- combn(length(faclevs)-1, 2)
}
else{
	combs <- matrix(1:length(faclevs), ncol=1)
}}

tmp.A <- matrix(0, nrow=length(quantseq), ncol=length(b))
A.list <- list()
k <- 1
for(i in 1:nrow(inds)){
	A.list[[k]] <- tmp.A
	A.list[[k]][,inds[i,1]] <- 1
	A.list[[k]][,inds[i,2]] <- quantseq
	k <- k+1
}
if(nrow(inds) > 1){
for(i in 1:ncol(combs)){
	A.list[[k]] <- tmp.A
	A.list[[k]][,inds[combs[1,i], 1]] <- -1
	A.list[[k]][,inds[combs[2,i], 1]] <- 1
	A.list[[k]][,inds[combs[1,i], 2]] <- -quantseq
	A.list[[k]][,inds[combs[2,i], 2]] <- quantseq
	k <- k+1
}
}

effs <- lapply(A.list, function(x)x%*%b)
se.effs <- lapply(A.list, function(x)sqrt(diag(x %*% vcov(obj)%*%t(x))))
allcombs <- combn(length(faclevs), 2)
list.labs <- apply(rbind(labs[allcombs[2,]], 
	labs[allcombs[1,]]), 2, 
	function(x)paste(x, collapse=" - "))

names(A.list) <- list.labs
dat <- data.frame(
	fit = do.call("c", effs),
	se.fit = do.call("c", se.effs),
	x = rep(quantseq, length(A.list)),
	contrast = rep(names(A.list), each=n)
	)
level <- level + ((1-level)/2)
dat$lower <- dat$fit - qt(level, 	
	obj$df.residual)*dat$se.fit
dat$upper <- dat$fit + qt(level, 
	obj$df.residual)*dat$se.fit
res <- dat
if(onlySig){
	sigs <- do.call(rbind, 
		by(dat[,c("lower", "upper")], 
		list(dat$contrast), function(x)
		c(max(x[,1]), min(x[,2]))))
	notsig <- which(sigs[,1] < 0 & sigs[,2] > 0)
	res <- dat[-which(dat$contrast %in% names(notsig)), ]
}
if(plot == FALSE){
	invisible(res)
}
else{
	rl <- range(c(res[, c("lower", "upper")]))
	p <- xyplot(fit ~ x | contrast, data=res, xlab = quantvar, ylab = "Predicted Difference", ylim = rl, 
		lower=res$lower, upper=res$upper, 
		prepanel = prepanel.ci, 
		panel = panel.ci, zl=TRUE)
	plot(p)
	return(p)
}
}

panel.ci <- function(x,y,subscripts,lower,upper,zl){
	panel.lines(x,y,col="black")
	panel.lines(x,lower[subscripts], col="black", lty=2)
	panel.lines(x,upper[subscripts], col="black", lty=2)
	if(zl)panel.abline(h=0, lty=3, col="gray50")
}

prepanel.ci <- function(x,y,subscripts, lower,upper){
    x2 <- as.numeric(x)
    list(xlim = range(x2, finite = TRUE), 
         ylim = range(c(lower[subscripts], upper[subscripts]), 
            finite = TRUE), 
         dx = diff(range(x2, finite = TRUE)), 
         dy = diff(range(c(lower[subscripts], upper[subscripts]), 
            finite = TRUE)))
}
panel.2cat <- function(x,y,subscripts,lower,upper){
	panel.points(x,y, pch=16, col="black")
	panel.arrows(x, lower[subscripts], x, upper[subscripts], code=3, angle=90, length=.2)
}	
crTest <- function(model, adjust.method="none",...){
	cl <- attr(terms(model), "dataClasses")
	cl <- cl[which(cl != "factor")]
    terms <- predictor.names(model)
	terms <- intersect(terms, names(cl))
    if (any(attr(terms(model), "order") > 1)) {
        stop("C+R plots not available for models with interactions.")
    }
	terms.list <- list()
	orders <- sapply(terms, function(x)df.terms(model, x))
	for(i in 1:length(terms)){
    tmp.x <- {if (df.terms(model, terms[i]) > 1) predict(model, type = "terms", term = terms[i])
    else model.matrix(model)[, terms[i]]}
	if(!is.null(colnames(tmp.x))){colnames(tmp.x) <- "x"}
	terms.list[[i]] <- data.frame(x=tmp.x, 
		y = residuals.glm(model, "partial")[,terms[i]])
}
lo.mods <- lapply(terms.list, function(z)loess(y ~ x, data=z,...))
lin.mods <- lapply(terms.list, function(z)lm(y ~ x, data=z))
lo.rss <- sapply(lo.mods, function(x)sum(residuals(x)^2))
lm.rss <- sapply(lin.mods, function(x)sum(residuals(x)^2))
num.df <- sapply(lo.mods, function(x)x$trace.hat) - sapply(lin.mods, function(x)x$rank)
denom.df <- sapply(terms.list, nrow) - sapply(lo.mods, function(x)x$trace.hat)
F.stats <- ((lm.rss-lo.rss)/num.df)/(lo.rss/denom.df)
pvals <- p.adjust(pf(F.stats, num.df, denom.df, lower.tail=F), method=adjust.method)
out <- data.frame(
	RSSp = sprintf("%.2f", lm.rss), 
	RSSnp = sprintf("%.2f", lo.rss),
	DFnum = sprintf("%.3f", num.df), 
	DFdenom = sprintf("%.3f", denom.df), 
	F = sprintf("%.3f", F.stats), 
	p = sprintf("%.3f", pvals))
rownames(out) <- terms
return(out)
}
crSpanTest <- 
function(model, spfromto, n=10, adjust.method = "none", adjust.type = c("both", "across", "within", "none")){
	span.seq <- seq(from=spfromto[1], to=spfromto[2], 
		length=n)
	out.list <- list()
	for(i in 1:length(span.seq)){
		out.list[[i]] <- crTest(model, adjust.method="none", 
			span = span.seq[i])
	}
	pvals <- sapply(out.list, function(x)
		as.numeric(as.character(x[,"p"])))
	if(adjust.type == "within"){
		pvals <- apply(pvals, 2, p.adjust, 	
			method=adjust.method)
	}
	if(adjust.type == "across"){
		pvals <- t(apply(pvals, 1, p.adjust, 
			method=adjust.method))
	}
	if(adjust.type == "both"){
		pvals <- matrix(p.adjust(c(pvals), 
			method=adjust.method), 
			nrow=nrow(pvals))
	}
	rownames(pvals) <- predictor.names(model)
	ret = list(x=span.seq, y=t(pvals))
}
scaleDataFrame <- 
function(data){
classes <- sapply(1:ncol(data), function(x)class(data[,x]))
dummies <- apply(data, 2, function(x)prod(x %in% c(0,1,NA)))
nn.ind <- union(which(dummies == 1), which(classes == "factor"))
num.ind <- setdiff(1:ncol(data), nn.ind)
if(length(num.ind) == 0){stop("No Numeric Variables in Data Frame")}
num.dat <- data[,num.ind]
nonnum.dat <- data[,nn.ind]
if(is.null(dimnames(num.dat))){
	num.dat <- data.frame(num.dat)
	names(num.dat) <- names(data)[num.ind]
}
if(is.null(dimnames(nonnum.dat))){
	nonnum.dat <- data.frame(nonnum.dat)
	names(nonnum.dat) <- names(data)[nn.ind]
}
newdat <- cbind(as.data.frame(scale(num.dat)), as.data.frame(nonnum.dat))
newdat
}
outXT <- function(obj, count=TRUE, prop.r = TRUE, prop.c = TRUE, prop.t = TRUE, 
	col.marg=TRUE, row.marg=TRUE, digits = 3, type = "word", file=NULL){
	if(!(type %in% c("word", "latex"))){stop("type must be one of 'word' or 'latex'")}
	tmp.list <- list()
	k <- 1
	if(count){tmp.list[[k]] <- matrix(sprintf("%.0f", c(t(obj$t))), ncol=nrow(obj$t)); k <- k+1}
	if(prop.r){tmp.list[[k]] <- matrix(sprintf(paste("%.", digits, "f", sep=""), 
		c(t(obj$prop.r))), ncol=nrow(obj$prop.r)); k <- k+1}
	if(prop.c){tmp.list[[k]] <- matrix(sprintf(paste("%.", digits, "f", sep=""), 
		c(t(obj$prop.c))), ncol=nrow(obj$prop.c)); k <- k+1}
	if(prop.t){tmp.list[[k]] <- matrix(sprintf(paste("%.", digits, "f", sep=""), 
		c(t(obj$prop.t))), ncol=nrow(obj$prop.t)); k <- k+1}
	out <- do.call(rbind, tmp.list)
	mat <- matrix(c(out), ncol=nrow(tmp.list[[1]]), byrow=T)
	rownames(mat) <- NULL
	rn <- rep("", length=nrow(mat))
	rn[seq(1,nrow(mat), by=length(tmp.list))] <- rownames(obj$t)
	mat <- cbind(rn, mat)
	colnames(mat) <- c("", colnames(obj$t))
	if(col.marg){mat <- rbind(mat, c("Total", as.character(apply(obj$t, 2, sum))))}
	if(row.marg){
		rmarg <- rep("", nrow(mat))
		rmarg[seq(1, length(rmarg)-as.numeric(col.marg), by=length(tmp.list))] <- as.character(apply(obj$t, 1, sum))
		mat <- cbind(mat, rmarg)
		colnames(mat)[ncol(mat)] <- "Total"
	}
	if(row.marg & col.marg){mat[nrow(mat), ncol(mat)] <- sum(c(obj$t))}
	sink(file=ifelse(is.null(file), "tmp_xt.txt", file), type="output", append=F)
	print(xtable(mat), include.rownames=F, include.colnames=T)
	sink(file=NULL)
	rl <- readLines(ifelse(is.null(file), "tmp_xt.txt", file))
	if(type == "word"){
		rl <- rl[-c(1:6,8)]
		rl <- rl[-c((length(rl)-3):length(rl))]
		rl <- gsub("\\\\", "", rl, fixed=T)
		rl <- gsub(" & ", ",", rl, fixed=T)
		if(!is.null(file)){writeLines(rl, con=file)}
	}
	return(noquote(rl))
}

glmChange2 <-
function (obj, varname, data, change=c("unit", "sd"), R=1500) 
{
    vars <- names(attr(terms(obj), "dataClasses"))[-1]
	vars <- gsub("poly\\((.*?),.*?\\)", "\\1", vars)
	vars <- gsub("bs\\((.*?),.*?\\)", "\\1", vars)
	vars <- gsub("log\\((.*?),.*?\\)", "\\1", vars)
    rn <- vars
    var.classes <- sapply(vars, function(x) class(data[[x]]))
	b <- mvrnorm(R, coef(obj), vcov(obj))
	change <- match.arg(change)
	delt <- switch(change, 
		unit=1, 
		sd = sd(data[[varname]], na.rm=T))
	if(is.numeric(data[[varname]])){
		d0 <- d1 <- data
		d0[[varname]] <- d0[[varname]]-(.5*delt)
		d1[[varname]] <- d1[[varname]]+(.5*delt)
		X0 <- model.matrix(obj, data=d0)
		X1 <- model.matrix(obj, data=d1)
		p0 <- family(obj)$linkinv(X0 %*% t(b))
		p1 <- family(obj)$linkinv(X1 %*% t(b))
		diff <- p1-p0
		eff <- colMeans(diff)
		res <- matrix(c(mean(eff), quantile(eff, c(.025,.975))), nrow=1)
		colnames(res) <- c("mean", "lower", "upper")
		rownames(res) <- varname
	}
	if(!is.numeric(data[[varname]]) & length(unique(na.omit(data[[varname]]))) == 2){
		l <- obj$xlevels[[varname]]
		X0 <- X1 <- model.matrix(obj)
		X0[,grep(paste(varname, l[2], sep=""), colnames(X0))] <- 0
		X1[,grep(paste(varname, l[2], sep=""), colnames(X0))] <- 1
		p0 <- family(obj)$linkinv(X0 %*% t(b))
		p1 <- family(obj)$linkinv(X1 %*% t(b))
		diff <- p1-p0
		eff <- colMeans(diff)
		res <- matrix(c(mean(eff), quantile(eff, c(.025,.975))), nrow=1)
		colnames(res) <- c("mean", "lower", "upper")
		rownames(res) <- varname
	}
	if(!is.numeric(data[[varname]]) & length(unique(na.omit(data[[varname]]))) > 2){
		l <- obj$xlevels[[varname]]
		X.list <- list()
		for(j in 1:length(l)){
			X.list[[j]] <- model.matrix(obj)
		}
		l1 <- paste(varname, l, sep="")
		X.list[[1]][,which(colnames(X.list[[1]]) %in% l1)] <- 0
		for(j in 2:length(l1)){
			X.list[[j]][, which(colnames(X.list[[1]]) == l1[j])] <- 1
		}
		combs <- combn(length(X.list), 2)
		d.list <- list()
		for(j in 1:ncol(combs)){
			d.list[[j]] <- family(obj)$linkinv(X.list[[combs[2,j]]] %*% t(b))-family(obj)$linkinv(X.list[[combs[1,j]]] %*% t(b))
		}
		eff <- sapply(d.list, colMeans)
		res <- apply(eff, 2, function(x)c(mean(x), quantile(x, c(.025,.975))))
		cl <- array(l[combs], dim=dim(combs))
		rownames(res) <- c("mean", "lower", "upper")
		colnames(res) <- apply(cl[c(2,1), ], 2, paste, collapse="-")
		res <- t(res)
	}
	return(res)
}
aveEffPlot <- function (obj, varname, data, R=1500, nvals=25, plot=TRUE,...) 
{
    vars <- names(attr(terms(obj), "dataClasses"))[-1]
	vars <- gsub("poly\\((.*?),.*?\\)", "\\1", vars)
	vars <- gsub("bs\\((.*?),.*?\\)", "\\1", vars)
	vars <- gsub("log\\((.*?),.*?\\)", "\\1", vars)
    rn <- vars
    var.classes <- sapply(vars, function(x) class(data[[x]]))
	b <- mvrnorm(R, coef(obj), vcov(obj))
	if(is.numeric(data[[varname]])){
		s <- seq(min(data[[varname]], na.rm=T), max(data[[varname]], na.rm=T), length=nvals)
		dat.list <- list()
		for(i in 1:length(s)){
			dat.list[[i]] <- data
			dat.list[[i]][[varname]] <- s[i]
		}
		mm <- lapply(dat.list, function(x)model.matrix(obj, data=x))
		probs <- lapply(mm, function(x)family(obj)$linkinv(x %*% t(b)))
		cmprobs <- sapply(probs, colMeans)
		ciprobs <- t(apply(cmprobs, 2, function(x)c(mean(x), quantile(x, c(.025,.975)))))
		colnames(ciprobs) <- c("mean", "lower", "upper")
		ciprobs <- cbind(s, ciprobs)
		tmp <- as.data.frame(ciprobs)
		if(plot){
			pl <- xyplot(mean ~ s, data=tmp,  xlab=varname, ylab="Predicted Value", ...,
				lower=tmp$lower, upper=tmp$upper, 
				prepanel = prepanel.ci, 
				panel = function(x,y, lower, upper){
					panel.lines(x,y, col="black", lty=1)
					panel.lines(x, lower, col="black", lty=2)
					panel.lines(x, upper, col="black", lty=2)
			})
			return(pl)
		}
		else{
			return(tmp)
		}
	}
	if(!is.numeric(data[[varname]])){
		l <- obj$xlevels[[varname]]
		dat.list <- list()
		for(j in 1:length(l)){
			dat.list[[j]] <- model.matrix(obj)
		}
		l1 <- paste(varname, l, sep="")
		dat.list[[1]][,which(colnames(dat.list[[1]]) %in% l1)] <- 0
		for(j in 2:length(l1)){
			dat.list[[j]][, which(colnames(dat.list[[1]]) == l1[j])] <- 1
		}
		probs <- lapply(dat.list, function(x)family(obj)$linkinv(x %*% t(b)))
		cmprobs <- sapply(probs, colMeans)
		ciprobs <- t(apply(cmprobs, 2, function(x)c(mean(x), quantile(x, c(.025,.975)))))
		colnames(ciprobs) <- c("mean", "lower", "upper")
		tmp <- as.data.frame(ciprobs)
		tmp$s <- factor(1:length(l), labels=l)
		if(plot){
			pl <- xyplot(mean ~ s, data=tmp, xlab="", ylab="Predicted Value", ...,
				lower=tmp$lower, upper=tmp$upper, 
				prepanel = prepanel.ci, 
				scales=list(x=list(at=1:length(l), labels=l)), 
				panel = function(x,y, lower, upper){
					panel.points(x,y, col="black", lty=1, pch=16)
					panel.segments(x, lower, x, upper, lty=1, col="black")
			})
			return(pl)
		}
		else{
			return(tmp)
		}
	}
}