[R] Probit Analysis and Interval Calculations for different LD50s

[Colleen Kenney]

I am encountering a problem with the calculation of Fieller and Delta
Method confidence intervals when performing probit analysis on
simulated data; my code is included below.  I am testing 5 dose
groups, with log doses (-0.2, -0.1, 0, 0.1, 0.2) and (1.8, 1.9, 2,
2.1, 2.2) so that the log(LD50) are 0 and 2, respectively.  However,
while I get the coverage as seen in the literature for the log doses
surrounding 0, I get very wide intervals when log(LD50)=2, with
everything else remaining constant.  Can anyone help please?


  nd=100
  N=10000
  m=5
  alpha=0.05

   x<-c(-0.2, -0.1, 0, 0.1, 0.2)

  logLD50<-0
  slope<-10

  for (i in 1:N){
    dose1[i]<-sum(rbinom(nd, 1, pnorm((x[1]-logLD50)*slope)))
    dose2[i]<-sum(rbinom(nd, 1, pnorm((x[2]-logLD50)*slope)))
    dose3[i]<-sum(rbinom(nd, 1, pnorm((x[3]-logLD50)*slope)))
    dose4[i]<-sum(rbinom(nd, 1, pnorm((x[4]-logLD50)*slope)))
    dose5[i]<-sum(rbinom(nd, 1, pnorm((x[5]-logLD50)*slope)))
}


  ld50<-function(b) -b[1]/b[2]


  for (i in 1:N){

  pw<-data.frame(x=x, n=rep(nd, m), y=c(dose1[i], dose2[i], dose3[i],
dose4[i], dose5[i]))
  pw$Ymat<-cbind(pw$y, nd-pw$y)
  pwp.1<-glm(Ymat~x, family=binomial(link=probit), data=pw)
  pwp<-summary(pwp.1)
  iter[i]<-pwp.1$iter
  ld[i]<-ld50(coef(pwp.1))
  a[i]<-coef(pwp.1)[1]
  b[i]<-coef(pwp.1)[2]

  nu11<-pwp$cov.unscaled[1,1]
  nu12<-pwp$cov.unscaled[1,2]
  nu22[i]<-pwp$cov.unscaled[2,2]
  mse[i]<- nu11/b^2+nu22*a^2/b^4-2*nu12*a/(b^3)
  s.ab<-sqrt(nu11/b^2+nu22*a^2/b^4-2*nu12*a/(b^3))
  z.alpha<-qnorm(1-alpha/2)
  g[i]<-z.alpha^2*nu22/b[i]^2
  fl.lower[i]<-ld[i]+g[i]/(1-g[i])*(ld[i]-nu12/nu22)-z.alpha/(b[i]*(1-g[i]))*sqrt(nu11-2*ld[i]*nu12+ld[i]^2*nu22-g[i]*(nu11-nu12^2/nu22))
 #Fieller interval
  fl.upper[i]<-ld[i]+g[i]/(1-g[i])*(ld[i]-nu12/nu22)+z.alpha/(b[i]*(1-g[i]))*sqrt(nu11-2*ld[i]*nu12+ld[i]^2*nu22-g[i]*(nu11-nu12^2/nu22))
  ci.lower[i]<-ld[i]-z.alpha*s.ab  #delta method interval
  ci.upper[i]<-ld[i]+z.alpha*s.ab

  }