family.R

#  File src/library/stats/R/family.R
#  Part of the R package, https://www.R-project.org
#
#  Copyright (C) 1995-2012 The R Core Team
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  https://www.R-project.org/Licenses/

family <- function(object, ...) UseMethod("family")

print.family <- function(x, ...)
{
    cat("\nFamily:", x$family, "\n")
    cat("Link function:", x$link, "\n\n")
    invisible(x)
}

power <- function(lambda = 1)
{
    if(!is.numeric(lambda) || is.na(lambda))
        stop("invalid argument 'lambda'")
    if(lambda <= 0) return(make.link("log"))
    if(lambda == 1) return(make.link("identity"))
    linkfun <- function(mu) mu^lambda
    linkinv <- function(eta)
        pmax(eta^(1/lambda), .Machine$double.eps)
    mu.eta <- function(eta)
        pmax((1/lambda) * eta^(1/lambda - 1), .Machine$double.eps)
    valideta <- function(eta) all(is.finite(eta)) && all(eta>0)
    link <- paste0("mu^", round(lambda, 3))
    structure(list(linkfun = linkfun, linkinv = linkinv,
                   mu.eta = mu.eta, valideta = valideta, name = link),
              class="link-glm")
}

## Written by Simon Davies Dec 1995
## Modified by Thomas Lumley 26 Apr 97
## added valideta(eta) function..
make.link <- function (link)
{
    switch(link,
           "logit" = {
               linkfun <- function(mu) .Call(C_logit_link, mu)
               linkinv <- function(eta) .Call(C_logit_linkinv, eta)
               mu.eta <- function(eta) .Call(C_logit_mu_eta, eta)
               valideta <- function(eta) TRUE
           },
           "probit" = {
               linkfun <- function(mu) qnorm(mu)
               linkinv <- function(eta) {
                   thresh <- - qnorm(.Machine$double.eps)
                   eta <- pmin(pmax(eta, -thresh), thresh)
                   pnorm(eta)
               }
               mu.eta <- function(eta)
                   pmax(dnorm(eta),.Machine$double.eps)
               valideta <- function(eta) TRUE
           },
           "cauchit" = {
               linkfun <- function(mu) qcauchy(mu)
               linkinv <- function(eta) {
                   thresh <- -qcauchy(.Machine$double.eps)
                   eta <- pmin(pmax(eta, -thresh), thresh)
                   pcauchy(eta)
               }
               mu.eta <- function(eta)
                   pmax(dcauchy(eta), .Machine$double.eps)
               valideta <- function(eta) TRUE
           },
           "cloglog" = {
               linkfun <- function(mu) log(-log(1 - mu))
               linkinv <- function(eta)
                   pmax(pmin(-expm1(-exp(eta)), 1 - .Machine$double.eps),
                        .Machine$double.eps)
               mu.eta <- function(eta) {
                   eta <- pmin(eta, 700)
                   pmax(exp(eta) * exp(-exp(eta)), .Machine$double.eps)
               }
               valideta <- function(eta) TRUE
           },
           "identity" = {
               linkfun <- function(mu) mu
               linkinv <- function(eta) eta
               mu.eta <- function(eta) rep.int(1, length(eta))
               valideta <- function(eta) TRUE
           },
           "log" = {
               linkfun <- function(mu) log(mu)
               linkinv <- function(eta)
                   pmax(exp(eta), .Machine$double.eps)
               mu.eta <- function(eta)
                   pmax(exp(eta), .Machine$double.eps)
               valideta <- function(eta) TRUE
           },
           "sqrt" = {
               linkfun <- function(mu) sqrt(mu)
               linkinv <- function(eta) eta^2
               mu.eta <- function(eta) 2 * eta
               valideta <- function(eta) all(is.finite(eta)) && all(eta>0)
           },
           "1/mu^2" = {
               linkfun <- function(mu) 1/mu^2
               linkinv <- function(eta) 1/sqrt(eta)
               mu.eta <- function(eta) -1/(2 * eta^1.5)
               valideta <- function(eta) all(is.finite(eta)) && all(eta>0)
           },
           "inverse" = {
               linkfun <- function(mu) 1/mu
               linkinv <- function(eta) 1/eta
               mu.eta <- function(eta) -1/(eta^2)
               valideta <- function(eta) all(is.finite(eta)) && all(eta != 0)
           },
           ## else :
           stop(gettextf("%s link not recognised", sQuote(link)),
                domain = NA)
           )# end switch(.)
    environment(linkfun) <- environment(linkinv) <- environment(mu.eta) <-
        environment(valideta) <- asNamespace("stats")
    structure(list(linkfun = linkfun, linkinv = linkinv,
                   mu.eta = mu.eta, valideta = valideta, name = link),
              class="link-glm")
}

poisson <- function (link = "log")
{
    linktemp <- substitute(link)
    if (!is.character(linktemp)) linktemp <- deparse(linktemp)
    okLinks <- c("log", "identity", "sqrt")
    if (linktemp %in% okLinks)
	stats <- make.link(linktemp)
    else if (is.character(link)) {
        stats <- make.link(link)
        linktemp <- link
    } else {
        ## what else shall we allow?  At least objects of class link-glm.
        if(inherits(link, "link-glm")) {
            stats <- link
            if(!is.null(stats$name)) linktemp <- stats$name
        } else {
            stop(gettextf('link "%s" not available for poisson family; available links are %s',
			  linktemp, paste(sQuote(okLinks), collapse =", ")),
		 domain = NA)
        }
    }
    variance <- function(mu) mu
    validmu <- function(mu) all(is.finite(mu)) && all(mu>0)
    dev.resids <- function(y, mu, wt)
    { ## faster than  2 * wt * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu))
	r <- mu*wt
	p <- which(y > 0)
	r[p] <- (wt * (y*log(y/mu) - (y - mu)))[p]
	2*r
    }
    aic <- function(y, n, mu, wt, dev) -2*sum(dpois(y, mu, log=TRUE)*wt)
    initialize <- expression({
	if (any(y < 0))
	    stop("negative values not allowed for the 'Poisson' family")
	n <- rep.int(1, nobs)
	mustart <- y + 0.1
    })
    simfun <- function(object, nsim) {
        ## A Poisson GLM has dispersion fixed at 1, so prior weights
        ## do not have a simple unambiguous interpretation:
        ## they might be frequency weights or indicate averages.
        wts <- object$prior.weights
        if (any(wts != 1)) warning("ignoring prior weights")
        ftd <- fitted(object)
        rpois(nsim*length(ftd), ftd)
    }
    structure(list(family = "poisson",
		   link = linktemp,
		   linkfun = stats$linkfun,
		   linkinv = stats$linkinv,
		   variance = variance,
		   dev.resids = dev.resids,
		   aic = aic,
		   mu.eta = stats$mu.eta,
		   initialize = initialize,
		   validmu = validmu,
		   valideta = stats$valideta,
                   simulate = simfun),
	      class = "family")
}

quasipoisson <- function (link = "log")
{
    linktemp <- substitute(link)
    if (!is.character(linktemp)) linktemp <- deparse(linktemp)
    okLinks <- c("log", "identity", "sqrt")
    if (linktemp %in% okLinks)
        stats <- make.link(linktemp)
    else if (is.character(link)) {
        stats <- make.link(link)
        linktemp <- link
    } else {
        ## what else shall we allow?  At least objects of class link-glm.
        if(inherits(link, "link-glm")) {
            stats <- link
            if(!is.null(stats$name)) linktemp <- stats$name
        } else {
	    stop(gettextf('link "%s" not available for quasipoisson family; available links are %s',
			  linktemp, paste(sQuote(okLinks), collapse =", ")),
		 domain = NA)
        }
    }
    variance <- function(mu) mu
    validmu <- function(mu) all(is.finite(mu)) && all(mu>0)
    dev.resids <- function(y, mu, wt)
    { ## faster than  2 * wt * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu))
	r <- mu*wt
	p <- which(y > 0)
	r[p] <- (wt * (y*log(y/mu) - (y - mu)))[p]
	2*r
    }
    aic <- function(y, n, mu, wt, dev) NA
    initialize <- expression({
	if (any(y < 0))
	    stop("negative values not allowed for the 'quasiPoisson' family")
	n <- rep.int(1, nobs)
	mustart <- y + 0.1
    })
    structure(list(family = "quasipoisson",
		   link = linktemp,
		   linkfun = stats$linkfun,
		   linkinv = stats$linkinv,
		   variance = variance,
		   dev.resids = dev.resids,
		   aic = aic,
		   mu.eta = stats$mu.eta,
		   initialize = initialize,
		   validmu = validmu,
		   valideta = stats$valideta),
	      class = "family")
}

gaussian <- function (link = "identity")
{
    linktemp <- substitute(link)
    if (!is.character(linktemp)) linktemp <- deparse(linktemp)
    okLinks <- c("inverse", "log", "identity")
    if (linktemp %in% okLinks)
	stats <- make.link(linktemp)
    else if (is.character(link)) {
        stats <- make.link(link)
        linktemp <- link
    } else {
        ## what else shall we allow?  At least objects of class link-glm.
        if(inherits(link, "link-glm")) {
            stats <- link
            if(!is.null(stats$name)) linktemp <- stats$name
        } else {
	    stop(gettextf('link "%s" not available for gaussian family; available links are %s',
			  linktemp, paste(sQuote(okLinks), collapse =", ")),
		 domain = NA)
        }
    }
    structure(list(family = "gaussian",
		   link = linktemp,
		   linkfun = stats$linkfun,
		   linkinv = stats$linkinv,
		   variance = function(mu) rep.int(1, length(mu)),
		   dev.resids = function(y, mu, wt) wt * ((y - mu)^2),
		   aic =	function(y, n, mu, wt, dev) {
                       nobs <- length(y)
                       nobs*(log(dev/nobs*2*pi)+1)+2 - sum(log(wt))
                   },
		   mu.eta = stats$mu.eta,
		   initialize = expression({
		       n <- rep.int(1, nobs)
                       if(is.null(etastart) && is.null(start) &&
                          is.null(mustart) &&
                          ((family$link == "inverse" && any(y == 0)) ||
                          (family$link == "log" && any(y <= 0))))
                           stop("cannot find valid starting values: please specify some")

		       mustart <- y }),
		   validmu = function(mu) TRUE,
		   valideta = stats$valideta
		   ),
	      class = "family")
}

binomial <- function (link = "logit")
{
    linktemp <- substitute(link)
    if (!is.character(linktemp)) linktemp <- deparse(linktemp)
    okLinks <- c("logit", "probit", "cloglog", "cauchit", "log")
    if (linktemp %in% okLinks)
        stats <- make.link(linktemp)
    else if (is.character(link)) {
        stats <- make.link(link)
        linktemp <- link
    } else {
        ## what else shall we allow?  At least objects of class link-glm.
        if(inherits(link, "link-glm")) {
            stats <- link
            if(!is.null(stats$name)) linktemp <- stats$name
        } else {
	    stop(gettextf('link "%s" not available for binomial family; available links are %s',
			  linktemp, paste(sQuote(okLinks), collapse =", ")),
	     domain = NA)
        }
    }
    variance <- function(mu) mu * (1 - mu)
    validmu <- function(mu) all(is.finite(mu)) && all(mu>0 &mu<1)
    dev.resids <- function(y, mu, wt) .Call(C_binomial_dev_resids, y, mu, wt)
    aic <- function(y, n, mu, wt, dev) {
        m <- if(any(n > 1)) n else wt
	-2*sum(ifelse(m > 0, (wt/m), 0)*
               dbinom(round(m*y), round(m), mu, log=TRUE))
    }
    initialize <- expression({
	if (NCOL(y) == 1) {
	    ## allow factors as responses
	    ## added BDR 29/5/98
	    if (is.factor(y)) y <- y != levels(y)[1L]
	    n <- rep.int(1, nobs)
            ## anything, e.g. NA/NaN, for cases with zero weight is OK.
            y[weights == 0] <- 0
	    if (any(y < 0 | y > 1))
		stop("y values must be 0 <= y <= 1")
            mustart <- (weights * y + 0.5)/(weights + 1)
            m <- weights * y
            if(any(abs(m - round(m)) > 1e-3))
                warning("non-integer #successes in a binomial glm!")
	}
	else if (NCOL(y) == 2) {
            if(any(abs(y - round(y)) > 1e-3))
                warning("non-integer counts in a binomial glm!")
	    n <- y[, 1] + y[, 2]
	    y <- ifelse(n == 0, 0, y[, 1]/n)
	    weights <- weights * n
            mustart <- (n * y + 0.5)/(n + 1)
	}
	else stop("for the 'binomial' family, y must be a vector of 0 and 1\'s\nor a 2 column matrix where col 1 is no. successes and col 2 is no. failures")
    })
    simfun <- function(object, nsim) {
        ftd <- fitted(object)
        n <- length(ftd)
        ntot <- n*nsim
        wts <- object$prior.weights
        if (any(wts %% 1 != 0))
            stop("cannot simulate from non-integer prior.weights")
        ## Try to fathom out if the original data were
        ## proportions, a factor or a two-column matrix
        if (!is.null(m <- object$model)) {
            y <- model.response(m)
            if(is.factor(y)) {
                ## ignote weights
                yy <- factor(1+rbinom(ntot, size = 1, prob = ftd),
                             labels = levels(y))
                split(yy, rep(seq_len(nsim), each = n))
            } else if(is.matrix(y) && ncol(y) == 2) {
                yy <- vector("list", nsim)
                for (i in seq_len(nsim)) {
                    Y <- rbinom(n, size = wts, prob = ftd)
                    YY <- cbind(Y, wts - Y)
                    colnames(YY) <- colnames(y)
                    yy[[i]] <- YY
                }
                yy
            } else
            rbinom(ntot, size = wts, prob = ftd)/wts
        } else rbinom(ntot, size = wts, prob = ftd)/wts
    }
    structure(list(family = "binomial",
		   link = linktemp,
		   linkfun = stats$linkfun,
		   linkinv = stats$linkinv,
		   variance = variance,
		   dev.resids = dev.resids,
		   aic = aic,
		   mu.eta = stats$mu.eta,
		   initialize = initialize,
		   validmu = validmu,
		   valideta = stats$valideta,
                   simulate = simfun),
	      class = "family")
}

quasibinomial <- function (link = "logit")
{
    linktemp <- substitute(link)
    if (!is.character(linktemp)) linktemp <- deparse(linktemp)
    okLinks <- c("logit", "probit", "cloglog", "cauchit", "log")
    if (linktemp %in% okLinks)
        stats <- make.link(linktemp)
    else if (is.character(link)) {
        stats <- make.link(link)
        linktemp <- link
    } else {
        ## what else shall we allow?  At least objects of class link-glm.
        if(inherits(link, "link-glm")) {
            stats <- link
            if(!is.null(stats$name)) linktemp <- stats$name
        } else {
	    stop(gettextf('link "%s" not available for quasibinomial family; available links are %s',
			  linktemp, paste(sQuote(okLinks), collapse =", ")),
	     domain = NA)
        }
    }
    variance <- function(mu) mu * (1 - mu)
    validmu <- function(mu) all(is.finite(mu)) && all(mu>0 &mu<1)
    dev.resids <- function(y, mu, wt) .Call(C_binomial_dev_resids, y, mu, wt)
    aic <- function(y, n, mu, wt, dev) NA
    initialize <- expression({
	if (NCOL(y) == 1) {
	    if (is.factor(y)) y <- y != levels(y)[1L]
	    n <- rep.int(1, nobs)
	    if (any(y < 0 | y > 1))
		stop("y values must be 0 <= y <= 1")
            mustart <- (weights * y + 0.5)/(weights + 1)
	}
	else if (NCOL(y) == 2) {
	    n <- y[, 1] + y[, 2]
	    y <- ifelse(n == 0, 0, y[, 1]/n)
	    weights <- weights * n
            mustart <- (n * y + 0.5)/(n + 1)
	}
	else stop("for the 'quasibinomial' family, y must be a vector of 0 and 1\'s\nor a 2 column matrix where col 1 is no. successes and col 2 is no. failures")
    })
    structure(list(family = "quasibinomial",
		   link = linktemp,
		   linkfun = stats$linkfun,
		   linkinv = stats$linkinv,
		   variance = variance,
		   dev.resids = dev.resids,
		   aic = aic,
		   mu.eta = stats$mu.eta,
		   initialize = initialize,
		   validmu = validmu,
		   valideta = stats$valideta),
	      class = "family")
}

Gamma <- function (link = "inverse")
{
    linktemp <- substitute(link)
    if (!is.character(linktemp)) linktemp <- deparse(linktemp)
    okLinks <- c("inverse", "log", "identity")
    if (linktemp %in% okLinks)
	stats <- make.link(linktemp)
    else if(is.character(link)) stats <- make.link(link)
    else {
        ## what else shall we allow?  At least objects of class link-glm.
        if(inherits(link, "link-glm")) {
            stats <- link
            if(!is.null(stats$name)) linktemp <- stats$name
        } else {
	    stop(gettextf('link "%s" not available for gamma family; available links are %s',
			  linktemp, paste(sQuote(okLinks), collapse =", ")),
	     domain = NA)
        }
    }
    variance <- function(mu) mu^2
    validmu <- function(mu) all(is.finite(mu)) && all(mu>0)
    dev.resids <- function(y, mu, wt)
	-2 * wt * (log(ifelse(y == 0, 1, y/mu)) - (y - mu)/mu)
    aic <- function(y, n, mu, wt, dev){
	n <- sum(wt)
	disp <- dev/n
	-2*sum(dgamma(y, 1/disp, scale=mu*disp, log=TRUE)*wt) + 2
    }
    initialize <- expression({
	if (any(y <= 0))
	    stop("non-positive values not allowed for the 'gamma' family")
	n <- rep.int(1, nobs)
	mustart <- y
    })
    simfun <- function(object, nsim) {
        wts <- object$prior.weights
        if (any(wts != 1)) message("using weights as shape parameters")
        ftd <- fitted(object)
        shape <- MASS::gamma.shape(object)$alpha * wts
        rgamma(nsim*length(ftd), shape = shape, rate = shape/ftd)
    }
    structure(list(family = "Gamma",
		   link = linktemp,
		   linkfun = stats$linkfun,
		   linkinv = stats$linkinv,
		   variance = variance,
		   dev.resids = dev.resids,
		   aic = aic,
		   mu.eta = stats$mu.eta,
		   initialize = initialize,
		   validmu = validmu,
		   valideta = stats$valideta,
                   simulate = simfun),
	      class = "family")
}

inverse.gaussian <- function(link = "1/mu^2")
{
    linktemp <- substitute(link)
    if (!is.character(linktemp)) linktemp <- deparse(linktemp)
    okLinks <- c("inverse", "log", "identity", "1/mu^2")
    if (linktemp %in% okLinks)
	stats <- make.link(linktemp)
    else if (is.character(link)) {
        stats <- make.link(link)
        linktemp <- link
    } else {
        ## what else shall we allow?  At least objects of class link-glm.
        if(inherits(link, "link-glm")) {
            stats <- link
            if(!is.null(stats$name)) linktemp <- stats$name
        } else {
	    stop(gettextf('link "%s" not available for inverse.gaussian family; available links are %s',
			  linktemp, paste(sQuote(okLinks), collapse =", ")),
                 domain = NA)
        }
    }
    variance <- function(mu) mu^3
    dev.resids <- function(y, mu, wt)  wt*((y - mu)^2)/(y*mu^2)
    aic <- function(y, n, mu, wt, dev)
	sum(wt)*(log(dev/sum(wt)*2*pi)+1)+3*sum(log(y)*wt)+2
    initialize <- expression({
	if(any(y <= 0))
	    stop("positive values only are allowed for the 'inverse.gaussian' family")
	n <- rep.int(1, nobs)
	mustart <- y
    })
    validmu <- function(mu) TRUE
    simfun <- function(object, nsim) {
        if(!requireNamespace("SuppDists", quietly = TRUE))
            stop("need CRAN package 'SuppDists' for simulation from the 'inverse.gaussian' family")
        wts <- object$prior.weights
        if (any(wts != 1)) message("using weights as inverse variances")
        ftd <- fitted(object)
        SuppDists::rinvGauss(nsim * length(ftd), nu = ftd,
                             lambda = wts/summary(object)$dispersion)
    }

    structure(list(family = "inverse.gaussian",
		   link = linktemp,
		   linkfun = stats$linkfun,
		   linkinv = stats$linkinv,
		   variance = variance,
		   dev.resids = dev.resids,
		   aic = aic,
		   mu.eta = stats$mu.eta,
		   initialize = initialize,
		   validmu = validmu,
		   valideta = stats$valideta,
                   simulate = simfun),
	      class = "family")
}

quasi <- function (link = "identity", variance = "constant")
{
    linktemp <- substitute(link)
    if (!is.character(linktemp)) linktemp <- deparse(linktemp)
    if (linktemp %in% c("logit", "probit", "cloglog", "identity",
                        "inverse", "log", "1/mu^2", "sqrt"))
        stats <- make.link(linktemp)
    else if (is.character(link)) {
        stats <- make.link(link)
        linktemp <- link
    } else {
        stats <- link
        linktemp <- if(!is.null(stats$name)) stats$name else deparse(linktemp)
    }
    vtemp <- substitute(variance)
    if (!is.character(vtemp)) vtemp <- deparse(vtemp)
    variance_nm <- vtemp
    switch(vtemp,
           "constant" = {
               varfun <- function(mu) rep.int(1, length(mu))
               dev.resids <- function(y, mu, wt) wt * ((y - mu)^2)
                   validmu <- function(mu) TRUE
               initialize <- expression({n <- rep.int(1, nobs); mustart <- y})
           },
           "mu(1-mu)" = {
               varfun <- function(mu) mu * (1 - mu)
               validmu <- function(mu) all(mu>0) && all(mu<1)
               dev.resids <- function(y, mu, wt) .Call(C_binomial_dev_resids, y, mu, wt)
               initialize <- expression({n <- rep.int(1, nobs)
                                         mustart <- pmax(0.001, pmin(0.999, y))})
           },
           "mu" = {
               varfun <- function(mu) mu
               validmu <- function(mu) all(mu>0)
               dev.resids <- function(y, mu, wt)
                   2 * wt * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu))
               ## 0.1 fudge here matches poisson: S has 1/6.
               initialize <- expression({n <- rep.int(1, nobs)
                                         mustart <- y + 0.1 * (y == 0)})
           },
           "mu^2" = {
               varfun <- function(mu) mu^2
               validmu <- function(mu) all(mu>0)
               dev.resids <- function(y, mu, wt)
		   pmax(-2 * wt * (log(ifelse(y == 0, 1, y)/mu) - (y - mu)/mu), 0)
               initialize <- expression({n <- rep.int(1, nobs)
                                         mustart <- y + 0.1 * (y == 0)})
           },
           "mu^3" = {
               varfun <- function(mu) mu^3
               validmu <- function(mu) all(mu>0)
               dev.resids <- function(y, mu, wt)
                   wt * ((y - mu)^2)/(y * mu^2)
               initialize <- expression({n <- rep.int(1, nobs)
                                         mustart <- y + 0.1 * (y == 0)})
           },
           variance_nm <- NA
           )# end switch(.)

    if(is.na(variance_nm)) {
        if(is.character(variance))
            stop(gettextf('\'variance\' "%s" is invalid: possible values are "mu(1-mu)", "mu", "mu^2", "mu^3" and "constant"', variance_nm), domain = NA)
        ## so we really meant the object.
        varfun <- variance$varfun
        validmu <- variance$validmu
        dev.resids <- variance$dev.resids
        initialize <- variance$initialize
        variance_nm <- variance$name
    }
    aic <- function(y, n, mu, wt, dev) NA
    structure(list(family = "quasi",
		   link = linktemp,
		   linkfun = stats$linkfun,
		   linkinv = stats$linkinv,
		   variance = varfun,
		   dev.resids = dev.resids,
		   aic = aic,
		   mu.eta = stats$mu.eta,
		   initialize = initialize,
		   validmu = validmu,
		   valideta = stats$valideta,
                   ## character form of the var fun is needed for gee
                   varfun = variance_nm),
	      class = "family")
}