[R] se.fit in predict.glm

[(Ted Harding)]

On 28-Apr-04 Peter Dalgaard wrote:
> Hmm, that doesn't actually hold up mathematically... You cannot just
> take upper/lower limits of both parameters and combine them.

Yes, excuse me. I had misunderstood what 'confint' does -- I misled
myself into thinking that it gave be the 'a' and 'b' for the
lower and the upper curve!

> I don't seem to get anything that drastic. Things look somewhat better
> if you use the link-scale estimates, but the response-scale curves are
> not hopeless. You do have to notice that these are pointwise CIs so
> having multiple curves straying outside is not necessarily a problem.
> 
> Just as a sanity check, did your plots look anything like the below:

Yes, that's useful -- and it also brought to the surface the other
error I'd made whereby a slip of eyesight caused me to transcribe
a correlation of -0.91 as 0.91, so I was simulating from the
wrong distribution ...

> y <- c(0,0,0,0,0,0,0,0,0,0,2,0,0,1,0,2,4,5,8,8,9,4,10)
> x <- 1:23
> d <- data.frame(x,y,n=12)
> m1 <- glm(cbind(y,n-y) ~ x,data=d,family=binomial(probit))
> confint(m1) # +-2SE approximation
> library(MASS)
> confint(m1) # profiling
> x1 <- with(predict(m1,se.fit=TRUE,type="response"),
>           fit+outer(se.fit,c(l=-2,u=2)))
> x2 <- with(predict(m1,se.fit=TRUE,type="link"),
>           fit+outer(se.fit,c(l=-2,u=2)))
> x2 <- pnorm(x2)
> 
> ab <- mvrnorm(20,coef(m1),vcov(m1))
> matplot(x2,type="l",col="black",lwd=3,lty=1)
> matlines(x1,type="l",col="red",lwd=3,lty=1)
> matlines(pnorm(t(ab%*%rbind(1,x))))

Yes, fairly similar when I do it right. Thanks for the above
code. I'd also overlooked 'vcov' etc and was copying SEs and
correlations from the output of 'summary.glm'.

The main lesson: one shouldn't stay up too late at these things.
Cheers,
Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 167 1972
Date: 28-Apr-04                                       Time: 23:54:20
------------------------------ XFMail ------------------------------