rm(list=ls())

## useful functions
expit <- function(x) 1/(1+exp(-x))
logit <- function(x) log(x/(1-x))

## 
## EXAMPLE 1: failure of wald confidence intervals
## 

num.success <- 29
num.total <- 31
  
## using likelihood
curve(dbinom(num.success, num.total, prob=x, log=TRUE),
      from=0.8, to=1, las=1)

## find max
f.binom <- function(x)
  dbinom(num.success, num.total, prob=x, log=TRUE)

mle <- optimize(f.binom, interval=c(0, 1), maximum=TRUE)$maximum

f.binom.shifted <- function(x)
  f.binom(x)-f.binom(mle)+1.92
f.lower.ci <- uniroot(f.binom.shifted, lower=0, upper=mle)$root
f.upper.ci <- uniroot(f.binom.shifted, lower=mle, upper=1)$root
ci.ll <- c(f.lower.ci, f.upper.ci)
ci.ll

## using GLM

## first, make a data.frame
ww <- data.frame(choice=c(rep(1,num.success),
                          rep(0,num.total-num.success)))


library(lme4)
out <- glm(choice ~ 1, family=binomial(link='logit'), data=ww)
coefs <- coef(summary(out))

## calculate intercept on our original scale
int <- expit(coefs['(Intercept)','Estimate'])
int.se <- coefs['(Intercept)','Std. Error']

## compare to mle
int
mle

## using the formula for "normal approximation" of the 95% confidence
## interval
int + c(qnorm(0.025),qnorm(0.975))*sqrt((1/nrow(ww))*int*(1-int))
## not good - goes outside the range of possible values!

## compare to our likelihood-based interval
ci.ll

## using R's confint function gives us the same as we calclulated with
## our likelihood based interval
expit(confint(out))

## 
## EXAMPLE 2: plotting a glm manually
## 

## simulate some data
set.seed(1)
## uniformly distributed predictor (could be body size, for example)
body.size <- runif(100, min=0, max=1)
range(body.size)

## poisson distributed response (could be fecundity, for example)
##
## let's assume an effect size of 2
##
## we are assuming a linear relationship in 'linear space', which
## means an exponential relationship in 'count space' - to go from
## 'count space', we use the log-link.  You might think - why not
## simulate our data in 'count' space - but what does it mean to have
## an effect size or 'slope' of 2 in count space?
intercept <- 0.5
effect.size <- 2
## calculate 'predictor'
eta <- intercept + effect.size*body.size

fecundity <- rpois(length(body.size), lambda=exp(eta))
dd <- data.frame(body.size=body.size, fecundity=fecundity)
head(dd)

## --------------------------------------------------
## plot the raw data
plot(x=dd$body.size, y=dd$fecundity,
     xlab='Body size', ylab='Fecundity',
     col='red', pch=16, las=1)
## --------------------------------------------------

## --------------------------------------------------
## fit glm
out <- glm(fecundity ~ body.size, data=dd,
           family=poisson(link='log'))
summary(out)
## model estimated our effect size exactly (note, the model is
## estimating our coefficients in 'linear space', which is how we
## specified them

## extract coefficients
coefs <- coef(summary(out))
b0 <- coefs['(Intercept)','Estimate']
b1 <- coefs['body.size','Estimate']

## add the best fit line
curve(exp(b0 + b1*x),
      from=min(dd$body.size), to=max(dd$body.size),
      col='blue', add=TRUE)
## --------------------------------------------------