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

Prof Brian Ripley ripley at stats.ox.ac.uk
Tue Apr 18 18:25:49 CEST 2006


Is this something that is fairly widely useful?  If so, it is something 
that we could add to the stats package.  (We ought in any case to make it 
easier to extend glm families, and this could be an example.  If so, I'd 
like a reference for its use.)

There is a good reason why logit is implmented as

> make.link("logit")
$linkfun
function (mu)
.Call("logit_link", mu, PACKAGE = "stats")
<environment: 01EB4150>

$linkinv
function (eta)
.Call("logit_linkinv", eta, PACKAGE = "stats")
<environment: 01EB4150>

$mu.eta
function (eta)
.Call("logit_mu_eta", eta, PACKAGE = "stats")
<environment: 01EB4150>

$valideta
function (eta)
TRUE
<environment: 01EB4150>

which is that these functions protect against underflow/overflow of exp, 
and why I suggested looking at make.link("probit").

You can write your linkfun and linkinv more accurately using plogis and 
qlogis.


On Mon, 17 Apr 2006, Mark Herzog wrote:

> I was a little hesitant to post to everyone until I figured out why
> there is a discrepancy in the intercept estimates when compared to the
> same model run in SAS vs. R.  Everything else comes out correctly,
> including the other coefficient estimates... so perhaps it is just the
> numerical method used.  I think glm in R is using IWLS, and SAS is using ML.
>
> If anyone has another idea feel free to let me know.  Here is my
> modified binomial family, now called logexposure.  Since there is no
> other choice for the link for this new logexposure "family", I have just
> embedded the make.link inside the logexposure family.
>
> Good luck and I'd appreciate any comments.
>
> logexposure<-function (link="logit",ExposureDays) {
>     logexp<-make.link("logit")
>     logexp$linkfun<-function(mu,expos=ExposureDays) {
>         log((mu^(1/expos))/(1-mu^(1/expos)))
>     }
>     logexp$linkinv<-function(eta,expos=ExposureDays) {
>         (exp(eta)/(1+exp(eta)))^expos
>     }
>     logexp$mu.eta<-function(eta,expos=ExposureDays) {
>         (expos*exp(eta*expos)*(1+exp(eta))^(-expos-1))
>     }

...

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595




More information about the R-help mailing list