R/methods.R

In daniel-gerhard/mcprofile: Testing Generalized Linear Hypotheses for Generalized Linear Model Parameters by Profile Deviance

expit <- function(x) 1/(1 + exp(-x))

plot.mcpCI <-
function(x, ...){
  Estimate <- x$estimate
  CM <- x$CM
  grp <- factor(rownames(CM), levels=rownames(CM)[nrow(CM):1])
  cidat <- data.frame(grp, Estimate)
  switch(x$alternative,
         two.sided = {
           lower <- x$confint$lower
           upper <- x$confint$upper
         }, greater = {
           lower <- x$confint$lower
           upper <- Inf
         }, lower = {
           lower <- x$confint$lower
           upper <- Inf
         })
  lower[is.na(lower)] <- -Inf
  upper[is.na(upper)] <- Inf
  cidat$lower <- lower
  cidat$upper <- upper
  ggplot(cidat, aes(x=Estimate, y=grp, xmin=lower, xmax=upper)) + geom_errorbarh(height=0.3) + geom_point() + ylab("") + xlab("")
}

plot.mcprofile <- function(x, ...){
  sdlist <- x$srdp  
  CM <- x$CM
  nr <- sapply(sdlist, nrow)
  grp <- factor(rep(rownames(CM), nr), levels=rownames(CM))
  b1 <- unlist(lapply(sdlist, function(x) x[,1]))
  z1 <- unlist(lapply(sdlist, function(x) x[,2]))
  sdd <- data.frame(b=b1, z=z1, grp)

  spl <- lapply(sdlist, function(x){
    x <- na.omit(x)
    sp <- try(interpSpline(x[,1], x[,2]))
    bcoord <- seq(min(x[,1]), max(x[,1]), length=200)
    preds <- if (class(sp)[1] == "try-error") y <- rep(NA, 200) else predict(sp, x=bcoord)$y
    data.frame(b=bcoord, z=preds)
  })
  pgrp <- factor(rep(rownames(CM), each=200), levels=rownames(CM))
  b <- unlist(lapply(spl, function(x) x[,1]))
  z <- unlist(lapply(spl, function(x) x[,2]))
  spd <- data.frame(b, z, grp=pgrp)
  
  ggplot(spd, aes(x=b, y=z)) + geom_line() + facet_wrap(~grp, scales = "free_x") + geom_hline(aes(yintercept=0), linetype=2) + geom_point(data=sdd, shape="|")
}

print.mcpCI <-
function(x, ...){
  cat("\n   mcprofile - Confidence Intervals \n\n")
  cat("level: \t\t", x$level,"\n")
  cat("adjustment:\t", x$adjust,"\n\n")
  print(data.frame(Estimate=x$estimate, x$confint), digits=3)
  cat("\n")
}
print.mcprofile <-
function(x, ...){
  cat("\n   Multiple Contrast Profiles\n\n")
  est <- x$CM %*% coefficients(x$object)
  vest <- x$CM %*% vcov(x$object) %*% t(x$CM)
  sdest <- sqrt(diag(vest))
  if (is.null(x$adjestimates)){
    print(data.frame(Estimate=est, Std.err=sdest), digits=3)
  } else {
    print(data.frame(Estimate=est, Std.err=sdest, adjEstimate=x$adjestimates), digits=3)
  }
  cat("\n")
}
print.mcpSummary <-
function(x, ...){
  cat("\n   mcprofile - Multiple Testing\n\n")
  cat("Adjustment:\t", x$adjust, "\n")
  cat("Margin:  \t", x$margin, "\n")
  cat("Alternative:\t", x$alternative, "\n\n")
  dat <- data.frame(round(x$CM %*% x$estimate,2), round(x$statistic,2), round(x$p.values,3))
  pname <- switch(x$alternative, less = paste("Pr(<", ifelse(is.null(x$df), "z", "t"), ")", sep = ""), greater = paste("Pr(>", ifelse(is.null(x$df), "z", "t"), ")", sep = ""), two.sided = paste("Pr(>|", ifelse(is.null(x$df), "z", "t"), "|)", sep = ""))
  names(dat) <- c("Estimate", "Statistic", pname)
  alt <- switch(x$alternative, two.sided = "==", less = ">=", greater = "<=")
  rownames(dat) <- paste(rownames(dat), alt, x$margin)
  error <- attr(x$p.values, "error")
  if (!is.null(error) && error > .Machine$double.eps) {
    sig <- which.min(abs(1/error - (10^(1:10))))
    sig <- 1/(10^sig)
  } else {
    sig <- .Machine$double.eps
  }
  printCoefmat(dat, digits = max(3, getOption("digits") - 3), has.Pvalue = TRUE, P.values = TRUE, eps.Pvalue = sig)
  cat("\n")
}