Nothing
#####################################################################################
## Author: Daniel Sabanés Bové [daniel *.* sabanesbove *a*t* ifspm *.* uzh *.* ch]
## Project: hypergsplines
##
## Time-stamp: <[getFamily.R] by DSB Mit 09/03/2011 16:53 (CET)>
##
## Description:
## Helper for glmModelData which extracts an S3 family object with additional elements
## from an option.
##
## History:
## 09/03/2011 copy from glmBfp
#####################################################################################
##' Helper function for glmModelData: Extracts an S3 family object
##'
##' Extracts an S3 family object, which (at least) the usual elements
##' "family", "link", "linkfun", "linkinv", "variance", "mu.eta", "dev.resids",
##' plus the additional elements:
##' "phi" which includes just the dispersion parameter to be used,
##' and "simulate" which can generate random variates for a given linear predictor and weight vector
##' (this function of course also uses the "phi" value)
##'
##' @param family the family argument passed to \code{\link{glmModelData}}
##' @param phi the dispersion argument passed to \code{\link{glmModelData}}
##'
##' @return The returned family object also includes a custom \sQuote{init} function,
##' which takes the response vector (or matrix) (\sQuote{y}) and the corresponding weight vector
##' (\sQuote{weights}), processes them to response vector \sQuote{y} and possibly altered weights
##' \sQuote{weights}, and includes starting values \sQuote{mustart} for the IWLS algorithm.
##' For example, here the binomial special case of two-column-response matrix is treated exactly in
##' the same way as in \code{\link{glm}}.
##'
##' @keywords internal
getFamily <- function(family,
phi)
{
## check that dispersion is positive etc.
phi <- as.numeric(phi)
stopifnot(identical(length(phi),
1L),
phi > 0)
## then process the family argument:
## convert character to function
if (is.character(family))
family <- get(family, mode = "function")
## call function to get object
if (is.function(family))
family <- family()
## check that family has correct form
stopifnot(inherits(family, "family"),
all(c("family", "link", "linkfun", "linkinv", "variance", "mu.eta", "dev.resids")
%in% names(family)))
## check here if family and link are supported.
## (we only want to use links which always map into the valid mu space)
## first check the family,
## and along the way, we force the dispersion parameter
## and also choose which links are supported for that family.
if(identical(family$family, "binomial"))
{
okLinks <- c("logit", "cauchit", "probit", "cloglog")
family$phi <- 1
family$simulate <- function(eta, weights)
{
## check if weights are non-negative integers
stopifnot(all((weights %% 1L) == 0L),
all(weights >= 0L))
## then simulate one observation for each lin pred value
rbinom(n=length(eta),
size=weights,
prob=family$linkinv(eta))
}
}
else if(identical(family$family, "gaussian"))
{
okLinks <- c("log", "identity")
family$phi <- phi
family$simulate <- function(eta, weights)
{
rnorm(n=length(eta),
mean=family$linkinv(eta),
sd=sqrt(phi / weights))
}
}
else if(identical(family$family, "poisson"))
{
okLinks <- c("log")
family$phi <- 1
family$simulate <- function(eta, weights)
{
## here we do not use the weights at all!
rpois(n=length(eta),
lambda=family$linkinv(eta))
}
}
else
{
stop(simpleError(paste("family",
family$family,
"not supported")))
}
## warn if phi has changed
if(family$phi != phi)
warning(simpleWarning(paste("dispersion was changed from",
phi,
"to",
family$phi)))
## then check the link
if(! (family$link %in% okLinks))
{
stop(simpleError(paste("link",
family$link,
"not supported for family",
family$family)))
}
## now attach a nicer initialize function, instead of just the expression e.g. in
## binomial()$initialize
family$init <- function(y, # the response vector or matrix
weights, # the weights vector
nobs=length(weights)) # number of observations
{
## create dummy variables which are sometimes needed in the
## initializer expression
etastart <- start <- mustart <- NULL
## evaluate the initializer expression, which modifies y and weights
## and computes mustart values
eval(family$initialize)
## return these
return(list(y=y,
weights=weights,
linPredStart=family$linkfun(mustart)))
## use start values on the linear predictor scale for ease and clarity in C++.
}
return(family)
}
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.