[R] Adding the complementary log-link to binomial() and make.link()

[Ben Bolker]

Ken Knoblauch <ken.knoblauch <at> inserm.fr> writes:

> 
> Roland Deutsch <roland.deutsch <at> tuwien.ac.at> writes:
> > in my research I frequently work with binomial 
> response models, which 
> > are of course part of the generalized linear 
> models. While I do use 
> > common link functions such as the logit, probit 
> and cloglog, I often 
> > have the need of invoking the lesser-known 
> Complementary Log link 
> > (Walter W. Piegorsch, 1992, "Complementary Log 
> Regression for 
> > Generalized Linear Models ", The American 
> Statistician , Vol. 46, No. 2, 
> > pp. 94-99 ) which is based on the exponential 
> distribution.
> > 
> > Before the release of R 3.0, I simply could 
> do so by adding the 
> > following lines to the long switch command 
> in make.link(...):
> > 
> > clog = {
> >          linkfun <- function(mu) qexp(mu)
> >          linkinv <- function(eta) pmax(.Machine$double.eps,pexp(eta))
> >          mu.eta <- function(eta) pmax(dexp(eta), .Machine$double.eps)
> >          valideta <- function(eta) all(eta > 0)
> >      },
> > 
> > and then add "clog" to the okLinks vector in 
> binomial(). However, I 
> 
> I wouldn't mess with make.link.  Why don't you just
> define your own custom link, for example following the
> example under help(family)?  This is what I did in
> the psyphy package and I just checked and
> they all still work under the current version of R.
> 
> > Thanks a lot,
> > Roland Deutsch
> 

This seems to work:

linkinfo <- list(link="clog",
                 linkfun=function(mu) qexp(mu),
                 linkinv=function(eta) pmax(.Machine$double.eps,pexp(eta)),
                 mu.eta=function(eta) pmax(dexp(eta), .Machine$double.eps),
                 valideta=function(eta) all(eta > 0))
binomClog <- binomial()
for (i in names(linkinfo))
    binomClog[[i]] <- linkinfo[[i]]

set.seed(101)
d <- data.frame(x=runif(1000))
d$y <- rbinom(1000,prob=binomClog$linkinv(1+2*d$x),size=1)
library(lattice)
xyplot(y~x,data=d,type=c("p","smooth"))
g1 <- glm(y~x,data=d,family=binomClog)
library(MASS)
confint(g1)