Nothing
#' Enrich objects of class \code{\link{family}}
#'
#' Enrich objects of class \code{\link{family}} with family-specific
#' mathematical functions
#'
#' @param object an object of class \code{\link{family}}
#' @param with a character vector with enrichment options for \code{object}
#' @param ... extra arguments to be passed to the \code{compute_*}
#' functions
#'
#' @details
#'
#' \code{\link[stats]{family}} objects specify characteristics of the
#' models used by functions such as \code{\link[stats]{glm}}. The
#' families implemented in the \code{stats} package include
#' \code{\link[stats]{binomial}}, \code{\link[stats]{gaussian}},
#' \code{\link[stats]{Gamma}}, \code{\link[stats]{inverse.gaussian}},
#' and \code{\link[stats]{poisson}}, which are all special cases of
#' the exponential family of distributions that have probability mass
#' or density function of the form \deqn{f(y; \theta, \phi) =
#' \exp\left\{\frac{y\theta - b(\theta) - c_1(y)}{\phi/m} -
#' \frac{1}{2}a\left(-\frac{m}{\phi}\right) + c_2(y)\right\} \quad y
#' \in Y \subset \Re\,, \theta \in \Theta \subset \Re\, , \phi >
#' 0}{f(y, theta, phi) = exp((y * theta - b(theta) - c_1(y))/(phi/m) -
#' a(-m/phi)/2 + c_2(y))} where \eqn{m > 0}{m > 0} is an observation
#' weight, and \eqn{a(.)}{a(.)}, \eqn{b(.)}{b(.)},
#' \eqn{c_1(.)}{c_1(.)} and \eqn{c_2(.)}{c_2(.)} are sufficiently
#' smooth, real-valued functions.
#'
#' The current implementation of \code{\link[stats]{family}} objects
#' includes the variance function (\code{variance}), the deviance
#' residuals (\code{dev.resids}), and the Akaike information criterion
#' (\code{aic}). See, also \code{\link{family}}.
#'
#' The \code{enrich} method can further enrich exponential
#' \code{\link{family}} distributions with \eqn{\theta}{theta} in
#' terms of \eqn{\mu}{mu} (\code{theta}), the functions
#' \eqn{b(\theta)}{b(theta)} (\code{bfun}), \eqn{c_1(y)}{c_1(y)}
#' (\code{c1fun}), \eqn{c_2(y)}{c_2(y)} (\code{c2fun}),
#' \eqn{a(\zeta)}{a(zeta)} (\code{fun}), the first two derivatives of
#' \eqn{V(\mu)}{V(mu)} (\code{d1variance} and \code{d2variance},
#' respectively), and the first four derivatives of
#' \eqn{a(\zeta)}{a(zeta)} (\code{d1afun}, \code{d2afun},
#' \code{d3afun}, \code{d4afun}, respectively).
#'
#' Corresponding enrichment options are also avaialble for
#' \code{\link[stats]{quasibinomial}},
#' \code{\link[stats]{quasipoisson}} and \code{\link[gnm]{wedderburn}}
#' families.
#'
#' The \code{\link[stats]{quasi}} families are enriched with
#' \code{d1variance} and \code{d2variance}.
#'
#' See \code{\link{enrich.link-glm}} for the enrichment of
#' \code{\link[=make.link]{link-glm}} objects.
#'
#' @return The object \code{object} of class \code{\link{family}} with
#' extra components. \code{get_enrichment_options.family()}
#' returns the components and their descriptions.
#'
#' @seealso \code{\link{enrich.link-glm}}, \code{\link{make.link}}
#'
#' @export
#' @examples
#'
#' ## An example from ?glm to illustrate that things still work with
#' ## enriched families
#' counts <- c(18,17,15,20,10,20,25,13,12)
#' outcome <- gl(3,1,9)
#' treatment <- gl(3,3)
#' print(d.AD <- data.frame(treatment, outcome, counts))
#' glm.D93 <- glm(counts ~ outcome + treatment, family = enrich(poisson()))
#' anova(glm.D93)
#' summary(glm.D93)
`enrich.family` <- function(object, with = "all", ...) {
if (is.null(with)) {
return(object)
}
dots <- list(...)
enrichment_options <- get_enrichment_options(object, option = with)
component <- unlist(enrichment_options$component)
compute <- unlist(enrichment_options$compute_function)
for (j in seq.int(length(component))) {
ccall <- call(compute[j], object = object)
for (nam in names(dots)) {
ccall[[nam]] <- dots[[nam]]
}
object[[component[j]]] <- eval(ccall)
}
if (is.null(attr(object, "enriched"))) {
attr(object, "enriched") <- TRUE
classes <- class(object)
class(object) <- c(paste0("enriched_", classes[1]), classes)
}
object
}
#' Available options for the enrichment objects of class family
#'
#' @param object the object to be enriched
#' @param option a character vector listing the options for enriching
#' the object
#' @param all_options if \code{TRUE} then output a data frame with the
#' available enrichment options, their descriptions, the names of
#' the components that each option results in, and the names of
#' the corresponding \code{compute_*} functions.
#' @return an object of class \code{enrichment_options}
#'
#' @details A check is being made whether the requested option is
#' available. No check is being made on whether the functions that
#' produce the components exist.
#' @examples
#' \dontrun{
#' get_enrichment_options.family(option = "all")
#' get_enrichment_options.family(all_options = TRUE)
#' }
#' @export
`get_enrichment_options.family` <- function(object, option, all_options = missing(option)) {
## List the enrichment options that you would like to make
## available for objects of class
out <- list()
out$option <- c('theta', 'bfun', 'c1fun', 'c2fun', 'd1variance', 'd2variance', 'afun', 'd1afun', 'd2afun', 'd3afun', 'd4fun', 'variance derivatives', 'function a derivatives')
## Provide the descriptions of the enrichment options
out$description <- c('natural parameter', 'cumulant transform', 'c1 function', 'c2 function', '1st derivative of the variance function', '2nd derivative of the variance function', 'a function', '1st derivative of the a function', '2nd derivative of the a function', '3rd derivative of the a function', '4th derivative of the a function', '1st and 2nd derivative of the variance function', '1st, 2nd and 3rd derivative of the a function')
## Add all as an option
out$option <- c(out$option, 'all')
out$description <- c(out$description, 'all available options')
out$component <- list('theta', 'bfun', 'c1fun', 'c2fun', 'd1variance', 'd2variance', 'afun', 'd1afun', 'd2afun', 'd3afun', 'd4afun', c('d1variance', 'd2variance'), c('d1afun', 'd2afun', 'd3afun'))
out$component[[length(out$component) + 1]] <- unique(unlist(out$component))
names(out$component) <- names(out$description) <- out$option
out$compute_function <- lapply(out$component, function(z) paste0('compute_', z))
class(out) <- 'enrichment_options'
if (all_options) {
return(out)
}
invalid_options <- !(option %in% out$option)
if (any(invalid_options)) {
stop(gettextf('some options have not been implemented: %s', paste0('"', paste(option[invalid_options], collapse = ', '), '"')))
}
out <- list(option = option,
description = out$description[option],
component = out$component[option],
compute_function = out$compute_function[option])
class(out) <- 'enrichment_options'
out
}
`compute_theta.family` <- function(object, ...) {
family <- object$family
switch(family,
"poisson" = function(mu) {
log(mu)
},
"quasipoisson" = function(mu) {
log(mu)
},
"gaussian" = function(mu) {
mu
},
"binomial" = function(mu) {
log(mu/(1 - mu))
},
"quasibinomial" = function(mu) {
log(mu/(1 - mu))
},
"wedderburn" = function(mu) {
log(mu/(1 - mu))
},
"Gamma" = function(mu) {
-1/mu
},
"inverse.gaussian" = function(mu) {
-1/2*mu^2
})
}
`compute_theta` <- function(object, ...) {
UseMethod('compute_theta')
}
`compute_bfun.family` <- function(object, ...) {
family <- object$family
switch(family,
"poisson" = function(theta) {
exp(theta)
},
"quasipoisson" = function(theta) {
exp(theta)
},
"gaussian" = function(theta) {
theta^2/2
},
"binomial" = function(theta) {
log(1 + exp(theta))
},
"quasibinomial" = function(theta) {
log(1 + exp(theta))
},
"wedderburn" = function(theta) {
log(1 + exp(theta))
},
"Gamma" = function(theta) {
-log(-theta)
},
"inverse.gaussian" = function(theta) {
-sqrt(-2*theta)
})
}
`compute_bfun` <- function(object, ...) {
UseMethod('compute_bfun')
}
`compute_c1fun.family` <- function(object, ...) {
family <- object$family
switch(family,
"poisson" = function(y) {
0
},
"quasipoisson" = function(y) {
0
},
"gaussian" = function(y) {
y^2/2
},
"binomial" = function(y) {
0
},
"quasibinomial" = function(y) {
0
},
"wedderburn" = function(y) {
0
},
"Gamma" = function(y) {
-log(y)
},
"inverse.gaussian" = function(y) {
1/(2*y)
})
}
`compute_c1fun` <- function(object, ...) {
UseMethod('compute_c1fun')
}
`compute_c2fun.family` <- function(object, ...) {
family <- object$family
switch(family,
"poisson" = function(y) {
-log(gamma(y+1))
},
"quasipoisson" = function(y) {
-log(gamma(y+1))
},
"gaussian" = function(y) {
-log(2*pi)/2
},
"binomial" = function(y, m) {
log(choose(m, m * y))
},
"quasibinomial" = function(y, m) {
log(choose(m, m * y))
},
"wedderburn" = function(y, m) {
log(choose(m, m * y))
},
"Gamma" = function(y) {
-log(y)
},
"inverse.gaussian" = function(y) {
-0.5 * log(2 * pi * y^3)
})
}
`compute_c2fun` <- function(object, ...) {
UseMethod('compute_c2fun')
}
`compute_d1variance.family` <- function(object, ...) {
family <- object$family
switch(family,
"poisson" = function(mu) {
rep(1, length(mu))
},
"quasipoisson" = function(mu) {
rep(1, length(mu))
},
"gaussian" = function(mu) {
rep(0, length(mu))
},
"binomial" = function(mu) {
1 - 2 * mu
},
"quasibinomial" = function(mu) {
1 - 2 * mu
},
"wedderburn" = function(mu) {
2 * mu * (1 - mu) * (1 - 2 * mu)
},
"Gamma" = function(mu) {
2 * mu
},
"inverse.gaussian" = function(mu) {
3 * mu^2
},
"quasi" = switch(object$varfun,
"constant" = function(mu) {
rep(0, length(mu))
},
"mu(1-mu)" = function(mu) {
1 - 2 * mu
},
"mu" = function(mu) {
rep(1, length(mu))
},
"mu^2" = function(mu) {
2 * mu
},
"mu^3" = function(mu) {
3 * mu^2
},
stop(sQuote(object$varfun), " variance function not supported")),
stop(sQuote(family), " family not supported"))
}
`compute_d1variance` <- function(object, ...) {
UseMethod('compute_d1variance')
}
`compute_d2variance.family` <- function(object, ...) {
family <- object$family
switch(family,
"poisson" = function(mu) {
rep(0, length(mu))
},
"quasipoisson" = function(mu) {
rep(0, length(mu))
},
"gaussian" = function(mu) {
rep(0, length(mu))
},
"binomial" = function(mu) {
rep(-2, length(mu))
},
"quasibinomial" = function(mu) {
rep(-2, length(mu))
},
"wedderburn" = function(mu) {
2 - 12 * mu * ( 1- mu)
},
"Gamma" = function(mu) {
rep(2, length(mu))
},
"inverse.gaussian" = function(mu) {
6 * mu
},
"quasi" = switch(object$varfun,
"constant" = function(mu) {
rep(0, length(mu))
},
"mu(1-mu)" = function(mu) {
rep(-2, length(mu))
},
"mu" = function(mu) {
rep(0, length(mu))
},
"mu^2" = function(mu) {
rep(2, length(mu))
},
"mu^3" = function(mu) {
6 * mu
},
stop(sQuote(object$varfun), " variance function not supported")),
stop(sQuote(family), " family not supported"))
}
`compute_d2variance` <- function(object, ...) {
UseMethod('compute_d2variance')
}
`compute_afun.family` <- function(object, ...) {
family <- object$family
switch(family,
"gaussian" = function(zeta) {
-log(-zeta)
},
"Gamma" = function(zeta) {
2*log(gamma(-zeta)) + 2*zeta*log(-zeta)
},
"inverse.gaussian" = function(zeta) {
-log(-zeta)
},
"poisson" = function(zeta) {
0
},
"quasipoisson" = function(zeta) {
0
},
"binomial" = function(zeta) {
0
},
"quasibinomial" = function(zeta) {
0
},
"wedderburn" = function(zeta) {
0
})
}
`compute_afun` <- function(object, ...) {
UseMethod('compute_afun')
}
`compute_d1afun.family` <- function(object, ...) {
family <- object$family
switch(family,
"gaussian" = function(zeta) {
-1/zeta
},
"Gamma" = function(zeta) {
-2*psigamma(-zeta, 0) + 2*log(-zeta) + 2
## this is the expectation of dev.resids + 2, because of the way dev.resids is implemented
},
"inverse.gaussian" = function(zeta) {
-1/zeta
},
"poisson" = function(zeta) {
0
},
"quasipoisson" = function(zeta) {
0
},
"binomial" = function(zeta) {
0
},
"quasibinomial" = function(zeta) {
0
},
"wedderburn" = function(zeta) {
0
})
}
`compute_d1afun` <- function(object, ...) {
UseMethod('compute_d1afun')
}
`compute_d2afun.family` <- function(object, ...) {
family <- object$family
switch(family,
"gaussian" = function(zeta) {
1/zeta^2
},
"Gamma" = function(zeta) {
2*psigamma(-zeta, 1) + 2/zeta
},
"inverse.gaussian" = function(zeta) {
1/zeta^2
},
"poisson" = function(zeta) {
0
},
"quasipoisson" = function(zeta) {
0
},
"binomial" = function(zeta) {
0
},
"quasibinomial" = function(zeta) {
0
},
"wedderburn" = function(zeta) {
0
})
}
`compute_d2afun` <- function(object, ...) {
UseMethod('compute_d2afun')
}
`compute_d3afun.family` <- function(object, ...) {
family <- object$family
switch(family,
"gaussian" = function(zeta) {
-2/zeta^3
},
"Gamma" = function(zeta) {
-2*psigamma(-zeta, 2) - 2/zeta^2
},
"inverse.gaussian" = function(zeta) {
-2/zeta^3
},
"poisson" = function(zeta) {
0
},
"quasipoisson" = function(zeta) {
0
},
"binomial" = function(zeta) {
0
},
"quasibinomial" = function(zeta) {
0
},
"wedderburn" = function(zeta) {
0
})
}
`compute_d3afun` <- function(object, ...) {
UseMethod('compute_d3afun')
}
`compute_d4afun.family` <- function(object, ...) {
family <- object$family
switch(family,
"gaussian" = function(zeta) {
6/zeta^4
},
"Gamma" = function(zeta) {
2 * psigamma(-zeta, 3) + 4/zeta^3
},
"inverse.gaussian" = function(zeta) {
6/zeta^4
},
"poisson" = function(zeta) {
0
},
"quasipoisson" = function(zeta) {
0
},
"binomial" = function(zeta) {
0
},
"quasibinomial" = function(zeta) {
0
},
"wedderburn" = function(zeta) {
0
})
}
`compute_d4afun` <- function(object, ...) {
UseMethod('compute_d4afun')
}
## ## Call that produced the initial enrichwith template for the current script:
## create_enrichwith_skeleton(class = "family", option = c("d1variance",
## "d2variance", "d1afun", "d2afun", "d3afun", "variance derivatives",
## "function a derivatives"), description = c("1st derivative of the variance function",
## "2nd derivative of the variance function", "1st derivative of the a function",
## "2nd derivative of the a function", "3rd derivative of the a function",
## "1st and 2nd derivative of the variance function", "1st, 2nd and 3rd derivative of the a function"),
## component = list("d1variance", "d2variance", "d1afun", "d2afun",
## "d3afun", c("d1variance", "d2variance"), c("d1afun",
## "d2afun", "d3afun")), path = "~/Downloads", attempt_rename = TRUE)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.