modifications {brglm}

Additive Modifications to the Binomial Responses and Totals for Use within ‘brglm.fit’

Description

Get, test and set the functions that calculate the additive modifications to the responses and totals in binomial-response GLMs, for the application of bias-reduction either via modified scores or via maximum penalized likelihood (where penalization is by Jeffreys invariant prior).

Usage

modifications(family, pl = FALSE)

Arguments

family

a family object of the form binomial(link = "link"), where "link" can be one of "logit", "probit", "cloglog" and "cauchit". The usual ways of giving the family name are supported (see family).

pl

logical determining whether the function returned corresponds to modifications for the penalized maximum likelihood approach or for the modified-scores approach to bias-reduction. Default value is FALSE.

Details

The function returned from modifications accepts the argument p which are the binomial probabilities and returns a list with components ar and at, which are the link-dependent parts of the additive modifications to the actual responses and totals, respectively.

Since the resulting function is used in brglm.fit, for efficiency reasons no check is made for p >= 0 | p <= 1, for length(at) == length(p) and for length(ap) == length(p).

Construction of custom pseudo-data representations

If y^* are the pseudo-responses (pseudo-counts) and m^* are the pseudo-totals then we call the pair (y^*, m^*) a pseudo-data representation. Both the modified-scores approach and the maximum penalized likelihood have a common property:

there exists (y^*, m^*) such that if we replace the actual data (y, m) with (y^*, m^*) in the expression for the ordinary scores (first derivatives of the likelihood) of a binomial-response GLM, then we end-up either with the modified-scores or with the derivatives of the penalized likelihood (see Kosmidis, 2007, Chapter 5).

Let \mu be the mean of the binomial response y (i.e. \mu=mp, where p is the binomial probability corresponding to the count y). Also, let d and d' denote the first and the second derivatives, respectively, of \mu with respect to the linear predictor \eta of the model. All the above are viewed as functions of p. The pseudo-data representations have the generic form

pseudo-response :	y^*=y + h a_r(p)
pseudo-totals :	m^*=m + h a_t(p),

where h is the leverage corresponding to y. The general expressions for a_r(p) ("r" for "response") and a_t(p) ("t" for "totals") are:

modified-scores approach

a_r(p) = d'(p)/(2w(p))

a_t(p) = 0,

maximum penalized likelihood approach

a_r(p) = d'(p)/w(p) + p - 0.5

a_t(p) = 0.

For supplying (y^*, m^*) in glm.fit (as is done by brglm.fit), an essential requirement for the pseudo-data representation is that it should mimic the behaviour of the original responses and totals, i.e. 0 \le y^* \le m^*. Since h \in [0, 1], the requirement translates to 0 \le a_r(p) \le a_t(p) for every p \in (0, 1). However, the above definitions of a_r(p) and a_t(p) do not necessarily respect this requirement.

On the other hand, the pair (a_r(p), a_t(p)) is not unique in the sense that for a given link function and once the link-specific structure of the pair has been extrapolated, there is a class of equivalent pairs that can be resulted following only the following two rules:

• add and subtract the same quantity from either a_r(p) or a_t(p).

• if a quantity is to be moved from a_r(p) to a_t(p) it first has to be divided by -p.

For example, in the case of penalized maximum likelihood, the pairs (d'(p)/w(p) + p - 0.5 , 0) and (d'(p)/w(p) + p , 0.5/p) are equivalent, in the sense that if the corresponding pseudo-data representations are substituted in the ordinary scores both return the same expression.

So, in order to construct a pseudo-data representation that corresponds to a user-specified link function and has the property 0 \le a_r(p) \le a_t(p) for every p \in (0, 1), one merely has to pursue a simple algebraic calculation on the initial pair (a_r(p), a_t(p)) using only the two aforementioned rules until an appropriate pair is resulted. There is always a pair!

Once the pair has been found the following steps should be followed.

1. For a user-specified link function the user has to write a modification function with name "br.custom.family" or "pml.custom.family" for pl=FALSE or pl=TRUE, respectively. The function should take as argument the probabilities p and return a list of two vectors with same length as p and with names c("ar", "at"). The result corresponds to the pair (a_r(p), a_t(p)).

2. Check if the custom-made modifications function is appropriate. This can be done via the function checkModifications which has arguments fun (the function to be tested) and Length with default value Length=100. Length is to be used when the user-specified link function takes as argument a vector of values (e.g. the logexp link in ?family). Then the value of Length should be the length of that vector.

3. Put the function in the search patch so that modifications can find it.

4. brglm can now be used with the custom family as glm would be used.

Note

The user could also deviate from modified-scores and maximum penalized likelihood and experiment with implemented (or not) links, e.g. probit, constructing his own pseudo-data representations of the aforementioned general form. This could be done by changing the link name, e.g. by

probitt <- make.link(probit) ; probitt$name <- "probitt"

and then setting a custom br.custom.family that does not necessarily depend on the probit link. Then, brglm could be used with pl=FALSE.

A further generalization would be to completely remove the hat value h in the generic expression of the pseudo-data representation and have general additive modifications that depend on p. To do this divide both ar and at by pmax(get("hats",parent.frame()),.Machine\$double.eps) within the custom modification function (see also Examples).

Author(s)

Ioannis Kosmidis, [e-mail]

References

Kosmidis I. and Firth D. (2021). Jeffreys-prior penalty, finiteness and shrinkage in binomial-response generalized linear models. Biometrika, 108, 71–82.

Kosmidis, I. (2007). Bias reduction in exponential family nonlinear models. PhD Thesis, Department of Statistics, University of Warwick.

See Also

brglm, brglm.fit

Examples

## Begin Example 1
## logistic exposure model, following the Example in ?family. See,
## Shaffer, T.  2004. Auk 121(2): 526-540.
# Definition of the link function
logexp <- function(days = 1) {
  linkfun <- function(mu) qlogis(mu^(1/days))
  linkinv <- function(eta) plogis(eta)^days
  mu.eta <- function(eta) days * plogis(eta)^(days-1) *
        binomial()$mu.eta(eta)
  valideta <- function(eta) TRUE
  link <- paste("logexp(", days, ")", sep="")
  structure(list(linkfun = linkfun, linkinv = linkinv,
    mu.eta = mu.eta, valideta = valideta, name = link),
    class = "link-glm")
}
# Here d(p) = days * p * ( 1 - p^(1/days) )
#      d'(p) = (days - (days+1) * p^(1/days)) * d(p)
#      w(p) = days^2 * p * (1-p^(1/days))^2 / (1-p)
# Initial modifications, as given from the general expressions above:
br.custom.family <- function(p) {
  etas <- binomial(logexp(.days))$linkfun(p)
  # the link function argument `.days' will be detected by lexical
  # scoping. So, make sure that the link-function inputted arguments
  # have unusual names, like `.days' and that
  # the link function enters `brglm' as
  # `family=binomial(logexp(.days))'.
  list(ar = 0.5*(1-p)-0.5*(1-p)*exp(etas)/.days,
       at = 0*p/p) # so that to fix the length of at
}
.days <-3
# `.days' could be a vector as well but then it should have the same
# length as the number of observations (`length(.days)' should be
# equal to `length(p)'). In this case, `checkModifications' should
# have argument `Length=length(.days)'.
#
# Check:
## Not run: checkModifications(br.custom.family)
# OOOPS error message... the condition is not satisfied
#
# After some trivial algebra using the two allowed operations, we
# get new modifications:
br.custom.family <- function(p) {
  etas <- binomial(logexp(.days))$linkfun(p)
  list(ar=0.5*p/p, # so that to fix the length of ar
       at=0.5+exp(etas)*(1-p)/(2*p*.days))
}
# Check:
checkModifications(br.custom.family)
# It is OK.
# Now,
modifications(binomial(logexp(.days)))
# works.
# Notice that for `.days <- 1', `logexp(.days)' is the `logit' link
# model and `a_r=0.5', `a_t=1'.
# In action:
library(MASS)
example(birthwt)
m.glm <- glm(formula = low ~ ., family = binomial, data = bwt)
.days <- bwt$age
m.glm.logexp <- update(m.glm,family=binomial(logexp(.days)))
m.brglm.logexp <- brglm(formula = low ~ ., family =
binomial(logexp(.days)), data = bwt)
# The fit for the `logexp' link via maximum likelihood
m.glm.logexp
# and the fit for the `logexp' link via modified scores
m.brglm.logexp
## End Example
## Begin Example 2
## Another possible use of brglm.fit:
## Deviating from bias reducing modified-scores:
## Add 1/2 to the response of a probit model.
y <- c(1,2,3,4)
totals <- c(5,5,5,5)
x1 <- c(1,0,1,0)
x2 <- c(1,1,0,0)
my.probit <- make.link("probit")
my.probit$name <- "my.probit"
br.custom.family <- function(p) {
   h <- pmax(get("hats",parent.frame()),.Machine$double.eps)
   list(ar=0.5/h,at=1/h)
}
m1 <- brglm(y/totals~x1+x2,weights=totals,family=binomial(my.probit))
m2 <- glm((y+0.5)/(totals+1)~x1+x2,weights=totals+1,family=binomial(probit))
# m1 and m2 should be the same.
# End example
# Begin example 3: Maximum penalized likelihood for logistic regression, 
# with the penalty being a powerof the Jeffreys prior (`.const` below)
# Setup a custom logit link
mylogit <- make.link("logit")
mylogit$name <- "mylogit"
## Set-up the custom family
br.custom.family <- function(p) {
     list(ar = .const * p/p, at = 2 * .const * p/p)
}
data("lizards")
## The reduced-bias fit is
.const <- 1/2
brglm(cbind(grahami, opalinus) ~ height + diameter +
          light + time, family = binomial(mylogit), data=lizards)
## which is the same as what brglm does by default for logistic regression
brglm(cbind(grahami, opalinus) ~ height + diameter +
          light + time, family = binomial(logit), data=lizards)
## Stronger penalization (e.g. 5/2) can be achieved by 
.const <- 5/2
brglm(cbind(grahami, opalinus) ~ height + diameter +
          light + time, family = binomial(mylogit), data=lizards)
# End example