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

[Ben Bolker]

Mark Herzog <mherzog <at> prbo.org> writes:

> 
> 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.
> 

  Hmmm.  When I do this with mle() I get the same answers
as you do with the logexposure family.  I'm stumped too.
Following on from your example ...

library(stats4)

mlefun <- function(data) {
  with(data,
       function(p1,p2,p3) {
         p <- c(p1,p2,p3) ## kluge to put list back into vector
         mu <- model.matrix(~parastat+patsize)%*%p
         prob <- plogis(mu)^expos
         -sum(dbinom(survive,size=trials,prob=prob,log=TRUE))
       })
}

m0 <- mlefun(nestdata)
m1 <- mle(minuslogl=m0,start=list(p1=2.7,p2=1,p3=-1))

 m1

Call:
mle(minuslogl = m0, start = list(p1 = 2.7, p2 = 1, p3 = -1))

Coefficients:
       p1        p2        p3
 2.746707  1.034936 -1.084389        



c1 <- confint(m1)
Profiling...
         2.5 %     97.5 %
p1  2.12739190  3.4903797
p2  0.01279248  2.0646978
p3 -2.12163290 -0.1156173

s1 <-  summary(m1)

Maximum likelihood estimation

Call:
mle(minuslogl = m0, start = list(p1 = 2.7, p2 = 1, p3 = -1))

Coefficients:
    Estimate Std. Error
p1  2.746707  0.3442989
p2  1.034936  0.5201348
p3 -1.084389  0.5093672

coef(s1)[,"Estimate"]+1.96*outer(coef(s1)[,"Std. Error"],c(-1,1))
          [,1]        [,2]
p1  2.07188136  3.42153311
p2  0.01547139  2.05439969
p3 -2.08274893 -0.08602966


  Standard errors are slightly different between glm()
and mle() [whether summary() or confint()].


  Ben Bolker