# [R] second try; writing user-defined GLM link function

Mark Herzog mherzog at prbo.org
Wed Apr 19 08:41:59 CEST 2006

```Excellent!

So to provide the final summary for everyone's sake, based on Dr.
Ripley's revisions, and now a revised logexposure "family" function to
use with the improved link function (until 2.4.0 makes it simpler) :

# Logistical Exposure Link Function
# See Shaffer, T.  2004. A unifying approach to analyzing nest success.
# Auk 121(2): 526-540.

logexp <- function(days = 1)
{
linkfun <- function(mu) qlogis(mu^(1/days))
linkinv <- function(eta) plogis(eta)^days
mu.eta <- function(eta)
days*.Call("logit_mu_eta", eta,
PACKAGE = "stats")*plogis(eta)^(days-1)
valideta <- function(eta) TRUE
link <- paste("logexp(", days, ")", sep="")
mu.eta = mu.eta, valideta = valideta, name = link),
}

# Modified binomial family function (that allows logexp link function)

variance <- function(mu) mu * (1 - mu)
validmu <- function(mu) all(mu > 0) && all(mu < 1)
dev.resids <- function(y, mu, wt) .Call("binomial_dev_resids",
y, mu, wt, PACKAGE = "stats")
aic <- function(y, n, mu, wt, dev) {
m <- if (any(n > 1))
n
else wt
-2 * sum(ifelse(m > 0, (wt/m), 0) * dbinom(round(m *
y), round(m), mu, log = TRUE))
}
initialize <- expression({
if (NCOL(y) == 1) {
if (is.factor(y)) y <- y != levels(y)[1]
n <- rep.int(1, nobs)
if (any(y < 0 | y > 1)) stop("y values must be 0 <= y <= 1")
mustart <- (weights * y + 0.5)/(weights + 1)
m <- weights * y
if (any(abs(m - round(m)) > 0.001))
warning("non-integer successes in a binomial glm!")
} else if (NCOL(y) == 2) {
if (any(abs(y - round(y)) > 0.001))
warning("non-integer counts in a binomial glm!")
n <- y[, 1] + y[, 2]
y <- ifelse(n == 0, 0, y[, 1]/n)
weights <- weights * n
mustart <- (n * y + 0.5)/(n + 1)
} else stop("for the binomial family,",
" y must be a vector of 0 and 1's\n",
"or a 2 column matrix where col 1 is",
" no. successes and col 2 is no. failures")
})
dev.resids=dev.resids, aic=aic,
mu.eta=logexp(ExposureDays)\$mu.eta, initialize=initialize,
validmu=validmu, valideta=logexp\$valideta), class = "family")
}

#Example
chat.glm.logexp<-glm(survive/trials~parastat+nest_ht*patsize,family=logexposure(ExposureDays=nestdata\$expos),data=nestdata)
# if you have MASS installed
library(MASS)
chat.step<-stepAIC(chat.glm.logexp,scope=list(upper=~parastat+nest_ht*patsize,lower=~1))
chat.step\$anova
summary(chat.step)

Mark Herzog, Ph.D.
Program Leader, San Francisco Bay Research
Wetland Division, PRBO Conservation Science
4990 Shoreline Highway 1
Stinson Beach, CA 94970
(415) 893-7677 x308
mherzog at prbo.org

Prof Brian Ripley wrote:
> A couple more comments:
>
> The null deviance is wrong here, as the code assumes that the link
> function maps constant vectors to constant vectors, which it does not
> here.  You can circumvent that by setting an offset.
>
> Even setting dispersion = 1 I get slightly different se's.
>
> Here's a more robust way to specify the link:
>
> logexp <- function(days = 1)
> {
>     linkfun <- function(mu) qlogis(mu^(1/days))
>     linkinv <- function(eta) plogis(eta)^days
>     mu.eta <- function(eta)
>       days*.Call("logit_mu_eta", eta, PACKAGE =
> "stats")*plogis(eta)^(days-1)
>     valideta <- function(eta) TRUE
>     link <- paste("logexp(", days, ")", sep="")
>                    mu.eta = mu.eta, valideta = valideta, name = link),
>               class = "link-glm")
> }
>
> and in 2.4.0 you will able to use binomial(logexp(nestdata\$exposure)).
>

```