# Read grouped beetle data from the web:
data="http://www.stat.ufl.edu/~aa/glm/data/Beetles2.dat"
beetle = read.table(data, header = T)

# Fit logistic regression with logdose as covariate:
fit=glm(cbind(dead,n-dead)~logdose,family=binomial,data=beetle)
summary(fit)

# Compute 95% Wald confidence intervals of the betas from scratch:
low=fit$coef-qnorm(0.975)*sqrt(diag(summary(fit)$cov.scaled))
up=fit$coef+qnorm(0.975)*sqrt(diag(summary(fit)$cov.scaled))
cbind(low,up)

# Compute 95% Wald confidence intervals of the betas:
ci.wald=confint.default(fit)
ci.wald

# Compute confidence intervals of the betas based 
# on the profile likelihood:
ci.lik=confint(fit)
ci.lik      

# Compute estimated odds ratio for 10% increase in dose  
# [i.e. an increase in log-dose equal to log10(1.1)=0.0414]
exp(log10(1.1)*fit$coef[2])

# Compute confidence interval of odds ratio for 10% 
# increase in dose (based on profile likelihood)
exp(log10(1.1)*ci.lik[2,])

# Compute linear predictor with SE's
fit.lin=predict(fit,se.fit=T)
fit.lin$fit
fit.lin$se.fit

# Compute confidence intervals for probabilities by
# transforming CI's for the linear predictor
logistic.func = function (x) {exp(x)/(1 + exp(x))}
low=logistic.func(fit.lin$fit-1.96*fit.lin$se.fit)
up=logistic.func(fit.lin$fit+1.96*fit.lin$se.fit)
cbind(low,up)

# Compute estimated probabilities with SE's (from the delta method)
fit.prob=predict(fit,type='response',se.fit=T)
fit.prob$fit
fit.prob$se.fit

# Compute confidence intervals for probabilities 
# based on SE's from the delta method (and normal approximation) 
cbind(fit.prob$fit-1.96*fit.prob$se.fit,fit.prob$fit+1.96*fit.prob$se.fit)
#Better to transform CI's for the linear predictor!