source("http://sites.williams.edu/bklingen/files/2012/09/plotModels.r")
dose <- c(0,1,4,12,24)
M1 <- list(family=binomial())
M2 <- list(model=pow(dose,0.5), family=binomial())
M3 <- list(model=pow(dose,0), family=binomial())
M4 <- list(model=pow(dose,-0.5), family=binomial())
M5 <- list(model=pow(dose,-1), family=binomial())
M6 <- list(family=binomial(link=log))
M7 <- list(model=expo(dose,2,scale=max(dose)), family=binomial(link="identity"))
M8 <- list(model=pow(dose,c(1,2),dmax=14), family=binomial())
M9 <- list(model=pow(dose,c(0,-1),dmax=8), family=binomial())
M10 <- list(model=pow(dose,c(0,1),dmax=8), family=binomial())
models <- list(M1,M2,M3,M4,M5,M6,M7,M8,M9,M10)
plotModels(dose, models, low=0.3, high=0.65) #plots candidate dose-resp. models for given placebo effect (low) and maximal efficacy (high)

#Inference with candidate models (adj. P-values, MED,...):
source("http://sites.williams.edu/bklingen/files/2012/09/perm_minP_GLM.r")
y <- c(38,52,67,59,58) # successes
n <- c(100,102,98,99,94) # sample sizes
resp <- cbind(y,n-y)
dr <- permT(dose,resp,models,perms=500,trace=50,clinRel=0.15,alpha=0.025)
summary(dr)
hist(dr, xlim=c(0,0.1), breaks=200) # plots the histogram of the minimum P-value
plot(dr) # plots the fitted model with the smallest adj. P-value, see Figure 4
plot(dr, which.model=c("M5","M9","M10","M4"), se.fit=FALSE) # similar to Figure 4