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

Ben Bolker bolker at ufl.edu
Mon Apr 17 19:13:23 CEST 2006

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

```