# tryCtest.r                  August, November 2009
 #
 # error handling in a simulation study
 #
  n <- 5              # dosage levels
  x <- c(0,1,2,3,4)   # doses
  m <- rep(10,5)      # ten patients at each dosage level
 #
 # ****************************************************************
 # function for computing estimators -- y's as input
  estims <- function(y) { # compute all estimators
     p0 <- sum(y)/sum(m)
     ll0 <- log(p0/(1-p0))
 # first logistic regression    
     ym <- cbind(y,m-y)    
     lstyxm <- list(ym,x,m)
     flaglg <- 0
     logreg <- tryCatch(glm(ym ~ x,family=binomial(link="logit") ),
         warning=function(e){})
     if(is.null(logreg)) { # if warning, flag and rerun
           flaglg <- 1
           print(c(1000,y))    
           logreg <- glm( ym ~ x,family=binomial(link="logit")) 
                         } 
     betalg <- c(coef(summary(logreg))[,(1:2)]) 
 # second minimize unlikelihood
     llikez <- function(z) { # creates unloglikelihood function     
      function(beta) {
     gam <- beta[1] + beta[2]*x
      -sum(z*gam-m*log(1+exp(gam)))} }   # - log-likelhood
     mle <- tryCatch(nlm(llikez(y),c(p0,0),print.level=0,hessian=TRUE),
         error = function(e){}  ) 
 #o
     betaml <- rep(NA,4) 
     flagml <- 1
     if(!is.null(mle)) {   # if failed, flag and set param ests
                   flagml <- 0
                   betaml <- c(mle$estimate,
                      sqrt(diag(chol2inv(chol(mle$hessian)))))
                          } else print(c(2000,y)) 
 # third weighted nonlinear regression
     vwls <- function(z,be0,be1) { # residual function to be minimized    
     p <- exp(be0+be1*x)/(1+exp(be0+be1*x))
     vwls <- (z-m*p)/sqrt( m*p*(1-p)) }   # weighted residuals
     nlsreg <- tryCatch(nls( ~vwls(y,be0,be1),
         start=list(be0=log(sum(y)/(sum(m)-sum(y))),be1=0),trace=FALSE),
         error = function(e){}  ) 
 #
     betanl <- rep(NaN,4)
     flagnl <- 1
     if(!is.null(nlsreg)) {   # if failed, flag and set param ests
                   flagnl <- 0
                   betanl <- c(coef(summary(nlsreg))[,(1:2)]) 
                          } else print(c(3000,y))
 # done with estimation
 # put two pairs of estimators together with warnings
    estims <- c(betalg,betaml,betanl,as.numeric(logreg$converged),
      flaglg,mle$code,flagml,flagnl )  
    }  # end of estimation
 # *********************************************************************
 #
 # set seed
   set.seed(5151917 + 4)
 # set parameters
  nrep <- 20
  beta0 <- -1          # true beta0
  beta0 <- 0          # true beta0
  beta1 <- 2          # true beta1
 # generate data
  p <- exp(beta0+beta1*x)/(1+exp(beta0+beta1*x))
  ally <- matrix(rbinom(n*nrep,m,p),n,nrep)
  dim(ally)
  t(ally[,(1:3)])
   ests <- apply(ally,2,estims)
 # do we get any errors thrown?
   c(nrow(ests),ncol(ests))
 # estimators in pairs: glm, nlm, nls
   t(ests)[,c(1,2,5,6,9,10)]
 # flags  1       2       3        4     5
 #  logreg$conv,flaglg,mle$code,flagml,flagnl   
   t(ests)[,(13:17)]
 # done
  rm(list=ls())
  q()