## this is a hack and could be made more object oriented!
#' @export
coeff.extract <- function (object) UseMethod("coeff.extract", object)
#' @export
coeff.extract.default <- function (object) {
## works for lm, glm, coxph, clm
tmp <- data.frame(beta = coef(object), se = sqrt(diag(vcov(object))))
colnames(tmp) <- c("Estimate", "Std. Error")
return(tmp)
}
#' @export
n.extract <- function (object) UseMethod("n.extract", object)
#' @export
n.extract.default <- function (object) {
## works for coxph, clm
return(n = tryCatch(object[["n"]], error = function(e) return(NA)))
}
#' @export
n.extract.lm <- function (object) {
## works for lm, glm; coerces logicals to 0/1 for sum
return(n = tryCatch(sum(!is.na(residuals(object))), error = function(e) return(NA)))
}
#' @export
test.extract <- function(object, partial = NULL) UseMethod("test.extract", object)
#' @export
test.extract.htest <- function(object, partial = NULL) {
## object with class "htest" is returned by chisq.test() and fisher.test()
stopifnot(is.null(partial))
ts <- if (object$method == "Fisher's Exact Test for Count Data") {
if (!is.null(object$or)) {
c(pval = object$p.value, effect = log(object$or), df = 1)
} else {
c(pval = object$p.value, effect = NA, df = NA) # !!!
}
} else if (object$method == "Pearson's Chi-squared test") {
c(pval = object$p.value, effect = NA, df = object$parameter)
} else if (object$method == "Welch Two Sample t-test") {
c(pval = object$p.value, effect = object$estimate[2] - object$estimate[1], df = 1)
} else if (object$method == "Pearson's product-moment correlation") {
c(pval = object$p.value, effect = object$estimate, df = 1)
} else {
c(pval = object$p.value, effect = NA, df = NA) # !!!
}
names(ts) <- c("pval", "effect", "df"); return(ts)
}
#' @export
test.extract.survdiff <- function(object, partial = NULL) {
stopifnot(is.null(partial))
df = length(object$n) - 1 # there must be a better way...
ts <- c(pval = pchisq(object$chisq, df = df, lower.tail = FALSE), effect = NA, df = df)
names(ts) <- c("pval", "effect", "df"); return(ts)
}
#' @export
test.extract.lm <- function(object, partial = NULL) {
objcoef <- coef(object)
if (is.null(partial) && length(objcoef) == 2 && names(objcoef)[1] == "(Intercept)") partial <- names(objcoef)[2]
ts <- if (!is.null(partial)) {
tryCatch({
objeffect <- objcoef[partial]
stopifnot(!is.na(objeffect))
objdrop <- drop1(object, partial, test = "Chisq")
c(pval = objdrop[["Pr(>Chi)"]][2],
effect = objeffect,
df = objdrop[["Df"]][2])},
error = function(e) return(c(pval = 1, effect = NA, df = 0)))
} else {
tryCatch({
objf <- summary(object)$fstatistic
c(pval = pf(objf["value"], objf["numdf"], objf["dendf"], lower.tail = FALSE),
effect = NA, df = objf["numdf"])},
error = function(e) return(c(pval = 1, effect = NA, df = 0)))
}
names(ts) <- c("pval", "effect", "df"); return(ts)
}
#' @export
test.extract.glm <- function(object, partial = NULL) {
objcoef <- coef(object)
if (is.null(partial) && length(objcoef) == 2 && names(objcoef)[1] == "(Intercept)") partial <- names(objcoef)[2]
ts <- if (!is.null(partial)) {
## same as for lm
tryCatch({
objeffect <- objcoef[partial]
stopifnot(!is.na(objeffect))
objdrop <- drop1(object, partial, test = "Chisq")
c(pval = objdrop[["Pr(>Chi)"]][2],
effect = objeffect,
df = objdrop[["Df"]][2])},
error = function(e) return(c(pval = 1, effect = NA, df = 0)))
} else {
tryCatch({
df <- object$df.null - object$df.residual
c(pval = pchisq(object$null.deviance - object$deviance, df = df, lower.tail = FALSE),
effect = NA, df = df)},
error = function(e) return(c(pval = 1, effect = NA, df = 0)))
}
names(ts) <- c("pval", "effect", "df"); return(ts)
}
#' @export
test.extract.coxph <- function(object, partial = NULL) {
objcoef <- coef(object)
if (is.null(partial) && length(objcoef) == 1) partial <- names(objcoef)[1]
ts <- if (!is.null(partial)) {
## same as for lm
objeffect <- objcoef[partial]
tryCatch({
stopifnot(!is.na(objeffect))
## suppressWarnings needed when object has only one term, which is dropped
objdrop <- suppressWarnings(drop1(object, partial, test = "Chisq"))
c(pval = objdrop[["Pr(>Chi)"]][2],
effect = objeffect,
df = objdrop[["Df"]][2])},
error = function(e) return(c(pval = 1, effect = NA, df = 0)))
} else {
## compute LRT for all terms
tryCatch({
df <- sum(!is.na(object$coefficients)) # there must be a better way...
c(pval = pchisq(2*(object$loglik[2] - object$loglik[1]), df = df, lower.tail = FALSE),
effect = NA, df = df)},
error = function(e) return(c(pval = 1, effect = NA, df = 0)))
}
names(ts) <- c("pval", "effect", "df"); return(ts)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.