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

[Mark Herzog]

It's probable that Roger's pdf that he provided is much more appropriate 
(or perhaps may have a good reference).  However, if interested, the 
specific example used for my code, is from the following paper (and is 
now widely being used in many avian nesting studies, but I leave to 
others to decide if useful in a much broader context):

Shaffer, T.  2004. A unifying approach to analyzing nest success. Auk 
121(2): 526-540.

In this paper Shaffer develops the link, but only cites 
McCullugh/Nelder, and Hosmer/Lemeshow (which I don't think discuss any 
other links for the binomial family besides those already available in R)
On his website ( 
http://www.npwrc.usgs.gov/resource/birds/nestsurv/nestsurv.htm ), 
Shaffer has some code examples in SAS (with the data embedded in his 
example code).  My previous email, with code and example, provides a 
direct link to a datafile with that same data.  I am fairly certain he 
would be willing to offer that data for use, if desired.

Mark
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:
> 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))
>>     }
> 
> ...
>