Ddose <- function(theta, data = y)
{
	opt <- maxdose(data = data)
	-2*(ldose(opt, data = data) - ldose(theta, data = data))
}

plotdose2 <- function(alpha=seq(-7,-3,.1),beta=seq(2.3,4,.05),data=y){
  #does perspective plot of Ddose likelihood surface
  n1 <- length(alpha)
  n2 <- length(beta)
  z <- matrix(0,n1,n2)
  for(i in 1:n1){for(j in 1:n2){
    z[i,j] <-  Ddose(c(alpha[i],beta[j]),data=data)}}
persp(alpha,beta,z,xlab="alpha",ylab="beta",zlab="Rel.Like.")
list(alpha=alpha,beta=beta,z=z)}

> zz <- plotdose2()
> contour(zz$alpha, zz$beta, zz$z, levels = c(-log(0.1),-log(0.05), -log(.01)))

Profile Likelihood

maxdose2_function(alpint=-5,beta=3,data=y){
  # Finds profile vlaues
xx_nlminb(alpint,ldose2,beta=beta,data=data)
xx$parameters}

ldose2_function(alpha,beta=3,data=y){
# Data has 3 columns dose, n, success. Gives negative loglikelihood

prob_exp(alpha+beta*data[,1])/(1+exp(alpha+beta*data[,1]))
like_sum(data[,3]*log(prob)+(data[,2]-data[,3])*log(1-prob))
-like}

Ddose2 <- function(beta, data = y)
{
	opt <- c(maxdose2(beta=beta,data = data),beta)
        opt1 <- maxdose(data=data)
	-2*(ldose(opt1, data = data) - ldose(opt, data = data))
}

plotdose2a <- function(beta=seq(2.3,4,.05),data=y){
  #does plot of Ddose2 likelihood curve
  n2 <- length(beta)
  z <- rep(0,n2)
  for(j in 1:n2){
    z[j] <-  Ddose2(beta[j],data=data)}
plot(beta,z,xlab="beta",ylab="D2",type='l')
  abline(h=-log(.1))
  abline(h=-log(.05))
  abline(h=-log(.01))
list(beta=beta,z=z)}

> plotdose2a()