[R] 95% confidence intercal with glm

Sam Sam_Smith at me.com
Wed Sep 29 16:04:56 CEST 2010


Hi Ben and list,

Sorry to be a pain! I have followed your code, and modified it -

> pp <- predict(model1,se.fit=TRUE, type = "response")
>> etaframe <-
> + with(pp,cbind(fit,lower=fit-1.96*se.fit,upper=fit+1.96*se.fit))
>> pframe <- plogis(etaframe) 
>> pframe

My response variable is 0 or 1, the probabilities are all high above my cut-off point of 0.445, i think this may have something to do with 

> you may need to multiply by the binomial N as
> appropriate.

However how do i multiply if my response is 0 or 1?

Furthermore, will the approach above work for a ordinal response?

Thanks

Sam

On 29 Sep 2010, at 13:44, Ben Bolker wrote:

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1


## from ?glm
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
d.AD <- data.frame(treatment, outcome, counts)
glm.D93 <- glm(counts ~ outcome + treatment, family=poisson,
             data=d.AD)

## predict on 'link' or 'linear predictor' scale, with SEs
pp <- predict(glm.D93,se.fit=TRUE)
etaframe <-
with(pp,cbind(fit,lower=fit-1.96*se.fit,upper=fit+1.96*se.fit))
pframe <- exp(etaframe)  ## inverse-link
## picture
pframe <- as.data.frame(cbind(obs=d.AD$counts,pframe))
library(plotrix)
plot(pframe$obs,ylim=c(5,30))
with(pframe,plotCI(1:9,fit,li=lower,ui=upper,col=2,add=TRUE))

If you're using a binomial model you need 'plogis' rather than 'exp'
as your
inverse link, and you may need to multiply by the binomial N as
appropriate.

On 10-09-29 06:07 AM, Sam wrote:
> I am looking to do the same but am a bit confused
> 
>> and apply the inverse link function for your model.
> 
> i understand up to this point and i understand what this means,
> however i am unsure why it needs to be done and how you do it - i.e
> i use family="binomial" is this wrong if i use this method to
> calculate 95% CI?
> 
> Thanks
> 
> Sam On 28 Sep 2010, at 14:50, Ben Bolker wrote:
> 
> zozio32 <remy.pascal <at> gmail.com> writes:
> 
>> 
>> 
>> Hi
>> 
>> I had to use a glm instead of my basic lm on some data due to
>> unconstant variance.
>> 
>> now, when I plot the model over the data, how can I easily get
>> the 95% confidence interval that sormally coming from:
>> 
>>> yv <- predict(modelVar,list(aveLength=xv),int="c")
>>> matlines(xv,yv,lty=c(1,2,2))
>> 
>> There is no "interval" argument to pass to the predict function
>> when using a glm, so I was wondering if I had to use an other
>> function
>> 
> 
> You need to use predict with se=TRUE; construct the confidence
> intervals by computing predicted values +- 1.96 times the standard
> error returned; and apply the inverse link function for your
> model.
> 
> If heteroscedasticity is your main problem, and not a specific
> (known) non-normal distribution, you might consider using the gls
> function from the nlme package with an appropriate 'weights'
> argument.
> 
> ______________________________________________ R-help at r-project.org
> mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do
> read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iEYEARECAAYFAkyjNKQACgkQc5UpGjwzenNW9gCgjfXin/9dI2y0rnk9wZWu97P8
TVgAn3WsG2ATva5WpLZXS1PXUTMRTTMi
=gvAc
-----END PGP SIGNATURE-----

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



More information about the R-help mailing list