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

[Mark Herzog]

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="")
     structure(list(linkfun = linkfun, linkinv = linkinv,
	mu.eta = mu.eta, valideta = valideta, name = link),
         class = "link-glm")
}

# Modified binomial family function (that allows logexp link function)
logexposure<-function (link="logexp",ExposureDays) {

     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")
     })
     structure(list(family="binomial", link=logexp(ExposureDays),
	linkfun=logexp(ExposureDays)$linkfun,
	linkinv=logexp(ExposureDays)$linkinv, variance=variance,
	dev.resids=dev.resids, aic=aic,
	mu.eta=logexp(ExposureDays)$mu.eta, initialize=initialize,
	validmu=validmu, valideta=logexp$valideta), class = "family")
}

#Example
nestdata<-read.table("http://www.branta.org/ShafferChatNestData/chat.txt")
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="")
>     structure(list(linkfun = linkfun, linkinv = linkinv,
>                    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)).
>