R/simu_results_summary_paraNLME.R

In unkyunglee/HDChangePoint: Estimating disease onset from change points of markers measured with error

#' Summary results of the parametric NLME, when the true data are generated by the logistic model and the parametric NLME assumes the correct logistic  model.
#'
#' @param nsim number of simulation runs.
#' @param n number of sample size.
#' @param res a nsim-length of list for simulation results.
#' @param theta1 a true scale parameter, which determines steepness of the logistic model.
#' @param theta3 a true parameter, which determines the maximum (asymptote) of logstic model.
#' @param beta  a (p-1)-length of true coefficient vector corresponding to subject specific covariates.
#' @param beta0 true intercept of the log-normal model for the inflection point.
#' @param time.length number of data points at which predictors are required for each individual longitudinal trajectory. This time point for graphs to be plotted.
#' @param file a character string of file name for a creating plot in .eps file extention.
#'
#' @return A list of the summarized simulation results including
#' \itemize{
#'        \item{res.table}{a data frame of the absolute biases, estimated standard deviations, average of the estimated standard errors and 95\% coverage probabilities for the fixed effects (\code{beta0}, \code{beta}).}
#'        \item{Rel.glob.logT.est}{a data frame of the relative average performance of inflection points for all subjects including the relative absolute biases, the relative empirical standard errors, the relative bootstrap standard deviations and 95\% boostrap confidence intervals.}
#' }
#' @importFrom graphics plot
#' @importFrom graphics legend
#' @importFrom graphics lines
#' @importFrom graphics mtext
#' @importFrom graphics par
#' @importFrom graphics segments
#' @importFrom grDevices postscript
#' @importFrom grDevices dev.off
#'
#' @export

logit.results<-function(nsim=100, n=80, res=results, theta1=6, theta3=1, time.length=20, beta=0.1, beta0=0.5, file="nlme_logist"){


  #############################################
  ### Fixed effects                       #####
  ###    #theta1                          #####
  ###    #theta2 =beta*W_cov+ranom.effect #####
  #############################################

  #######################
  # TRUE Fixed Effects ##
  #######################

  true.fixed.eff<-c(theta1,theta3, beta, beta0)



  ##############
  ## estimates #
  ##############

  ## main fixed effects
  estimated.theta1<-do.call("rbind",lapply(res,function(x) x$est.theta1))
  estimated.theta3<-do.call("rbind",lapply(res,function(x) x$est.theta3))

  ## baseline effects
  estimated.beta<-do.call("rbind",lapply(res,function(x) x$est.beta))
  estimated.beta0<-do.call("rbind",lapply(res,function(x) x$est.beta0))


  ## average of fixed effects
  mean.est.fixed.eff<-apply(cbind(estimated.theta1, estimated.theta3,  estimated.beta, estimated.beta0),2, mean)

  ## standard deviation of fixed effects
  sd.est.fixed.eff<-apply(cbind(estimated.theta1, estimated.theta3,  estimated.beta, estimated.beta0),2, sd)



  ####################################
  ### bias of fixed estimates
  ####################################
  abs.bias.fixed.eff<-abs(mean.est.fixed.eff-true.fixed.eff)


  #############################################
  ## standard errors of fixed effects
  #############################################
  estimated.str.theta1<-do.call("rbind",lapply(res,function(x) x$est.str.theta1))
  estimated.str.theta3<-do.call("rbind",lapply(res,function(x) x$est.str.theta3))

  estimated.str.beta0<-do.call("rbind",lapply(res,function(x) x$est.str.beta0))
  estimated.str.beta<-do.call("rbind",lapply(res,function(x) x$est.str.beta))


  est.avr.str.theta1<-mean(estimated.str.theta1)
  est.avr.str.theta3<-mean(estimated.str.theta3)
  est.avr.str.beta0<-mean(estimated.str.beta0)
  est.avr.str.beta<-mean(estimated.str.beta)


  ########################################
  ## coverage probabilities for theta1
  ########################################
  estimated.lowercp.theta1<-do.call("rbind",lapply(res,function(x) x$cp.lower.theta1))
  estimated.uppercp.theta1<-do.call("rbind",lapply(res,function(x) x$cp.upper.theta1))

  ind1theta1<- do.call("rbind",lapply(estimated.uppercp.theta1,function(x) 1*(x-(theta1)>0)))
  ind2theta1<- do.call("rbind",lapply(estimated.lowercp.theta1,function(x) 1*(x-(theta1)<0)))
  indtheta1<-ind1theta1*ind2theta1

  cptheta1<-mean(indtheta1)

  ########################################
  ## coverage probabilities for theta3 ###
  ########################################
  estimated.lowercp.theta3<-do.call("rbind",lapply(res,function(x) x$cp.lower.theta3))
  estimated.uppercp.theta3<-do.call("rbind",lapply(res,function(x) x$cp.upper.theta3))


  ind1theta3<- do.call("rbind",lapply(estimated.uppercp.theta3,function(x) 1*(x-(theta3)>0)))
  ind2theta3<- do.call("rbind",lapply(estimated.lowercp.theta3,function(x) 1*(x-(theta3)<0)))
  indtheta3<-ind1theta3*ind2theta3

  cptheta3<-mean(indtheta3)

  ########################################
  ## coverage probabilities for beta
  ########################################
  estimated.lowercp.beta0<-do.call("rbind",lapply(res,function(x) x$cp.lower.beta0))
  estimated.uppercp.beta0<-do.call("rbind",lapply(res,function(x) x$cp.upper.beta0))


  ind1b0<- do.call("rbind",lapply(estimated.uppercp.beta0,function(x) 1*(x-(beta0)>0)))
  ind2b0<- do.call("rbind",lapply(estimated.lowercp.beta0,function(x) 1*(x-(beta0)<0)))
  indb0<-ind1b0*ind2b0

  cpb0<-mean(indb0)


  ########################################
  ## coverage probabilities for beta
  ########################################
  estimated.lowercp.beta<-do.call("rbind",lapply(res,function(x) x$cp.lower.beta))
  estimated.uppercp.beta<-do.call("rbind",lapply(res,function(x) x$cp.upper.beta))


  ind1b<- do.call("rbind",lapply(estimated.uppercp.beta,function(x) 1*(x-(beta)>0)))
  ind2b<- do.call("rbind",lapply(estimated.lowercp.beta,function(x) 1*(x-(beta)<0)))
  indb<-ind1b*ind2b

  cpb<-mean(indb)

  ######################
  ### summarization ####
  ######################

  ## theta1
  estimates.theta1<-data.frame( abs.bias.fixed.eff["theta1"], sd.est.fixed.eff["theta1"],est.avr.str.theta1, cptheta1)
  colnames(estimates.theta1)<-c("Bias", "ESD", "ASE", "CP")

  ## theta3
  estimates.theta3<-data.frame( abs.bias.fixed.eff["theta3"], sd.est.fixed.eff["theta3"],est.avr.str.theta3, cptheta3)
  colnames(estimates.theta3)<-c("Bias", "ESD", "ASE", "CP")

  ## beta0
  estimates.beta0<-data.frame( abs.bias.fixed.eff["theta2.(Intercept)"], sd.est.fixed.eff["theta2.(Intercept)"],  est.avr.str.beta0,  cpb0)
  colnames(estimates.beta0)<-c("Bias", "ESD", "ASE", "CP")


  ## beta
  estimates.beta<-data.frame( abs.bias.fixed.eff["theta2.W.cov"], sd.est.fixed.eff["theta2.W.cov"],  est.avr.str.beta,  cpb)
  colnames(estimates.beta)<-c("Bias", "ESD", "ASE", "CP")

  #estimates.theta1,estimates.theta3,"theta1", "theta3",
  fixed.est<-rbind(estimates.beta0, estimates.beta)
  rownames(fixed.est)<-c("beta0", "beta")

  fixed.est<-round(fixed.est, 2)

  ######################################
  ## within individual scale parameter #
  ######################################
  estimated.sigma.eps<-do.call("rbind",lapply(res,function(x) x$sig.eps))
  mean.est.sigma.eps<-mean(estimated.sigma.eps)
  sd.est.sigma.eps<-sd(estimated.sigma.eps)


  #############################################
  ### Random effect where                   ###
  ###    #theta2 =beta*W_cov+ranom.effect   ###
  #############################################

  ## random effect for Each ID (columwise) based on 100 simulation
  estimated.random.eff<-do.call("rbind",lapply(res,function(x) x$est.rand.ef))
  avr.est.random.eff<-apply(estimated.random.eff,2, mean)


  #######################################################
  ## Inflection Point                               #####
  ## est.logT is the same as fixed effect theta2    #####
  ## est.logT=\hat(beta)*W_cov+\hat(ranom.effect)   #####
  #######################################################

  ###############
  # True logT  ##
  ###############

  TRUE.logT<-do.call("rbind",lapply(res,function(x) x$true.logT))
  tr.logT<-apply(TRUE.logT,2, mean)
  avr.tr.logT<-mean(tr.logT)

  #############
  # Estimates #
  #############

  estimated.logT<-do.call("rbind",lapply(res,function(x) x$est.logT))
  ind.abs.bias<-abs(estimated.logT- TRUE.logT) ## individual biases for each subject
  avr.ind.abs.bias<-apply(ind.abs.bias, 1, mean) ## average ind.abs.bias across subjects

  ##############
  ## bias logT #
  ##############
  avr.sim.ind.abs.bias<-mean(avr.ind.abs.bias)


  ##############
  ## sd logT #
  ##############
  #avr.est.logT<-apply(estimated.logT,2, mean)
  sd.est.logT<-apply(estimated.logT,2, sd) # sd of logT over simulation runs
  avr.sd.est.logT<-mean(sd.est.logT)       # avr of sd of logT across subjects

  ## bootstrap str

  boot.ind.sd<-do.call("rbind",lapply(res,function(x) x$ind.sd))
  avr.boot.ind.sd<-apply(boot.ind.sd,2, mean) # average of bootstrap sd over simulation runs
  avr.boot.str<-mean(avr.boot.ind.sd) # average of bootstrap sd across subjects


  # sd(results[[1]]$boot.est.logT[1,])
  boot.est.cp<-do.call("rbind",lapply(res,function(x) x$cp.boot))
  boot.ind.cp<-array(0, dim=c(n, 2, nsim), dimnames=list(paste("id",1:n,sep=""),paste(c("95%lower", "95%upper"))
                                                         ,paste("sim",1:nsim,sep="")))

  for(s in 1:nsim){
    # s<-1
    boot.ind.cp[,, s]<-boot.est.cp[(1+(s-1)*80):(s*80),]
  }

  avr.boot.ind.cp<-array(0, dim=c(n, 2), dimnames=list(paste("id",1:n,sep=""),paste(c("95%lower", "95%upper"))))

  for (id in 1:n){
    avr.boot.ind.cp[id,]<-apply(boot.ind.cp[id,,1:nsim],1,mean)
  }



  ###################################
  ## average of inflection points ###
  ###################################
  avr.est.logT<-apply(estimated.logT, 2, mean)
  glob.logT<-mean(avr.est.logT)

  ## global bootstrap 95% coverage probability
  ind.logT.cp<-cbind(avr.est.logT, avr.boot.ind.cp)
  glob.logT.cp<-apply(ind.logT.cp, 2, mean)
  #################################
  ### global logT summarization ###
  #################################


  glob.logT.est<-cbind(avr.sim.ind.abs.bias,  avr.sd.est.logT, avr.boot.str)
  colnames(glob.logT.est)<-c("Bias", "ESD", "ASE")
  rownames(glob.logT.est)<-c("est")
  glob.logT.est<-round(glob.logT.est, 3)

  Rel.glob.logT.est<-cbind(avr.sim.ind.abs.bias/avr.tr.logT,  avr.sd.est.logT/avr.tr.logT, avr.boot.str/avr.tr.logT, rbind(glob.logT.cp[2:3]))
  colnames(Rel.glob.logT.est)<-c("Rel.Bias", "Rel.ESD", "Rel.Boot.sd", "Boot.95%.cp.lower", "Boot.95%.cp.upper")
  rownames(Rel.glob.logT.est)<-c("est")
  Rel.glob.logT.est<-round(Rel.glob.logT.est, 2)


  ######################################
  ###### plot for estimated omega ######
  ######################################

  #num.interp<-45
  newl<-time.length


  true.tms<-array(0, dim=c(nsim, newl, n),
                  dimnames = list(paste("sim",1:nsim,sep=""),paste("visit",1:newl,sep=""), paste("id",1:id,sep="")
                  ))
  subj.logS<-true.tms
  est.tms<-true.tms

  ind.avr.true.tms<-array(0, dim=c(newl, n),
                          dimnames = list(paste("visit",1:newl,sep=""), paste("id",1:id,sep="")
                          ))

  ind.avr.true.xx<-ind.avr.true.tms
  ind.avr.est.tms<-ind.avr.true.tms

  for (id in 1:n){
    for (sim in 1:nsim){

      true.tms[sim,,id]<-res[[sim]]$new.pred.data[, "ture.pseudo.tms", id]
      est.tms[sim,,id]<-res[[sim]]$new.pred.data[, "est.pseudo.tms", id]
      subj.logS[sim,,id]<-res[[sim]]$new.pred.data[, "pseudo.logS", id]
    }
    ind.avr.true.tms[, id]<- apply(true.tms[,, id], 2, mean)
    ind.avr.est.tms[, id]<- apply(est.tms[,, id], 2, mean)
    ind.avr.true.xx[, id]<- apply(subj.logS[,, id], 2, mean)
  }


  mean.true.tms<-apply(ind.avr.true.tms,1,mean)
  mean.est.tms<-apply(ind.avr.est.tms,1,mean)
  mean.true.xx<-apply(ind.avr.true.xx,1,mean)




  ################################################
  # Plot of Average of Longitudinal trajectories #
  ################################################

  postscript(paste(file,"average_longitudinal_trajectory.eps",sep=""))
  par(mar=c(6.1,5.1,2.5,1.1))

  ## Create a blank plot

  ymin<-min(mean.true.tms)-0.1
  ymax<-max(mean.true.tms)+0.5
  xmin<-min(mean.true.xx)-0.3
  xmax<-max(mean.true.xx)+0.3

  #average of trajectories of all subjects
  plot( mean.true.xx,      mean.true.tms,  ylim=c(ymin,  ymax), xlim=c(xmin, xmax),xlab="x", ylab=expression("paraNLME"(omega(x))),
        cex.lab=1.3, cex.axis=1.0,
        main="Average of Individual Trajectories",  font=1.2,
        col="black", lty=1, lwd=2, type="l") # sub='(a)', cex.sub=1.3,

  lines( mean.true.xx,  mean.est.tms, col="red", lty=2, lwd=2)

  segments( glob.logT, 0,  glob.logT,median(mean.est.tms)+0.5,  col="blue", lty=3, lwd=2)

  legend( xmin, ymax,c("True Trajectory", "Estimated Trajectory"), col=c("black", "red"),
          lty=c(1,2),bty="n",cex=1.3,lwd=rep(2,2))
  mtext('(a)', outer=F,side=1,line=4.5, cex=1.3)

  dev.off()



  ###########################################
  ## First and Second Derivative Estimation #
  ###########################################


  ### True First Derivative : all simulations are the same at true values
  true.first.deriv<-apply(res[[1]]$true.first.deriv,1,mean)
  true.second.deriv<-apply(res[[1]]$true.second.deriv,1,mean)


  ## predicted values for omega at the population and subject levels

  subj.first.deriv<-array(0, dim=c(nsim, newl,n),
                          dimnames = list(paste("nsim",1:nsim,sep=""),paste("visit",1:newl,sep="")
                                          ,paste("ID",1:n,sep="")))

  subj.second.deriv<-subj.first.deriv
  mean.subj.first.deriv<-array(0, dim=c(newl,n),
                               dimnames = list(paste("visit",1:newl,sep="")
                                               ,paste("ID",1:n,sep="")))

  mean.subj.second.deriv<-mean.subj.first.deriv


  for (id in 1:n){

    subj.first.deriv[,,id]<-do.call("rbind",lapply(res,function(x) x$est.first.deriv[,id]))
    mean.subj.first.deriv[,id]<-apply(subj.first.deriv[,,id],2,mean)

    subj.second.deriv[,,id]<-do.call("rbind",lapply(res,function(x) x$est.second.deriv[,id]))
    mean.subj.second.deriv[,id]<-apply(subj.second.deriv[,,id],2,mean)

  }

  est.first.deriv<-apply(mean.subj.first.deriv,1,mean)
  est.second.deriv<-apply(mean.subj.second.deriv,1,mean)

  ## array list to store predicted tms, true tms, and visits
  logit.true.first<-true.first.deriv
  logit.est.first<-est.first.deriv

  ## array list to store predicted tms, true tms, and visits
  logit.true.second<-true.second.deriv
  logit.est.second<-est.second.deriv



  ####################################
  # Plot of Second Derivatives of Omega
  ####################################

  postscript(paste(file,"second_derivative.eps",sep=""))
  #pdf(file="avr_traj_logist_second.pdf")
  par(mar=c(6.1,5.1,2.5,1.1))
  #average of trajectories of all subjects
  ymin<-min(logit.est.second)-0.1
  ymax<-max(logit.est.second)+1.8

  plot(mean.true.xx, logit.true.second, ylim=c(ymin, ymax), xlim=c(xmin, xmax),
       xlab="x", ylab=expression("paraNLME" (partialdiff^{2}~omega/ partialdiff~x^2)), col="black", lty=1, type="l",
       main="Average Second Derivative of Trajectories", cex.lab=1.3, cex.axis=1.0,
       lwd=2,  font=1.3)
  lines(mean.true.xx, logit.est.second, col="red", lty=2, lwd=2)
  legend( xmin, ymax,c("True Trajectory", "Estimated Trajectory"), col=c("black", "red"),
          lty=c(1,2),bty="n",cex=1.3,lwd=rep(2,2))
  mtext('(b)', outer=F,side=1,line=4.5, cex=1.3)


  dev.off()

  return(list(fixed.est=fixed.est, #glob.logT=glob.logT,
              Rel.glob.logT.est=Rel.glob.logT.est



  ))

}


#' Summary results of the parametric NLME, when the true data are generated by the arctangent model and the parametric NLME assumes the correct arctangent  model.
#'
#' @param nsim number of simulation runs.
#' @param n number of sample size.
#' @param res a nsim-length of list for simulation results.
#' @param time.length number of data points at which predictors are required for each individual longitudinal trajectory. This time point for graphs to be plotted.
#' @param true.gam1 a true scale parameter for the arctangent model.
#' @param true.gam3 a true parameter, which determines the steepness of the arctangent model.
#' @param true.gam4 a true parameter, which determines the minimum and maximum of the arctangent model.
#' @param beta a (p-1)-length of true coefficient vector corresponding to subject specific covariates.
#' @param beta0 true intercept of the log-normal model for the inflection point.
#' @param file a character string of file name for a creating plot in .eps file extention.
#'
#'
#'
#' @return A list of the summarized simulation results including
#'\itemize{
#'        \item{res.table}{a data frame of the absolute biases, estimated standard deviations, average of the estimated standard errors and 95\% coverage probabilities for the fixed effects (\code{beta0}, \code{beta}).}
#'        \item{Rel.glob.logT.est}{a data frame of the relative average performance of inflection points for all subjects including the relative absolute biases, the relative empirical standard errors, the relative bootstrap standard deviations and 95\% boostrap confidence intervals.}
#'}
#'
#' @importFrom graphics plot
#' @importFrom graphics legend
#' @importFrom graphics lines
#' @importFrom graphics mtext
#' @importFrom graphics par
#' @importFrom graphics segments
#' @importFrom grDevices postscript
#' @importFrom grDevices dev.off
#' @export


arctan.results<-function(nsim=100, n=80, res=results, time.length=20, true.gam1=(2.45/pi), true.gam3=(pi/1.1), true.gam4=0.8,  beta=0.1, beta0=2,  file="nlme_arctan_"){

  #############################################
  ### Fixed effects                       #####
  ###    #theta1                          #####
  ###    #theta2 =beta*W_cov+ranom.effect #####
  #############################################

  ############
  ## estimates
  ############

  ## main fixed effects
  estimated.gam1<-do.call("rbind",lapply(res,function(x) x$est.gam1))
  estimated.gam3<-do.call("rbind",lapply(res,function(x) x$est.gam3))
  estimated.gam4<-do.call("rbind",lapply(res,function(x) x$est.gam4))

  ## baseline effects
  estimated.beta<-do.call("rbind",lapply(res,function(x) x$est.beta))
  estimated.beta0<-do.call("rbind",lapply(res,function(x) x$est.beta0))


  ## average of fixed effects
  mean.est.fixed.eff<-apply(cbind(estimated.gam1, estimated.gam3,estimated.gam4,estimated.beta, estimated.beta0),2, mean)

  ## standard deviation of fixed effects
  sd.est.fixed.eff<-apply(cbind(estimated.gam1, estimated.gam3,  estimated.gam4, estimated.beta, estimated.beta0),2, sd)

  ##
  true.fixed.eff<-c(true.gam1, true.gam3, true.gam4, beta, beta0)


  ####################################
  ### bias of fixed estimates
  ####################################
  abs.bias.fixed.eff<-abs(mean.est.fixed.eff-true.fixed.eff)


  #############################################
  ## standard errors of fixed effects
  #############################################
  estimated.str.gam1<-do.call("rbind",lapply(res,function(x) x$est.str.gam1))
  estimated.str.gam3<-do.call("rbind",lapply(res,function(x) x$est.str.gam3))
  estimated.str.gam4<-do.call("rbind",lapply(res,function(x) x$est.str.gam4))

  estimated.str.beta0<-do.call("rbind",lapply(res,function(x) x$est.str.beta0))
  estimated.str.beta<-do.call("rbind",lapply(res,function(x) x$est.str.beta))


  est.avr.str.gam1<-mean(estimated.str.gam1)
  est.avr.str.gam3<-mean(estimated.str.gam3)
  est.avr.str.gam4<-mean(estimated.str.gam4)
  est.avr.str.beta0<-mean(estimated.str.beta0)
  est.avr.str.beta<-mean(estimated.str.beta)

  ########################################
  ## coverage probabilities for theta1
  ########################################
  estimated.lowercp.gam1<-do.call("rbind",lapply(res,function(x) x$cp.lower.gam1))
  estimated.uppercp.gam1<-do.call("rbind",lapply(res,function(x) x$cp.upper.gam1))



  ind1gam1<- do.call("rbind",lapply(estimated.uppercp.gam1,function(x) 1*(x-(true.gam1)>0)))
  ind2gam1<- do.call("rbind",lapply(estimated.lowercp.gam1,function(x) 1*(x-(true.gam1)<0)))
  indgam1<-ind1gam1*ind2gam1

  cpgam1<-mean(indgam1)

  ########################################
  ## coverage probabilities for theta3 ###
  ########################################
  estimated.lowercp.gam3<-do.call("rbind",lapply(res,function(x) x$cp.lower.gam3))
  estimated.uppercp.gam3<-do.call("rbind",lapply(res,function(x) x$cp.upper.gam3))

  ind1gam3<- do.call("rbind",lapply(estimated.uppercp.gam3,function(x) 1*(x-(true.gam3)>0)))
  ind2gam3<- do.call("rbind",lapply(estimated.lowercp.gam3,function(x) 1*(x-(true.gam3)<0)))
  indgam3<-ind1gam3*ind2gam3

  cpgam3<-mean(indgam3)

  ########################################
  ## coverage probabilities for theta3 ###
  ########################################
  estimated.lowercp.gam4<-do.call("rbind",lapply(res,function(x) x$cp.lower.gam4))
  estimated.uppercp.gam4<-do.call("rbind",lapply(res,function(x) x$cp.upper.gam4))


  ind1gam4<- do.call("rbind",lapply(estimated.uppercp.gam4,function(x) 1*(x-(true.gam4)>0)))
  ind2gam4<- do.call("rbind",lapply(estimated.lowercp.gam4,function(x) 1*(x-(true.gam4)<0)))
  indgam4<-ind1gam4*ind2gam4

  cpgam4<-mean(indgam4)


  ########################################
  ## coverage probabilities for beta
  ########################################
  estimated.lowercp.beta0<-do.call("rbind",lapply(res,function(x) x$cp.lower.beta0))
  estimated.uppercp.beta0<-do.call("rbind",lapply(res,function(x) x$cp.upper.beta0))

  ind1b0<- do.call("rbind",lapply(estimated.uppercp.beta0,function(x) 1*(x-(beta0)>0)))
  ind2b0<- do.call("rbind",lapply(estimated.lowercp.beta0,function(x) 1*(x-(beta0)<0)))
  indb0<-ind1b0*ind2b0

  cpb0<-mean(indb0)



  ########################################
  ## coverage probabilities for beta
  ########################################
  estimated.lowercp.beta<-do.call("rbind",lapply(res,function(x) x$cp.lower.beta))
  estimated.uppercp.beta<-do.call("rbind",lapply(res,function(x) x$cp.upper.beta))

  ind1b<- do.call("rbind",lapply(estimated.uppercp.beta,function(x) 1*(x-(beta)>0)))
  ind2b<- do.call("rbind",lapply(estimated.lowercp.beta,function(x) 1*(x-(beta)<0)))
  indb<-ind1b*ind2b

  cpb<-mean(indb)



  ######################
  ### summarization ####
  ######################

  ## theta1
  estimates.gam1<-data.frame( abs.bias.fixed.eff["gam1"], sd.est.fixed.eff["gam1"],est.avr.str.gam1, cpgam1)
  colnames(estimates.gam1)<-c("abs.bias", "sd", "str", "cp")

  ## theta3
  estimates.gam3<-data.frame( abs.bias.fixed.eff["gam3"], sd.est.fixed.eff["gam3"],est.avr.str.gam3, cpgam3)
  colnames(estimates.gam3)<-c( "abs.bias", "sd", "str", "cp")

  ## theta4
  estimates.gam4<-data.frame(  abs.bias.fixed.eff["gam4"], sd.est.fixed.eff["gam4"],est.avr.str.gam4, cpgam4)
  colnames(estimates.gam4)<-c("abs.bias", "sd", "str", "cp")


  ## beta0
  estimates.beta0<-data.frame( abs.bias.fixed.eff["gam2.(Intercept)"], sd.est.fixed.eff["gam2.(Intercept)"],  est.avr.str.beta0,  cpb0)
  colnames(estimates.beta0)<-c("abs.bias", "sd", "str", "cp")


  ## beta
  estimates.beta<-data.frame(  abs.bias.fixed.eff["gam2.W.cov"], sd.est.fixed.eff["gam2.W.cov"],  est.avr.str.beta,  cpb)
  colnames(estimates.beta)<-c("abs.bias", "sd", "str", "cp")

  #fixed.est<-rbind(estimates.gam1,estimates.gam3, estimates.gam4, estimates.beta0, estimates.beta)
  #rownames(fixed.est)<-c("gam1", "gam3","gam4", "beta0", "beta")

  fixed.est<-rbind(estimates.beta0, estimates.beta)
  rownames(fixed.est)<-c( "beta0", "beta")
  fixed.est<-round(fixed.est, 2)


  ######################################
  ## within individual scale parameter #
  ######################################
  estimated.sigma.eps<-do.call("rbind",lapply(res,function(x) x$sig.eps))
  mean.est.sigma.eps<-mean(estimated.sigma.eps)
  sd.est.sigma.eps<-sd(estimated.sigma.eps)



  #############################################
  ### Random effect where                   ###
  ###    #theta2 =beta*W_cov+ranom.effect   ###
  #############################################

  ## random effect for Each ID (columwise) based on 100 simulation
  estimated.random.eff<-do.call("rbind",lapply(res,function(x) x$est.rand.ef))
  avr.est.random.eff<-apply(estimated.random.eff,2, mean)


  #######################################################
  ## Inflection Point                               #####
  ## est.logT is the same as fixed effect theta2    #####
  ## est.logT=\hat(beta)*W_cov+\hat(ranom.effect)   #####
  #######################################################

  ###############
  # True logT  ##
  ###############

  TRUE.logT<-do.call("rbind",lapply(res,function(x) x$true.logT))
  tr.logT<-apply(TRUE.logT,2, mean)
  avr.tr.logT<-mean(tr.logT)

  #############
  # Estimates #
  #############

  estimated.logT<-do.call("rbind",lapply(res,function(x) x$est.logT))
  ind.abs.bias<-abs(estimated.logT-TRUE.logT) ## individual biases for each subject
  avr.ind.abs.bias<-apply(ind.abs.bias, 1, mean) ## average ind.abs.bias across subjects


  ##############
  ## bias logT #
  ##############
  avr.sim.ind.abs.bias<-mean(avr.ind.abs.bias)

  ##############
  ## sd logT #
  ##############
  #avr.est.logT<-apply(estimated.logT,2, mean)
  sd.est.logT<-apply(estimated.logT,2, sd) # sd of logT over simulation runs
  avr.sd.est.logT<-mean(sd.est.logT)       # avr of sd of logT across subjects

  #bias.logT<-abs(avr.est.logT-tr.logT)   ## mean(bias.logT)

  ## bootstrap str
  #boot.est.logT=boot.est.logT,  ind.sd=ind.sd, cp.boot=cp.boot
  boot.ind.sd<-do.call("rbind",lapply(res,function(x) x$ind.sd))
  avr.boot.ind.sd<-apply(boot.ind.sd,2, mean) # average of bootstrap sd over simulation runs
  avr.boot.str<-mean(avr.boot.ind.sd) # average of bootstrap sd across subjects



  # sd(results[[1]]$boot.est.logT[1,])
  boot.est.cp<-do.call("rbind",lapply(res,function(x) x$cp.boot))
  boot.ind.cp<-array(0, dim=c(n, 2, nsim), dimnames=list(paste("id",1:n,sep=""),paste(c("95%lower", "95%upper"))
                                                         ,paste("sim",1:nsim,sep="")))

  for(s in 1:nsim){
    # s<-1
    boot.ind.cp[,, s]<-boot.est.cp[(1+(s-1)*80):(s*80),]
  }

  avr.boot.ind.cp<-array(0, dim=c(n, 2), dimnames=list(paste("id",1:n,sep=""),paste(c("95%lower", "95%upper"))))

  for (id in 1:n){
    avr.boot.ind.cp[id,]<-apply(boot.ind.cp[id,,1:nsim],1,mean)
  }



  ###################################
  ## average of inflection points ###
  ###################################
  avr.est.logT<-apply(estimated.logT, 2, mean)
  glob.logT<-mean(avr.est.logT)
  ## global bootstrap 95% coverage probability
  ind.logT.cp<-cbind(avr.est.logT, avr.boot.ind.cp)
  glob.logT.cp<-apply(ind.logT.cp, 2, mean)
  #################################
  ### global logT summarization ###
  #################################


  glob.logT.est<-cbind(avr.sim.ind.abs.bias,  avr.sd.est.logT, avr.boot.str)
  #colnames(glob.logT.est)<-c("Bias", "ESD", "ASE")
  #rownames(glob.logT.est)<-c("est")
  #glob.logT.est<-round(glob.logT.est, 3)

  Rel.glob.logT.est<-cbind(avr.sim.ind.abs.bias/avr.tr.logT,  avr.sd.est.logT/avr.tr.logT, avr.boot.str/avr.tr.logT, rbind(glob.logT.cp[2:3]))
  colnames(Rel.glob.logT.est)<-c("Rel.Bias", "Rel.ESD", "Rel.Boot.sd", "Boot.95%.cp.lower", "Boot.95%.cp.upper")
  rownames(Rel.glob.logT.est)<-c("est")
  Rel.glob.logT.est<-round(Rel.glob.logT.est, 2)


  ######################################
  ###### plot for estimated omega ######
  ######################################


  newl<-time.length



  true.tms<-array(0, dim=c(nsim, newl, n),
                  dimnames = list(paste("sim",1:nsim,sep=""),paste("visit",1:newl,sep=""), paste("id",1:id,sep="")
                  ))
  subj.logS<-true.tms
  est.tms<-true.tms

  ind.avr.true.tms<-array(0, dim=c(newl, n),
                          dimnames = list(paste("visit",1:newl,sep=""), paste("id",1:id,sep="")
                          ))

  ind.avr.true.xx<-ind.avr.true.tms
  ind.avr.est.tms<-ind.avr.true.tms

  for (id in 1:n){
    for (sim in 1:nsim){

      true.tms[sim,,id]<-res[[sim]]$new.pred.data[, "ture.pseudo.tms", id]
      est.tms[sim,,id]<-res[[sim]]$new.pred.data[, "est.pseudo.tms", id]
      subj.logS[sim,,id]<-res[[sim]]$new.pred.data[, "pseudo.logS", id]
    }
    ind.avr.true.tms[, id]<- apply(true.tms[,, id], 2, mean)
    ind.avr.est.tms[, id]<- apply(est.tms[,, id], 2, mean)
    ind.avr.true.xx[, id]<- apply(subj.logS[,, id], 2, mean)
  }


  mean.true.tms<-apply(ind.avr.true.tms,1,mean)
  mean.est.tms<-apply(ind.avr.est.tms,1,mean)
  mean.true.xx<-apply(ind.avr.true.xx,1,mean)




  ################################################
  # Plot of Average of Longitudinal trajectories #
  ################################################

  par(mar=c(6.1,5.1,2.5,1.1))

  ## Create a blank plot
  ymin<-min(mean.true.tms)-0.1
  ymax<-max(mean.true.tms)+0.5
  xmin<-min(mean.true.xx)-0.3
  xmax<-max(mean.true.xx)+0.3

  #average of trajectories of all subjects
  plot( mean.true.xx,      mean.true.tms,  ylim=c(ymin,  ymax), xlim=c(xmin, xmax), xlab="x", ylab=expression("paraNLME"(omega(x))),
        cex.lab=1.3, cex.axis=1.0,
        main="Average of Individual Trajectories",  font=1.2,
        col="black", lty=1, lwd=2, type="l") # sub='(a)', cex.sub=1.3,

  lines( mean.true.xx,  mean.true.tms, col="red", lty=2, lwd=2)

  segments(glob.logT, 0, glob.logT,median(mean.est.tms)+0.5,  col="blue", lty=3, lwd=2)

  legend( xmin, ymax, c("True Trajectory", "Estimated Trajectory"), col=c("black", "red"),
          lty=c(1,2),bty="n",cex=1.3,lwd=rep(2,2))
  mtext('(e)', outer=F,side=1,line=4.5, cex=1.3)

  dev.off()




  ###########################################
  ## First and Second Derivative Estimation #
  ###########################################


  ### True First Derivative : all simulations are the same at true values
  true.first.deriv<-apply(res[[1]]$true.first.deriv,1,mean)
  true.second.deriv<-apply(res[[1]]$true.second.deriv,1,mean)


  ## predicted values for omega at the population and subject levels

  subj.first.deriv<-array(0, dim=c(nsim, newl,n),
                          dimnames = list(paste("nsim",1:nsim,sep=""),paste("visit",1:newl,sep="")
                                          ,paste("ID",1:n,sep="")))

  subj.second.deriv<-subj.first.deriv
  mean.subj.first.deriv<-array(0, dim=c(newl,n),
                               dimnames = list(paste("visit",1:newl,sep="")
                                               ,paste("ID",1:n,sep="")))

  mean.subj.second.deriv<-mean.subj.first.deriv


  for (id in 1:n){

    subj.first.deriv[,,id]<-do.call("rbind",lapply(res,function(x) x$est.first.deriv[,id]))
    mean.subj.first.deriv[,id]<-apply(subj.first.deriv[,,id],2,mean)

    subj.second.deriv[,,id]<-do.call("rbind",lapply(res,function(x) x$est.second.deriv[,id]))
    mean.subj.second.deriv[,id]<-apply(subj.second.deriv[,,id],2,mean)

  }

  est.first.deriv<-apply(mean.subj.first.deriv,1,mean)
  est.second.deriv<-apply(mean.subj.second.deriv,1,mean)

  ## array list to store predicted tms, true tms, and visits
  arctan.true.first<-true.first.deriv
  arctan.est.first<-est.first.deriv

  ## array list to store predicted tms, true tms, and visits
  arctan.true.second<-true.second.deriv
  arctan.est.second<-est.second.deriv






  ####################################
  # Plot of Second Derivatives of Omega
  ####################################

  postscript(paste(file,"second_derivative.eps",sep=""))
  #pdf(file="arctan_true_second.pdf")
  par(mar=c(6.1,5.1,2.5,1.1))

  ymin<-min(logit.est.second)-0.1
  ymax<-max(logit.est.second)+2.8
  #average of trajectories of all subjects
  plot(mean.true.xx, logit.true.second, ylim=c(ymin, ymax), xlim=c(xmin, xmax),
       xlab="x",  ylab=expression("paraNLME" (partialdiff^{2}~omega/ partialdiff~x^2)), col="black", lty=1, type="l",
       main="Average Second Derivative of Trajectories", cex.lab=1.3, cex.axis=1.0,
       lwd=2,  font=1.3)
  lines(mean.true.xx, logit.est.second, col="red", lty=2, lwd=2)
  legend( xmin, ymax,c("True Trajectory", "Estimated Trajectory"), col=c("black", "red"),
          lty=c(1,2),bty="n",cex=1.3,lwd=rep(2,2))
  mtext('(f)', outer=F,side=1,line=4.5, cex=1.3)

  dev.off()




  return(list(fixed.est=fixed.est,
              Rel.glob.logT.est=Rel.glob.logT.est
  ))
}


#' Summary results of the parametric NLME, when the true data are generated by the logistic model, but the parametric NLME assumes the incorrect arctangent model.
#'
#' @param nsim number of simulation runs.
#' @param n number of sample size.
#' @param res a nsim-length of list for simulation results.
#' @param beta a (p-1)-length of true coefficient vector corresponding to subject specific covariates.
#' @param beta0 true intercept of the log-normal model for the inflection point.
#' @param time.length number of data points at which predictors are required for each individual longitudinal trajectory. This time point for graphs to be plotted.
#' @param file a character string of file name for a creating plot in .eps file extention.
#'
#' @return A list of the summarized simulation results including
#' \itemize{
#'        \item{res.table}{a data frame of the absolute biases, estimated standard deviations, average of the estimated standard errors and 95\% coverage probabilities for the fixed effects (\code{beta0}, \code{beta}).}
#'        \item{Rel.glob.logT.est}{a data frame of the relative average performance of inflection points for all subjects including the relative absolute biases, the relative empirical standard errors, the relative bootstrap standard deviations and 95\% boostrap confidence intervals.}
#'}
#'
#' @importFrom graphics plot
#' @importFrom graphics legend
#' @importFrom graphics lines
#' @importFrom graphics mtext
#' @importFrom graphics par
#' @importFrom graphics segments
#' @importFrom grDevices postscript
#' @importFrom grDevices dev.off
#' @export

missp.arctan.results<-function(nsim=100, n=80, res=results, time.length=time.length, beta=0.1, beta0=0.5, file="nlme_arctan_missp_"){


  #############################################
  ### Fixed effects                       #####
  ###    #theta1                          #####
  ###    #theta2 =beta*W_cov+ranom.effect #####
  #############################################

  ############
  ## estimates
  ############

  estimated.beta<-do.call("rbind",lapply(res,function(x) x$est.beta))
  estimated.beta0<-do.call("rbind",lapply(res,function(x) x$est.beta0))


  ## average of fixed effects
  mean.est.fixed.eff<-apply(cbind(estimated.beta, estimated.beta0),2, mean) #estimated.gam1, estimated.gam3,estimated.gam4,

  ## standard deviation of fixed effects
  sd.est.fixed.eff<-apply(cbind( estimated.beta, estimated.beta0),2, sd) #estimated.gam1, estimated.gam3,  estimated.gam4,

  ##
  true.fixed.eff<-c(beta, beta0) #gam1,gam3, gam4,


  ####################################
  ### bias of fixed estimates
  ####################################
  abs.bias.fixed.eff<-abs(mean.est.fixed.eff-true.fixed.eff)


  #############################################
  ## standard errors of fixed effects
  #############################################

  estimated.str.beta0<-do.call("rbind",lapply(res,function(x) x$est.str.beta0))
  estimated.str.beta<-do.call("rbind",lapply(res,function(x) x$est.str.beta))

  est.avr.str.beta0<-mean(estimated.str.beta0)
  est.avr.str.beta<-mean(estimated.str.beta)


  ########################################
  ## coverage probabilities for beta
  ########################################
  estimated.lowercp.beta0<-do.call("rbind",lapply(res,function(x) x$cp.lower.beta0))
  estimated.uppercp.beta0<-do.call("rbind",lapply(res,function(x) x$cp.upper.beta0))


  ind1b0<- do.call("rbind",lapply(estimated.uppercp.beta0,function(x) 1*(x-(beta0)>0)))
  ind2b0<- do.call("rbind",lapply(estimated.lowercp.beta0,function(x) 1*(x-(beta0)<0)))
  indb0<-ind1b0*ind2b0

  cpb0<-mean(indb0)



  ########################################
  ## coverage probabilities for beta
  ########################################
  estimated.lowercp.beta<-do.call("rbind",lapply(res,function(x) x$cp.lower.beta))
  estimated.uppercp.beta<-do.call("rbind",lapply(res,function(x) x$cp.upper.beta))

  ind1b<- do.call("rbind",lapply(estimated.uppercp.beta,function(x) 1*(x-(beta)>0)))
  ind2b<- do.call("rbind",lapply(estimated.lowercp.beta,function(x) 1*(x-(beta)<0)))
  indb<-ind1b*ind2b

  cpb<-mean(indb)



  ######################
  ### summarization ####
  ######################

  ## beta0
  #"est": mean.est.fixed.eff["gam2.(Intercept)"],
  estimates.beta0<-data.frame( abs.bias.fixed.eff["gam2.(Intercept)"], sd.est.fixed.eff["gam2.(Intercept)"],  est.avr.str.beta0,  cpb0)
  colnames(estimates.beta0)<-c( "BIAS", "ESD", "ASE", "CP")


  ## beta
  #"est": mean.est.fixed.eff["gam2.W.cov"],
  estimates.beta<-data.frame(  abs.bias.fixed.eff["gam2.W.cov"], sd.est.fixed.eff["gam2.W.cov"],  est.avr.str.beta,  cpb)
  colnames(estimates.beta)<-c( "BIAS", "ESD", "ASE", "CP")

  fixed.est<-rbind(estimates.beta0, estimates.beta)
  rownames(fixed.est)<-c("beta0", "beta") #"gam1", "gam3","gam4",
  fixed.est<-round(fixed.est, 2)

  ######################################
  ## within individual scale parameter #
  ######################################
  estimated.sigma.eps<-do.call("rbind",lapply(res,function(x) x$sig.eps))
  mean.est.sigma.eps<-mean(estimated.sigma.eps)
  sd.est.sigma.eps<-sd(estimated.sigma.eps)


  #############################################
  ### Random effect where                   ###
  ###    #theta2 =beta*W_cov+ranom.effect   ###
  #############################################

  ## random effect for Each ID (columwise) based on 100 simulation
  estimated.random.eff<-do.call("rbind",lapply(res,function(x) x$est.rand.ef))
  avr.est.random.eff<-apply(estimated.random.eff,2, mean)


  #######################################################
  ## Inflection Point                               #####
  ## est.logT is the same as fixed effect theta2    #####
  ## est.logT=\hat(beta)*W_cov+\hat(ranom.effect)   #####
  #######################################################

  ###############
  # True logT  ##
  ###############

  TRUE.logT<-do.call("rbind",lapply(res,function(x) x$true.logT))
  tr.logT<-apply(TRUE.logT,2, mean)
  avr.tr.logT<-mean(tr.logT)

  #############
  # Estimates #
  #############

  estimated.logT<-do.call("rbind",lapply(res,function(x) x$est.logT))
  ind.abs.bias<-abs(estimated.logT-TRUE.logT) ## individual biases for each subject
  avr.ind.abs.bias<-apply(ind.abs.bias, 1, mean) ## average ind.abs.bias across subjects

  ##############
  ## bias logT #
  ##############
  avr.sim.ind.abs.bias<-mean(avr.ind.abs.bias)

  ##############
  ## sd logT #
  ##############
  #avr.est.logT<-apply(estimated.logT,2, mean)
  sd.est.logT<-apply(estimated.logT,2, sd) # sd of logT over simulation runs
  avr.sd.est.logT<-mean(sd.est.logT)       # avr of sd of logT across subjects



  ## bootstrap str

  boot.ind.sd<-do.call("rbind",lapply(res,function(x) x$ind.sd))
  avr.boot.ind.sd<-apply(boot.ind.sd,2, mean) # average of bootstrap sd over simulation runs
  avr.boot.str<-mean(avr.boot.ind.sd) # average of bootstrap sd across subjects


  # sd(results[[1]]$boot.est.logT[1,])
  boot.est.cp<-do.call("rbind",lapply(res,function(x) x$cp.boot))
  boot.ind.cp<-array(0, dim=c(n, 2, nsim), dimnames=list(paste("id",1:n,sep=""),paste(c("95%lower", "95%upper"))
                                                         ,paste("sim",1:nsim,sep="")))

  for(s in 1:nsim){
    # s<-1
    boot.ind.cp[,, s]<-boot.est.cp[(1+(s-1)*80):(s*80),]
  }

  avr.boot.ind.cp<-array(0, dim=c(n, 2), dimnames=list(paste("id",1:n,sep=""),paste(c("95%lower", "95%upper"))))

  for (id in 1:n){
    avr.boot.ind.cp[id,]<-apply(boot.ind.cp[id,,1:nsim],1,mean)
  }



  ###################################
  ## average of inflection points ###
  ###################################
  avr.est.logT<-apply(estimated.logT, 2, mean)
  glob.logT<-mean(avr.est.logT)


  ## global bootstrap 95% coverage probability
  ind.logT.cp<-cbind(avr.est.logT, avr.boot.ind.cp)
  glob.logT.cp<-apply(ind.logT.cp, 2, mean)
  #################################
  ### global logT summarization ###
  #################################


  glob.logT.est<-cbind(avr.sim.ind.abs.bias,  avr.sd.est.logT, avr.boot.str)
  #colnames(glob.logT.est)<-c("Bias", "ESD", "ASE")
  #rownames(glob.logT.est)<-c("est")
  #glob.logT.est<-round(glob.logT.est, 3)

  Rel.glob.logT.est<-cbind(avr.sim.ind.abs.bias/avr.tr.logT,  avr.sd.est.logT/avr.tr.logT, avr.boot.str/avr.tr.logT, rbind(glob.logT.cp[2:3]))
  colnames(Rel.glob.logT.est)<-c("Rel.Bias", "Rel.ESD", "Rel.Boot.sd", "Boot.95%.cp.lower", "Boot.95%.cp.upper")
  rownames(Rel.glob.logT.est)<-c("est")
  Rel.glob.logT.est<-round(Rel.glob.logT.est, 2)



  ######################################
  ###### plot for estimated omega ######
  ######################################


  newl<-time.length



  true.tms<-array(0, dim=c(nsim, newl, n),
                  dimnames = list(paste("sim",1:nsim,sep=""),paste("visit",1:newl,sep=""), paste("id",1:n,sep="")
                  ))
  subj.logS<-true.tms
  est.tms<-true.tms

  ind.avr.true.tms<-array(0, dim=c(newl, n),
                          dimnames = list(paste("visit",1:newl,sep=""), paste("id",1:n,sep="")
                          ))

  ind.avr.true.xx<-ind.avr.true.tms
  ind.avr.est.tms<-ind.avr.true.tms

  for (id in 1:n){
    for (sim in 1:nsim){

      true.tms[sim,,id]<-res[[sim]]$new.pred.data[, "ture.pseudo.tms", id]
      est.tms[sim,,id]<-res[[sim]]$new.pred.data[, "est.pseudo.tms", id]
      subj.logS[sim,,id]<-res[[sim]]$new.pred.data[, "pseudo.logS", id]
    }
    ind.avr.true.tms[, id]<- apply(true.tms[,, id], 2, mean)
    ind.avr.est.tms[, id]<- apply(est.tms[,, id], 2, mean)
    ind.avr.true.xx[, id]<- apply(subj.logS[,, id], 2, mean)
  }


  mean.true.tms<-apply(ind.avr.true.tms,1,mean)
  mean.est.tms<-apply(ind.avr.est.tms,1,mean)
  mean.true.xx<-apply(ind.avr.true.xx,1,mean)




  ################################################
  # Plot of Average of Longitudinal trajectories #
  ################################################

  postscript(paste(file,"average_longitudinal_trajectory.eps",sep=""))
  #pdf(file="missp_arctan_true.pdf")
  par(mar=c(6.1,5.1,2.5,1.1))

  ## Create a blank plot

  #average of trajectories of all subjects

  ymin<-min(mean.true.tms)-0.1
  ymax<-max(mean.true.tms)+0.5
  xmin<-min(mean.true.xx)-0.3
  xmax<-max(mean.true.xx)+0.3

  plot( mean.true.xx,      mean.true.tms, ylim=c(ymin,  ymax), xlim=c(xmin, xmax), xlab="x", ylab=expression("paraNLME"(omega(x))),
        cex.lab=1.3, cex.axis=1.0,
        main="Average of Individual Trajectories",  font=1.2,
        col="black", lty=1, lwd=2, type="l") # sub='(a)', cex.sub=1.3,

  lines( mean.true.xx,  mean.est.tms, col="blue", lty=3, lwd=2)

  segments( glob.logT, 0,  glob.logT,median(mean.est.tms)+0.5,  col="orange", lty=4, lwd=2)

  legend(xmin, ymax,c("True Trajectory", "Estimated Trajectory"), col=c("black", "blue"),
         lty=c(1,3),bty="n",cex=1.3,lwd=rep(2,2))
  mtext('(a)', outer=F,side=1,line=4.5, cex=1.3)

  dev.off()




  ###########################################
  ## First and Second Derivative Estimation #
  ###########################################


  ### True First Derivative : all simulations are the same at true values
  true.first.deriv<-apply(res[[1]]$true.first.deriv,1,mean)
  true.second.deriv<-apply(res[[1]]$true.second.deriv,1,mean)


  ## predicted values for omega at the population and subject levels

  subj.first.deriv<-array(0, dim=c(nsim, newl,n),
                          dimnames = list(paste("nsim",1:nsim,sep=""),paste("visit",1:newl,sep="")
                                          ,paste("ID",1:n,sep="")))

  subj.second.deriv<-subj.first.deriv
  mean.subj.first.deriv<-array(0, dim=c(newl,n),
                               dimnames = list(paste("visit",1:newl,sep="")
                                               ,paste("ID",1:n,sep="")))

  mean.subj.second.deriv<-mean.subj.first.deriv

  #res[[1]]$est.first.deriv
  for (id in 1:n){

    subj.first.deriv[,,id]<-do.call("rbind",lapply(res,function(x) x$est.first.deriv[,id]))
    mean.subj.first.deriv[,id]<-apply(subj.first.deriv[,,id],2,mean)

    subj.second.deriv[,,id]<-do.call("rbind",lapply(res,function(x) x$est.second.deriv[,id]))
    mean.subj.second.deriv[,id]<-apply(subj.second.deriv[,,id],2,mean)

  }

  est.first.deriv<-apply(mean.subj.first.deriv,1,mean)
  est.second.deriv<-apply(mean.subj.second.deriv,1,mean)

  ## array list to store predicted tms, true tms, and visits
  logit.true.first<-true.first.deriv
  logit.est.first<-est.first.deriv

  ## array list to store predicted tms, true tms, and visits
  logit.true.second<-true.second.deriv
  logit.est.second<-est.second.deriv


  ####################################
  # Plot of Second Derivatives of Omega
  ####################################

  postscript(paste(file,"second_derivative.eps",sep=""))
  #pdf(file="missp_arctan_second.pdf")
  par(mar=c(6.1,5.1,2.5,1.1))

  ymin<-min(logit.est.second)-0.1
  ymax<-max(logit.est.second)+1.3
  #average of trajectories of all subjects
  plot(mean.true.xx, logit.true.second,ylim=c(ymin, ymax), xlim=c(xmin, xmax),
       xlab="x", ylab=expression("paraNLME" (partialdiff^{2}~omega/ partialdiff~x^2)), col="black", lty=1, type="l",
       main="Average Second Derivative of Trajectories", cex.lab=1.3, cex.axis=1.0,
       lwd=2,  font=1.3)
  lines(mean.true.xx, logit.est.second, col="blue", lty=3, lwd=2)
  mtext('(b)', outer=F,side=1,line=4.5, cex=1.3)

  dev.off()



  return(list( fixed.est=fixed.est, Rel.glob.logT.est=Rel.glob.logT.est
               #glob.logT.est=glob.logT.est


  ))

}


#' Summary results of the parametric NLME, when the true data are generated by the arctangent model, but the parametric NLME assumes the incorrect logistic model.
#'
#' @param nsim number of simulation runs.
#' @param n number of sample size.
#' @param res a nsim-length of list for simulation results.
#' @param beta a (p-1)-length of true coefficient vector corresponding to subject specific covariates.
#' @param beta0 true intercept of the log-normal model for the inflection point.
#' @param time.length number of data points at which predictors are required for each individual longitudinal trajectory. This time point for graphs to be plotted.
#' @param file a character string of file name for a creating plot in .eps file extention.
#'
#' @return A list of the summarized simulation results including
#' \itemize{
#'        \item{res.table}{a data frame of the absolute biases, estimated standard deviations, average of the estimated standard errors and 95\% coverage probabilities for the fixed effects (\code{beta0}, \code{beta}).}
#'        \item{Rel.glob.logT.est}{a data frame of the relative average performance of inflection points for all subjects including the relative absolute biases, the relative empirical standard errors, the relative bootstrap standard deviations and 95\% boostrap confidence intervals.}
#' }
#'
#' @importFrom graphics plot
#' @importFrom graphics legend
#' @importFrom graphics lines
#' @importFrom graphics mtext
#' @importFrom graphics par
#' @importFrom graphics segments
#' @importFrom grDevices postscript
#' @importFrom grDevices dev.off
#' @export


missp.logist.results<-function(nsim=100, n=80, res=results, time.length=20, beta=0.1, beta0=2, file="nlme_logisic"){


  #############################################
  ### Fixed effects                       #####
  ###    #theta1                          #####
  ###    #theta2 =beta*W_cov+ranom.effect #####
  #############################################

  #######################
  # TRUE Fixed Effects ##
  #######################

  true.fixed.eff<-c( beta, beta0)



  ##############
  ## estimates #
  ##############

  ## baseline effects
  estimated.beta<-do.call("rbind",lapply(res,function(x) x$est.beta))
  estimated.beta0<-do.call("rbind",lapply(res,function(x) x$est.beta0))


  ## average of fixed effects
  mean.est.fixed.eff<-apply(cbind( estimated.beta, estimated.beta0),2, mean)

  ## standard deviation of fixed effects
  sd.est.fixed.eff<-apply(cbind(  estimated.beta, estimated.beta0),2, sd)



  ####################################
  ### bias of fixed estimates
  ####################################
  abs.bias.fixed.eff<-abs(mean.est.fixed.eff-true.fixed.eff)


  #############################################
  ## standard errors of fixed effects
  #############################################

  estimated.str.beta0<-do.call("rbind",lapply(res,function(x) x$est.str.beta0))
  estimated.str.beta<-do.call("rbind",lapply(res,function(x) x$est.str.beta))



  est.avr.str.beta0<-mean(estimated.str.beta0)
  est.avr.str.beta<-mean(estimated.str.beta)


  ########################################
  ## coverage probabilities for theta1
  ########################################


  ########################################
  ## coverage probabilities for beta
  ########################################
  estimated.lowercp.beta0<-do.call("rbind",lapply(res,function(x) x$cp.lower.beta0))
  estimated.uppercp.beta0<-do.call("rbind",lapply(res,function(x) x$cp.upper.beta0))


  ind1b0<- do.call("rbind",lapply(estimated.uppercp.beta0,function(x) 1*(x-(beta0)>0)))
  ind2b0<- do.call("rbind",lapply(estimated.lowercp.beta0,function(x) 1*(x-(beta0)<0)))
  indb0<-ind1b0*ind2b0

  cpb0<-mean(indb0)


  ########################################
  ## coverage probabilities for beta
  ########################################
  estimated.lowercp.beta<-do.call("rbind",lapply(res,function(x) x$cp.lower.beta))
  estimated.uppercp.beta<-do.call("rbind",lapply(res,function(x) x$cp.upper.beta))

  ind1b<- do.call("rbind",lapply(estimated.uppercp.beta,function(x) 1*(x-(beta)>0)))
  ind2b<- do.call("rbind",lapply(estimated.lowercp.beta,function(x) 1*(x-(beta)<0)))
  indb<-ind1b*ind2b

  cpb<-mean(indb)

  ######################
  ### summarization ####
  ######################

  ## beta0
  estimates.beta0<-data.frame( abs.bias.fixed.eff["theta2.(Intercept)"], sd.est.fixed.eff["theta2.(Intercept)"],  est.avr.str.beta0,  cpb0)
  colnames(estimates.beta0)<-c( "abs.bias", "sd", "str", "cp")


  ## beta
  estimates.beta<-data.frame(  abs.bias.fixed.eff["theta2.W.cov"], sd.est.fixed.eff["theta2.W.cov"],  est.avr.str.beta,  cpb)
  colnames(estimates.beta)<-c( "abs.bias", "sd", "str", "cp")

  fixed.est<-rbind( estimates.beta0, estimates.beta) #estimates.theta1,estimates.theta3,
  rownames(fixed.est)<-c( "beta0", "beta") #"theta1", "theta3",

  fixed.est<-round(fixed.est,2)

  ######################################
  ## within individual scale parameter #
  ######################################
  estimated.sigma.eps<-do.call("rbind",lapply(res,function(x) x$sig.eps))
  mean.est.sigma.eps<-mean(estimated.sigma.eps)
  sd.est.sigma.eps<-sd(estimated.sigma.eps)


  #############################################
  ### Random effect where                   ###
  ###    #theta2 =beta*W_cov+ranom.effect   ###
  #############################################

  ## random effect for Each ID (columwise) based on 100 simulation
  estimated.random.eff<-do.call("rbind",lapply(res,function(x) x$est.rand.ef))
  avr.est.random.eff<-apply(estimated.random.eff,2, mean)


  #######################################################
  ## Inflection Point                               #####
  ## est.logT is the same as fixed effect theta2    #####
  ## est.logT=\hat(beta)*W_cov+\hat(ranom.effect)   #####
  #######################################################

  ###############
  # True logT  ##
  ###############

  TRUE.logT<-do.call("rbind",lapply(res,function(x) x$true.logT))
  tr.logT<-apply(TRUE.logT,2, mean)
  avr.tr.logT<-mean(tr.logT)

  #############
  # Estimates #
  #############

  estimated.logT<-do.call("rbind",lapply(res,function(x) x$est.logT))
  ind.abs.bias<-abs(estimated.logT-TRUE.logT) ## individual biases for each subject
  avr.ind.abs.bias<-apply(ind.abs.bias, 1, mean) ## average ind.abs.bias across subjects


  ##############
  ## bias logT #
  ##############
  avr.sim.ind.abs.bias<-mean(avr.ind.abs.bias)

  ##############
  ## sd logT #
  ##############
  #avr.est.logT<-apply(estimated.logT,2, mean)
  sd.est.logT<-apply(estimated.logT,2, sd) # sd of logT over simulation runs
  avr.sd.est.logT<-mean(sd.est.logT)       # avr of sd of logT across subjects


  ## bootstrap str
  #boot.est.logT=boot.est.logT,  ind.sd=ind.sd, cp.boot=cp.boot
  boot.ind.sd<-do.call("rbind",lapply(res,function(x) x$ind.sd))
  avr.boot.ind.sd<-apply(boot.ind.sd,2, mean) # average of bootstrap sd over simulation runs
  avr.boot.str<-mean(avr.boot.ind.sd) # average of bootstrap sd across subjects


  # sd(results[[1]]$boot.est.logT[1,])
  boot.est.cp<-do.call("rbind",lapply(res,function(x) x$cp.boot))
  boot.ind.cp<-array(0, dim=c(n, 2, nsim), dimnames=list(paste("id",1:n,sep=""),paste(c("95%lower", "95%upper"))
                                                         ,paste("sim",1:nsim,sep="")))

  for(s in 1:nsim){
    # s<-1
    boot.ind.cp[,, s]<-boot.est.cp[(1+(s-1)*80):(s*80),]
  }

  avr.boot.ind.cp<-array(0, dim=c(n, 2), dimnames=list(paste("id",1:n,sep=""),paste(c("95%lower", "95%upper"))))

  for (id in 1:n){
    avr.boot.ind.cp[id,]<-apply(boot.ind.cp[id,,1:nsim],1,mean)
  }



  ###################################
  ## average of inflection points ###
  ###################################
  avr.est.logT<-apply(estimated.logT, 2, mean)
  glob.logT<-mean(avr.est.logT)

  ## global bootstrap 95% coverage probability
  ind.logT.cp<-cbind(avr.est.logT, avr.boot.ind.cp)
  glob.logT.cp<-apply(ind.logT.cp, 2, mean)
  #################################
  ### global logT summarization ###
  #################################


  glob.logT.est<-cbind(avr.sim.ind.abs.bias,  avr.sd.est.logT, avr.boot.str)
  #colnames(glob.logT.est)<-c("Bias", "ESD", "ASE")
  #rownames(glob.logT.est)<-c("est")
  #glob.logT.est<-round(glob.logT.est, 3)

  Rel.glob.logT.est<-cbind(avr.sim.ind.abs.bias/avr.tr.logT,  avr.sd.est.logT/avr.tr.logT, avr.boot.str/avr.tr.logT, rbind(glob.logT.cp[2:3]))
  colnames(Rel.glob.logT.est)<-c("Rel.Bias", "Rel.ESD", "Rel.Boot.sd", "Boot.95%.cp.lower", "Boot.95%.cp.upper")
  rownames(Rel.glob.logT.est)<-c("est")
  Rel.glob.logT.est<-round(Rel.glob.logT.est, 2)


  ######################################
  ###### plot for estimated omega ######
  ######################################


  newl<-time.length



  true.tms<-array(0, dim=c(nsim, newl, n),
                  dimnames = list(paste("sim",1:nsim,sep=""),paste("visit",1:newl,sep=""), paste("id",1:id,sep="")
                  ))
  subj.logS<-true.tms
  est.tms<-true.tms

  ind.avr.true.tms<-array(0, dim=c(newl, n),
                          dimnames = list(paste("visit",1:newl,sep=""), paste("id",1:id,sep="")
                          ))

  ind.avr.true.xx<-ind.avr.true.tms
  ind.avr.est.tms<-ind.avr.true.tms

  for (id in 1:n){
    for (sim in 1:nsim){

      true.tms[sim,,id]<-res[[sim]]$new.pred.data[, "ture.pseudo.tms", id]
      est.tms[sim,,id]<-res[[sim]]$new.pred.data[, "est.pseudo.tms", id]
      subj.logS[sim,,id]<-res[[sim]]$new.pred.data[, "pseudo.logS", id]
    }
    ind.avr.true.tms[, id]<- apply(true.tms[,, id], 2, mean)
    ind.avr.est.tms[, id]<- apply(est.tms[,, id], 2, mean)
    ind.avr.true.xx[, id]<- apply(subj.logS[,, id], 2, mean)
  }


  mean.true.tms<-apply(ind.avr.true.tms,1,mean)
  mean.est.tms<-apply(ind.avr.est.tms,1,mean)
  mean.true.xx<-apply(ind.avr.true.xx,1,mean)




  ################################################
  # Plot of Average of Longitudinal trajectories #
  ################################################

  postscript(paste(file,"average_longitudinal_trajectory.eps",sep=""))
  par(mar=c(6.1,5.1,2.5,1.1))

  ## Create a blank plot
  ymax<-max(mean.true.tms)+0.5
  xmin<-min(mean.true.xx)-0.3
  xmax<-max(mean.true.xx)+0.3
  ymin<-min(mean.true.tms)-0.1

  #average of trajectories of all subjects

  plot( mean.true.xx, mean.true.tms, ylim=c(ymin,  ymax), xlim=c(xmin, xmax), xlab="x", ylab=expression("paraNLME"(omega(x))),
        cex.lab=1.3, cex.axis=1.0,
        main="Average of Individual Trajectories",  font=1.2,
        col="black", lty=1, lwd=2, type="l") # sub='(a)', cex.sub=1.3,

  lines( mean.true.xx,  mean.est.tms, col="blue", lty=3, lwd=2)

  segments(glob.logT, 0, glob.logT, median(mean.est.tms)+0.5,  col="orange", lty=4, lwd=2)

  legend(xmin, ymax,c("True Trajectory", "Estimated Trajectory"), col=c("black", "blue"),
         lty=c(1,3),bty="n",cex=1.3,lwd=rep(2,2))
  mtext('(e)', outer=F,side=1,line=4.5, cex=1.3)

  dev.off()



  ###########################################
  ## First and Second Derivative Estimation #
  ###########################################


  ### True First Derivative : all simulations are the same at true values
  true.first.deriv<-apply(res[[1]]$true.first.deriv,1,mean)
  true.second.deriv<-apply(res[[1]]$true.second.deriv,1,mean)


  ## predicted values for omega at the population and subject levels

  subj.first.deriv<-array(0, dim=c(nsim, newl,n),
                          dimnames = list(paste("nsim",1:nsim,sep=""),paste("visit",1:newl,sep="")
                                          ,paste("ID",1:n,sep="")))

  subj.second.deriv<-subj.first.deriv
  mean.subj.first.deriv<-array(0, dim=c(newl,n),
                               dimnames = list(paste("visit",1:newl,sep="")
                                               ,paste("ID",1:n,sep="")))

  mean.subj.second.deriv<-mean.subj.first.deriv


  for (id in 1:n){

    subj.first.deriv[,,id]<-do.call("rbind",lapply(res,function(x) x$est.first.deriv[,id]))
    mean.subj.first.deriv[,id]<-apply(subj.first.deriv[,,id],2,mean)

    subj.second.deriv[,,id]<-do.call("rbind",lapply(res,function(x) x$est.second.deriv[,id]))
    mean.subj.second.deriv[,id]<-apply(subj.second.deriv[,,id],2,mean)

  }

  est.first.deriv<-apply(mean.subj.first.deriv,1,mean)
  est.second.deriv<-apply(mean.subj.second.deriv,1,mean)

  ## array list to store predicted tms, true tms, and visits
  arctan.true.first<-true.first.deriv
  logit.est.first<-est.first.deriv

  ## array list to store predicted tms, true tms, and visits
  arctan.true.second<-true.second.deriv
  logit.est.second<-est.second.deriv


  ####################################
  # Plot of Second Derivatives of Omega
  ####################################

  postscript(paste(file,"second_derivative.eps",sep=""))
  #pdf(file="arctan_missp_logsit_sec.pdf")
  par(mar=c(6.1,5.1,2.5,1.1))

  ymin<-min(logit.true.second)-0.8
  ymax<-max(logit.true.second)+2.3
  #average of trajectories of all subjects
  plot(mean.true.xx, logit.true.second, ylim=c(ymin, ymax), xlim=c(xmin, xmax),
       xlab="x", ylab=expression("paraNLME" (partialdiff^{2}~omega/ partialdiff~x^2)), col="black", lty=1, type="l",
       main="Average Second Derivative of Trajectories", cex.lab=1.3, cex.axis=1.0,
       lwd=2,  font=1.3)
  lines(mean.true.xx, logit.est.second, col="blue", lty=3, lwd=2)
  mtext('(f)', outer=F,side=1,line=4.5, cex=1.3)

  dev.off()



  return(list(fixed.est=fixed.est,
              Rel.glob.logT.est=Rel.glob.logT.est
              #glob.logT.est=glob.logT.est
  ))

}