R/collectedfunctions.R

Defines functions cov_fn cov_plots gof_plots vpc_plots lhcut lhfactor stackvar lhmutate lhmerge_wide_data lhwide lhlong findiff lhjoin lhorder txt install.pack lhtemplate lhdoc lhcbind lhloess loadpack m6 addvar addtime rt2tad nmid diftm format_time tab chclass one nodup duprow hl1cpt hl2cpt hl3cpt cround1 dup1 dup2 bround geom geocv cv se cilow ciup nmiss lhcontab roundbatch indiv.tab

Documented in addtime addvar bround chclass cilow ciup cov_fn cov_plots cround1 cv diftm dup1 dup2 duprow findiff format_time geocv geom gof_plots hl1cpt hl2cpt hl3cpt indiv.tab install.pack lhcbind lhcontab lhcut lhdoc lhfactor lhjoin lhloess lhlong lhmerge_wide_data lhmutate lhorder lhtemplate lhwide loadpack m6 nmid nmiss nodup one roundbatch rt2tad se stackvar tab txt vpc_plots

#' Covariate plots function
#'
#' @param data Data frame. Use the codes below as template 
#' @keywords cov_plots()
#' @export
#' @examples
cov_fn<-function(...){
library(grid)
require(gpairs)
gpairs<-function (x, upper.pars = list(scatter = "points", conditional = "barcode", 
                                       mosaic = "mosaic"), lower.pars = list(scatter = "points", 
                                                                             conditional = "boxplot", mosaic = "mosaic"), diagonal = "default", 
                  outer.margins = list(bottom = unit(2, "lines"), left = unit(2, 
                                                                              "lines"), top = unit(2, "lines"), right = unit(2, "lines")), 
                  xylim = NULL, outer.labels = NULL, outer.rot = c(0, 90), 
                  gap = 0.05, buffer = 0.02, reorder = NULL, cluster.pars = NULL, 
                  stat.pars = NULL, scatter.pars = NULL, bwplot.pars = NULL, 
                  stripplot.pars = NULL, barcode.pars = NULL, mosaic.pars = NULL, 
                  axis.pars = NULL, diag.pars = NULL, whatis = FALSE) 
{
  if (!is.data.frame(x)) {
    if (is.matrix(x)) 
      x <- as.data.frame(x)
    else stop("What did you give me? You might want to use Excel. (Only one column in argument to gpairs.\n\n")
  }
  zc <- function(x) length(unique(x)) <= 1
  if (any(sapply(x, zc), na.rm = TRUE)) {
    warning(paste(sum(sapply(x, zc), na.rm = TRUE), "columns with less than two distinct values eliminated"))
    x <- x[, !(sapply(x, zc))]
  }
  if (!is.null(lower.pars) & !is.list(lower.pars)) {
    warning("lower.pars is not a list, proceed with caution.")
  }
  if (!is.null(upper.pars) & !is.list(upper.pars)) {
    warning("upper.pars is not a list, proceed with caution.")
  }
  if (!is.null(reorder)) {
    if (pmatch(reorder, "cluster", nomatch = FALSE)) {
      if (is.null(cluster.pars)) {
        cluster.pars <- list(dist.method = "euclidean", 
                             hclust.method = "complete")
      }
      x.num <- as.matrix(as.data.frame(lapply(x, as.numeric)))
      x.clust <- hclust(dist(t(x.num), method = cluster.pars$dist.method), 
                        method = cluster.pars$hclust.method)
      x <- x[, x.clust$order]
    }
  }
  if (is.null(lower.pars$scatter.pars)) {
    lower.pars$scatter.pars <- "points"
  }
  if (is.null(lower.pars$conditional)) {
    lower.pars$conditional <- "boxplot"
  }
  if (is.null(lower.pars$mosaic)) {
    lower.pars$mosaic <- "mosaic"
  }
  if (is.null(upper.pars$scatter.pars)) {
    upper.pars$scatter.pars <- "points"
  }
  if (is.null(upper.pars$conditional)) {
    upper.pars$conditional <- "barcode"
  }
  if (is.null(upper.pars$mosaic)) {
    upper.pars$mosaic <- "mosaic"
  }
  if (!is.list(outer.margins)) {
    if (length(outer.margins) == 4) {
      if (is.unit(outer.margins[1])) {
        outer.margins <- list(bottom = outer.margins[1], 
                              left = outer.margins[2], top = outer.margins[3], 
                              right = outer.margins[4])
      } else {
        outer.margins <- list(bottom = unit(outer.margins[1], 
                                            "lines"), left = unit(outer.margins[2], "lines"), 
                              top = unit(outer.margins[3], "lines"), right = unit(outer.margins[4], 
                                                                                  "lines"))
      }
    } else {
      stop("outer.margins are not valid.")
    }
  }
  if (is.null(outer.labels)) {
    outer.labels$top <- rep(FALSE, ncol(x))
    outer.labels$top[seq(2, ncol(x), by = 2)] <- TRUE
    outer.labels$left <- rep(FALSE, ncol(x))
    outer.labels$left[seq(2, ncol(x), by = 2)] <- TRUE
    outer.labels$right <- !outer.labels$left
    outer.labels$bottom <- !outer.labels$top
  } else {
    if (pmatch(as.character(outer.labels), "all", nomatch = FALSE)) {
      all.labeling <- TRUE
    } else if (pmatch(as.character(outer.labels), "none", nomatch = FALSE)) {
      all.labeling <- FALSE
    } else {
      stop("argument to outer.labels not understood\n")
    }
    outer.labels <- NULL
    outer.labels$top <- rep(all.labeling, ncol(x))
    outer.labels$left <- rep(all.labeling, ncol(x))
    outer.labels$bottom <- rep(all.labeling, ncol(x))
    outer.labels$right <- rep(all.labeling, ncol(x))
  }
  if (is.null(stat.pars$fontsize)) {
    stat.pars$fontsize <- 7
  }
  if (is.null(stat.pars$signif)) {
    stat.pars$signif <- 0.05
  }
  if (is.null(stat.pars$verbose)) {
    stat.pars$verbose <- FALSE
  }
  if (is.null(stat.pars$use.color)) {
    stat.pars$use.color <- TRUE
  }
  if (is.null(stat.pars$missing)) {
    stat.pars$missing <- "missing"
  }
  if (is.null(stat.pars$just)) {
    stat.pars$just <- "centre"
  }
  if (is.null(scatter.pars$pch)) {
    scatter.pars$pch <- 1
  }
  if (is.null(scatter.pars$size)) {
    scatter.pars$size <- unit(0.25, "char")
  }
  if (is.null(scatter.pars$col)) {
    scatter.pars$col <- "black"
  }
  if (is.null(scatter.pars$plotpoints)) {
    scatter.pars$plotpoints <- TRUE
  }
  if (is.null(axis.pars$n.ticks)) {
    axis.pars$n.ticks <- 5
  }
  if (is.null(axis.pars$fontsize)) {
    axis.pars$fontsize <- 9
  }
  if (axis.pars$n.ticks < 3) {
    axis.pars$n.ticks <- 3
    warning("Fewer than 3 axis ticks might cause problems.")
  }
  if (is.null(diag.pars$fontsize)) {
    diag.pars$fontsize <- 9
  }
  if (is.null(diag.pars$show.hist)) {
    diag.pars$show.hist <- TRUE
  }
  if (is.null(diag.pars$hist.color)) {
    diag.pars$hist.color <- "black"
  }
  if (is.null(stripplot.pars$pch)) {
    stripplot.pars$pch <- 1
  }
  if (is.null(stripplot.pars$size)) {
    stripplot.pars$size <- unit(0.5, "char")
  }
  if (is.null(stripplot.pars$col)) {
    stripplot.pars$col <- "black"
  }
  if (is.null(stripplot.pars$jitter)) {
    stripplot.pars$jitter <- FALSE
  }
  if (is.null(barcode.pars$nint)) {
    barcode.pars$nint <- 0
  }
  if (is.null(barcode.pars$ptsize)) {
    barcode.pars$ptsize <- unit(0.25, "char")
  }
  if (is.null(barcode.pars$ptpch)) {
    barcode.pars$ptpch <- 1
  }
  if (is.null(barcode.pars$bcspace)) {
    barcode.pars$bcspace <- NULL
  }
  if (is.null(barcode.pars$use.points)) {
    barcode.pars$use.points <- FALSE
  }
  if (is.null(mosaic.pars$gp_labels)) {
    mosaic.pars$gp_labels <- gpar(fontsize = 9)
  }
  if (is.null(mosaic.pars$gp_args)) {
    mosaic.pars$gp_args <- list()
  }
  draw.axis <- function(x, y, axis.pars, xpos, ypos, cat.labels = NULL, 
                        horiz = NULL, xlim = NULL, ylim = NULL) {
    x <- as.numeric(x)
    y <- as.numeric(y)
    if (is.null(xlim)) {
      px <- pretty(x, axis.pars$n.ticks)
      px <- px[px > min(x, na.rm = TRUE) & px < max(x, 
                                                    na.rm = TRUE)]
    } else {
      px <- pretty(xlim, axis.pars$n.ticks)
      px <- px[px > min(xlim, na.rm = TRUE) & px < max(xlim, 
                                                       na.rm = TRUE)]
    }
    if (is.null(ylim)) {
      py <- pretty(y, axis.pars$n.ticks)
      py <- py[py > min(y, na.rm = TRUE) & py < max(y, 
                                                    na.rm = TRUE)]
    } else {
      py <- pretty(ylim, axis.pars$n.ticks)
      py <- py[py > min(ylim, na.rm = TRUE) & py < max(ylim, 
                                                       na.rm = TRUE)]
    }
    k <- length(cat.labels)
    if (!is.null(xpos)) {
      if (!is.null(cat.labels) && !horiz) {
        grid.text(cat.labels, x = unit(1:k, "native"), 
                  y = unit(rep(1 * (1 - xpos), k), "npc") + unit(rep(-1 * 
                                                                       xpos + 1 * (1 - xpos), k), "lines"), rot = outer.rot[1], 
                  gp = gpar(fontsize = axis.pars$fontsize))
      } else grid.xaxis(at = px, gp = gpar(fontsize = axis.pars$fontsize), 
                        main = xpos)
    }
    if (!is.null(ypos)) {
      if (!is.null(cat.labels) && horiz) {
        grid.text(cat.labels, y = unit(1:k, "native"), 
                  x = unit(rep(1 * (1 - ypos), k), "npc") + unit(rep(-1 * 
                                                                       ypos + 1 * (1 - ypos), k), "lines"), rot = outer.rot[2], 
                  gp = gpar(fontsize = axis.pars$fontsize))
      }else grid.yaxis(at = py, gp = gpar(fontsize = axis.pars$fontsize), 
                       main = ypos)
    }
  }
  qq.panel <- function(x, y, scatter.pars, axis.pars, xpos, 
                       ypos, xlim, ylim) {
    pushViewport(viewport(xscale = xlim, yscale = ylim))
    draw.axis(x, y, axis.pars, xpos, ypos, NULL, NULL, xlim, 
              ylim)
    popViewport(1)
    pushViewport(viewport(xscale = xlim, yscale = ylim, clip = TRUE))
    grid.rect(gp = gpar(fill = scatter.pars$frame.fill, col = scatter.pars$border.col))
    x <- sort(x)
    y <- sort(y)
    grid.lines(unit(x, "native"), unit(y, "native"))
    popViewport(1)
  }
  scatterplot.panel <- function(x, y, type, scatter.pars, axis.pars, 
                                xpos, ypos, xylim) {
    if (is.null(xylim)) {
      xlim <- range(x, na.rm = TRUE) + c(-buffer * (max(x, 
                                                        na.rm = TRUE) - min(x, na.rm = TRUE)), buffer * 
                                           (max(x, na.rm = TRUE) - min(x, na.rm = TRUE)))
      ylim <- range(y, na.rm = TRUE) + c(-buffer * (max(y, 
                                                        na.rm = TRUE) - min(y, na.rm = TRUE)), buffer * 
                                           (max(y, na.rm = TRUE) - min(y, na.rm = TRUE)))
    }    else {
      xlim <- xylim
      ylim <- xylim
    }
    pushViewport(viewport(xscale = xlim, yscale = ylim))
    draw.axis(x, y, axis.pars, xpos, ypos, NULL, NULL, xlim, 
              ylim)
    popViewport(1)
    pushViewport(viewport(xscale = xlim, yscale = ylim, clip = TRUE))
    grid.rect(gp = gpar(fill = scatter.pars$frame.fill, col = scatter.pars$border.col))
    if (scatter.pars$plotpoints & (type == "points" || type == 
                                   "lm" || type == "ci" || type == "symlm" || type == 
                                   "loess")) {
      grid.points(x, y, pch = scatter.pars$pch, size = scatter.pars$size, 
                  gp = gpar(col = scatter.pars$col))
    }
    if (type == "lm") {
      xy.lm <- lm(y ~ x)
      panel.abline(xy.lm$coef[1], xy.lm$coef[2], col = "red", 
                   lwd = 2)
    }
    if (type == "ci") {
      xy.lm <- lm(y ~ x)
      xy <- data.frame(x = seq(min(x, na.rm = TRUE), max(x, 
                                                         na.rm = TRUE), length.out = 20))
      yhat <- predict(xy.lm, newdata = xy, interval = "confidence")
      ci <- data.frame(lower = yhat[, "lwr"], upper = yhat[, 
                                                           "upr"])
      grid.lines(x = c(xy$x), y = c(ci$lower), default.units = "native")
      grid.lines(x = c(xy$x), y = c(ci$upper), default.units = "native")
      grid.polygon(x = c(xy$x, xy$x[length(xy$x):1]), y = c(ci$lower, 
                                                            ci$upper[length(ci$upper):1]), gp = gpar(fill = "grey"), 
                   default.units = "native")
    }
    if (type == "loess") {
      junk <- try(panel.loess(x, y, color = "red", span = 1))
      if (class(junk) == "try-error") 
        warning("An error in loess occurred and was ignored; no line was plotted.")
    }
    if (type == "symlm") {
      pcs <- try(prcomp(cbind(x, y)))
      if (class(pcs) == "try-error") 
        warning("An error in symlm occurred and was ignored; no line was plotted.")
      else {
        slope <- abs(pcs$rotation[1, 2]/pcs$rotation[1, 
                                                     1])
        if (cor(x, y) < 0) 
          slope <- -1 * slope
        panel.abline(pcs$center[2] - slope * pcs$center[1], 
                     slope, col = "blue")
      }
    }
    if (type == "corrgram") {
      pear.test <- cor.test(x, y, method = "pearson", alternative = "two.sided")
      corr <- format(pear.test$estimate, digits = 2)
      if (as.numeric(corr) > 0) {
        panel.fill(col = hsv(h = 0.5, s = abs(as.numeric(corr)), 
                             v = 1), border = hsv(h = 0.5, s = abs(as.numeric(corr)), 
                                                  v = 1))
        grid.lines(x = unit(c(0, 1), "npc"), y = unit(c(0, 
                                                        1), "npc"), gp = gpar(col = "white", lwd = 2))
      } else {
        panel.fill(col = hsv(h = 0, s = abs(as.numeric(corr)), 
                             v = 1), border = hsv(h = 0, s = abs(as.numeric(corr)), 
                                                  v = 1))
        grid.lines(x = unit(c(0, 1), "npc"), y = unit(c(1, 
                                                        0), "npc"), gp = gpar(col = "white", lwd = 2))
      }
    }
    if (type == "qqplot") {
      qq.panel(x, y, scatter.pars, axis.pars, xpos, ypos, 
               xlim, ylim)
    }
    if (type == "stats") {
      complete.obs <- nrow(na.omit(cbind(x, y)))
      missing <- length(x) - complete.obs
      pear.test <- cor.test(x, y, method = "pearson", alternative = "two.sided")
      corr <- sprintf("%03.2f", pear.test$estimate)
      rho.test <- cor.test(x, y, method = "spearman", alternative = "two.sided")
      tau.test <- cor.test(x, y, method = "kendall", alternative = "two.sided")
      rho <- sprintf("%03.2f", rho.test$estimate)
      tau <- sprintf("%03.2f", tau.test$estimate)
      xy.lm <- lm(y ~ x)
      r2 <- sprintf("%03.2f", summary(xy.lm)$r.squared)
      p <- sprintf("%06.4f", pf(q = as.numeric(summary(xy.lm)$fstatistic)[1], 
                                df1 = as.numeric(summary(lm(xy.lm))$fstatistic)[2], 
                                df2 = as.numeric(summary(lm(xy.lm))$fstatistic)[3], 
                                lower.tail = FALSE))
      bonfp <- stat.pars$signif/(N * (N - 1))/2
      sig <- 1
      sigrho <- NULL
      sigtau <- NULL
      sigcor <- NULL
      sigp <- NULL
      if (pear.test$p.value < bonfp) {
        sig <- sig + 1
        sigcor <- "*"
      }
      if (rho.test$p.value < bonfp) {
        sig <- sig + 1
        sigrho <- "*"
      }
      if (tau.test$p.value < bonfp) {
        sig <- sig + 1
        sigtau <- "*"
      }
      if (as.numeric(p) < bonfp) {
        sig <- sig + 1
        sigp <- "*"
      }
      if (mean(as.numeric(rho), as.numeric(tau), as.numeric(corr)) > 
          0) {
        text.color <- "black"
        if (sig == 1) 
          box.color <- 0.5
        else if (sig > 1 && sig < 5) 
          box.color <- 0.75
        else if (sig == 5) 
          box.color <- 1
      } else if (mean(as.numeric(rho), as.numeric(tau), as.numeric(corr)) < 
                 0) {
        text.color <- "white"
        if (sig == 1) 
          box.color <- 0.5
        else if (sig > 1 && sig < 5) 
          box.color <- 0.25
        else if (sig == 5) 
          box.color <- 0
      }
      if (!stat.pars$use.color) {
        panel.fill(col = grey(box.color), border = grey(box.color))
      } else {
        text.color <- "black"
        if (as.numeric(corr) > 0) {
          panel.fill(col = hsv(h = 0.5, s = abs(as.numeric(corr)), 
                               v = 1), border = hsv(h = 0.5, s = abs(as.numeric(corr)), 
                                                    v = 1))
        } else {
          panel.fill(col = hsv(h = 0, s = abs(as.numeric(corr)), 
                               v = 1), border = hsv(h = 0, s = abs(as.numeric(corr)), 
                                                    v = 1))
        }
      }
      if (!is.na(stat.pars$verbose)) {
        if (stat.pars$verbose == TRUE) {
          # grid.text(bquote(rho == .(rho) * .(sigrho)), 
          #          x = 0.5, y = 0.9, just = stat.pars$just, 
          #       gp = gpar(fontsize = stat.pars$fontsize, 
          #                col = text.color))
          # grid.text(bquote(tau == .(tau) * .(sigtau)), 
          #          x = 0.5, y = 0.7, just = stat.pars$just, 
          #         gp = gpar(fontsize = stat.pars$fontsize, 
          #                  col = text.color))
          grid.text(paste("r=", corr, sigcor, sep = ""), 
                    x = 0.5, y = 0.5, just = stat.pars$just, 
                    gp = gpar(fontsize = stat.pars$fontsize, 
                              col = text.color))
          grid.text(paste("p=", p, sigp, sep = ""), x = 0.5, 
                    y = 0.3, just = stat.pars$just, gp = gpar(fontsize = stat.pars$fontsize, 
                                                              col = text.color))
          if (missing > 0) 
            grid.text(paste(missing, stat.pars$missing), 
                      x = 0.5, y = 0.1, just = stat.pars$just, 
                      gp = gpar(fontsize = stat.pars$fontsize, 
                                col = "red"))
        } else {
          
          grid.text(paste(corr, sigcor, sep = ""), x = 0.5, 
                    y = 0.7, just = stat.pars$just, gp = gpar(fontsize = stat.pars$fontsize, 
                                                              col = text.color))
          if (missing > 0) 
            grid.text(paste(missing, "missing"), x = 0.5, 
                      y = 0.3, just = stat.pars$just, gp = gpar(fontsize = stat.pars$fontsize, 
                                                                col = text.color))
        }
      }
    }
    popViewport(1)
  }
  mosaic.panel <- function(x, y, mosaic.pars, axis.pars, xpos, 
                           ypos) {
    if (!is.null(xpos) & !is.null(ypos)) {
      strucplot(table(y, x), margins = c(0, 0, 0, 0), newpage = FALSE, 
                pop = FALSE, keep_aspect_ratio = FALSE, shade = mosaic.pars$shade, 
                legend = FALSE, gp = mosaic.pars$gp, gp_args = mosaic.pars$gp_args, 
                labeling_args = list(tl_labels = c(xpos, !ypos), 
                                     gp_labels = mosaic.pars$gp_labels, varnames = c(FALSE, 
                                                                                     FALSE), rot_labels = c(outer.rot, outer.rot)))
    }  else {
      if (is.null(xpos) & is.null(ypos)) {
        strucplot(table(y, x), margins = c(0, 0, 0, 0), 
                  shade = mosaic.pars$shade, legend = FALSE, 
                  gp = mosaic.pars$gp, gp_args = mosaic.pars$gp_args, 
                  newpage = FALSE, pop = FALSE, keep_aspect_ratio = FALSE, 
                  labeling = NULL)
      }  else {
        if (is.null(xpos)) {
          strucplot(table(y, x), margins = c(0, 0, 0, 
                                             0), newpage = FALSE, pop = FALSE, keep_aspect_ratio = FALSE, 
                    shade = mosaic.pars$shade, legend = FALSE, 
                    gp = mosaic.pars$gp, gp_args = mosaic.pars$gp_args, 
                    labeling_args = list(labels = c(TRUE, FALSE), 
                                         tl_labels = c(ypos, FALSE), gp_labels = mosaic.pars$gp_labels, 
                                         varnames = c(FALSE, FALSE), rot_labels = c(outer.rot, 
                                                                                    outer.rot)))
        } else {
          strucplot(table(y, x), margins = c(0, 0, 0, 
                                             0), newpage = FALSE, pop = FALSE, keep_aspect_ratio = FALSE, 
                    shade = mosaic.pars$shade, legend = FALSE, 
                    gp = mosaic.pars$gp, gp_args = mosaic.pars$gp_args, 
                    labeling_args = list(labels = c(FALSE, TRUE), 
                                         tl_labels = c(FALSE, !xpos), gp_labels = mosaic.pars$gp_labels, 
                                         varnames = c(FALSE, FALSE), rot_labels = c(outer.rot, 
                                                                                    outer.rot)))
        }
      }
    }
  }
  boxplot.panel <- function(x, y, type, axis.pars, xpos, ypos, 
                            xylim) {
    xlim <- NULL
    ylim <- NULL
    old.color <- trellis.par.get("box.rectangle")$col
    trellis.par.set(name = "box.rectangle", value = list(col = "black"))
    trellis.par.set(name = "box.umbrella", value = list(col = "black"))
    trellis.par.set(name = "box.dot", value = list(col = "black"))
    trellis.par.set(name = "plot.symbol", value = list(col = "black"))
    if (is.factor(x)) {
      cat.labels <- levels(x)
      k <- length(levels(x))
      cat.var <- as.numeric(x)
      cont.var <- y
      horiz <- FALSE
    } else {
      cat.labels <- levels(y)
      k <- length(levels(y))
      cat.labels <- cat.labels[k:1]
      cat.var <- k + 1 - as.numeric(y)
      cont.var <- x
      horiz <- TRUE
    }
    if (horiz) {
      if (is.null(xylim)) {
        xlim <- range(cont.var, na.rm = TRUE) + c(-buffer * 
                                                    (max(cont.var, na.rm = TRUE) - min(cont.var, 
                                                                                       na.rm = TRUE)), buffer * (max(cont.var, na.rm = TRUE) - 
                                                                                                                   min(cont.var, na.rm = TRUE)))
      } else {
        xlim <- xylim
      }
      pushViewport(viewport(xscale = xlim, yscale = c(0.5, 
                                                      max(cat.var, na.rm = TRUE) + 0.5)))
      if (is.null(ypos)) 
        cat.labels <- NULL
      draw.axis(cont.var, cat.var, axis.pars, xpos, ypos, 
                cat.labels, horiz, xlim, ylim)
      popViewport(1)
      pushViewport(viewport(xscale = xlim, yscale = c(0.5, 
                                                      max(cat.var, na.rm = TRUE) + 0.5), clip = TRUE))
      if (type == "boxplot") 
        panel.bwplot(cont.var, cat.var, horizontal = horiz, 
                     col = "black", pch = "|", gp = gpar(box.umbrella = list(col = "black")))
      if (type == "stripplot") 
        panel.stripplot(cont.var, cat.var, horizontal = horiz, 
                        jitter.data = stripplot.pars$jitter, col = stripplot.pars$col, 
                        cex = stripplot.pars$size, pch = stripplot.pars$pch)
    }else {
      if (is.null(xylim)) {
        ylim <- range(cont.var, na.rm = TRUE) + c(-buffer * 
                                                    (max(cont.var, na.rm = TRUE) - min(cont.var, 
                                                                                       na.rm = TRUE)), buffer * (max(cont.var, na.rm = TRUE) - 
                                                                                                                   min(cont.var, na.rm = TRUE)))
      }else {
        ylim <- xylim
      }
      pushViewport(viewport(yscale = ylim, xscale = c(0.5, 
                                                      max(cat.var, na.rm = TRUE) + 0.5)))
      if (is.null(xpos)) 
        cat.labels <- NULL
      draw.axis(cat.var, cont.var, axis.pars, xpos, ypos, 
                cat.labels, horiz, xlim, ylim)
      popViewport(1)
      pushViewport(viewport(yscale = ylim, xscale = c(0.5, 
                                                      max(cat.var, na.rm = TRUE) + 0.5), clip = TRUE))
      if (type == "boxplot") 
        panel.bwplot(cat.var, cont.var, horizontal = horiz, 
                     col = "black", pch = "|", gp = gpar(box.umbrella = list(col = "black")))
      if (type == "stripplot") 
        panel.stripplot(cat.var, cont.var, horizontal = horiz, 
                        jitter.data = stripplot.pars$jitter, col = stripplot.pars$col, 
                        cex = stripplot.pars$size, pch = stripplot.pars$pch)
    }
    grid.rect(gp = gpar(fill = NULL))
    popViewport(1)
    trellis.par.set(name = "box.rectangle", value = list(col = old.color))
    trellis.par.set(name = "box.umbrella", value = list(col = old.color))
    trellis.par.set(name = "box.dot", value = list(col = old.color))
    trellis.par.set(name = "plot.symbol", value = list(col = old.color))
  }
  diag.panel <- function(x, varname, diag.pars, axis.pars, 
                         xpos, ypos, xylim) {
    x <- x[!is.na(x)]
    if (is.null(xylim)) {
      xlim <- range(as.numeric(x), na.rm = TRUE) + c(-buffer * 
                                                       (max(as.numeric(x), na.rm = TRUE) - min(as.numeric(x), 
                                                                                               na.rm = TRUE)), buffer * (max(as.numeric(x), 
                                                                                                                             na.rm = TRUE) - min(as.numeric(x), na.rm = TRUE)))
    }else {
      xlim <- xylim
    }
    ylim <- xlim
    pushViewport(viewport(xscale = xlim, yscale = ylim))
    draw.axis(as.numeric(x), as.numeric(x), axis.pars, xpos, 
              ypos, NULL, NULL, xlim, ylim)
    popViewport(1)
    pushViewport(viewport(xscale = xlim, yscale = ylim, clip = TRUE))
    if (!diag.pars$show.hist) {
      grid.rect()
      grid.text(varname, 0.5, 0.5, gp = gpar(fontsize = diag.pars$fontsize, 
                                             fontface = 2))
    }
    popViewport(1)
    if (diag.pars$show.hist) {
      if (!is.factor(x)) {
        pushViewport(viewport(xscale = xlim, yscale = c(0, 
                                                        100), clip = TRUE))
        panel.histogram(as.numeric(x), breaks = NULL, 
                        type = "percent", col = diag.pars$hist.color)
      }else {
        pushViewport(viewport(xscale = c(min(as.numeric(x), 
                                             na.rm = TRUE) - 1, max(as.numeric(x), na.rm = TRUE) + 
                                           1), yscale = c(0, 100), clip = TRUE))
        panel.barchart(1:length(table(x)), 100 * table(x)/sum(table(x)), 
                       horizontal = FALSE, col = diag.pars$hist.color)
      }
      grid.text(varname, 0.5, 0.85, gp = gpar(fontsize = diag.pars$fontsize))
      popViewport(1)
    }
  }
  grid.newpage()
  N <- ncol(x)
  vp.main <- viewport(x = outer.margins$bottom, y = outer.margins$left, 
                      width = unit(1, "npc") - outer.margins$right - outer.margins$left, 
                      height = unit(1, "npc") - outer.margins$top - outer.margins$bottom, 
                      just = c("left", "bottom"), name = "main", clip = "off")
  pushViewport(vp.main)
  for (i in 1:N) {
    for (j in 1:N) {
      if (diagonal == "default") 
        labelj <- j
      else labelj <- N - j + 1
      x[is.infinite(x[, i]), i] <- NA
      x[is.infinite(x[, j]), j] <- NA
      vp <- viewport(x = (labelj - 1)/N, y = 1 - i/N, width = 1/N, 
                     height = 1/N, just = c("left", "bottom"), name = as.character(i * 
                                                                                     N + j))
      pushViewport(vp)
      vp.in <- viewport(x = 0.5, y = 0.5, width = 1 - gap, 
                        height = 1 - gap, just = c("center", "center"), 
                        name = paste("IN", as.character(i * N + j)))
      pushViewport(vp.in)
      xpos <- NULL
      if (i == 1 && outer.labels$top[j]) {
        xpos <- FALSE
      }
      if (i == N && outer.labels$bottom[j]) {
        xpos <- TRUE
      }
      ypos <- NULL
      if (j == N && outer.labels$right[i]) {
        ypos <- FALSE
      }
      if (j == 1 && outer.labels$left[i]) {
        ypos <- TRUE
      }
      if (!is.null(ypos) & diagonal != "default") {
        ypos <- !ypos
      }
      if (i == j) {
        diag.panel(x[, i], names(x)[i], diag.pars, axis.pars, 
                   xpos, ypos, xylim)
      }else {
        if (is.factor(x[, i]) + is.factor(x[, j]) == 
            1) {
          if (i < j & upper.pars$conditional != "barcode") 
            boxplot.panel(x[, j], x[, i], upper.pars$conditional, 
                          axis.pars, xpos, ypos, xylim)
          if (i > j & lower.pars$conditional != "barcode") 
            boxplot.panel(x[, j], x[, i], lower.pars$conditional, 
                          axis.pars, xpos, ypos, xylim)
          if (i < j & upper.pars$conditional == "barcode") {
            if (is.factor(x[, i])) {
              barcode(split(x[, j], x[, i])[length(levels(x[, 
                                                            i])):1], horizontal = TRUE, xlim = xylim, 
                      labelloc = ypos, axisloc = xpos, labelouter = TRUE, 
                      newpage = FALSE, fontsize = axis.pars$fontsize, 
                      buffer = buffer, nint = barcode.pars$nint, 
                      ptsize = barcode.pars$ptsize, ptpch = barcode.pars$ptpch, 
                      bcspace = barcode.pars$bcspace, use.points = barcode.pars$use.points)
            } else {
              if (!is.null(ypos)) 
                ypos <- !ypos
              barcode(split(x[, i], x[, j])[length(levels(x[, 
                                                            j])):1], horizontal = FALSE, xlim = xylim, 
                      labelloc = xpos, axisloc = ypos, labelouter = TRUE, 
                      newpage = FALSE, fontsize = axis.pars$fontsize, 
                      buffer = buffer, nint = barcode.pars$nint, 
                      ptsize = barcode.pars$ptsize, ptpch = barcode.pars$ptpch, 
                      bcspace = barcode.pars$bcspace, use.points = barcode.pars$use.points)
            }
          }
          if (i > j & lower.pars$conditional == "barcode") {
            if (is.factor(x[, i])) {
              barcode(split(x[, j], x[, i])[length(levels(x[, 
                                                            i])):1], horizontal = TRUE, xlim = xylim, 
                      labelloc = ypos, axisloc = xpos, labelouter = TRUE, 
                      newpage = FALSE, fontsize = axis.pars$fontsize, 
                      buffer = buffer, nint = barcode.pars$nint, 
                      ptsize = barcode.pars$ptsize, ptpch = barcode.pars$ptpch, 
                      bcspace = barcode.pars$bcspace, use.points = barcode.pars$use.points)
            } else {
              if (!is.null(ypos)) 
                ypos <- !ypos
              barcode(split(x[, i], x[, j])[length(levels(x[, 
                                                            j])):1], horizontal = FALSE, xlim = xylim, 
                      labelloc = xpos, axisloc = ypos, labelouter = TRUE, 
                      newpage = FALSE, fontsize = axis.pars$fontsize, 
                      buffer = buffer, nint = barcode.pars$nint, 
                      ptsize = barcode.pars$ptsize, ptpch = barcode.pars$ptpch, 
                      bcspace = barcode.pars$bcspace, use.points = barcode.pars$use.points)
            }
          }
        }
        if (is.factor(x[, i]) + is.factor(x[, j]) == 
            0) {
          if (i < j) 
            type <- upper.pars$scatter
          else type <- lower.pars$scatter
          scatterplot.panel(x[, j], x[, i], type, scatter.pars, 
                            axis.pars, xpos, ypos, xylim)
        }
        if (is.factor(x[, i]) + is.factor(x[, j]) == 
            2) {
          if (i < j) 
            mosaic.panel(x[, j], x[, i], mosaic.pars, 
                         axis.pars, xpos, ypos)
          else mosaic.panel(x[, j], x[, i], mosaic.pars, 
                            axis.pars, xpos, ypos)
        }
      }
      popViewport(1)
      upViewport()
    }
  }
  popViewport()
  if (whatis) 
    whatis(x)
}

gpar <- function(...) {
  gp <- validGP(list(...))
  class(gp) <- "gpar"
  gp
}

is.gpar <- function(x) {
  inherits(x, "gpar")
}

print.gpar <- function(x, ...) {
  print(unclass(x), ...)
  invisible(x)
}

validGP <- function(gpars) {
  # Check a (non-NULL) gpar is not of length 0
  check.length <- function(gparname) {
    if (length(gpars[[gparname]]) == 0)
      stop(gettextf("'gpar' element '%s' must not be length 0", gparname),
           domain = NA)
  }
  # Check a gpar is numeric and not NULL
  numnotnull <- function(gparname) {
    if (!is.na(match(gparname, names(gpars)))) {
      if (is.null(gpars[[gparname]]))
        gpars[[gparname]] <<- NULL
      else {
        check.length(gparname)
        gpars[[gparname]] <<- as.numeric(gpars[[gparname]])
      }
    }
  }
  # fontsize, lineheight, cex, lwd should be numeric and not NULL
  numnotnull("fontsize")
  numnotnull("lineheight")
  numnotnull("cex")
  numnotnull("lwd")
  numnotnull("lex")
  # gamma defunct in 2.7.0
  if ("gamma" %in% names(gpars)) {
    warning("'gamma' 'gpar' element is defunct")
    gpars$gamma <- NULL
  }
  numnotnull("alpha")
  # col and fill are converted in C code
  # BUT still want to check length > 0
  if (!is.na(match("col", names(gpars)))) {
    if (is.null(gpars$col))
      gpars$col <- NULL
    else
      check.length("col")
  }
  if (!is.na(match("fill", names(gpars)))) {
    if (is.null(gpars$fill))
      gpars$fill <- NULL
    else
      check.length("fill")
  }
  # lty converted in C code
  # BUT still want to check for NULL and check length > 0
  if (!is.na(match("lty", names(gpars)))) {
    if (is.null(gpars$lty))
      gpars$lty <- NULL
    else
      check.length("lty")
  }
  if (!is.na(match("lineend", names(gpars)))) {
    if (is.null(gpars$lineend))
      gpars$lineend <- NULL
    else
      check.length("lineend")
  }
  if (!is.na(match("linejoin", names(gpars)))) {
    if (is.null(gpars$linejoin))
      gpars$linejoin <- NULL
    else
      check.length("linejoin")
  }
  # linemitre should be larger than 1
  numnotnull("linemitre")
  if (!is.na(match("linemitre", names(gpars)))) {
    if (any(gpars$linemitre < 1))
      stop("invalid 'linemitre' value")
  }
  # alpha should be 0 to 1
  if (!is.na(match("alpha", names(gpars)))) {
    if (any(gpars$alpha < 0 || gpars$alpha > 1))
      stop("invalid 'alpha' value")
  }
  # font should be integer and not NULL
  if (!is.na(match("font", names(gpars)))) {
    if (is.null(gpars$font))
      gpars$font <- NULL
    else {
      check.length("font")
      gpars$font <- as.integer(gpars$font)
    }
  }
  # fontfamily should be character
  if (!is.na(match("fontfamily", names(gpars)))) {
    if (is.null(gpars$fontfamily))
      gpars$fontfamily <- NULL
    else {
      check.length("fontfamily")
      gpars$fontfamily <- as.character(gpars$fontfamily)
    }
  }
  # fontface can be character or integer;  map character to integer
  # store value in font
  # Illegal to specify both font and fontface
  if (!is.na(match("fontface", names(gpars)))) {
    if (!is.na(match("font", names(gpars))))
      stop("must specify only one of 'font' and 'fontface'")
    gpars$font <-
      if (is.null(gpars$fontface)) NULL # remove it
    else {
      check.length("fontface")
      if (is.numeric(gpars$fontface))
        as.integer(gpars$fontface)
      else
        vapply(as.character(gpars$fontface),
               function(ch) # returns integer
                 switch(ch,
                        plain = 1L,
                        bold  = 2L,
                        italic=, oblique = 3L,
                        bold.italic = 4L,
                        symbol= 5L,
                        # These are Hershey variants
                        cyrillic=5L,
                        cyrillic.oblique=6L,
                        EUC   = 7L,
                        stop("invalid fontface ", ch)), 0L)
    }
  }
  gpars
}

# Method for subsetting "gpar" objects
`[.gpar` <- function(x, index, ...) {
  if (length(x) == 0)
    return(gpar())
  maxn <- do.call("max", lapply(x, length))
  newgp <- lapply(x, rep, length.out=maxn)
  newgp <- lapply(X = newgp, FUN = "[", index, ...)
  class(newgp) <- "gpar"
  newgp
}

# possible gpar names
# The order must match the GP_* values in grid.h
.grid.gpar.names <- c("fill", "col", "gamma", "lty", "lwd", "cex",
                      "fontsize", "lineheight", "font", "fontfamily",
                      "alpha", "lineend", "linejoin", "linemitre",
                      "lex",
                      # Keep fontface at the end because it is never
                      # used in C code (it gets mapped to font)
                      "fontface")

set.gpar <- function(gp) {
  if (!is.gpar(gp))
    stop("argument must be a 'gpar' object")
  temp <- grid.Call(L_getGPar)
  # gamma defunct in 2.7.0
  if ("gamma" %in% names(gp)) {
    warning("'gamma' 'gpar' element is defunct")
    gp$gamma <- NULL
  }
  # Special case "cex" (make it cumulative)
  if (match("cex", names(gp), nomatch=0L))
    tempcex <- temp$cex * gp$cex
  else
    tempcex <- temp$cex
  # Special case "alpha" (make it cumulative)
  if (match("alpha", names(gp), nomatch=0L))
    tempalpha <- temp$alpha * gp$alpha
  else
    tempalpha <- temp$alpha
  # Special case "lex" (make it cumulative)
  if (match("lex", names(gp), nomatch=0L))
    templex <- temp$lex * gp$lex
  else
    templex <- temp$lex
  # All other gpars
  temp[names(gp)] <- gp
  temp$cex <- tempcex
  temp$alpha <- tempalpha
  temp$lex <- templex
  # Do this as a .Call.graphics to get it onto the base display list
  grid.Call.graphics(L_setGPar, temp)
}

get.gpar <- function(names=NULL) {
  if (is.null(names)) {
    result <- grid.Call(L_getGPar)
    # drop gamma
    result$gamma <- NULL
  } else {
    if (!is.character(names) ||
        !all(names %in% .grid.gpar.names))
      stop("must specify only valid 'gpar' names")
    # gamma deprecated
    if ("gamma" %in% names) {
      warning("'gamma' 'gpar' element is defunct")
      names <- names[-match("gamma", names)]
    }
    result <- unclass(grid.Call(L_getGPar))[names]
  }
  class(result) <- "gpar"
  result
}

# When editing a gp slot, only update the specified gpars
# Assume gp is NULL or a gpar
# assume newgp is a gpar (and not NULL)
mod.gpar <- function(gp, newgp) {
  if (is.null(gp))
    gp <- newgp
  else
    gp[names(newgp)] <- newgp
  gp
}
}


#' template for Covariate plots
#'
#' @param data Data frame. Use the codes below as template 
#' @keywords cov_plots()
#' @export
#' @examples
cov_plots<-function(...){
  library(reshape2)
  library(plyr)
  require(PCSmisc)
  require(foreign)
  require(Hmisc)
  library(sas7bdat)
  library(lhtool)
  library(lhwordtool)
  library(gridExtra)
  library(ggplot2)
  library(gridExtra)
  require(reshape)
  library(ggpubr)
  cov_fn()
  library(dplyr)
  dir<-"//certara.com/sites/S02-Cary/Consulting/Projects/projects/"
  cov1<-read.csv(file.path(dir,"dataset.csv"))
  cov1<-cov1[cov1$mdv1==0,]

  cov<-cov1[!duplicated(cov1$subject),]

  r<-read.csv("Eta.csv")

  head(r)
  unique(r$Name)
  names(r)<-tolower(names(r))
  
  head(cov)
  #ID for for merging
  id<-"id"
  #Keep ETA
  keta<-c("ndur","ntlag","nv","ncl")
  #Keep Categorical COV
  head(r)
  #cov<-cov[cov$id%in%r$id,]
  
  cca<-c("study","sex","agecat","country","bmic", "ethnic","egfrcatn","phasen","oster","health","hepcat","corti","race","ethnic1","patient")
  
  #KEEP Cont COV
  cco<-c("wt","bsa","age","bmi" ,"alb","alt","ast","bil","dose","uacr","bvas","egfr")
  
  cov<-chclass(cov,cca,"char")
  
  #ORDER BASELINE CATEGORICAL
  one(cov,"race")
  cov$race<-as.character(cov$race)
  dat1<-reflag(cov,"sex",c("Female", "Male"),c("Female","Male"))
  dat1<-reflag(dat1,"agecat",c(">=18 & <50", ">=50 & <65",">=65 & <75", ">=75"))
  dat1<-reflag(dat1,"bmic",c("0","1","Missing"),c("<30 kg/m^2", ">=30 kg/m^2","Missing"))
  dat1<-reflag(dat1,"ethnic1",c("Hispanic or Latino","Not Hispanic or Latino","Unknown","Not reported"),c("Hispanic or Latino","Not Hispanic or Latino","Unknown","Unknown"))
  
  dat1<-reflag(dat1,"race",c(" White","White"," Black Or African American","Black"," Asian","Asian"," Other"," White_ Native Hwaiian/Pacific Islander"),c("White","White","Black","Black","Asian","Asian","Other","Other"))
  
  dat1<-reflag(dat1,"patient",c("other", "AAV","IgAN"))
  
  dat1<-reflag(dat1,"hepcat",c("Healthy","Mild","Moderate","Others"))
  
  dat1<-reflag(dat1,"ethnic",c( "non-Japanese","Japanese", "Unknown"))
  dat1<-reflag(dat1,"health",c(0,3,1,2),c("Other","Other","AAV","IgAN"))
  dat1<-reflag(dat1,"egfrcatn",c("0","1", "2", "3" ),c("Normal","Mild", "Moderate", "Severe" ))
  dat1<-reflag(dat1,"phasen",c("1","2" ),
               c("Phase 1","Phase 2"))
  dat1<-reflag(dat1,"oster",c("0","1"),c("No","Yes"))
  dat1<-reflag(dat1,"corti",c("No","No with AAV","Yes"))
  dat1<-chclass(dat1,cco,"num")
  
  dat2<-lhcut(dat1[!is.na(dat1$bmi),],"bmi",c(20,30))
  dat3<-dat1[is.na(dat1$bmi),];dat3$bmicat<-"Missing"
  dat1<-rbind(dat2,dat3)
  

  #GROUP CONTINUOUS COVARIATES
  ###########RUN THE FOLLOWING LINES AS IS FOR MOST OF THE TIMES#############################
  cov<-chclass(cov,cco,"num")
  names(r)<-tolower(names(r))
  
  
  eta<-r[,c(id,keta)]
  cateta<-join(eta,dat1[,c(id,cca)],type="left")
  nrow(eta);nrow(cateta)
  head(cateta)
  
  
  cat1<-lhlong(cateta,names(cateta[,(length(keta)+2):ncol(cateta)]))
  
  head(cat1)
  
  names(cat1)[names(cat1)=="variable"]<-"Covariate"
  names(cat1)[names(cat1)=="value"]<-"Categorical"
  cat1<-chclass(cat1,c("Covariate","Categorical"),"char")
  cat1<-addvar(cat1,c("Covariate","Categorical"),keta[1],"length(x)","yes","count")
  cat1$Cat1<-paste0(cat1$Categorical,"\n (n=",cat1$count,")")
  
  cat1<-lhlong(cat1,keta)
  head(cat1)
  cat1<-chclass(cat1,c("Covariate","Categorical","variable"),"char")
  unique(cat1$Categorical)
  head(cat1)
  cat1$variable<-factor(cat1$variable,levels=keta)
  
  catnum<-addvar(nodup(cat1,c("Covariate","Categorical"),"var"),"Covariate","Categorical","length(x)","no","catnumber")
  
  for(i in cca[cca%in%catnum$Covariate[catnum$catnumber>0]]){
    head(cateta)
    dcat<-cat1[cat1$Covariate%in%i,]
    ord<-sort(unique(cateta[,i]))
    label<-nodup(dcat,c("Categorical","Cat1"),"var")
    label$Categorical<-factor(label$Categorical,levels=ord)
    lablel<-label$Cat1[order(label$Categorical)]
    dcat$Cat1<-factor(dcat$Cat1,levels=lablel)

    p<-ggplot(dcat,aes(x=Cat1,y=value))+
      geom_boxplot(outlier.shape = NA)+
      geom_jitter(position=position_jitter(0.1),col="grey")+
      geom_hline(yintercept=0,linetype=2, color="red",size=1)+
      ylab("Individual Random Effect")+xlab("")+
      facet_wrap(~variable,scale="free")+
      theme_bw()+
      theme(axis.text.x = element_text(angle =45, hjust =0.5, vjust = 0.5))
    ggsave(paste0("1_",i,"_boxplot.png"),p,width=5.7,height=5.7)
  }
  
  
  ########GPAIRS##########
  coneta<-join(eta,cov[,c("id",cco)],type="left")

  names(coneta)<-tolower(names(coneta))
  cco<-tolower(cco)
  cco<-cco
  coneta<-chclass(coneta,cco,"num")
  
  congr<-NULL
  congr[[1]]<-c("wt","age","bsa","bmi")
  congr[[2]]<-c("alb","alt","ast","bil","dose")
  congr[[3]]<-c("uacr","bvas","egfr")

  for(i in 1:length(congr)){
    png(file=paste0("scatter",i,".png"),width=768,height=768,pointsize = 16)
    print(gpairs(x=coneta[,c(keta,congr[[i]])], 
                 upper.pars = list(conditional = 'boxplot', scatter = 'loess'),
                 lower.pars = list(scatter = 'stats',conditional = "barcode"),
                 diag.pars = list(fontsize = 9, show.hist = TRUE, hist.color = "gray"),
                 stat.pars =list(fontsize = 11, signif =F, verbose =T, use.color = TRUE, missing = 'missing', just = 'centre'),
                 scatter.pars = list(pch = 20)))
    dev.off()
  }

 
  # Distribution of ETA 
  p1<-ggplot(r,aes(x=ntlag))+  geom_density(fill="royalblue3",col=NA,alpha=0.3)+
    geom_histogram(aes(x=ntlag,y=..density..),fill=NA,col="black")+
    geom_vline(xintercept=0,col="red",size=1.5,linetype = "dashed")
  
  p2<-ggplot(r,aes(x=ndur))+  geom_density(fill="royalblue3",col=NA,alpha=0.3,bin=60)+
    geom_histogram(aes(x=ndur,y=..density..),fill=NA,col="black")+
    geom_vline(xintercept=0,col="red",size=1.5,linetype = "dashed")
  
  p3<-ggplot(r,aes(x=ncl))+  geom_density(fill="royalblue3",col=NA,alpha=0.3,bin=60)+
    geom_histogram(aes(x=ncl,y=..density..),fill=NA,col="black")+
    geom_vline(xintercept=0,col="red",size=1.5,linetype = "dashed")
  
  p4<-ggplot(r,aes(x=nv))+  geom_density(fill="royalblue3",col=NA,alpha=0.3,bin=60)+
    geom_histogram(aes(x=nv,y=..density..),fill=NA,col="black")+
    geom_vline(xintercept=0,col="red",size=1.5,linetype = "dashed")
  
  library(gridExtra)
  ggsave("ETA_CWRES_distribution.png",
         ggarrange(p1,p2,p3,p4),dpi = 300, width =12, height = 8,units = c("in"))  

}
  
  
#' template for GOF
#'
#' @param data Data frame. Use the codes below as template 
#' @keywords gof_plots()
#' @export
#' @examples 

gof_plots<-function(...){
  rm(list=ls())
  library(reshape2)
  library(plyr)
  require(PCSmisc)
  require(foreign)
  require(Hmisc)
  library(sas7bdat)
  library(lhtool)
  library(ggplot2)
  library(lhwordtool)
  library(gridExtra)
  library(ggpubr)
    ###########
  #install.packages("qqplotr")
  ###########
  # Functions
  qqplot.cwres <- function(dat, ...) {
    ylim <-c(-10, 10)
    xlim <-c(-4, 4)
    xlab <- "Quantiles of Standard Normal"
    ylab <- "Conditional Weighted Residuals"
    with(dat, qqnorm(CWRES, ylim=ylim, xlim=xlim, xlab=xlab, ylab=ylab, ...))
    with(dat, qqline(CWRES))
    abline(a=0,b=1,col="red") #add the standard normal qqplot
  }
  
  histogram.cwres <- function(dat, ...) {
    ylim <-c(0, 0.55)
    xlab <- "Conditional Weighted Residuals"
    with(dat, hist(CWRES, ylim=ylim, xlab=xlab, main="", freq=FALSE, ...))
    abline(v=0, lty=2, lwd=3, col="gray")
    xs <- seq(-10, 140, len=100)
    lines(xs, dnorm(xs), col="gray", lwd=3)
  }
  #------------------------------------------------------------------------------

  dir<-"//certara.com/sites/S02-Cary/Consulting/Projects/projects/"
  dat<-read.csv(file.path(dir,"dataset.csv"))
  cov<-dat[!duplicated(dat$subject),]#read.csv(file.path(path,"Individual of Covariates_ximab.csv"))  
  r<-read.csv("Residuals.csv")
  dat<-dat[dat$mdv==0,]
  unique(r$IVAR-dat$rtime)
  nrow(dat)
  nrow(r)
  
  
  r$cwres<-r$CWRES
  # CWRES distribution
  png(file = "dist_QQ_CWRES.png", width = 8, height = 6, units = 'in', res = 300)
  par(mfrow = c(1, 2))
  qqplot.cwres(r)
  histogram.cwres(r)
  dev.off()
  
  IPRED<-"Individual Predicted Concentration (ng/mL)"
  PRED<-"Population Predicted Concentration (ng/mL)"
  DV<-"Observed Concentration (ng/mL)"
  TAD<-"Time After Dose (h)"
  RTIME<-"Time After First Dose (h)"
  conc<-"Concentration (ng/mL)"
  IVAR<-"Time After First Dose (h)"
  CWRES<-"Conditional Weighted Residuals"
  
  r$IPRED<-exp(r$IPRED);r$DV<-exp(r$DV)
  r$PRED<-exp(r$PRED)
  
  dem<-nodup(dat,"subject","all")

  r1<-join(r,dem[,c("id","study","subject")],type="left")
  
  outl<-r1[abs(r1$CWRES)>4,c("study","subject","IVAR","TAD","DV","PRED","IPRED","CWRES")]
  outl<-outl[order(outl$TAD),]
  
  exp<-dat[dat$subject%in%outl$subject[outl$PRED>1],c("subject","rtime","TAD","dv1")]
  nrow(outl)/nrow(dat)*100
  
  exp<-left_join(exp,lhmutate(outl[outl$PRED>1,c("subject","IVAR","DV","CWRES")],c("DV=outlier","IVAR=rtime")))
  write.csv(outl,"Outliers.csv")
  

  limx<-limy<-range(c(r1$IPRED,r1$DV))
  limx2<-limy2<-range(c(r1$PRED,r1$DV))
  limx1<-limy1<-c(0.001,10000)
  
  r1$STUDYID<-as.character(r1$study)

  
  library(RColorBrewer)
  cols<-brewer.pal(n =6, name = 'Paired')
  cols<-c("#A6CEE3", "#1F78B4", "#B2DF8A", "#33A02C", "#FB9A99", "#E31A1C", "black", "#FF7F00")
  #cols<-c("grey","green","red","black","blue")
  
  p1<-ggplot(r1,aes(x=IPRED,y=DV))+
    geom_point(aes(col=STUDYID))+
    scale_x_continuous(limits=limy,labels = function(x) format(x, scientific =F))+
    scale_y_continuous(limits=limx,labels = function(x) format(x, scientific =F))+
    xlab(IPRED)+ylab(DV)+
    geom_abline(slope=1,size=1)+
    geom_hline(aes(yintercept=0,col="Identity"))+
    geom_smooth(aes(col="LOESS"),se=F,span=1,size=1)+
    
    scale_colour_manual(name="",values=cols, 
                        guide = "legend") + 
    theme_bw()+
    guides(colour=guide_legend(override.aes=list(linetype=c(rep(0,6),1,1),shape=c(rep(16,6),NA,NA))))+
    theme(axis.text.y = element_text(angle = 90, hjust =0.5, vjust = 0.5),legend.position=c(0.8,0.2),legend.box.margin=NULL, legend.background = element_rect(fill="transparent"))
  
  
  p2<-ggplot(r1,aes(x=IPRED,y=DV))+
    #geom_point(col="grey")+
    geom_point(aes(col=factor(STUDYID)))+
    scale_y_log10(limits =limx,breaks = scales::trans_breaks("log10", function(x) 10^x),
                  labels = scales::trans_format("log10", scales::math_format(10^.x)))+
    scale_x_log10(limits =limx , breaks = scales::trans_breaks("log10", function(x) 10^x),
                  labels = scales::trans_format("log10", scales::math_format(10^.x)))+
    annotation_logticks(sides = "lb") +
    geom_abline(slope=1)+
    xlab(IPRED)+ylab(DV)+
    geom_abline(slope=1,size=1)+
    geom_hline(aes(yintercept=0,col="Identity"))+
    geom_smooth(aes(col="LOESS"),se=F,span=1,size=1)+
    
    scale_colour_manual(name="",values=cols, 
                        guide = "legend") + 
    theme_bw()+
    guides(colour=guide_legend(override.aes=list(linetype=c(rep(0,6),1,1),shape=c(rep(16,6),NA,NA))))+
    theme(axis.text.y = element_text(angle = 90, hjust =0.5, vjust = 0.5),legend.position=c(0.8,0.2),legend.box.margin=NULL, legend.background = element_rect(fill="transparent"))
  
  p1a<-ggplot(r1,aes(x=PRED,y=DV))+
    #geom_point(col="grey")+
    geom_point(aes(col=factor(STUDYID)))+
    scale_x_continuous(limits=limy2,labels = function(x) format(x, scientific =F))+
    scale_y_continuous(limits=limx2,labels = function(x) format(x, scientific =F))+
    xlab(PRED)+ylab(DV)+
    geom_abline(slope=1,size=1)+
    geom_hline(aes(yintercept=0,col="Identity"))+
    geom_smooth(aes(col="LOESS"),se=F,span=1,size=1)+
    
    scale_colour_manual(name="",values=cols, 
                        guide = "legend") + 
    theme_bw()+
    guides(colour=guide_legend(override.aes=list(linetype=c(rep(0,6),1,1),shape=c(rep(16,6),NA,NA))))+
    theme(axis.text.y = element_text(angle = 90, hjust =0.5, vjust = 0.5),legend.position=c(0.8,0.2),legend.box.margin=NULL, legend.background = element_rect(fill="transparent"))
  
  
  p2a<-ggplot(r1,aes(x=PRED,y=DV))+
    #geom_point(ae(col="grey")+
    geom_point(aes(col=factor(STUDYID)))+
    scale_y_log10(limits =limx2,breaks = scales::trans_breaks("log10", function(x) 10^x),
                  labels = scales::trans_format("log10", scales::math_format(10^.x)))+
    scale_x_log10(limits =limx2 , breaks = scales::trans_breaks("log10", function(x) 10^x),
                  labels = scales::trans_format("log10", scales::math_format(10^.x)))+
    annotation_logticks(sides = "lb") +
    geom_abline(slope=1)+
    xlab(PRED)+ylab(DV)+
    geom_abline(slope=1,size=1)+
    geom_hline(aes(yintercept=0,col="Identity"))+
    geom_smooth(aes(col="LOESS"),se=F,span=1,size=1)+
    
    scale_colour_manual(name="",values=cols, 
                        guide = "legend") + 
    theme_bw()+
    guides(colour=guide_legend(override.aes=list(linetype=c(rep(0,6),1,1),shape=c(rep(16,6),NA,NA))))+
    theme(axis.text.y = element_text(angle = 90, hjust =0.5, vjust = 0.5),legend.position=c(0.8,0.2),legend.box.margin=NULL, legend.background = element_rect(fill="transparent"))
  
  ggsave("GOF.png",
         ggarrange(p1,p1a,p2,p2a, ncol=2, nrow=2, common.legend = TRUE, legend="bottom"),dpi = 300, width =12, height = 8,units = c("in"))
  
  #######################################################################################
  #cols<-brewer.pal(n =9, name = 'Paired')
  cols<-c("#A6CEE3", "#1F78B4", "#B2DF8A", "#33A02C", "#FB9A99", "#E31A1C", "black", "#FF7F00","#CAB2D6")
  
  p3a<-ggplot(r1,aes(x=TAD,y=DV))+
    geom_point(aes(color=STUDYID))+
    #geom_point(aes(col=factor(STUDYID)))+
    xlab(TAD)+ylab(conc)+
    #geom_abline(slope=1)+
    scale_y_log10(limits =limx,breaks = scales::trans_breaks("log10", function(x) 10^x),labels = scales::trans_format("log10", scales::math_format(10^.x)))+
    
    scale_x_continuous(breaks =seq(0,1500,7*24),
                       labels =seq(0,1500,7*24))+
    annotation_logticks(sides = "l") +
    geom_smooth(aes(y=DV,col="LOESS_OBS"),se=F)+
    geom_smooth(aes(y=IPRED,col="LOESS_IPRED"),se=F)+
    geom_smooth(aes(y=PRED,col="LOESS_PRED"),se=F)+
    scale_colour_manual(name="",values=cols, 
                        guide = "legend") + 
    theme_bw()+
    guides(colour=guide_legend(override.aes=list(linetype=c(rep(0,6),1,1,1),shape=c(rep(16,6),NA,NA,NA))))+
    theme(axis.text.y = element_text(angle = 90, hjust =0.5, vjust = 0.5),legend.position=c(0.8,0.2),legend.box.margin=NULL, legend.background = element_rect(fill="transparent"))
  
  p3b<-ggplot(r1,aes(x=TAD,y=DV))+
    geom_point(aes(color=STUDYID))+
    #geom_point(aes(col=factor(STUDYID)))+
    xlab(TAD)+ylab("Concentration (ng/mL)")+
    #geom_abline(slope=1)+
    scale_x_continuous(breaks =seq(0,1500,7*24),
                       labels =seq(0,1500,7*24))+
    #scale_y_continuous(limits=limx,breaks =seq(0,10^6,10^5),labels = seq(0,10))+
    geom_smooth(aes(y=DV,col="LOESS_OBS"),se=F)+
    geom_smooth(aes(y=IPRED,col="LOESS_IPRED"),se=F)+
    geom_smooth(aes(y=PRED,col="LOESS_PRED"),se=F)+
    scale_colour_manual(name="",values=cols, 
                        guide = "legend") + 
    theme_bw()+
    guides(colour=guide_legend(override.aes=list(linetype=c(rep(0,6),1,1,1),shape=c(rep(16,6),NA,NA,NA))))+
    theme(axis.text.y = element_text(angle = 0, hjust =0.5, vjust = 0.5),legend.position=c(0.8,0.2),legend.box.margin=NULL, legend.background = element_rect(fill="transparent"))
  
  
  p3c<-ggplot(r1,aes(x=TAD,y=DV))+
    geom_point(aes(color=STUDYID))+
    #geom_point(aes(col=factor(STUDYID)))+
    xlab(TAD)+ylab(conc)+
    #geom_abline(slope=1)+
    scale_y_log10(limits =limx,breaks = scales::trans_breaks("log10", function(x) 10^x),labels = scales::trans_format("log10", scales::math_format(10^.x)))+
    scale_x_log10()+
    #scale_x_continuous()+
    annotation_logticks(sides = "l") +
    geom_smooth(aes(y=DV,col="LOESS_OBS"),se=F)+
    geom_smooth(aes(y=IPRED,col="LOESS_IPRED"),se=F)+
    geom_smooth(aes(y=PRED,col="LOESS_PRED"),se=F)+
    scale_colour_manual(name="",values=cols, 
                        guide = "legend") +
    theme_bw()+
    guides(colour=guide_legend(override.aes=list(linetype=c(rep(0,6),1,1,1),shape=c(rep(16,6),NA,NA,NA))))+
    theme(axis.text.y = element_text(angle = 90, hjust =0.5, vjust = 0.5),legend.position=c(0.8,0.2),legend.box.margin=NULL, legend.background = element_rect(fill="transparent"))
  
  ggsave("profile.png",
         ggarrange(p3a,p3b,p3c, ncol=2, nrow=2, common.legend = TRUE, legend="bottom"),dpi = 300, width =12, height = 8,units = c("in"))
  
  ########################################################################################
  
  #cols<-c("grey","black","blue")
  cols<-c("#A6CEE3", "#1F78B4", "#B2DF8A", "#33A02C", "#FB9A99", "#E31A1C", "black", "#FF7F00")
  p4<-ggplot(r1,aes(x=PRED,y=CWRES))+
    geom_point(aes(col=factor(STUDYID)))+
    xlab(PRED)+ylab(CWRES)+
    geom_hline(aes(yintercept=0, col="Line_zero")) +
    geom_smooth(aes(col="LOESS"),se=F)+
    scale_x_continuous(labels = function(x) format(x, scientific =F))+
    scale_y_continuous(limits=c(-6,6),breaks=seq(-6,6,2))+
    theme_bw()+
    scale_colour_manual(name="Observed",values=cols, 
                        guide = "legend") + 
    guides(colour=guide_legend(override.aes=list(linetype=c(rep(0,6),1,1),shape=c(rep(16,6),NA,NA))))+
    theme(legend.position=c(0.8,0.2),legend.box.margin=NULL, legend.background = element_rect(fill="transparent"))
  
  p4a<-ggplot(r1,aes(x=IVAR,y=CWRES))+
    geom_point(aes(col=factor(STUDYID)))+
    xlab(IVAR)+ylab(CWRES)+
    geom_hline(aes(yintercept=0, col="Line_zero")) +
    geom_smooth(aes(col="LOESS"),se=F,span=1)+
    # scale_x_continuous(labels = function(x) format(x, scientific =F))+
    scale_x_continuous(breaks =seq(0,1500,7*24),
                       labels =seq(0,1500,7*24))+
    scale_y_continuous(limits=c(-6,6),breaks=seq(-6,6,2))+
    theme_bw()+
    scale_colour_manual(name="Observed",values=cols, 
                        guide = "legend") + 
    guides(colour=guide_legend(override.aes=list(linetype=c(rep(0,6),1,1),shape=c(rep(16,6),NA,NA))))+
    theme(legend.position=c(0.8,0.2),legend.box.margin=NULL, legend.background = element_rect(fill="transparent"))
  
  p4b<-ggplot(r1,aes(x=TAD,y=CWRES))+
    geom_point(aes(col=factor(STUDYID)))+
    xlab(TAD)+ylab(CWRES)+
    geom_hline(aes(yintercept=0, col="Line_zero")) +
    geom_smooth(aes(col="LOESS"),se=F,span=1)+
    # scale_x_continuous(labels = function(x) format(x, scientific =F))+
    scale_x_continuous(breaks =seq(0,1500,7*24),
                       labels =seq(0,1500,7*24))+
    scale_y_continuous(limits=c(-6,6),breaks=seq(-6,6,2))+
    theme_bw()+
    scale_colour_manual(name="Observed",values=cols, 
                        guide = "legend") + 
    guides(colour=guide_legend(override.aes=list(linetype=c(rep(0,6),1,1),shape=c(rep(16,6),NA,NA))))+
    theme(legend.position=c(0.8,0.2),legend.box.margin=NULL, legend.background = element_rect(fill="transparent"))

  ggsave("CWRES.png",
         ggarrange(p4,p4a,p4b, ncol=1, nrow=3, common.legend = TRUE, legend="bottom"),dpi = 300, width =12, height = 8,units = c("in"))
  
}

#' template for VPC
#'
#' @param data Data frame Use the codes below as template
#' @keywords vpc_plots()
#' @export
#' @examples 

vpc_plots<-function(...){
  
  range(output$SIM.bin)
  out1$SIM.bin<-as.numeric(as.character(out1$SIM.bin))
  
  p<-ggplot(out1, aes(x=SIM.bin,group=OBS.variable)) +
    geom_ribbon(aes(ymin=SIM.low, ymax=SIM.up, fill=SIM.variable, col=SIM.variable, group=SIM.variable), alpha=0.1,col=NA) +
    geom_line(aes(y=SIM.med, col=SIM.variable, group=SIM.variable), size=1) +
    geom_line(aes(y=OBS.value, linetype=SIM.variable), size=1) +
    scale_y_log10()+
    scale_x_continuous(breaks=seq(0,100,12))+
    scale_colour_manual(
      name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
      breaks=c("qt05", "qt50", "qt95"),
      values=c("red", "blue", "red"),
      labels=c("5%", "50%", "95%")) +
    scale_fill_manual(
      name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
      breaks=c("qt05", "qt50", "qt95"),
      values=c("red", "blue", "red"),
      labels=c("5%", "50%", "95%")) +
    scale_linetype_manual(
      name="Observed Percentiles\n(black lines)",
      breaks=c("qt05", "qt50", "qt95"),
      values=c("dotted", "solid", "dashed"),
      labels=c("5%", "50%", "95%")) +
    guides(
      fill=guide_legend(order=2),
      colour=guide_legend(order=2),
      linetype=guide_legend(order=1)) +
    theme(legend.position="top",
          legend.key.width=grid::unit(2, "cm")) +
    labs(x="Time (h)", y="Concentration (ng/mL)")
  p+geom_point(aes(y=DV),col="grey")
  
  blq<-out1[,c("SIM.bin","blqo","blqs")]
  blq$SIM.bin<-as.numeric(as.character(blq$SIM.bin))
  
  ##EXAMPLE BLQ PLOT
  ggplot(blq,aes(x=SIM.bin,y=blqo))+
    geom_point(aes(col="Observed"))+
    geom_line(aes(col="Observed"))+
    geom_point(aes(x=SIM.bin,y=blqs,col="Predicted"))+
    geom_line(aes(x=SIM.bin,y=blqs,col="Predicted"))
   }


#' Compute VPC stats 
#'
#' @param obs.data Observed data. Type vpc_plots forbasic R code to be used with the output 
#' @param sim.data Sim data
#' @param bin Binning. This should be created in both datasets
#' @param prob Quantile probability.  
#' @param sort Sorting as vector 
#' @param dv Concentration or response.Same name in both datasets   
#' @param tad Time after dose. In obs.data
#' @param rtime Time after first dose. In obs.data
#' @param blq LLOQ. Could be single value or vector. If vector, both datasets should have the same variable name. This will compute the % of BLQ for each bin.
#' @param replicate Replicate in sim.data only
#' @param pred.corr Pred-corrected. ex: C("PRED","lin") use the column PRED and corrected in linear scale or log for log scale
#' @keywords lhvpc_stat()
#' @export
#' @examples 

lhvpc_stat<-function (obs.data = obs, sim.data = sim, bin = "bin", prob = c(0.05, 
                                                                            0.5, 0.95), sort = NULL, dv = "DV", tad = "TAD", rtime = "IVAR", 
                      blq = NULL, replicate = "REPLICATE", pred.corr = NULL) 
{
  library(reshape)
  if (!is.null(pred.corr)) {
    medpred <- median(obs.data[, pred.corr[1]])
    if (pred.corr[2] == "lin") {
      obs.data[, dv] <- obs.data[, dv] * medpred/obs.data[, 
                                                          pred.corr[1]]
    } else {
      obs.data[, dv] <- obs.data[, dv] + medpred - obs.data[, 
                                                            pred.corr[1]]
    }
  } else {
    obs.data[, dv] <- obs.data[, dv]
  }
  var <- NULL
  for (i in 1:length(prob)) {
    namqt <- paste0("qt", prob[i] * 100)
    exp<-paste0("quantile(x,",prob[i],")")
    var <- c(var, namqt)
    if (prob[i] == prob[1]) {
      obs1 <- addvar(obs.data, c(sort, "bin"), dv, exp, 
                     "no", namqt)
    } else {
      obs1 <- dplyr::left_join(obs1, addvar(obs.data, c(sort, 
                                                        "bin"), dv,exp, "no", namqt))
    }
  }
  obs1 <- lhlong(obs1, var)
  if (!is.null(pred.corr)) {
    medpred <- median(sim.data[, pred.corr[1]])
    if (pred.corr[2] == "lin") {
      sim.data[, dv] <- sim.data[, dv] * medpred/sim.data[, 
                                                          pred.corr[1]]
    } else {
      sim.data[, dv] <- sim.data[, dv] + medpred - sim.data[, 
                                                            pred.corr[1]]
    }
  } else {
    sim.data[, dv] <- sim.data[, dv]
  }
  var <- NULL
  for (i in 1:length(prob)) {
    namqt <- paste0("qt", prob[i] * 100)
    exp<-paste0("quantile(x,",prob[i],")")
    var <- c(var, namqt)
    if (prob[i] == prob[1]) {
      s1 <- addvar(sim.data, c(sort, bin, replicate), dv, 
                   exp, "no", namqt)
    }else {
      s1 <- dplyr::left_join(s1, addvar(sim.data, c(sort, 
                                                    bin, replicate), dv, exp, "no", 
                                        namqt))
    }
  }
  s1 <- lhlong(s1, var)
  namvar <- c("low", "med", "up")
  for (j in 1:3) {
    exp<-paste0("quantile(x,",prob[j],")")
    if (j == 1) {
      s2 <- addvar(s1, c(bin, "variable"), "value", exp, 
                   "no", namvar[1])
    } else {
      s2 <- dplyr::left_join(s2, addvar(s1, c(bin, "variable"), 
                                        "value", exp, "no", namvar[j]))
    }
  }
  if (!is.null(blq)) {
    if (is.numeric(blq)) {
      obs.data$blq <- obs.data[, dv] <= blq
      sim.data$blq <- sim.data[, dv] <= blq
    }
    else {
      obs.data$blq <- obs.data[, dv] <= obs.data[, blq]
      sim.data$blq <- sim.data[, dv] <= sim.data[, blq]
    }
    blqo <- addvar(obs.data, bin, "blq", "sum(x)", "no", 
                   "blqo")
    blqo <- dplyr::left_join(blqo, addvar(obs.data, bin, 
                                          "blq", "length(x)", "no", "nblqo"))
    blqo$blqo <- blqo$blqo/blqo$nblqo * 100
    blqs <- addvar(sim.data, bin, "blq", "sum(x)", "no", 
                   "blqs")
    blqs <- dplyr::left_join(blqs, addvar(sim.data, bin, 
                                          "blq", "length(x)", "no", "nblqs"))
    blqs$blqs <- blqs$blqs/blqs$nblqs * 100
  } else {
    blqo <- NULL
    blqs <- NULL
  }
  names(obs1) <- paste0("OBS.", names(obs1))
  names(s2) <- paste0("SIM.", names(s2))
  output <- lhcbind(obs1, s2)
  out1 <- dplyr::left_join(lhmutate(obs.data[, c(sort, "bin", 
                                                 dv, tad, rtime)], "bin=SIM.bin"), output)
  if (!is.null(blqo)) {
    out1 <- dplyr::left_join(lhmutate(blqo, "bin=SIM.bin"), 
                             out1)
    out1 <- dplyr::left_join(lhmutate(blqs, "bin=SIM.bin"), 
                             out1)
  } else {
    out1 <- out1
  }
  out1
}

#' Look for keyword across dataset
#'
#' @param data Data frame 
#' @param var Variable to be cut
#' @param file.or.path File path ("xxx/zzz/") or data frame (c("lb","dm")) 
#' @param fpattern File extention: ex. "csv", all csv file in folder will be used  
#' @param filename specify full file name. set fpattern to Null
#' @keywords lhseek()
#' @export
#' @examples 




lhseek<-function (var, file.or.path, fpattern = NULL, filename = NULL)
  
{
  
  if (!is.null(fpattern)) {
    
    list <- dir(file.or.path)[grep(fpattern, dir(file.or.path))]
    
  }
  
  if (!is.null(filename)) {
    
    list = filename
    
  }
  
  if (is.null(filename) & is.null(fpattern)) {
    
    list = file.or.path
    
  }
  
  out <- NULL
  
  for (i in tolower(var)) {
    
    for (j in list) {
      
      b <- function(x) {
        
      }
      
      if (!is.null(filename) | !is.null(fpattern)) {
        
        d <- file.or.path
        
        xt <- substring(j, nchar(j) - 2)
        
        if(xt=="dat"){xt="sas"}else{xt=xt}
        
        rt <- c("read.csv", "haven::read_xpt", "haven::read_sas")
        
        rt <- rt[grep(xt, rt)]
        
        txt <- paste0(rt, "(file.path(d,j))")
        
      } else {
        
        txt <- j
        
      }
      
      body(b) <- parse(text = txt)
      
      inp <- as.data.frame(b())
      
      for (z in names(inp)) {
        
        if (length(tolower(z)[grep(i, tolower(z))]) >
            
            0) {
          
          xx <- paste0(j, ":", z)
          
        } else {
          
          xx = NULL
          
        }
        
        if (length(unique(inp[, z][grep(i, tolower(inp[,
                                                       
                                                       z]))])) > 0) {
          
          out1 <- paste0(xx, j, ";", z, unique(inp[,
                                                   
                                                   z][grep(i, inp[, z])]), ":", i)
          
        } else {
          
          out1 <- c(xx, NULL)
          
        }
        
        out <- c(out, out1)
        
      }
      
      print(out)
      
    }
    
  }
  
}





#' Cut values and create category
#'
#' @param data Data frame 
#' @param var Variable to be cut
#' @param breaks break points 
#' @param labels category name. If fancy, the categories will be created according to the break points 
#' @param right If false, right value will be exclusive
#'  @param newvar vector name of the categorical. If default, the var with suffix"cat" will be used as default name
#' @keywords lhcut()
#' @export
#' @examples 

lhcut<-function(data,var,breaks,labels="fancy",right=F,newvar="default"){
  brk=c(min(data[,var]),breaks,max(data[,var])^2)
  if(newvar=="default"){nvar=paste0(var,"cat")}else{nvar=newvar}
  if(labels=="fancy"){
    lab1=c(paste0("<",breaks))
    lab2<-c(paste0(">=",breaks))
    lab3<-c(paste0("<=",breaks))
    lab4<-c(paste0(">",breaks))
    lab11<-NULL;lab22<-NULL
    for(i in 1:(length(breaks)-1)){
      lab11<-c(lab11,paste(lab2[i],"&",lab1[i+1]))
      lab22<-c(lab22,paste(lab4[i],"&",lab3[i+1]))
    }
    if(right){
      labels1<-c(lab3[1],lab22,lab4[length(lab4)])}else{ 
        labels1<-c(lab1[1],lab11,lab2[length(lab2)])
      } }else{labels1=labels}
  
  data[,nvar]<-cut(data[,var],breaks=brk,labels=labels1,right=right)
  print(addvar(data,nvar,var,"range(x)","no"))
  data}


#' Change factor level of a variable using matched level of another variable in the same dataset
#'
#' @param data Data frame 
#' @param leader lead Variable to be used for factor level of the follower variable
#' @param follower follower Variable  
#' @keywords lhfactor()
#' @export
#' @examples 

lhfactor<-function(data=test,leader="AGEcat",follower="catt"){
  lab<-nodup(data,c(leader,follower),"var");lab<-lab[order(lab[,leader]),follower]
  data<-reflag(data,follower,lab)
}




#' Combine variables in the same column
#'
#' @param data Data frame 1 and 2 with long vectors and values. Note: no duplicated sorting vector allowed 
#' @param combine.var Variable name, c(var1,var2). var1 will be stacked over var 2
#' @keywords stackvar()
#' @export
#' @examples
stackvar<-function(data,combine.var=c("xxx","variable")){
  z$dum<-seq(nrow(z))
  z1<-z
  z1$dum<-z1$dum-1
  z1<-nodup(z1,combine.var[1],"all")
  keep<-c(combine.var[1],combine.var[2],"dum")
  z1[,combine.var[2]]<-z1[,combine.var[1]];z1[,!names(z1)%in%keep]<-""
  z<-rbind(z,z1[,names(z)]); z<-z[order(z$dum),]
  z[,c(combine.var[1],"dum")]<-NULL
  z}



#' mutate variable names
#'
#' @param data Data frame 1 and 2 with long vectors and values. Note: no duplicated sorting vector allowed 
#' @param mutate Vector to be mutated ex. "xxx=yyy" for renaming xxx as yyy
#' @keywords lhmutate()
#' @export
#' @examples
lhmutate<-function(data,mutate){
  keep<-sub("=.*","",mutate)%in%names(data)
  imp<-sub(".*=","",mutate)[keep]
  bimp<-sub("=.*","",mutate)[keep]
  print(c("Not found:",sub("=.*","",mutate)[!sub("=.*","",mutate)%in%names(data)]))
  
  for(i in 1:length(bimp)){
    names(data)[names(data)==bimp[i]]<-imp[i]
  }
  data
}



#' merge and make long to wide data frame
#'
#' @param dat1,dat2 Data frame 1 and 2 with long vectors and values. Note: no duplicated sorting vector allowed 
#' @param long.vector Vector containing wide variable names
#' @param value Column containing values associated with each wide variable. 
#' @keywords lhmerge_long_data()
#' @export
#' @examples
lhmerge_long_data<-function (dat1=dd1, dat2=dd2, long.vector1="PARAMCD", long.vector2="PARAMCD", value1="AVAL", value2="CHG") 
{
  names(dat1)[names(dat1) == long.vector1] <- "x"
  names(dat1)[names(dat1) == value1] <- "y"
  names(dat2)[names(dat2) == long.vector2] <- "x"
  names(dat2)[names(dat2) == value2] <- "y"
  dat2$x[dat2$x == dat1$x] <- paste0(dat2$x[dat2$x == dat1$x], 
                                     "dat2")
  datall <- lhrbind(dat1, dat2)
  head(datall)
  if(nrow(dup1(datall, names(datall)[!names(datall)%in%"y"],"all"))>0){
    print("duplicated data")
    output<-dup1(datall, names(datall)[!names(datall)%in%"y"],"all")
  }else{output <- lhwide(datall, "y", "x")}
}

#' merge wide data frames
#'
#' @param dat1,dat2 Data frame 1 and 2 with long vectors and values. Note: no duplicated sorting vector allowed 
#' @param by1,by2 Vectors containing matched values from each data frames. 
#' @keywords lhmerge_long_data()
#' @export
#' @examples
lhmerge_wide_data<-function(dat1,dat2,by1,by2){
  d1<-chclass(dat1,names(dat1),"char");d2<-chclass(dat2,names(dat2),"char")
  d1<-lhlong(d1,names(d1)[!names(d1)%in%by1]);d2<-lhlong(d2,names(d2)[!names(d2)%in%by2])
  head(d1)
  d1<-chclass(d1,names(d1),"char");d2<-chclass(d2,names(d2),"char")
  d2$variable<-as.character(d2$variable);d1$variable<-as.character(d1$variable)
  d2$variable[d2$variable%in%d1$variable]<-paste0(d2$variable[d2$variable%in%d1$variable],"d2")
  datall<-lhrbind(d1,d2)
  head(datall)
  unique(datall$variable)
  test<-nodup(datall,c(by1,by2,"variable"),"all")
  test<-test[order(test$USUBJID,as.Date(test$date)),]
  output<-lhwide(test,"value","variable")
}


#' Reshape wide
#'
#' @param data Dataset
#' @param wide.data Name of vector containing data to be dcasted
#'  @param wide.vector Name of vector to be reshape as heading 
#' @param data Dataset 
#' @keywords lhwide()
#' @export
#' @examples
#' lhwide()
lhwide<-function(data,wide.data,wide.vector){
  b <- function(x) {}
  x1<-paste(paste(names(data)[!names(data)%in%c(wide.data,wide.vector)],collapse="+"),"~",wide.vector)
  body(b) <- parse(text = x1)
  z1<-dcast(data,b())}

#' Reshape long
#'
#' @param data Dataset
#' @param long.vector List of vector to be melted
#'
#' @keywords lhlong()
#' @export
#' @examples
#' lhlong()
lhlong<-function(data,long.vector){
  z1<-melt(data,names(data)[!names(data)%in%long.vector])
}


#' find different values between two datasets
#'
#' @param dat1,dat2 Dataset 1 and 2" 
#' @keywords findiff()
#' @export
#' @examples
#' findiff()

findiff<-function(dat1,dat2){
 # stopifnot(nrow(dat1)==nrow(dat2))
  dum1a<-""
  nm1<-"dat1"
  dum2a<-""
  nm2<-"dat2"
  for(i in 1:length(names(dat1))){
    nm1<-paste(nm1,names(dat1)[i],sep="/")
    dum1a<-paste(dum1a,dat1[,names(dat1)[i]],sep="/")
    nm2<-paste(nm2,names(dat2)[i],sep="/")
    dum2a<-paste(dum2a,dat2[,names(dat2)[i]],sep="/")
  }
  a<-setdiff(dum1a,dum2a)
  b<-setdiff(dum2a,dum1a)
  out<-data.frame(nm1=unique(a));names(out)<-nm1
  out1<-data.frame(nm2=unique(b));names(out1)<-nm2
  row1<-data.frame(N1=length(dum1a),N2=length(dum2a))
  out3<-lhcbind(out,out1)
  out3<-lhcbind(out3,row1)
  out3 }


#' lhjoin funtion
#'Join two datasets and print report of joining procedure
#' @param dat1,by1 Data frame 1 and variables to be matched. If NULL, match="all"
#' @param dat2,by2 Data frame 2 and variables to be matched. If by1=NULL then by2=NULL then match="all"
#' @param type could be "full", "left","right" or "inner"
#' @keywords lhjoin()
#' @export
#' @examples
#' lhjoin()
lhjoin<-function(dat1,by1=NULL,dat2,by2=NULL,type="full"){
  invar<-intersect(names(dat1),names(dat2))
  if(is.null(by1)){
    by1=invar}else{
      by1=by1}
  if(is.null(by2)){
    by2=invar
  }else{
    by2=by2
    names(dat2)[names(dat2)%in%invar]<-paste0("df2_",names(dat2)[names(dat2)%in%invar])
  }
  
  if(length(by1)>1){
    dat1[,"dum"]<-dat1[,by1[1]]
    for(i in 2:length(by1)){
      dat1[,"dum"]<-paste(dat1[,"dum"],dat1[,by1[i]],sep="-")
    }}else{dat1[,"dum"]<-dat1[,by1[1]]}
  
  by2[!by2%in%by1]<-paste0("df2_",by2[!by2%in%by1])
  
  if(length(by2)>1){
    dat2[,"dum"]<-dat2[,by2[1]]
    for(i in 2:length(by2)){
      dat2[,"dum"]<-paste(dat2[,"dum"],dat2[,by2[i]],sep="-")
    }}else{dat2[,"dum"]<-dat2[,by2[1]]}
  
  dat<-plyr::join(dat1,dat2,by="dum",type=type)
  
  report<-data.frame(nrow_data1=nrow(dat1),
                     nrow_data2=nrow(dat2),
                     nrow_joint=nrow(dat))
  for(c in 1:length(by1)){    
    x<-data.frame(z=setdiff(dat1[,by1[c]],dat2[,by2[c]]))
    names(x)<-paste0(by1[c],"_not_in_data2")
    y<-data.frame(z=setdiff(dat2[,by2[c]],dat1[,by1[c]]))
    names(y)<-paste0(by1[c],"_not_in_data1")
    zz<-lhcbind(x,y)
    report<-lhcbind(report,zz)
  }
  print(report)
  dat}



#' lhorder funtion
#'
#' Make simple table. Use data frame created by addvar2
#' @param dat Datframe
#' @param var Order by variables. ex: ":Trt,:Agegr"
#' @keywords lhorder()
#' @export
#' @examples
#' lhorder()

lhorder<-function(dat,var){
  data<-dat
  x<-paste0("data[order(",gsub(":","data$",var),"),]")
  b<- function(x) {}
  body(b) <- parse(text = x)
  data<-b()
}



#' lhtab funtion
#'
#' Make simple table. Use data frame created by addvar2
#' @param data Datframe
#' @param vh Vertical and horizontal headers. ex: "Trt+Agegr~Param"
#' @param value Values. Example: c("mean","SD","Mean (CV)")
#' @param ord Order variables. ex: ":Trt,:Agegr"
#' @param save.name Save table as word document. Enter the file name: ex "test.docx"
#' @param output output="csv" for csv output, else output will be in FlexTable format
#' @keywords lhtab()
#' @export
#' @examples
#' lhtab()

lhtab<-function (data, vh, value, ord = NULL, save.name = NULL, output = "csv") 
{
  library(reshape)
  b <- function(x) {
  }
  body(b) <- parse(text = vh)
  
  #vh<-"group~label+ss+test"
  v <- gsub("+", ":", sub("~.*", "", vh), fixed = T)
  v <- unlist(strsplit(v, ":")) 
  h <- gsub("+", ":", sub(".*~", "", vh), fixed = T)
  list(gsub(":", ",", h))
  h <- unlist(strsplit(h, ":"))
  data$dum<-""
  for(i in h){
    if(i==h[1]){
      data$dum<-data[,i]}else{data$dum<-paste(data$dum,data[,i],sep="_")}
  }  
  
  
  hd<-nodup(data,c("dum",v,h),"var")
  w<-NULL
  for(uu in value){
    data[,"stats"]<-uu
    w1 <- reshape(data[,c(v,"dum",uu,"stats")], 
                  timevar ="dum",
                  idvar =c(v,"stats"),
                  direction = "wide")
    for(u in names(data[,c(v,"dum",uu,"stats")])){
      rm<-paste0(u,".")
      names(w1)<-gsub(rm,"",names(w1),fixed = T)
    }
    w<-rbind(w,w1)}
  hw<-NULL
  
  for(d in h){
    hd[,"stats"]<-"stats"
    hw1 <- reshape(hd[,c("dum",v,d,"stats")], 
                   timevar ="dum",
                   idvar =c(v,"stats"),
                   direction = "wide")
    
    for(u in names(hd[,c("dum",v,d,"stats")])){
      rm<-paste0(u,".")
      names(hw1)<-gsub(rm,"",names(hw1),fixed = T)
    }
    hw1<-nodup(hw1,names(hw1)[!names(hw1)%in%c(v,"stats")],"all")
    hw<-rbind(hw,hw1)
  }
  hw<-hw[,unique(names(hw))]
  for(vv in v){
    hw[,vv]<-vv 
  }
  
  setdiff(names(w),names(hw))
  
  hw1<-rbind(hw,w)
  head(hw1,10)
  
  
  if (!is.null(ord)) {
    y <- paste0(ord, ",:stats")
  }else {
    y <- ":stats"
  }
  stor<-c("stats",value)
  hw1[,"stats"]<-factor(hw1[,"stats"],level=stor)
  head(hw1)
  
  hw1 <- lhorder(hw1,y)
  
  hw1 <- chclass(hw1, names(hw1), "char")
  tab <- ReporteRs::FlexTable(hw1, header.columns = FALSE)
  
  for (y in c(v,"stats")) {
    tab = ReporteRs::spanFlexTableRows(tab, j = y, runs = as.character(hw1[, 
                                                                           y]))
  }
  t4 <- hw1
  colnames(t4) <- NULL
  rownames(t4) <- NULL
  for (z in 1:length(h)) {
    tab = ReporteRs::spanFlexTableColumns(tab, i = z, runs = paste(t4[z, 
                                                                      ]))
  }
  tab[1:length(h), ] = ReporteRs::textProperties(font.weight = "bold")
  tab[, names(hw1)] = ReporteRs::parCenter()
  if (!is.null(save.name)) {
    doc <- ReporteRs::docx()
    doc <- ReporteRs::addFlexTable(doc, tab)
    ReporteRs::writeDoc(doc, save.name)
    ReporteRs::writeDoc(doc, save.name)
  }
  if (output != "csv") {
    res <- tab
  }
  else {
    res <- hw1
  }
  res
}




#' txt funtion
#'
#' Clone expression function for adding special formats and symbol to plots.
#' @param c text. Example: c("Concentration mg L","-1::s"," AUC::u"," Delta::i","moles::e"). s=subscript, u=underline, Delta= capital greek delta letter, i= italic, e=superscript
#' @keywords txt()
#' @export
#' @examples
#' txt()

txt<-function(c){
  z1<-""
  for(j in 1:length(c)){
    if(length(grep("::",sub(".*::","::", c[j])))==0){z=c[j]}else{
      if(length(grep(":e",sub(".*:e",":e", c[j])))!=0){
        z=paste0("^{",gsub(sub(".*::", "::",c[j]),"",c[j]),"}")}
      if(length(grep(":s",sub(".*:s",":s", c[j])))!=0){
        z=paste0("[",gsub(sub(".*::", "::",c[j]),"",c[j]),"]")}
      if(length(grep(":u",sub(".*:u",":u", c[j])))!=0){
        z=paste0(" underline(",gsub(sub(".*::", "::",c[j]),"",c[j]),")")}
      if(length(grep(":i",sub(".*:i",":i", c[j])))!=0){
        z=paste0(" italic(",gsub(sub(".*::", "::",c[j]),"",c[j]),")")}}
    z1=paste0(z1,z)}
  z1=gsub(" ","~",z1)
  z1=paste0("expression(",z1,")")
  b <- function(x) {}
  body(b)<-parse(text=z1)
  text<-b()
}




#' install.pack
#'
#' To install require packages. Use ipak function to install desired packages.
#' @param packages pre-define packages list.
#' @keywords install.pack
#' @export
#' @examples
#' install.pack()

install.pack<-function(...){
  packages <- c("SASxport", "reshape", "Hmisc", "tidyr","ReporteRs","plyr","downloader")
  ipak(packages)}


#' lhtemplate
#'
#' To install require packages. Use ipak function to install desired packages.
#' @param repo leonpheng
#' @param user logo.style
#' @keywords lhtemplate
#' @export
#' @examples
#' lhtemplate()
lhtemplate <- function(){
  require(downloader)
  dir.create("c:/lhtemplate")
  wf<-"c:/lhtemplate"
  url <- sprintf("https://github.com/%s/%s/archive/master.zip", "leonpheng","logo.style")
  tmp <- tempfile(fileext = "styleapp.zip",tmpdir=wf)
  download(url, tmp)
  unzip(tmp,exdir=wf)

  #download_repo("logo.style","leonpheng")
  zipF<-paste0(wf,"/logo.style-master/styleapp.zip")
  unzip(zipF,exdir=wf)
  zipF<-paste0(wf,"/logo.style-master/logostyle.zip")
  unzip(zipF,exdir=wf)
  frm<-dir(wf)
  index<-c(grep("zip",frm))
  frm<-frm[index]
  for(i in frm){
    file.remove(paste(wf,i,sep="/"))
  }
  unlink(paste0(wf,"/logo.style-master"), recursive = T)
}

#' lhdoc
#'
#' Create doc for word document using this fonction or you can do it manually.
#' Type lhtext and copy the template to R workspace and start writing.
#' @param t template name see at the end of lhtext function.
#' @param toc.level maximimum toc level
#' @param template Word document template could be used for styles. Styles should be mapped in style.to.map. Template is also available at github: to load it, just run  lhtemp() once to download and store the templates in your PC at "c:lhtemplate. Note that the templates and logo are also used in xptdef package.

#' @param TOC Set to F if no TOC wanted
#' @param style.to.map Map the styles in template to be used. Ex: mypar is for footnote (font size)
#' @keywords lhdoc
#' @export
#' @examples
lhdoc<-function(toc.level=4,template="default",TOC=F,
        style.to.map="default",logo.certara=F){

  library(ReporteRs)
  library(flextable)
  library(dplyr)

  if(template=="default"){
  doc<-docx(template = "c:/lhtemplate/styleapp.docx", empty_template = TRUE)}else{doc<-docx(template =template, empty_template = TRUE)}
if("default"%in%style.to.map){
  doc <- map_title(doc, stylenames =c("Heading1",
                                     "Heading2", "Heading2", "fnt"))}else{doc <- map_title(doc, stylenames =style.to.map)}
  if(logo.certara){
    doc<-addImage(doc,"c:/lhtemplate/logo.jpg", par.properties = parProperties(text.align = "center"),width = 3.35, height = 1.6)
  }
  if(TOC){
    doc <-addTOC(doc,level_max =toc.level)
  }
  doc
  }


#' lhrbind
#'
#' r bind 2 data frames regardless number of columns.
#' @param dat1,dat2 data frames.
#'
#' @keywords lhrbind
#' @export
#' @examples
#' lhrbind()

lhrbind<-function (dat1, dat2, na.replace = NA, all.character = T) 
{
  dat1[, setdiff(names(dat2), names(dat1))] <- na.replace
  dat2[, setdiff(names(dat1), names(dat2))] <- na.replace
  if (all.character) {
    dat <- rbind(chclass(dat1, names(dat1), "char"), chclass(dat2, 
                                                             names(dat2), "char"))
    print("Warning: all vectors in new dataset are character")
    dat
  }
  else (dat <- rbind(dat1, dat2))
}

#' lhcbind
#'
#' C bind 2 data frames regardless number of row length.
#' @param dat1,dat2 data frames.
#'
#' @keywords lhcbind
#' @export
#' @examples
#' lhcbind()

lhcbind<-function(dat1,dat2){
  dat1=as.data.frame(dat1)
  dat2=as.data.frame(dat2)
  r1<-nrow(dat1)
  r2<-nrow(dat2)
  if(r1>r2){
    r3=as.data.frame(matrix(ncol=ncol(dat2),nrow=r1-r2,data=""))
    names(r3)<-names(dat2)
    r3=rbind(dat2,r3)
    dat=cbind(dat1,r3)
  }
  if(r1<r2){
    r3=as.data.frame(matrix(ncol=ncol(dat1),nrow=r2-r1,data=""))
    names(r3)<-names(dat1)
    r3=rbind(dat1,r3)
    dat=cbind(r3,dat2)
  }
  if(r1==r2){dat=cbind(dat1,dat2)}
  dat
}

#' lhloess
#'
#' Compute the LOESS data for ploting.
#' @param data data.
#' @param x Independent variable
#' @param y Dependent variable
#' @param by Sort by. Only one sorting variabele is accepted. If more than 1 variables, create a unique sorting using paste(var1,var2,etc)
#' @param span LOESS stiffness
#' @keywords lhloess
#' @export
#' @examples
#' lhloess()


lhloess<-function(data,x,y,by,span=1){
  library(plyr)
  data$x=data[,x]
  data$y=data[,y]
  data$by=data[,by]
  head(data)
  dat=NULL
  for(i in unique(data$by)){
    tmp<-data[data$by==i,c(x,"x","y")]
    head(tmp)
    tmp1<-with(tmp,unlist(predict(loess(y~x,tmp,span=span),x)))
    tmp$loess<-tmp1
    dat<-rbind(dat,tmp)
  }
  #data$x<-data$y<-data$by<-NULL
  data<-join(data,dat)
}



#' load.pack1
#'
#' To install require packages. Use ipak function to install desired packages.
#' @param packages pre-define packages list.
#' @keywords load.pack1
#' @export
#' @examples
#' load.pack1()

loadpack<-function(...){
  library("ReporteRs")
  require(dplyr)
  library(plyr)
  require(stats)     # To format summary
  require(PCSmisc)     #
     #This package is needed to add a dataset label,it can only be used  with the 32bit version of R. If dataset label is not needed, the package SASxport can be used
  require(tidyr)
}


#' Ben's function
#'
#' internal used
#' @param packages pre-define packages list.
#' @keywords blk.locf
#' @export
#' @examples
#' blk.locf2()

blk.locf2<-function (x, id, na.action = c("fill", "carry.back"), fill = NA)
{
  .checkID(id)
  if (length(x) != length(id)) {
    stop("Lengths of x and id must match")
  }
  na.action <- match.arg(na.action)
  ii <- ifelse(is.na(x), NA, .myseqalong(x))
  ii <- unlist(tapply(ii, id, function(y, na.action) {
    if (all(is.na(y)))
      return(y)
    z <- cumsum(!is.na(y))
    if (na.action == "carry.back") {
      z[z == 0] <- 1
    }
    ifelse(z == 0, NA, y[!is.na(y)][ifelse(z > 0, z, 1)])
  }, na.action = na.action, simplify = FALSE))
  y <- x
  y[!is.na(ii)] <- x[ii[!is.na(ii)]]
  y[is.na(ii)] <- fill
  y
}



#########
#' TAD from ADDL
#'
#' This function allows you to derive time after dose from ADDL.
#' @param data data frame
#' @param id ID vector
#' @param ii dose interval vector
#' @param addl additional dose vector
#' @param rtime relative time after first dose vector
#' @param evid EVID vector
#' @param dose amount adminstered (ex: AMT) vector
#' @param dose.expand If "yes", all dosing rows in ADDL will be outputed
#' @keywords tad
#' @export
#' @examples
#' tad_addl()

tad_addl<-function (data, id, ii, addl, 
                    rtime, evid, dose.expand = "yes") 
{
  data <- chclass(data, c(rtime, evid, addl, ii), "num")
  data[, addl][is.na(data[, addl])] <- 0
  data[, ii][is.na(data[, ii])] <- 0
  data[, "TAD"] <- data[, "tad"] <- NULL
  nam <- names(data)
  data <- data[order(data[, id], data[, rtime]), ]
  dose <- data[data[, evid] == 1, ]
  dose[, "TAD"] <- 0
  datp <- data[data[, evid] != 1, ]
  dat0 <- data[data[, evid] == 1, ]
  datr <- NULL
  for (i in 1:nrow(dat0)) {
    dat1 <- dat0[i, ]
    if (dat1[, addl] == 0) {
      dat2 <- dat1
    }    else {
      dat2 <- as.data.frame(matrix(ncol = ncol(dat1), nrow = dat1[, 
                                                                  addl]+1))
      names(dat2) <- names(dat1)
      dat2[, names(dat2)] <- dat1
      dat2$dum <- seq(0, nrow(dat2) - 1, 1)
      dat2[, rtime] <- dat2[, rtime] + (dat2[, ii] * dat2$dum)
      dat2$dum <- NULL
    }
    datr <- rbind(datr, dat2)
  }
  dup1(datr, names(datr), "all")
  datr <- nodup(datr, names(datr), "all")
  datr[, addl] <- datr[, ii] <- 0
  setdiff(names(datr), names(datp))
  datr$loc1 <- datr[, rtime]
  datr$lhdose <- "yes"
  datp$loc1 <- NA
  datp$lhdose <- "no"
  datp1 <- rbind(datp, datr)
  datp1 <- datp1[order(datp1[, id], datp1[, rtime]), ]
  head(datp1)
  datp1 <- locf2(datp1, id, "loc1")
  datp1$TAD <- datp1[, rtime] - datp1$loc1
  datp1$TAD[datp1$TAD < 0] <- 0
  range(datp1$TAD)
  if (dose.expand != "yes") {
    d1 <- datp1[datp1$lhdose == "no", ]
    d1$loc1 <- d1$lhdose <- NULL
    data <- rbind(d1, dose)
    data <- data[order(data[, id], data[, rtime]), ]
  }  else {
    data <- datp1[, !names(datp1) %in% c("loc1", "lhdose")]
  }
  data
}


###########
#' BLQ M6 Method
#'
#' This function allows you to create data with BLQ using M6 method.
#' @param data data frame
#' @param id ID vector
#' @param evid EVID vector
#' @keywords m6
#' @export
#' @examples
#' tad_addl()
m6<-function(data,id,evid,mdv,blq.flag,time,dv,lloq){
  dat<-data
  #id="id";time="RTIME";mdv="mdv";evid="evid";blq.flag="blqf";dv="dv";lloq=0.01
  dat$cum<-cumsum(dat[,evid])
  dat$cum1<-cumsum(dat[,blq.flag])

  good<-addvar(dat[dat[,evid]==0&dat[,mdv]==0,],c(id,"cum"),time,"max(x)","no","good")
  good1<-addvar(dat[dat[,time]>0&dat[,blq.flag]==1,],c(id,"cum1"),time,"min(x)","no","good1")
  good1[,time]<-good1$good1
  good

  m4<-plyr::join(dat,good)
  m4<-plyr::join(m4,good1)
  good2<-addvar(m4[m4$good<=m4$good1,],c(id,"cum"),"good1","min(x)","no","good2")
  good2[,time]<-good2$good2

  good3<-addvar(m4[m4$good>=m4$good1,],c(id,"cum"),"good1","max(x)","no","good3")
  good3[,time]<-good3$good3

  m4<-plyr::join(m4,good2)
  m4<-plyr::join(m4,good3)
  m4$dvm6<-m4[,dv]
  m4$mdvm6<-m4[,mdv]

m4$dvm6[m4[,time]==m4$good2|m4[,time]==m4$good3]<-lloq/2
m4$mdvm6[m4[,time]==m4$good2|m4[,time]==m4$good3]<-0
m4$cum<-m4$cum1<- m4$good<-m4$good1<-m4$good2<-m4$good3<-NULL
  m4
}


#-------------------------
#' Reflag variables
#'
#' This function allows you to change variable name (ex: "M" to "Male").
#' @param dat data frame
#' @param var Vector to be changed
#' @param orignal.flag Original names (ex:c("M","F"))
#' @param new.flag New names (ex:c("Male","Female"))
#' @param newvar Create new vector
#' @keywords reflag
#' @export
#' @examples
#' reflag(dat,var="SEX",c("M","F"),c("Male","Female"),"SEXCH"))
reflag<-function (dat, var, orignal.flag, new.flag=NULL,newvar=NULL,to.factor=T,missing=c("",".","NA",NA)) 
{
  if(is.null(new.flag)){
    new.flag=orignal.flag
  }else{new.flag}
  forgot<-setdiff(dat[,var],orignal.flag)
  forgot<-forgot[!forgot%in%missing]
  print(paste("forgot:",forgot))
  stopifnot(length(forgot)==0)
  dat[,var]<-as.character(dat[,var])
  dat[dat[,var]%in%missing,var]<-"missing or unknown"
  orignal.flag<-as.character(orignal.flag)
  new.flag<-as.character(new.flag)
  dat[,var]<-as.character(dat[,var])
  if(!is.null(newvar)){
    dat[,newvar]<-factor(dat[,var],levels=c(orignal.flag,"missing or unknown"),
                         labels=c(new.flag,"missing or unknown"))
    if(to.factor==F){
      dat[,newvar]<-as.character(dat[,newvar])
    }}else{dat[,var]<-factor(dat[,var],levels=c(orignal.flag,"missing or unknown"),
                             labels=c(new.flag,"missing or unknown"))
    if(to.factor==F){
      dat[,var]<-as.character(dat[,var])
    }}
  dat
}
#-------------------------
#' Derived 1 variable and 1 function
#'
#' This function allows you to add derived variable (ex: add mean value by ID).
#' @param dat data frame
#' @param sort sort derived variable by (ex:c("ID","SEX"))
#' @param var variable to be derived
#' @param fun deriving funtion ex:"mean(x)")
#' @param add.to.data if "yes" result will be appended to dat
#' @param name column name of derived variable
#' @keywords addvar
#' @export
#' @examples
#' addvar()
addvar<-function(dat,sort,var,fun,add.to.data="yes",name=NULL){
  library(plyr)
  d<-dat
  a<-fun
  if(is.null(name)){name=paste0(gsub("(x)",var,fun))}
  b<-function(x){}
  body(b)<-parse(text=a)

  if(length(sort)>1){dd<-(aggregate(d[,var],d[,sort],b))}else{dd<-(aggregate(d[,var],list(d[,sort]),b))}
  names(dd)<-c(sort,name)
  if(add.to.data=="yes"){out<-plyr::join(dat,dd,type="left")}else{out<-dd}}

#-------------------------
#' Derived more variables and functions
#'
#' This function allows you to add derived variable (ex: add mean value by ID).
#' @param dat data frame
#' @param sort sort derived variable by (ex:c("ID","SEX"))
#' @param var variable to be derived
#' @param fun deriving funtion ex:c("mean(x)=mean","length(x[is.na(x)])")
#'
#' @keywords addvar
#' @export
#' @examples
#' addvar()

addvar2<-function (dat, sort, var, fun, rounding = "sigfig(x,3)") 
{
  tmp1 <- NULL
  stn <- NULL
  for (z in 1:length(fun)) {
    fy = gsub("=", "", sub(".*=", "=", fun[z]))
    fx <- gsub(sub(".*=", "=", fun[z]), "", fun[z])
    tmp <- NULL
    stn <- c(stn, fy)
    for (v in var) {
      t <- addvar(dat = dat, sort = sort, var = v, fun = fx, 
                  add.to.data = "no", name = fy)
      t$var <- v
      tmp <- rbind(tmp, t)
      tmp[, fy] <- as.numeric(as.character(tmp[, fy]))
      rounding1 <- "round(x,10)"
      if (!fy %in% c("N", "n")) {
        a <- gsub("x", "tmp[,fy]", rounding1)
        b <- function(x) {
        }
        body(b) <- parse(text = a)
        tmp[, fy] <- b()
      }
      else {
        tmp
      }
    }
    if (z == 1) {
      tmp1 <- tmp
    }
    else {
      tmp1 <- join(tmp1, tmp)
    }
  }
  a <- rounding
  b <- function(x) {
  }
  body(b) <- parse(text = a)
  for (z in stn[!stn %in% c("N", "n")]) {
    for(x in 1:nrow(tmp1)){
      tmp1[,z] <- b(as.numeric(tmp1[, z]))
    }
  }
  tmp1
}


#######ADD TIME########
#-------------------------
#' Add time in hour to calendar date/time
#'
#' This function allows you to add derived variable (ex: add mean value by ID).
#' @param datetime date/time vector to be computed
#' @param timehour Time to be added in hour
#' @param format date and time format
#' @param tz Time zone (default="GMT")
#' @param add.to.data if "yes" result will be appended to dat
#' @param name column name of derived variable
#' @keywords addtime
#' @export
#' @examples
#' addtime()

addtime<-function(datetime,timehour,format="%Y-%m-%d %H:%M",tz="GMT"){
  output<-substring(strptime(datetime,format=format,tz=tz)+timehour*60*60,1,16)
  output}


###TAD calculation using elapsed RTIME#########
#-------------------------
#' Derive TAD from RTIME
#'
#' This function allows you to derive Time after dose from Time after first dose.
#' @param data data frame
#' @param id subject id
#' @param time time after first dose
#' @param evid evid
#' @keywords rt2tad
#' @export
#' @examples
#' rt2tad()
rt2tad<-function(data,id,time,evid){
  #id="uid";time="T";evid="evid"
  data$cumsum<-unlist(tapply(data[,evid],list(data[,id]),cumsum))
  nrow(data)
  data$time<-data[,time]
  d<-data[data[,evid]==1,c(id,time,"cumsum")]
  names(d)<-c(id,"time1","cumsum")
  data$sort<-seq(1,nrow(data),1)
  d<-nodup(d,c(id,"cumsum"),"all")
  data1<-plyr::join(data,d,type="left")
  data1<-chclass(data1,c("time","time1"),"num")
  data1$tad<-with(data1,time-time1)
  data1$ndose<-data1$cumsum
  data1<-data1[order(data1$sort),];data1$sort<-data1$time1<-data1$time<-data1$cumsum<-NULL
  data1
}
####create NMEM UNIQUE SUBJECT####
#-------------------------
#' NMID
#'
#' This function allows you to create NMID.
#' @param data data frame
#' @param id subject id
#' @param varname column name
#' @keywords nmid
#' @export
#' @examples
#' nmid()
nmid<-function(data,id,varname="NMID"){
  id=id;varname=varname
  data<-data
  data$ord<-seq(1,nrow(data),1)
  idat<-data.frame(id=unique(data[,id]),varname=seq(1,length(unique(data[,id])),1))
  names(idat)<-c(id,varname)
  data<-merge(data,idat)
  data<-data[order(data[,varname],data$ord),]
  data$ord<-NULL
  data
}
###Elapse Time###
#-------------------------
#' Compute delta using calendar date and time
#'
#' This function allows you to compute the delta (time1-time2).
#' @param tm1 data frame
#' @param tm2 subject id
#' @param form1 date/time format 1
#' @param form2 date/time format 2
#' @keywords diftm
#' @export
#' @examples
#' diftm()
diftm<-function(tm1,tm2,unit="hour",form1="%Y-%m-%d %H:%M",form2="%Y-%m-%d %H:%M",tz="GMT"){
dat<- as.numeric(difftime(strptime(tm1,format=form1,tz=tz),strptime(tm2,format=form2,tz=tz), units=unit))
 dat}

######change date and time format
#-------------------------
#' Reformat calendar date/time
#'
#' This function allows you to change the date/time format.
#' @param dttm original date/time
#' @param tm2 subject id
#' @param form1 date/time format 1 to be changed
#' @param form2 new date/time format
#' @keywords format_time
#' @export
#' @examples
#' format_time()

format_time<-function(dttm,form1,form2="%Y-%m-%d %H:%M",tz="GMT"){
  strftime(strptime(dttm,format=form1,tz=tz),format=form2,tz=tz)
}
####SUMMARY table
#-------------------------
#' Quick table
#'
#' This function allows you to create quick and dirty table.
#' @param data original date/time
#' @param var1 var1
#' @param var2 var2
#' @keywords tab
#' @export
#' @examples
#' tab()
tab<-function(data,var1,var2){
  data[,"var1"]<-data[,var1]
  data[,"var2"]<-data[,var2]
  xtabs(~var1+var2,data)
}
##Change class
#-------------------------
#' Change variable class
#'
#' This function allows you to change variable class ("num" or "char").
#' @param data data
#' @param var variable (ex:c("DV","MDV"))
#' @param class class ("char" or "num")
#' @keywords chclass
#' @export
#' @examples
#' chclass()
chclass<-function(data,var,class="char"){
   for(i in var){
    if (class=="num"){
      data[,i]<-as.numeric(as.character(data[,i]))}
    else {data[,i]<-as.character(data[,i])}
  }
  data
}
#print unique variable only
#-------------------------
#' one
#'
#' one.
#' @param data data
#' @param var variable
#' @keywords one
#' @export
#' @examples
#' one()
one<-function(data,var){
  for(i in var){
    print(i)
    print(data[!duplicated(data[,i]),i])
  }
}
# keep unique only
#-------------------------
#' No duplicate
#'
#' This function allows you to remove duplicates.
#' @param data data
#' @param var variable (ex:c("DV","MDV"))
#' @param all if all="all", all columns in data will be kept (ex:all=c("ID","DV"))
#' @keywords nodup
#' @export
#' @examples
#' nodup()
nodup<-function(data,var,all,item){
  if(all=="all"){d1<-data[!duplicated(data[,var]),names(data)]}else{
    if(all=="var"){d1<-data[!duplicated(data[,var]),var]}else{
      d1<-data[!duplicated(data[,var]),c(var,item)]}}
  d1
}
#Output duplicated row for checking or remove duplicates if remove is set to non-NULL
#-------------------------
#' Check duplicates
#'
#' This function allows you to check duplicates.
#' @param data data
#' @param var variable (ex:c("DV","MDV"))
#' @param remove if remove="yes", duplicates will be removed)
#' @keywords duprow
#' @export
#' @examples
#' duprow()
duprow<-function(data,var=NULL,remove=NULL){
  flag="flag"
  data[,flag]<-""
  if(is.null(var)){
    var=names(data)}
  for(i in 1:length(var)){
    data[,flag]<-paste(data[,flag],data[,var[i]],sep="")
  }
  if(is.null(remove)){
    data[duplicated(data[,"flag"]),]}
  else{data1<-data[!duplicated(data[,"flag"]),]
       data1[,"flag"]<-NULL
       data1}
}


########Compute Rtime and tad#########
#-------------------------
#' Derive TAD and RTIME
#'
#' This function allows you to derive TAD and RTIME from calendar date/time.
#' @param data data
#' @param id subject id
#' @param date date variable "%Y-%m-%d"
#' @param time time variable  "%H:%M")
#' @param EVID evid variable (evid>0 for dose)
#' @keywords tadRT
#' @export
#' @examples
#' tadRT()
tadRT<-function (data, id, date, time, EVID, tz = "UTC") 
{
  locf <- function(x) {
    good <- !is.na(x)
    positions <- seq(length(x))
    good.positions <- good * positions
    last.good.position <- cummax(good.positions)
    last.good.position[last.good.position == 0] <- NA
    x[last.good.position]
  }
  data$DTTM <- data$TAD <- data$RTIME <- NULL
  data$DTTM <- as.character(paste(data[, date], data[, time], 
                                  sep = " "))
  data <- chclass(data, c(date, time), "char")
  data$tadtm <- NA
  data <- data[order(data[, id], as.Date(data[, date]), data[, 
                                                             time]), ]
  head(data)
  dtm <- data[data[, EVID] > 0, ]
  rtime <- dtm[!duplicated(dtm[, id]), c(id, "DTTM")]
  names(rtime)[2] <- "FDDTM"
  nodose <- data[data[, EVID] == 0, ]
  dose <- data[data[, EVID] > 0, ]
  dose$tadtm <- as.character(dose$DTTM)
  data <- rbind(dose, nodose)
  data$tadtm <- as.character(data$tadtm)
  head(data)
  data$DTTM <- strftime(strptime(data$DTTM, format = "%Y-%m-%d %H:%M", 
                                 tz = tz), format = "%Y-%m-%d %H:%M", tz = tz)
  data <- data[order(data[, id], data$DTTM), ]
  data$WT1 <- unlist(tapply(data$tadtm, data[, id], locf))
  data$tadtm <- rev(locf(rev(data$WT1)))
  data <- data[order(data[, id], as.Date(data[, date]), data[, 
                                                             time]), ]
  head(data)
  data$DTTM <- strftime(strptime(data$DTTM, format = "%Y-%m-%d %H:%M", 
                                 tz = tz), format = "%Y-%m-%d %H:%M", tz = tz)
  data$tadtm <- strftime(strptime(data$tadtm, format = "%Y-%m-%d %H:%M", 
                                  tz = tz), format = "%Y-%m-%d %H:%M", tz = tz)
  data$TAD <- as.numeric(difftime(strptime(data$tadtm, format = "%Y-%m-%d %H:%M", 
                                           tz = tz), strptime(data$DTTM, format = "%Y-%m-%d %H:%M", 
                                                              tz = tz), units = "hour")) * (-1)
  data <- merge(data, rtime, all.x = T)
  data$RTIME <- as.numeric(difftime(strptime(data$DTTM, format = "%Y-%m-%d %H:%M", 
                                             tz = tz), strptime(data$FDDTM, format = "%Y-%m-%d %H:%M", 
                                                                tz = tz), units = "hour"))
  data$WT1 <- NULL
  data$tadtm <- NULL
  data$FDDTM <- NULL
  data <- data[order(data[, id], as.Date(data[, date]), data[, 
                                                             time]), ]
  data$RTIME <- round(data$RTIME, 4)
  data$TAD <- round(data$TAD, 4)
  data
}


#' LOCF and LOCB
#'
#' LOCF LOCB function
#' @param data data
#' @param var variable to locf
#' @param by sort variable
#' @param locb carry backward
#' @keywords locb2
#' @export
#' @examples
#' locf2()
locf2<-function (data, by, var, locb = T) 
{
  locf <- function(x) {
    good <- !is.na(x)
    positions <- seq(length(x))
    good.positions <- good * positions
    last.good.position <- cummax(good.positions)
    last.good.position[last.good.position == 0] <- NA
    x[last.good.position]
  }
  data$dumy<-seq(1,nrow(data))
  data[, var] <- unlist(tapply(data[, var], data[, by], 
                               locf))
  
  if (locb) {
    data <- data[order(data[,by],rev(data$dumy)),]
    data[, var] <- unlist(tapply(data[, var], data[, by], 
                                 locf))
    data <- data[order(data[,by],data$dumy),]
  }
  data$dumy<-NULL
  data
}


########
# 1 cpt
#' One compartment micro constants and HL
#'
#' This function allows you to derive TAD and RTIME from calendar date/time.
#' @param data data
#' @keywords hl1cpt
#' @export
#' @examples
#' hl1cpt()
hl1cpt<-function(data,cl,v,output){
all<-c("HL","k")
  ifelse(output=="all",output<-all,output)

  k<-data[,cl]/data[,v]
  HL<-log(2)/k
  datf<-data.frame(k=k,HL=HL)
  data[,output]<-datf[,output]
  data
  }

#df<-data.frame(id=1:5,cla=1:5/2,v=2:6*4,cl2=3:7/20,v2=3:7*100,cl3=3:7/10,v3=3:7*50)
#df<-hl3cpt(df,"cla","cl2","cl3","v","v2","v3","all")


#Two-compartment
#' Two compartment micro constants and HL
#'
#' This function allows you to derive micro constants and HL.
#' @param data data
#' @keywords hl2cpt
#' @export
#' @examples
#' hl2cpt(pkdat1,"cl","cl2","v","v2","all")
#
hl2cpt<-function(data,cl,cl2,v,v2,output){
  all<-c("HLa","HLb","alfa","beta","k","k12","k21")
  ifelse(output=="all",output<-all,output)
df<-data
  k<-df[,cl]/df[,v]
  k12<-df[,cl2]/df[,v]
  k21<-df[,cl2]/df[,v2]
beta1<-(1/2)*(k12+k21+k-(sqrt((k12+k21+k)^2-(4*k21*k))))
alfa<-k21*k/beta1
alfaHL<-log(2)/alfa    # to be verify with excel
betaHL<-log(2)/beta1    # to be verified with excel
datf<-data.frame(k=k,k12=k12,k21=k21,alfa=alfa,beta=beta1,HLa=alfaHL,HLb=betaHL)
  data[,output]<-datf[,output]
  data
}

# Three CPT
#' Three compartment micro constants and HL
#'
#' This function allows you to derive micro constants and HL.
#' @param data data
#' @keywords hl3cpt
#' @export
#' @examples
#' hl3cpt(pkdat1,"cl","cl2",,"cl3","v","v2","v3","all")
#
hl3cpt<-function(data,Cl,Cl2,Cl3,V,V2,V3,output){
  all<-c("HLa","HLb","HLg","A","B","C","alpha","beta","gama")
  ifelse(output=="all",output<-all,output)
  df<-data
  k<-df[,Cl]/df[,V]
  k12<-df[,Cl2]/df[,V]
  k21<-df[,V]*k12/df[,V2]
  k13<-df[,Cl3]/df[,V]
  k31<-df[,V]*k13/df[,V3]
  a0<-k*k21*k31
  a1<-(k*k31) + (k21*k31) + (k21*k13) + (k*k21) + (k31*k12)
  a2<-k + k12 + k13 + k21 + k31
  p<-a1 - (a2^2)/3
  q<-2*(a2^3)/27 - a1*a2/3 + a0
  r1<-sqrt((-1)*p^3/27)
  r2<-2*(r1^(1/3))
  phi<-acos(-1*q/(2*r1))/3
  gama<-(-1)*((cos(phi)*r2)-(a2/3))             # gama instead of alpha: Formula error found in Ref. Dubois A. et al., "Mathematical Expressions of the Pharmacokinetic and Pharmacodynamic Models implemented in the PFIM"
  alpha<-(-1)*(cos(phi+(2*pi/3))*r2-a2/3)       # alpha instead of beta: Formula error found in Ref. Dubois A. et al., "Mathematical Expressions of the Pharmacokinetic and Pharmacodynamic Models implemented in the PFIM"
  beta<-(-1)*(cos(phi+(4*pi/3))*r2-a2/3)        # beta instead of gamma: Formula error found in Ref. Dubois A. et al., "Mathematical Expressions of the Pharmacokinetic and Pharmacodynamic Models implemented in the PFIM"
  alfaHL<-log(2)/alpha
  betaHL<-log(2)/beta
  gamaHL<-log(2)/gama
  A=(1/df[,V])*((k21-alpha)/(alpha-beta))*((k31-alpha)/(alpha-gama))
  B=(1/df[,V])*((k21-beta)/(beta-alpha))*((k31-beta)/(beta-gama))
  C=(1/df[,V])*((k21-gama)/(gama-beta))*((k31-gama)/(gama-alpha))
  datf<-data.frame(HLa=alfaHL,HLb=betaHL,HLg=gamaHL,alpha=alpha,beta=beta,gama=gama,A=A,B=B,C=C)
  data[,output]<-datf[,output]
  data
}


#Round old method
#' Rounding as per Excel Internal use
#'
#' This function allows you to round value as per Excel method.
#' @keywords cround
#' @export
#' @examples
#' cround1()
cround1= function(x,n,asnum=T){
  vorz = sign(x)
  z = abs(x)*10^n
  z = z + 0.5
  z = trunc(z)
  z = z/10^n
ifelse(is.na(x),output<-NA,
  output<-sprintf(paste("%.",n,"f",sep=""),z*vorz))
if(asnum){output<-as.numeric(as.character(output))}else{
output}
output
}

#' Round up 
#'
#' This function allows you to round value as in Excel.
#' @param z Vector or single value to be rounded
#' @param y number of significant figure 
#' @keywords rounding
#' @export
#' @examples
#' cround()
cround<-function (z, y) 
{
  if(length(z)>1){
    output<-NULL
    for(i in 1:length(z)){
      output1<-cround1(as.numeric(z[i]),y)
      output<-rbind(output,output1)}
    output
  }else{output<-cround1(as.numeric(z),y)
  output
  }}


#sigfig Internal Use
#' Significant figure
#'
#' This function allows you to round value in significant figure.
#' @keywords sigfig
#' @export
#' @examples
#' sigfig1()
sigfig1<-function (x, y) 
{
  sround = function(x, n) {
    vorz = sign(x)
    z = abs(x) * 10^n
    z = z + 0.5
    z = trunc(z)
    z = z/10^n
    ifelse(is.na(x), sro <- NA, sro <- z * vorz)
    sro
  }
  nround <- ifelse(x == 0, y - 1, y - 1 - floor(log10(abs(x))))
  if (!is.na(x) & ceiling(log10(abs(x))) >= 3) {
    output <- as.character(cround(x, 0))
  }else {
    if (!is.na(x) & ceiling(log10(abs(x)))<3) {
      output <- sprintf(paste("%.", nround, "f", sep = ""), 
                        sround(x, nround))
    }else{
      output <- NA
    }
  }
  output
}

#Sigfig 
#' Significant figure
#'
#' This function allows you to round value in significant figure.
#' @param z Vector or single value to be rounded
#' @param y number of significant figure 
#' @keywords sigfig
#' @export
#' @examples
#' sigfig()
sigfig<-function (z, y) 
{
  if(length(z)>1){
    output<-NULL
    for(i in 1:length(z)){
      output1<-sigfig1(as.numeric(z[i]),y)
      output<-rbind(output,output1)}
    output
  }else{output<-sigfig1(as.numeric(z),y)
  output
  }}

#Only output unique duplicate item
#' Filter unique duplicated row
#'
#' This function allows you to filter duplicated rows but only show unique row
#' @param data data
#' @param data data
#' @param all display all columns (all="all")
#' @param select display selected variables only
#' @keywords dup1
#' @export
#' @examples
#' dup1()
dup1<-function(data,var,all,select){
  d1<-data[duplicated(data[,var]),]
  if(all=="all"){d1<-d1}else{
    if(all=="var"){d1<-d1[,var]}else{
      d1<-d1[,c(var,select)]}}
  d1
}

#Output duplicated items with all or partial variabes
#' Filter all duplicated rows
#'
#' This function allows you to filter duplicated rows but only show unique row
#' @param data data
#' @param data data
#' @param all display all columns (all="all")
#' @param select display selected variables only
#' @keywords dup2
#' @export
#' @examples
#' dup2()
dup2<-function(data,var,all,select){
  d1<-data
  d1$dum<-""
  for(i in var){
    d1$dum<-paste(d1$dum,d1[,i],sep="-")
  }
  dup<-d1[duplicated(d1$dum),"dum"]
  d1<-d1[d1$dum%in%dup,]
  if(all=="all"){d1<-d1[,names(data)]}else{
    if(all=="var"){d1<-d1[,var]}else{
      d1<-d1[,c(var,select)]}}
  d1
}

#TABLE FUNCTIONS###############
#' bround Table function
#' @param data data
#' @keywords bround
#' @export
#' @examples
#' bround()
bround<-function(data,var,rtype="sigfig",dec=3){
  data<-chclass(data,var,"num")
  for(i in var){
    data[is.na(data[,i]),i]<-9999999999999
    if(rtype=="sigfig"){data[,i]<-sigfig(data[,i],dec)}else{data[,i]<-cround(data[,i],dec)}
    data[data[,i]=="9999999999999",i]<-"NA"
  }
  data
}

#' geom Table function
#'
#' @param x data
#' @keywords geom
#' @export
#' @examples
#' geom()


geom <- function(x) {
  exp(mean(log(x[x > 0]), na.rm=TRUE))
}

#' geocv Table function
#' @param x data
#' @keywords geocv
#' @export
#' @examples
#' geocv()

geocv <- function(x) {
  100*sqrt(exp(var(log(x[x > 0]), na.rm=TRUE)) - 1)
}

#' cv Table function
#' @param x data
#' @keywords cv
#' @export
#' @examples
#' cv()

cv <- function(x) {
  abs(sd(x,na.rm=TRUE)/mean(x,na.rm=TRUE)*100)
}

#########
# mean<-"mean(x,na.rm=T)=mean"
# sd<-"sd(x,na.rm=T)=sd"
# cv<-"cv(x)=cv"
# qt05<-"quantile(x,0.05,na.rm=TRUE)=qt05"
# qt95<-"quantile(x,0.95,na.rm=TRUE)=qt05"

#per95<-function(x){quantile(x,percentile,na.rm=TRUE)}

#' se Table function
#' internal use.
#' @param x data
#' @keywords se
#' @export
#' @examples
#' se()

se<-function(x){sd(x,na.rm=TRUE)/(length(x))^0.5}

#' cilow Table function
#'
#' internal use
#' @param x data
#' @keywords generic
#' @export
#' @examples
#' cilow()

cilow<-function(x){mean(x,na.rm=TRUE)-((sd(x,na.rm=TRUE)/(length(x))^0.5)*qt(0.975,df=length(x)-1))}    #1.96)}

#' ciup Table function
#' internal use.
#' @param x data
#' @keywords ciup
#' @export
#' @examples
#' ciup()
ciup<-function(x){mean(x,na.rm=TRUE)+((sd(x,na.rm=TRUE)/(length(x))^0.5)*qt(0.975,df=length(x)-1))}

#' nmiss Table function
#' internal use.
#' @param x data
#' @keywords nmiss
#' @export
#' @examples
#' nmiss()
nmiss<-function(x){length(x[is.na(x)])}

#nmiss<-function(x){length(x[is.na(x)])}
#upci<-function(x){mean(x,na.rm=TRUE)+((sd(x,na.rm=TRUE)/(length(x#))^0.5)*qt(0.975,df=length(x)-1))}
#loci<-function(x){mean(x,na.rm=TRUE)-((sd(x,na.rm=TRUE)/(length(x#))^0.5)*qt(0.975,df=length(x)-1))}    #1.96)}
#se<-function(x){sd(x,na.rm=TRUE)/(length(x))^0.5}
#per95<-function(x){quantile(x,percentile,na.rm=TRUE)}
#per05<-function(x){quantile(x,1-percentile,na.rm=TRUE)}
#Gmean <- function(x) {
 # exp(mean(log(x[x > 0]), na.rm=TRUE))
#}
# Gcv <- function(x) {
#   100*sqrt(exp(var(log(x[x > 0]), na.rm=TRUE)) - 1)
# }
# 
# }
#####################TABLE STATS##################

#' conti<-function(input=input,var=var,by=by,round.type="sigfig",digit=3,quanti){
#'   statsfun()
#'   cv<-function(x){
#'     abs(sd(x,na.rm=TRUE)/mean(x,na.rm=TRUE)*100)
#'   }
#'   nmiss<-function(x){length(x[is.na(x)])}
#'   upci<-function(x){mean(x,na.rm=TRUE)+((sd(x,na.rm=TRUE)/(length(x))^0.5)*qt(0.975,df=length(x)-1))}
#'   loci<-function(x){mean(x,na.rm=TRUE)-((sd(x,na.rm=TRUE)/(length(x))^0.5)*qt(0.975,df=length(x)-1))}    #1.96)}
#'   se<-function(x){sd(x,na.rm=TRUE)/(length(x))^0.5}
#'   per95<-function(x){quantile(x,percentile,na.rm=TRUE)}
#'   per05<-function(x){quantile(x,1-percentile,na.rm=TRUE)}
#'   Gmean <- function(x) {
#'     exp(mean(log(x[x > 0]), na.rm=TRUE))
#'   }
#'   Gcv <- function(x) {
#'     100*sqrt(exp(var(log(x[x > 0]), na.rm=TRUE)) - 1)
#'   }
#' 
#'   input<-input[,c(var,by)]
#'   l<-stats::reshape(input,
#'              varying = var,
#'              v.names = "score",
#'              timevar = "var",
#'              times = var,
#'              direction = "long")
#'   l1<-l[!is.na(l$score),]
#'   l2<-l[is.na(l$score),]
#'   l1$qt<-quanti
#'   sum<-plyr::ddply(l1,c(by,"var"),summarise,
#'              n=length(score),
#'              #nmiss=nmiss(score),
#'              mean=mean(score,na.rm=T),
#'              cv=cv(score),
#'              sd=sd(score,na.rm=T),
#'              se=se(score),
#'              median=median(score,na.rm=T),
#'              min=min(score,na.rm=T,na.rm=T),
#'              max=max(score),
#'              lo_qt975=quantile(score,0.025,na.rm=T),
#'              hi_qt975=quantile(score,0.975,na.rm=T),
#'              lo_qt95=quantile(score,0.05,na.rm=T),
#'              hi_qt95=quantile(score,0.95,na.rm=T),
#'              lo_qtxx=quantile(score,1-unique(qt),na.rm=T),
#'              hi_qtxx=quantile(score,unique(qt),na.rm=T),
#'              lowCI=loci(score),
#'              HiCI=upci(score),
#'              GeoMean=Gmean(score),
#'              GeoCV=Gcv(score)
#'   )
#' 
#'   variable=c("mean","cv","sd","se","median","min","GeoMean", "GeoCV",
#'              "max","lo_qt975", "hi_qt975", "lo_qt95",  "hi_qt95","lo_qtxx","hi_qtxx","lowCI","HiCI")
#' 
#'   l<-reshape(sum,
#'              varying = variable,
#'              v.names = "value",
#'              timevar = "toround",
#'              times = variable,
#'              direction = "long")
#'   l<-l[!is.na(l$value),]
#'   if(round.type=="sigfig"){
#'     l$value<-sigfig(l$value,digit)}else{l$value<-cround(l$value,digit)}
#'   l$id<-NULL
#'   keep<-names(l)[!names(l)%in%c("toround","value")]
#'   w <- stats::reshape(l,
#'                timevar = "toround",
#'                idvar = keep,
#'                direction = "wide")
#'   names(w)<-gsub("value.","",names(w))
#'   if(nrow(l2)>0){
#'     sum1<-plyr::ddply(l2,c(by,"var"),summarise,
#'                 nmiss=nmiss(score))
#'     w<-plyr::join(w,sum1)
#'     w$nmiss[is.na(w$nmiss)]<-0}else{w$nmiss=0}
#'   w$score<-w$id<-NULL
#'   w
#' }


#' Funtion for Descriptove stats of continuous covariate
#'
#'
#' Descriptove stats of continuous covariates
#' @param data datset or data frame (ex:data=PKdatat)
#' @param var List of continuous covariates (ex:c("CRCL","WT"))
#' @param by  Stratification variable (ex: by="study")
#' @param colby Specify sorting variable to be displayed vertically. (ex: colby=by or colby="var")
#' @param rowby Specidy sorting variable horizontally. If by > 1, set rowby=by.
#' @param sumry Result summary format to be displayed. Type summary.set to see different options of summary sets.
#' @param round.type rounding method, sigfig by default. (ex:"round" for rounding)
#' @param digit Number of round digits or significant figures
#' @keywords con.tab
#' @export
#' @examples
#' con.tab(data=dat,var=c("WT","crcl","ALT"),by=c("study"),colby="var",rowby=by,sumry=set1)

lhcontab<-function(dat,var,by,format="by.var",sumry=1,round.type="sigfig",digit=3,quant=0.90){
  statsfun()
  cv<-function(x){
    abs(sd(x,na.rm=TRUE)/mean(x,na.rm=TRUE)*100)
  }
  nmiss<-function(x){length(x[is.na(x)])}
  upci<-function(x){mean(x,na.rm=TRUE)+((sd(x,na.rm=TRUE)/(length(x))^0.5)*qt(0.975,df=length(x)-1))}
  loci<-function(x){mean(x,na.rm=TRUE)-((sd(x,na.rm=TRUE)/(length(x))^0.5)*qt(0.975,df=length(x)-1))}    #1.96)}
  se<-function(x){sd(x,na.rm=TRUE)/(length(x))^0.5}
  per95<-function(x){quantile(x,percentile,na.rm=TRUE)}
  per05<-function(x){quantile(x,1-percentile,na.rm=TRUE)}
  Gmean <- function(x) {
    exp(mean(log(x[x > 0]), na.rm=TRUE))
  }
  Gcv <- function(x) {
    100*sqrt(exp(var(log(x[x > 0]), na.rm=TRUE)) - 1)
  }

  p1<-"(";p2<-")"
  b1<-"[";  b2<-"]"
  a<-"{"; a<-"}"
  d<-"-"; com<-","
  bsl<-"/";lb<-"\n"
  sp<-" "
  set=NULL
  set[[1]]=c("paste0(mean,sp,p1,cv,p2,lb,median,sp,b1,min,d,max,b2)","Mean(CV)Med[min-max]")
  set[[2]]=c("paste0(mean,sp,p1,sd,p2,lb,median,sp.b1,min,d,max,b2)","Mean(SD)Med[min-max]")
  set[[3]]=c("paste0(mean,sp,p1,cv,com,sp,n,p2,lb,median,sp,b1,min,d,max,b2)","Mean(CV,N) Med[min-max]")
  set[[4]]=c("paste0(mean,sp,p1,sd,com,sp,n,p2,lb,median,sp,b1,min,d,max,b2)","Mean(SD,N) Med[min-max]")
  set[[5]]=c("paste0(n,sp,p1,nmiss,p2,lb,mean,sp,p1,cv,p2,lb,median,sp,b1,min,d,max,b2)","N(Nmis)Mean(CV) Med[min-max]")
  set[[6]]=c("paste0(n,sp,p1,nmiss,p2,lb,mean,sp,p1,sd,p2,lb,median,sp,b1,min,d,max,b2)","N(Nmis)Mean(SD) Med[min-max]")
  set[[7]]=c("paste0(mean,sp,p1,cv,p2,lb,median,sp,b1,lowCI,d,HiCI,b2)","Mean(CV)Med[95CI]")
  set[[8]]= c("paste0(mean,sp,p1,cv,p2,lb,median,sp,b1,lo_qt95,d,hi_qt95,b2)","Mean(CV)Med[95PI]")
set[[9]]=c("paste0(mean,sp,p1,cv,p2,lb,median,sp,b1,lo_qt975,d,hi_qt975,b2)","Mean(CV)Med[97.5PI]")
set[[10]]=c("paste0(mean,sp,p1,cv,p2,lb,median,sp,b1,lo_qtxx,d,hi_qtxx,b2)","Mean(CV) Median[xxPI]")



  data<-dat[,c(var,by)]
    data<-data[,c(var,by)]
    data<-chclass(data,var,"num")
  #  lo_qtxx<-paste0("loqt",quant)
  #  hi_qtxx<-paste0("hiqt",quant)

    con1<-conti(data,var=var,by=by,round.type=round.type,digit=digit,quanti=quant)
    data1<-data
    head(data1)
    data1[,by]<-"Overall"
    con2<-conti(data1,var=var,by=by,round.type=round.type,digit=digit,quanti=quant)
    con1<-rbind(con1,con2)
   # names(con1)[names(con1)=="lo_qtxx"]<-lo_qtxx
  #  names(con1)[names(con1)=="hi_qtxx"]<-hi_qtxx
    grp<-sumry

    head(con1)
    if(is.null(grp)){w<-con1}else{
      groupstats<-set[[grp]][1]
      con1$summary<-with(con1,eval(parse(text=groupstats)))
      con1<-con1[,c(by,"var","summary")]
      if(format=="by.var"){
        w <- reshape(con1,
                            timevar = "var",
                            idvar = by,
                            direction = "wide")
        names(w)<-gsub("summary.","",names(w))
        nm<-names(w)
        w$stats<-set[[grp]][2]
        w<-w[,c("stats",nm)]
        }else{
            tx<-NULL
            for(i in var){
            ty<-t(con1[con1$var==i,c("var","summary")])
            tx<-rbind(tx,ty)}
            tx<-tx[row.names(tx)!="var",]
            row.names(tx)<-var
            dim(tx)
            con2<-chclass(con1,by,"char")
            byby<-nodup(con2,by,"var")
            t1<-t(byby)
            dim(t1)
            row.names(t1)<-by
            w<-as.data.frame(rbind(t1,tx))
            names(w)
            names(w)<-c(paste0(rep(set[[grp]][2],dim(w)[2]-1),1:dim(w)[2]-1))
            nm<-names(w)
           w$Covariate<-row.names(w)
           w<-w[,c("Covariate",nm)]
        }}
    row.names(w)<-NULL
    w
      }

################################################
#' roundbatch
#'
#' internal use
#' @keywords roundbatch
#' @export
#' @examples

roundbatch<-function(data,variable,toround,nb){
  head(data)
  data<-sum
  l<-stats::reshape(data,
             varying = variable,
             v.names = "value",
             timevar = "toround",
             times = variable,
             direction = "long")
  l<-l[!is.na(l$value),]
  if(toround=="sigfig"){
    l$value<-sigfig(l$value,nb)}else{l$value<-cround(l$value,nb)}
  l$id<-NULL
  keep<-names(l)[!names(l)%in%c("toround","value")]
  w <- stats::reshape(l,
               timevar = "toround",
               idvar = keep,
               direction = "wide")
  names(w)<-gsub("value.","",names(w))
  w
}


################COUNTS CATEGORICAL###############

#' Funtion for Descriptove Stats of Categorical Covariate
#'
#'
#' Descriptove stats of categorical covariates
#' cat.tab(data,var,by,colby="var",rowby=by)
#' @param data datset or data frame (ex:data=PKdatat)
#' @param var List of continuous covariates (ex:c("SEX","RACE"))
#' @param by  Stratification variable (ex: by="study")
#' @keywords cat.tab
#' @export
#' @examples
#' cat.tab(data=dat,var=c("SEX","RACE"),by=c("study"),colby="var",rowby=by)

lhcattab<-function (data, var, by)
  {
  rowby = by
  dat1 <- chclass(data[, c(var, by)], c(var, by), "char")

  tot <- stats::reshape(dat1, varying = var, v.names = "value",
                        timevar = "var", times = var, direction = "long")
  tot$id <- NULL
  tot1 <- addvar(tot, c(by, "var"), "var", "length(x)", "no",
                 "tot")
  tot2 <- addvar(tot, c(by, "var", "value"), "var", "length(x)",
                 "no", "subt")
  tot11 <- nodup(tot1, c(by), "all")
  names(tot11)[names(tot11) == "tot"] <- "N="
  tot12 <- addvar(tot, c("var", "value"), "var", "length(x)",
                  "no", "Overall")
  tot13 <- addvar(tot, c("var"), "var", "length(x)", "no",
                  "tot")
  tot12 <- plyr::join(tot12, tot13)
  tot12$Overall <- with(tot12, paste0(Overall, " (", sigfig(Overall/tot *
                                                              100, 3), "%)"))
  tot12$tot <- NULL
  tot11$"N=" <- paste0(tot11$"N=", " (", sigfig(tot11$"N="/max(tot13$tot) *
                                                  100, 3), "%)")
  tot4 <- plyr::join(tot2, tot1)
  tot4$summary <- with(tot4, paste0(subt, " (", sigfig(subt/tot *
                                                         100, 3), "%)"))
  tot3 <- tot4[, c(by, "var", "value", "summary")]

  tto<-addvar(tot4,c(rowby),"tot","max(x)","yes","all")
  tto[,c("var","value")]<-"all"
  tto$all<-paste0(tto$all," (100%)")

  tto0<-tto;tto0$tot<-tto0$subt<-tto0$summary<-NULL
  w0 <- stats::reshape(tto0, timevar =rowby, idvar = c("var",
                                                       "value"), direction = "wide")
  names(w0)<-gsub("all.","",names(w0))

  tto1<-tot4;tto1$tot<-tto1$subt<-NULL
  w <- stats::reshape(tto1, timevar =rowby, idvar = c("var",
                                                      "value"), direction = "wide")
  names(w)<-gsub("summary.","",names(w))

  w<-rbind(w,w0)

  w1<-addvar(tot4,c("var","value"),"subt","sum(x)","no","overall")
  w2<-addvar(tot4,c("var"),"subt","sum(x)","no","overall1")
  w1<-plyr::join(w1,w2)

  w1$overall<-paste0(w1$overall," (",sigfig(w1$overall/w1$overall1*100,3),"%)")
  w1a<-w1[1,]
  w1a$var<-w1a$value<-"all"
  w1a$overall<-paste0(w1a$overall1," (100%)")
  w1<-rbind(w1,w1a)

  w<-plyr::join(w,w1[,c("var",  "value","overall")])
  w<-w[order(w$var),]
  w
}



#' Individual table with descriptive statse
#'
#' Listing of individual data and descriptove stats
#' @param data datset or data frame (ex:data=PKdatat)
#' @param id unique identifier
#' @param by  Stratification variable (ex: by="study")
#' @param variables Specify sorting variable to be displayed vertically. (ex: colby=by or colby="var")
#' @param rtype rounding type. (sigfig by default)
#' @param dec round decimal or number of significant figures
#' @keywords ind.tab
#' @export
#' @examples
#' ind.tab(data=dat,id="NMID",by=c("study"))
indiv.tab<-function(data,id,by,variables,rtype="sigfig",dec=3){
  id<-id#
  data<-data[,c(id,by,variables)]#[!duplicated(data$id),]
  strat1<-by#c("phase")# mandatory
  convar<-variables #mandatory
  d1<-data[,c(id,strat1,convar)]
  d1<-chclass(d1,convar,"num")
  head(d1)

  t1<-NULL
  for(i in unique(d1[,strat1])){
    d0<-d1[d1[,strat1]%in%i,]
    l<-stats::reshape(d0,
               varying = c(convar),
               v.names = "score",
               timevar = "subj",
               times = c(convar),
               #new.row.names = 1:1000,
               direction = "long")
    head(l)
    l$id<-NULL
    str(l)
    st<-plyr::ddply(l,c(by,"subj"),summarise,
              N=round(length(score),0),
              Nmiss=round(length(score[is.na(score)]),0),
              Means=sigfig(mean(score,na.rm=T),3),
              SD=sigfig(sd(score,na.rm=T),3),
              cv=sigfig(cv(score),3),
              Median=sigfig(median(score,na.rm=T),3),
              Minimum=sigfig(min(score,na.rm=T),3),
              Maximum=sigfig(max(score,na.rm=T),3),
              GeoMean=sigfig(Gmean(score),3),
              GeoCV=sigfig(Gcv(score),3))
    keep<-names(st[,3:length(names(st))])
    l1<-stats::reshape(st,
                varying = c(keep),
                v.names = "Stats",
                timevar = "Results",
                times = c(keep),
                #new.row.names = 1:1000,
                direction = "long")
    l1$id<-NULL

    w<-stats::reshape(l1,
               timevar = "subj",
               idvar = c(strat1, "Results"),
               direction = "wide")
    names(w)<-gsub("Stats.","",names(w))
    head(d0)
    x1<-setdiff(names(d0),names(w))
    x2<-setdiff(names(w),names(d0))
    w[,x1]<-""
    d0[,x2]<-""
    d0<-d0[,c(id,strat1,x2,convar)]
    #d0<-chclass(d0,convar,"num")
    if(!is.null(rtype)){
      d0<-bround(d0,convar,rtype=rtype,dec=dec)}
    t<-rbind(d0,w)
    t<-t[,c(id,strat1,x2,convar)]
    t1<-rbind(t1,t)
  }
  t1
}

#' Calculate AUC Using the Trapezoidal Method
#'
#' Calculate the area under the curve (AUC) for each subject over the time interval for dv using the trapezoidal rule.
#' nca.analysis(data,n_lambda=3,id,time,dv,ss,partialAUC="no")
#' data
#' @param data.frame containing the data to use for the AUC calculation
#' @param time	chronologically ordered time variable present in data
#' @param id	variable in data defining subject level data
#' @param dv	dependent variable used to calculate AUC present in data
#' @keywords AUC
#' @export
#' @examples
#' AUC(data, time = 'TIME', id = 'ID', dv = 'DV')

AUC<-function (data, time = "TIME", id = "ID", dv = "DV")
{
  if (any(is.na(data[[id]])))
    warning("id contains NA")
  if (any(is.na(data[[time]])))
    warning("time contains NA")
  if (any(is.na(data[[dv]])))
    warning("dv contains NA")
  data <- data[order(data[[id]], -data[[time]]), ]
  nrec <- length(data[[time]])
  data$diff <- c(data[[time]][-nrec] - data[[time]][-1], 0)
  data$meanDV <- c((data[[dv]][-1] + data[[dv]][-nrec])/2,
                   0)
  data$dAUC <- data$diff * data$meanDV
  data <- data[order(data[[id]], data[[time]]), ]
  data <- data[duplicated(data[[id]]), ]
  AUC <- aggregate.data.frame(data$dAUC, by = list(data[[id]]),
                              FUN = sum)
  names(AUC) <- c(id, "AUC")
  return(AUC)
}

#' Derive Common NCA parameters using single and multiple profiles
#'
#' nca.cal()
#' @param data datset or data frame (ex:data=PKdatat)
#' @param id unique subject identifier
#' @param n_lambda  number of points for estimating the Lambda
#' @param time Sampling time after dose (TAD)
#' @param dv Concentration
#' @param dosing.regimen PK profile profile flag (ex: single, sd, md ,ss).As vector or numeric. Note that only one prefile per dose is to be analyzed.
#' @param label.for.sd Label to identify sd profile used in dosing.regimen. Need to compute the effective half-life
#' @param label.for.ss Label to identify ss profile used in dosing.regimen. Need to compute the effective half-life
#' @param tau Dose interval. As vector or numeric.
#' @param dose Dose administered. As vector or numeric. Required to estimate CL and Vss 
#' @param partialAUC Time interval for partial AUC. Ex: c(0,6,0,12,6,12) for AUC0-6, AUC0-12 and AUC6-12
#' @param partialConc Concentration at particular time (Ex:c(1,4) for concentration after 1 and 4 h)
#' @keywords nca.cal
#' @export
#' @examples 
#' p<-nca.cal(data=data, n_lambda = 3, id = "id", time = "time", dv = "dv", 
#' tau="ii",dose="dose",dosing.regimen = "ss",label.for.sd="single", partialAUC = NULL, partialConc = NULL)
#' 

nca.cal<-function (data, n_lambda = 3, id = "id", time = "time", dv = "dv",sort.by=NULL, ss.variable = "ss", sd_label_in_ss = NULL, ss_label_in_ss = NULL,tau = NULL, dose = NULL, partialAUC = NULL, partialConc = NULL) 
{
  dat <- data
  dat$time1 <- dat[, time]
  dat$time <- dat[, time]
  dat[, "tad"] <- NULL
  dat$id <- dat[, id]
  dat$dv <- dat[, dv]
  if (is.numeric(ss.variable)) {
    dat$ss <- ss.variable
  } else {
    dat$ss <- dat[,ss.variable]
  }
  if (is.numeric(dose)) {
    dat$dose <- dose
  } else {
    if (!is.null(dose)) {
      dat$dose <- dat[, dose]
    } else {
      dat$dose <- "dose required"
    }
  }
  if (is.numeric(tau) & !is.null(tau)) {
    dat$tau <- tau
  }else{
    if (!is.null(tau)) {
      dat$tau <- dat[, tau]
    }else {dat}
  }
  dat$idss <- paste(dat$id, dat$ss, sep = "-")
  if(!is.null(sort.by)){
    for(i in sort.by){
      dat$idss <- paste(dat$idss, dat[,i], sep = "-")  
    }}else{dat$idss<-dat$idss}
  
  dat <- dat[order(dat$idss, dat$time), ]
  dat <- addvar(dat, "idss", "time", "min(x)", "yes", "tad")
  dat$time <- dat$time - dat$tad
  idss <- nodup(dat, c("idss", "id", "ss", "dose",sort.by), "var")
  dat$dvtm <- dat$dv * dat$time
  datauc <- dat
  auclast <- AUC(datauc, time = time, id = "idss", dv = dv)
  names(auclast) <- c("idss", "AUClast")
  auclast <- plyr::join(auclast, idss)
  aucmlast <- AUC(datauc, time = time, id = "idss", dv = "dvtm")
  names(aucmlast) <- c("idss", "AUMClast")
  aucmlast <- plyr::join(aucmlast, idss)
  head(dat)
  dat$tad1 <- dat$time
  aucpart <- NULL
  if (!is.null(partialAUC)) {
    nauc <- length(partialAUC)/2
    for (z in seq(1, length(partialAUC), 2)) {
      tm1 <- partialAUC[z]
      tm2 <- partialAUC[z + 1]
      auc <- AUC(dat[dat[, "tad1"] >= tm1 & dat[, "tad1"] <= 
                       tm2, ], time = "tad1", id = "idss", dv = dv)
      names(auc) <- c("idss", paste0("AUC", tm1, "-", tm2))
      if (z == 1) {
        aucpart <- rbind(aucpart, auc)
      } else {
        aucpart[, paste0("AUC", tm1, "-", tm2)] <- auc[, 
                                                       2]
      }
    }
    aucpart <- join(aucpart, idss)
    aucpart$idss <- NULL
    aucpart
  } else {
    aucpart <- NULL
  }
  Cpart <- NULL
  if (!is.null(partialConc)) {
    nauc <- length(partialConc)
    for (z in 1:length(partialConc)) {
      tm1 <- partialConc[z]
      partc <- dat[dat[, "tad1"] == tm1, c("idss", "dv")]
      names(partc) <- c("idss", paste0("C", tm1))
      if (z == 1) {
        Cpart <- rbind(Cpart, partc)
      } else {
        Cpart[, paste0("C", tm1)] <- partc[, 2]
      }
    }
    Cpart
    Cpart <- join(Cpart, idss)
    Cpart$idss <- NULL
    Cpart
  }else {
    Cpart <- NULL
  }
  if (n_lambda != "no") {
    dat1 <- dat
    dat1$tmp <- seq(nrow(dat1))
    dat1 <- addvar(dat1, "idss", "tmp", "max(x)", "yes", 
                   "tmp2")
    head(dat1)
    dat1$tmp <- dat1$tmp2 - dat1$tmp
    dat1 <- dat1[dat1$tmp < n_lambda, ]
    str(dat1)
    dat1[, c("idss", "time", "dv")]
    test1 <- ddply(dat1[, c("idss", "time", "dv")], .(idss), 
                   summarize, interc = lm(log(dv) ~ time)$coef[1], Lambda = lm(log(dv) ~ time)$coef[2] * -1, R2 = summary(lm(log(dv) ~ time))$r.squared, HL = (log(2)/lm(log(dv) ~ time)$coef[2]) *  -1, that = max(time))
    test1$n_lambda <- n_lambda
    test1$Clast_hat <- with(test1, exp(-Lambda * that + interc))
    test1a <- ddply(dat1[, c("idss", "time", "dvtm")], .(idss), 
                    summarize, intercc = lm(log(dvtm) ~ time)$coef[1], 
                    Lambdac = lm(log(dvtm) ~ time)$coef[2] * -1, 
                    R2c = summary(lm(log(dvtm) ~ time))$r.squared, 
                    HLc = (log(2)/lm(log(dvtm) ~ time)$coef[2]) * -1, 
                    thatc = max(time))
    test1a$n_lambdac <- n_lambda
    test1a$Clast_hatc <- with(test1a, exp(-Lambdac * thatc + 
                                            intercc))
  }else {
    test1 <- NULL
  }
  if (TRUE %in% c(test1$HL < 0)) {
    test1$Warning.HL.Negative = ifelse(test1$HL, "yes", "")
  }
  max <- ddply(dat[, c("idss", "dv", "time", "time1")], .(idss), 
               summarize, Cmax = max(dv), Tmax = time1[dv == max(dv)], 
               Cmin = min(dv[time >= time[dv == max(dv)]]), Tlast = max(dat$time1), 
               Clast = dv[time == max(time)])
  maxa <- ddply(dat, .(idss), summarize, Clastc = dvtm[time == 
                                                         max(time)])
  head(dat)
  test <- plyr::join(max, idss)
  test <- plyr::join(test, maxa)
  test <- plyr::join(test, auclast)
  test <- plyr::join(test, aucmlast)
  if (n_lambda != "no") {
    test <- join(test, test1)
    test <- join(test, test1a)
    test$AUCinf_obs <- abs(as.numeric(as.character(test$AUClast)) + 
                             test$Clast/test$Lambda)
    test$AUMCinf_obs <- abs(as.numeric(as.character(test$AUMClast)) + 
                              test$Clastc/test$Lambdac)
    test$AUCinf_pred <- abs(as.numeric(as.character(test$AUClast)) + 
                              test$Clast_hat/test$Lambda)
    test$AUMCinf_pred <- abs(as.numeric(as.character(test$AUMClast)) + 
                               test$Clast_hatc/test$Lambdac)
    test$MRTlast <- test$AUMClast/test$AUClast
    test$MRTobs <- test$AUMCinf_obs/test$AUCinf_obs
    test$MRTpred <- test$AUMCinf_pred/test$AUCinf_pred
  }else {
    test$MRTlast <- test$AUMClast/test$AUClast
  }
  if (!is.null(Cpart)) {
    test <- plyr::join(test, Cpart)
  }
  if (!is.null(aucpart)) {
    test <- plyr::join(test, aucpart)
  }
  test$idss <- test$interc <- test$that <- NULL
  n <- names(test)
  n <- n[!n %in% c("id", "ss")]
  test <- test[, c("id", "ss", n)]
  if ("AUCinf_pred" %in% names(test) & is.numeric(test$dose)) {
    test$CLobs <- with(test, dose/AUCinf_obs)
    test$CLpred <- with(test, dose/AUCinf_pred)
    test$Vss_obs <- with(test, CLobs * AUMCinf_obs)
    test$Vss_pred <- with(test, CLpred * AUMCinf_pred)
    test <- test[, !names(test) %in% c("Lambdac", "R2c", 
                                       "HLc", "thatc", "n_lambdac", "Clast_hatc", "intercc")]
    names(test)[names(test) == "HL"] <- "HL_Lambda_z"
  }else {
    test$CLobs <- "Dose required"
    test$CLpred <- "Dose required"
    test$Vss_obs <- "Dose required"
    test$Vss_pred <- "Dose required"
  }
  if (!is.null(ss_label_in_ss) | !is.null(sd_label_in_ss)) {
    if (!is.null(ss_label_in_ss) & !is.null(sd_label_in_ss)) {
      head(dat)
      ss1 <- dat[dat$ss == ss_label_in_ss & dat$time <= dat$tau, 
                 ]
      single <- dat[dat$ss == sd_label_in_ss & dat$time <= 
                      dat$tau, ]
    }
    if (!is.null(sd_label_in_ss) & is.null(ss_label_in_ss)) {
      ss1 <- dat[dat$ss != sd_label_in_ss & dat$time <= dat$tau, 
                 ]
      single <- dat[dat$ss == sd_label_in_ss & dat$time <= 
                      dat$tau, ]
    }
    if (is.null(sd_label_in_ss) & !is.null(ss_label_in_ss)) {
      ss1 <- dat[dat$ss == ss_label_in_ss & dat$time <= dat$tau, 
                 ]
      single <- dat[dat$ss != ss_label_in_ss & dat$time <= 
                      dat$tau, ]
    }
    
    auctau <- AUC(ss1, time = "time", id = "idss", dv = "dv")
    names(auctau) <- c("idss", "AUCtau")
    aucsdtau <- AUC(single, time = "time", id = "idss", dv = "dv")
    names(aucsdtau) <- c("idss", "AUCsd_tau")
    test<-join(test,idss,type="left")
    #test$idss <- paste(test$id, test$ss, sep = "-")
    test <- join(test, aucsdtau, type = "left")
    test <- join(test, auctau, type = "left")
    ident <- nodup(test, c("id", "idss"), "var")
    a <- join(auctau, ident)
    a$idss <- NULL
    b <- join(aucsdtau, ident)
    b$idss <- NULL
    EHL <- join(a, b)
    EHL$Rc <- with(EHL, AUCtau/AUCsd_tau)
    tau1 <- nodup(single, c("id", "tau"), "var")
    EHL <- join(EHL, tau1)
    EHL$EHL <- with(EHL, log(2) * tau/(log(Rc/(Rc - 1))))
    test <- join(test, EHL[, c("id", "tau", "EHL")], type = "left")
  }
  if ("AUCinf_obs" %in% names(test)) {
    single <- dat[dat$time <= dat$tau, ]
    aucsdtau <- AUC(single, time = "time", id = "idss", dv = "dv")
    names(aucsdtau) <- c("idss", "AUCsd_tau")
    head(test)
    test$idss <- paste(test$id, test$ss, sep = "-")
    test <- join(test, aucsdtau, type = "left")
    test$Rc <- with(test, AUCinf_obs/AUCsd_tau)
    test$EHL <- with(test, log(2) * tau/(log(Rc/(Rc - 1))))
    test <- join(test, nodup(dat, c("idss", "tau"), "var"))
  } else {
    if (!c("tau") %in% names(dat)) {
      test$EHL <- "AUCinf or TAU is missing"
    }
  }
  keep <- names(test)[!names(test) %in% c("Lambdac", "R2c", 
                                          "HLc", "thatc", "n_lambdac", "Clast_hatc", "intercc", 
                                          "idss", "Rc")]
  test <- test[, keep]
}
###########################
leonpheng/lhtool documentation built on Dec. 16, 2019, 3:52 a.m.