# 
# Data:
#	   CC   CT   TT
# Cases     6    8   75
# Controls 10   66  163
#
#  Likelihood first
# 
x <- c(0,1,2)
y <- c(6,8,75)
z <- c(10,66,163)
logitmod <- glm(cbind(y,z)~x,family="binomial")
thetahat <- logitmod$coeff[2]		# Log odds ratio
V <- vcov(logitmod)[2,2]		# se^2
exp(thetahat) 
exp(thetahat-1.96*sqrt(V))
exp(thetahat+1.96*sqrt(V))
summary(logitmod)
# 
# Bayesian
#
source("http://faculty.washington.edu/jonno/BFDP.R")
pi1 <- c(1/100,1/1000,1/10000,1/100000) # Prior on alternative
Upper975 <- 1.5
W <- (log(Upper975)/1.96)^2 # 97.5 point of prior is log(1.5) so that we 
                            # believe with prior prob 0.95 that the odds
                            # corresponding to one T allele more,
			    # lies in (2/3,1.5)
BFcall <- BFDPfunV(thetahat,V,W,pi1)
r <- W/(V+W)
exp(r*thetahat)
exp(r*thetahat-1.96*sqrt(r*V))
exp(r*thetahat+1.96*sqrt(r*V))
cat("log odds, standard error and ratio: ",thetahat,sqrt(V),
     thetahat/sqrt(V),"\n")
# Posterior distribution is asymptotically lognormal with mean r x thetahat
# and variance r x thetahat
RRseq <- seq(.1,4,.01)
postseq <- dlnorm(RRseq,mean=r*thetahat,sd=sqrt(r*V))
plot(postseq~RRseq,type="n",xlab="Odds Ratio",ylab="Posterior")
priorseq <- dlnorm(RRseq,mean=0,sd=sqrt(W))
lines(priorseq~RRseq,lty=2)
lines(postseq~RRseq)
legend("topright",legend=c("Prior Distribution","Posterior Distribution"),
lty=2:1,bty="n")
cat("Bayes Factor of H0 over H1: ",BFcall$BF,"\n")
cat("Posterior probs of the null, given priors on the null of: ",1-pi1," are ",BFcall$BFDP,"\n")