# R/family.r In bigReg: Generalized Linear Models (GLM) for Large Data Sets

#### Documented in binomial_family_Gamma_gaussian_inverse.gaussian_poisson_quasi_quasibinomial_quasipoisson_

```#------------------------------------------------------------------------------------------------
# Family functions for blm and bglm big data functions
#------------------------------------------------------------------------------------------------

# The binomial family
#' @title binomial family function
#' @description binomial family function
#' @export
binomial_ <- function (link = "logit")
{
okLinks <- c("logit", "probit", "cloglog", "cauchit", "log")
}
else {
if (!is.null(stats\$name))
}
else {
domain = NA)
}
}
variance <- function(mu) mu * (1 - mu)
validmu <- function(mu) all(mu > 0) && all(mu < 1)
dev.resids <- binomial()\$dev.resids # function(y, mu, wt) .Call(stats:::C_binomial_dev_resids, y, mu, wt)
partialAIC <- function(y, n, mu, wt, dev)NA
finalAIC <- function(paic, n, dev)NA
initialize <- expression({
if (NCOL(y) == 1) {
if (is.factor(y)) y <- y != levels(y)[1L]
#n <- rep.int(1, nrow(y))
y[wghts == 0] <- 0
if (any(y < 0 | y > 1)) stop("y values must be 0 <= y <= 1")
mustart <- (wghts * y + 0.5)/(wghts + 1)
m <- wghts * y
if (any(abs(m - round(m)) > 0.001)) warning("non-integer #successes in a binomial glm!")
}
})
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")
if (!is.null(m <- object\$model)) {
y <- model.response(m)
if (is.factor(y)) {
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
}
partialAIC = partialAIC, finalAIC = finalAIC, mu.eta = stats\$mu.eta,
initialize = initialize, validmu = validmu, valideta = stats\$valideta,
simulate = simfun), class = "family")
}

# The poisson family
#' @title poisson family function
#' @description poisson family function
#' @export
poisson_ <- function (link = "log")
{
} else {
## what else shall we allow?  At least objects of class link-glm.
} else {
domain = NA)
}
}
variance <- 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))
partialAIC <- function(y, n, mu, wt, dev) -2*sum(dpois(y, mu, log=TRUE)*wt)
finalAIC <- function(paic, n, dev){
paic
}
initialize <- expression({
if (any(y < 0))
stop("negative values not allowed for the 'Poisson' family")
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",
variance = variance,
dev.resids = dev.resids,
partialAIC = partialAIC, finalAIC = finalAIC,
mu.eta = stats\$mu.eta,
initialize = initialize,
validmu = validmu,
valideta = stats\$valideta,
simulate = simfun),
class = "family")
}

# Gaussian family
#' @title gaussian family function
#' @description gaussian family function
#' @export
gaussian_ <- function (link = "identity")
{
} else {
## what else shall we allow?  At least objects of class link-glm.
} else {
domain = NA)
}
}
partialAIC <-  function(y, n, mu, wt, dev) {
sum(log(wt))
}
finalAIC <- function(paic, n, dev){
n*(log(dev/n*2*pi) + 1) + 2 - paic
}
structure(list(family = "gaussian",
variance = function(mu) rep.int(1, length(mu)),
dev.resids = function(y, mu, wt) wt * ((y - mu)^2),
partialAIC = partialAIC, finalAIC = finalAIC,
mu.eta = stats\$mu.eta,
initialize = expression({
n <- rep.int(1, nrow(y))
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")
}

#' @title quasipoisson family function
#' @description quasipoisson family function
#' @export
quasipoisson_ <- function (link = "log")
{
} else {
## what else shall we allow?  At least objects of class link-glm.
} else {
domain = NA)
}
}
variance <- 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))
partialAIC <- function(y, n, mu, wt, dev) NA
finalAIC <- function(paic, n, dev) NA
initialize <- expression({
if (any(y < 0))
stop("negative values not allowed for the 'Poisson' family")
mustart <- y + 0.1
})
structure(list(family = "quasipoisson",
variance = variance,
dev.resids = dev.resids,
partialAIC = partialAIC, finalAIC = finalAIC,
mu.eta = stats\$mu.eta,
initialize = initialize,
validmu = validmu,
valideta = stats\$valideta),
class = "family")
}

#' @title quasibinomial family function
#' @description quasibinomial family function
#' @export
quasibinomial_ <- function (link = "logit")
{
okLinks <- c("logit", "probit", "cloglog", "cauchit", "log")
} else {
## what else shall we allow?  At least objects of class link-glm.
} else {
domain = NA)
}
}
variance <- function(mu) mu * (1 - mu)
validmu <- function(mu) all(mu>0) && all(mu<1)
dev.resids <- function(y, mu, wt)
2 * wt * (y * log(ifelse(y == 0, 1, y/mu)) +
(1 - y) * log(ifelse(y == 1, 1, (1 - y)/(1 - mu))))
partialAIC <- function(y, n, mu, wt, dev) NA
finalAIC <- function(paic, n, dev) NA
initialize <- expression({
if (NCOL(y) == 1) {
if (is.factor(y)) y <- y != levels(y)[1L]
y[wghts == 0] <- 0
if (any(y < 0 | y > 1)) stop("y values must be 0 <= y <= 1")
mustart <- (wghts * y + 0.5)/(wghts + 1)
}
})
structure(list(family = "quasibinomial",
variance = variance,
dev.resids = dev.resids,
partialAIC = partialAIC, finalAIC = finalAIC,
mu.eta = stats\$mu.eta,
initialize = initialize,
validmu = validmu,
valideta = stats\$valideta),
class = "family")
}

#' @title Gamma family function
#' @description Gamma family function
#' @export
Gamma_ <- function (link = "inverse")
{
else {
## what else shall we allow?  At least objects of class link-glm.
} else {
domain = NA)
}
}
variance <- function(mu) mu^2
validmu <- function(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
}
partialAIC <- function(y, n, mu, wt, dev) sum(wt)
finalAIC <- function(paic, n, dev) NA
initialize <- expression({
if (any(y <= 0))
stop("non-positive values not allowed for the 'gamma' family")
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",
variance = variance,
dev.resids = dev.resids,
partialAIC = partialAIC,
finalAIC = finalAIC,
mu.eta = stats\$mu.eta,
initialize = initialize,
validmu = validmu,
valideta = stats\$valideta,
simulate = simfun),
class = "family")
}

#' @title inverse.gaussian family function
#' @description inverse.gaussian family function
#' @export
{
okLinks <- c("inverse", "log", "identity", "1/mu^2")
} else {
## what else shall we allow?  At least objects of class link-glm.
} else {
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
partialAIC <- function(y, n, mu, wt, dev) NA
finalAIC <- function(paic, n, dev) NA
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) NA

structure(list(family = "inverse.gaussian",
variance = variance,
dev.resids = dev.resids,
partialAIC = partialAIC, finalAIC = finalAIC,
mu.eta = stats\$mu.eta,
initialize = initialize,
validmu = validmu,
valideta = stats\$valideta,
simulate = simfun),
class = "family")
}

#' @title quasi family function
#' @description quasi family function
#' @param variance choice character
#' @export
quasi_ <- function(link = "identity", variance = "constant")
{
if (linktemp %in% c("logit", "probit", "cloglog", "identity",
"inverse", "log", "1/mu^2", "sqrt"))
} else {
}
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({
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)
2 * wt * (y * log(ifelse(y == 0, 1, y/mu)) +
(1 - y) * log(ifelse(y == 1, 1, (1 - y)/(1 - mu))))
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({
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({
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
partialAIC <- function(y, n, mu, wt, dev) NA
finalAIC <- function(paic, n, dev) NA
structure(list(family = "quasi",
variance = varfun,
dev.resids = dev.resids,
partialAIC = partialAIC, finalAIC = finalAIC,
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")
}

# The family function
#' @title family function
#' @description family function
#' @param distr distr character one of "binomial", "poisson", "gaussian", "quasipoisson",
#'              "quasibinomial", "Gamma", "inverse.gaussian", "quasi"
#' @export
{
switch(distr,