# 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). ```r 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 = 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 $\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. c("$0 \\le y*\n$", "$ \\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, ioannis.kosmidis@warwick.ac.uk ## 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 ```r ## 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 ```