>  # 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)
[1]  5 20
>   t(ally[,(1:3)])
     [,1] [,2] [,3] [,4] [,5]
[1,]    6    9   10   10   10
[2,]    7    8   10   10   10
[3,]    4    7    9    9   10
>    ests <- apply(ally,2,estims)
[1] 1000    7   10   10   10   10
[1] 3000    7   10   10   10   10
[1] 1000    3   10   10   10   10
[1] 3000    3   10   10   10   10
[1] 1000    3   10   10   10   10
[1] 3000    3   10   10   10   10
>  # do we get any errors thrown?
>    c(nrow(ests),ncol(ests))
[1] 17 20
>  # estimators in pairs: glm, nlm, nls
>    t(ests)[,c(1,2,5,6,9,10)]
            [,1]       [,2]        [,3]      [,4]        [,5]      [,6]
 [1,]  0.3539471  2.1358209  0.35394722  2.135820  0.37055080 1.9999717
 [2,]  0.6243082  1.4365662  0.62431127  1.436563  0.63981431 1.2082791
 [3,] -0.3338466  1.0997360 -0.33384666  1.099735 -0.29026580 1.0142629
 [4,]  0.1021951  0.8361676  0.10219573  0.836167  0.12680258 0.7731008
 [5,]  0.3539471  2.1358209  0.35394722  2.135820  0.37055080 1.9999717
 [6,] -0.4333046  2.7830642 -0.43330489  2.783065 -0.42181925 2.7032828
 [7,] -0.1266738  1.9321766 -0.12667369  1.932174 -0.09506039 1.7513829
 [8,]  0.3539471  2.1358209  0.35394722  2.135820  0.37055080 1.9999717
 [9,]  0.2427025  1.6783887  0.24270210  1.678389  0.27056555 1.4720464
[10,]  0.8472979 24.1275786  0.84729862 14.158359         NaN       NaN
[11,]  0.6243082  1.4365662  0.62431127  1.436563  0.56203273 1.1502644
[12,] -0.8472979 25.4248407 -0.84729774 19.890232         NaN       NaN
[13,]  0.1088975  1.4785771  0.10889695  1.478577  0.13209424 1.2441324
[14,]  0.3539471  2.1358209  0.35394722  2.135820  0.37055080 1.9999717
[15,] -0.6151895  1.9485861 -0.61519147  1.948587 -0.56590616 1.7452546
[16,]  0.3539471  2.1358209  0.35394722  2.135820  0.37055080 1.9999717
[17,] -0.5085853  2.2145797 -0.50858480  2.214579 -0.47629427 2.0609164
[18,] -0.8472979 25.4248407 -0.84729774 19.890232         NaN       NaN
[19,] -2.2118024  4.4381797 -2.21180204  4.438177 -2.20474606 4.4169224
[20,] -0.0368765  2.4477028 -0.03687927  2.447704 -0.02306711 2.3421351
>  # flags  1       2       3        4     5
>  #  logreg$conv,flaglg,mle$code,flagml,flagnl   
>    t(ests)[,(13:17)]
      [,1] [,2] [,3] [,4] [,5]
 [1,]    1    0    1    0    0
 [2,]    1    0    1    0    0
 [3,]    1    0    1    0    0
 [4,]    1    0    1    0    0
 [5,]    1    0    1    0    0
 [6,]    1    0    1    0    0
 [7,]    1    0    1    0    0
 [8,]    1    0    1    0    0
 [9,]    1    0    1    0    0
[10,]    1    1    2    0    1
[11,]    1    0    1    0    0
[12,]    1    1    1    0    1
[13,]    1    0    1    0    0
[14,]    1    0    1    0    0
[15,]    1    0    1    0    0
[16,]    1    0    1    0    0
[17,]    1    0    1    0    0
[18,]    1    1    1    0    1
[19,]    1    0    1    0    0
[20,]    1    0    1    0    0
>  # done
>   rm(list=ls())
>   q()