Additive Modifications to the Binomial Responses and Totals for Use within ‘brglm.fit’
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).
modifications(family, pl = FALSE)
family |
a family object of the form |
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
|
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).
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 μ be the mean of the binomial response y (i.e. μ = m p, 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 μ with respect to the linear predictor η 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 ≤ y*
≤ m*. Since h \in [0, 1], the requirement translates to
0 ≤ a_r(p) ≤ 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 ≤ a_r(p) ≤ 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.
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)).
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.
Put the function in the search patch so that
modifications can find it.
brglm can now be used with the custom family as
glm would be used.
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).
Ioannis Kosmidis, ioannis.kosmidis@warwick.ac.uk
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.
## 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 examplePlease choose more modern alternatives, such as Google Chrome or Mozilla Firefox.