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)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
