R/PlotTrajeR.R

In trajeR: Group Based Modeling Trajectory

###################################################################################
###################################################################################
# modification of plot's method for class trajectory
####################################################################################
###################################################################################
#################################################################################
# plot data
#################################################################################
plotTrajCensoredNormalLikelihood <- function(
  beta,
  sigma,
  theta,
  delta,
  plotcov,
  Y,
  A,
  X,
  mean,
  alpha,
  ng,
  n,
  Time,
  col,
  degre,
  ymin,
  ymax,
  method,
  ...
) {
  period <- length(Time)
  if (is.null(plotcov)) {
    plotcov <- matrix(rep(0, period), nrow = 1)
  }
  nbeta <- degre + 1
  ntheta <- length(theta) / ng
  betatmp <- beta
  j <- 1
  beta <- list()
  for (i in seq_along(nbeta)) {
    beta[[i]] <- betatmp[j:sum(nbeta[1:i])]
    j <- sum(nbeta[1:i]) + 1
  }
  if (any(is.na(delta))) {
    nw <- 0
    delta <- rep(list(0), ng)
  } else {
    nw <- length(delta) / ng
    deltatmp <- delta
    delta <- list()
    for (i in 1:ng) {
      delta[[i]] <- deltatmp[((i - 1) * nw + 1):(i * nw)]
    }
  }
  if (length(col) == 2 * ng) {
    cols1 <- col[1:ng]
    cols2 <- col[(ng + 1):(2 * ng)]
  } else {
    cols1 <- grDevices::gray.colors(
      ng,
      start = 0.3,
      end = 0.7,
      gamma = 2.2,
      alpha = 0.5
    )
    cols2 <- rep("black", ng)
  }
  pas <- seq(Time[1], Time[period], (Time[period] - Time[1]) / 100)
  if (is.null(Y) | is.null(A)) {
    vec <- c()
    for (i in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        sum(beta[[i]] * pas[s]**(0:(nbeta[i] - 1)))
      })
      tmp[tmp < ymin] <- ymin
      tmp[tmp > ymax] <- ymax
      vec <- cbind(vec, tmp)
    }
    graphics::matplot(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1,
      ...
    )
  } else {
    Ygr <- cbind(Y, rep(NA, n))
    if (is.null(X)) {
      X <- cbind(rep(1, n))
    }
    if (method == "L") {
      tmp <- sapply(1:ng, function(s) {
        piik(theta, i, s, ng, X) *
          gkCNORM_cpp(
            beta,
            sigma,
            i,
            s,
            nbeta,
            A,
            Y,
            ymin,
            ymax,
            TCOV = NULL,
            delta = NULL,
            nw = 0
          )
      })
    } else {
      tmp <- sapply(1:ng, function(s) {
        theta[s] *
          gkCNORM_cpp(
            beta,
            sigma,
            i,
            s,
            nbeta,
            A,
            Y,
            ymin,
            ymax,
            TCOV = NULL,
            delta = NULL,
            nw = 0
          )
      })
    }
    ncol <- which.max(tmp / sum(tmp))
    plot(
      A[1, ],
      Y[1, ],
      type = "b",
      ylim = c(min(Y, na.rm = T), max(Y, na.rm = T)),
      pch = 16,
      col = cols1[ncol],
      ...
    )
    Ygr[1, period + 1] <- ncol
    for (i in 2:n) {
      if (method == "L") {
        tmp <- sapply(1:ng, function(s) {
          piik(theta, i, s, ng, X) *
            gkCNORM_cpp(
              beta,
              sigma,
              i,
              s,
              nbeta,
              A,
              Y,
              ymin,
              ymax,
              TCOV = NULL,
              delta = NULL,
              nw = 0
            )
        })
      } else {
        tmp <- sapply(1:ng, function(s) {
          theta[s] *
            gkCNORM_cpp(
              beta,
              sigma,
              i,
              s,
              nbeta,
              A,
              Y,
              ymin,
              ymax,
              TCOV = NULL,
              delta = NULL,
              nw = 0
            )
        })
      }
      ncol <- which.max(tmp / sum(tmp))
      graphics::lines(A[i, ], Y[i, ], type = "b", pch = 16, col = cols1[ncol])
      Ygr[i, period + 1] <- ncol
    }
    vec <- c()
    for (k in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        sum(beta[[k]] * pas[s]**(0:(nbeta[k] - 1)))
      })
      tmp[tmp < ymin] <- ymin
      tmp[tmp > ymax] <- ymax
      vec <- cbind(vec, tmp)
    }
    graphics::matlines(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1
    )
    for (k in 1:ng) {
      if (mean == TRUE) {
        graphics::lines(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k],
          lty = 2
        )
        graphics::points(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k]
        )
      }
    }
  }
  if (is.matrix(plotcov) == FALSE) {
    plotcov <- matrix(plotcov, ncol = period * nw, byrow = TRUE)
  }
  if (any(plotcov != 0)) {
    pas <- seq(Time[1], Time[period], (Time[period] - Time[1]) / 100)
    for (k in 1:ng) {
      for (row in 1:nrow(plotcov)) {
        vec <- sapply(1:length(pas), function(s) {
          sum(
            beta[[k]] *
              pas[s]**(0:(nbeta[k] - 1)) +
              sum(
                delta[[k]] *
                  plotcov[
                    row,
                    seq(
                      from = floor(pas[s]),
                      to = floor(pas[s]) + (nw - 1) * period,
                      by = period
                    )
                  ]
              )
          )
        })
        graphics::lines(
          pas,
          vec,
          type = "l",
          pch = 16,
          col = cols2[k],
          lwd = 2,
          lty = 2
        )
      }
    }
  }
}
##########################################################################################
# plot LOGIT function
##########################################################################################
plotTrajLOGIT <- function(
  beta,
  theta,
  delta,
  Y,
  A,
  X,
  ng,
  n,
  Time,
  dec,
  col,
  degre,
  plotcov,
  mean,
  alpha,
  ...
) {
  period <- length(Time)
  Ygr <- cbind(Y, rep(NA, n))
  if (is.null(plotcov)) {
    plotcov <- matrix(rep(0, period), nrow = 1)
  }
  nbeta <- degre + 1
  ntheta <- length(theta) / ng
  betatmp <- beta
  j <- 1
  beta <- list()
  if (any(is.na(delta))) {
    nw <- 0
    delta <- rep(list(0), ng)
  } else {
    nw <- length(delta) / ng
  }
  for (i in 1:length(nbeta)) {
    beta[[i]] <- betatmp[j:sum(nbeta[1:i])]
    j <- sum(nbeta[1:i]) + 1
  }
  if (length(col) == 2 * ng) {
    cols1 <- col[1:ng]
    cols2 <- col[(ng + 1):(2 * ng)]
  } else {
    cols1 <- grDevices::gray.colors(
      ng,
      start = 0.3,
      end = 0.7,
      gamma = 2.2,
      alpha = 0.5
    )
    cols2 <- rep("black", ng)
  }
  pas <- seq(Time[1], Time[period], (Time[period] - Time[1]) / 100)
  if (is.null(Y) | is.null(A)) {
    vec <- c()
    for (i in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        exp(sum(beta[[i]] * pas[s]**(0:(nbeta[i] - 1))))
      })
      vec <- cbind(vec, tmp / (1 + tmp))
    }
    graphics::matplot(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1,
      ...
    )
  } else {
    if (is.null(X)) {
      X <- cbind(rep(1, n))
    }
    tmp <- sapply(1:ng, function(s) {
      piik(theta, 1, s, ng, X) *
        gkLOGIT_cpp(
          beta,
          i - 1,
          s - 1,
          nbeta,
          A,
          Y,
          TCOV = NULL,
          delta = NULL,
          nw = 0
        )
    })
    ncol <- which.max(tmp / sum(tmp))
    a <- dec * stats::runif(n * period, 0, 1) * 2 * 3.14159
    r <- dec * stats::runif(n * period, 0, 0.01)
    decx <- r * cos(a)
    decy <- r * sin(a)
    d <- max(decx)
    dy <- max(decy)
    xlim <- c(min(A) - d, max(A) + d)
    plot(
      A[1, ] + decx[1:period],
      Y[1, ] + decy[1:period],
      type = "b",
      ylim = c(-0.2, 1.2),
      xlim = xlim,
      pch = 16,
      col = cols1[ncol],
      ...
    )
    graphics::axis(side = 2, at = seq(0, 1, 0.25), labels = seq(0, 1, 0.25))
    graphics::rect(
      xleft = xlim[1] - d,
      ybottom = 1 - dy,
      xright = xlim[2] + d,
      ytop = 1 + dy,
      col = "gray91",
      border = NA
    )
    graphics::rect(
      xleft = xlim[1] - d,
      ybottom = -dy,
      xright = xlim[2] + d,
      ytop = dy,
      col = "gray91",
      border = NA
    )
    Ygr[1, period + 1] <- ncol
    for (i in 2:n) {
      tmp <- sapply(1:ng, function(s) {
        piik(theta, i, s, ng, X) *
          gkLOGIT_cpp(
            beta,
            i - 1,
            s - 1,
            nbeta,
            A,
            Y,
            TCOV = NULL,
            delta = NULL,
            nw = 0
          )
      })
      ncol <- which.max(tmp / sum(tmp))
      graphics::lines(
        A[i, ] + decx[((i - 1) * period + 1):(i * period)],
        Y[i, ] + decy[i],
        type = "b",
        pch = 16,
        col = grDevices::adjustcolor(cols1[ncol], alpha.f = alpha)
      )
      graphics::points(
        A[i, ] + decx[((i - 1) * period + 1):(i * period)],
        Y[i, ] + decy[i],
        pch = 16,
        col = cols1[ncol]
      )
      Ygr[i, period + 1] <- ncol
    }
    vec <- c()
    for (k in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        exp(sum(beta[[k]] * pas[s]**(0:(nbeta[k] - 1))))
      })
      vec <- cbind(vec, tmp / (1 + tmp))
    }
    graphics::matlines(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1
    )
    for (k in 1:ng) {
      if (mean == TRUE) {
        graphics::lines(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k],
          lty = 2
        )
        graphics::points(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k]
        )
      }
    }
  }
  if (is.matrix(plotcov) == FALSE) {
    plotcov <- matrix(plotcov, ncol = period)
  }
  if (any(plotcov != 0)) {
    pas <- seq(Time[1], Time[period], (Time[period] - Time[1]) / 100)
    for (k in 1:ng) {
      for (row in 1:nrow(plotcov)) {
        vec <- sapply(1:length(pas), function(s) {
          exp(sum(
            beta[[k]] *
              pas[s]**(0:(nbeta[k] - 1)) +
              sum(delta[[k]] * plotcov[row, floor(pas[s])])
          ))
        })
        vec <- vec / (1 + vec)
        graphics::lines(
          pas,
          vec,
          type = "l",
          pch = 16,
          col = cols2[k],
          lwd = 2,
          lty = 2
        )
      }
    }
  }
}
##########################################################################################
# plot POIS function
##########################################################################################
plotTrajPOIS <- function(
  beta,
  theta,
  delta,
  Y,
  A,
  X,
  TCOV,
  ng,
  n,
  Time,
  dec,
  col,
  degre,
  plotcov,
  mean,
  alpha,
  method,
  ...
) {
  period <- length(Time)
  Ygr <- cbind(Y, rep(NA, n))
  if (is.null(plotcov)) {
    plotcov <- matrix(rep(0, period), nrow = 1)
  }
  nbeta <- degre + 1
  ntheta <- length(theta) / ng
  betatmp <- beta
  j <- 1
  beta <- list()
  if (any(is.na(delta))) {
    nw <- 0
    delta <- rep(list(0), ng)
  } else {
    nw <- length(delta) / ng
  }
  for (i in 1:length(nbeta)) {
    beta[[i]] <- betatmp[j:sum(nbeta[1:i])]
    j <- sum(nbeta[1:i]) + 1
  }
  if (length(col) == 2 * ng) {
    cols1 <- col[1:ng]
    cols2 <- col[(ng + 1):(2 * ng)]
  } else {
    cols1 <- grDevices::gray.colors(
      ng,
      start = 0.3,
      end = 0.7,
      gamma = 2.2,
      alpha = 0.5
    )
    cols2 <- rep("black", ng)
  }
  pas <- seq(Time[1], Time[period], (Time[period] - Time[1]) / 100)
  if (is.null(Y) | is.null(A)) {
    vec <- c()
    for (i in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        exp(sum(beta[[i]] * pas[s]**(0:(nbeta[i] - 1))))
      })
      vec <- cbind(vec, tmp)
    }
    graphics::matplot(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1,
      ...
    )
  } else {
    if (is.null(X)) {
      X <- cbind(rep(1, n))
    }
    tmp <- sapply(1:ng, function(s) {
      piik(theta, 1, s, ng, X) *
        gkPois_cpp(
          beta,
          i - 1,
          s - 1,
          nbeta,
          A,
          Y,
          TCOV = NULL,
          delta = NULL,
          nw = 0
        )
    })
    ncol <- which.max(tmp / sum(tmp))
    a <- dec * stats::runif(n * period, 0, 1) * 2 * 3.14159
    r <- dec * stats::runif(n * period, 0, 0.01)
    decx <- r * cos(a)
    decy <- r * sin(a)
    d <- max(decx)
    dy <- max(decy)
    xlim <- c(min(A), max(A))
    ylim <- c(min(Y), max(Y))
    plot(
      A[1, ] + decx[1:period],
      Y[1, ] + decy[1:period],
      type = "b",
      ylim = c(min(Y, na.rm = T), max(Y, na.rm = T)),
      pch = 16,
      col = cols1[ncol],
      ...
    )
    Ygr[1, period + 1] <- ncol
    for (i in 2:n) {
      if (method == "L") {
        tmp <- sapply(1:ng, function(s) {
          piik(theta, i, s, ng, X) *
            gkPois_cpp(beta, i - 1, s - 1, nbeta, A, Y, TCOV, delta, nw)
        })
      } else {
        tmp <- sapply(1:ng, function(s) {
          theta[s] *
            gkPois_cpp(beta, i - 1, s - 1, nbeta, A, Y, TCOV, delta, nw)
        })
      }
      ncol <- which.max(tmp / sum(tmp))
      graphics::lines(
        A[i, ] + decx[((i - 1) * period + 1):(i * period)],
        Y[i, ] + decy[i],
        type = "b",
        pch = 16,
        col = grDevices::adjustcolor(cols1[ncol], alpha.f = alpha)
      )
      graphics::points(
        A[i, ] + decx[((i - 1) * period + 1):(i * period)],
        Y[i, ] + decy[i],
        pch = 16,
        col = cols1[ncol]
      )
      Ygr[i, period + 1] <- ncol
    }
    vec <- c()
    for (k in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        exp(sum(beta[[k]] * pas[s]**(0:(nbeta[k] - 1))))
      })
      vec <- cbind(vec, tmp)
    }
    graphics::matlines(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1
    )
    for (k in 1:ng) {
      if (mean == TRUE) {
        graphics::lines(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k],
          lty = 2
        )
        graphics::points(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k]
        )
      }
    }
  }
  if (is.matrix(plotcov) == FALSE) {
    plotcov <- matrix(plotcov, ncol = period)
  }
  if (any(plotcov != 0)) {
    pas <- seq(Time[1], Time[period], (Time[period] - Time[1]) / 100)
    for (k in 1:ng) {
      for (row in 1:nrow(plotcov)) {
        vec <- sapply(1:length(pas), function(s) {
          exp(sum(
            beta[[k]] *
              pas[s]**(0:(nbeta[k] - 1)) +
              sum(delta[[k]] * plotcov[row, floor(pas[s])])
          ))
        })
        graphics::lines(
          pas,
          vec,
          type = "l",
          pch = 16,
          col = cols2[k],
          lwd = 2,
          lty = 2
        )
      }
    }
  }
}
##########################################################################################
# plot ZIP function
##########################################################################################
plotTrajZIP <- function(
  beta,
  nu,
  theta,
  delta,
  Y,
  A,
  X,
  TCOV,
  ng,
  n,
  Time,
  dec,
  col,
  degre,
  plotcov,
  mean,
  alpha,
  method,
  degre.nu,
  period,
  ...
) {
  if (is.null(plotcov)) {
    plotcov <- matrix(rep(0, period), nrow = 1)
  }
  nbeta <- degre + 1
  nnu <- degre.nu + 1
  ntheta <- length(theta) / ng
  j <- 1
  betatmp <- beta
  beta <- list()
  for (i in 1:ng) {
    beta[[i]] <- betatmp[j:sum(nbeta[1:i])]
    j <- sum(nbeta[1:i]) + 1
  }
  j <- 1
  nutmp <- nu
  nu <- list()
  for (i in 1:ng) {
    nu[[i]] <- nutmp[j:sum(nnu[1:i])]
    j <- sum(nnu[1:i]) + 1
  }
  if (any(is.na(delta))) {
    nw <- 0
    delta <- rep(list(0), ng)
  } else {
    nw <- length(delta) / ng
    deltatmp <- delta
    ndeltacum <- cumsum(c(0, rep(nw, ng)))
    delta <- list()
    for (i in 1:ng) {
      delta[[i]] <- deltatmp[(ndeltacum[i] + 1):(ndeltacum[i + 1])]
    }
  }
  if (length(col) == 2 * ng) {
    cols1 <- col[1:ng]
    cols2 <- col[(ng + 1):(2 * ng)]
  } else {
    cols1 <- grDevices::gray.colors(
      ng,
      start = 0.3,
      end = 0.7,
      gamma = 2.2,
      alpha = 0.5
    )
    cols2 <- rep("black", ng)
  }
  pas <- seq(Time[1], Time[2], (Time[2] - Time[1]) / 100)
  if (is.null(Y) | is.null(A)) {
    vec <- c()
    for (k in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        exp(sum(beta[[k]] * pas[s]**(0:(nbeta[k] - 1))))
      })
      tmp2 <- sapply(1:length(pas), function(s) {
        exp(sum(nu[[k]] * pas[s]**(0:(nnu[k] - 1)))) /
          (1 + exp(sum(nu[[k]] * pas[s]**(0:(nnu[k] - 1)))))
      })
      vec <- cbind(vec, (1 - tmp2) * tmp)
    }
    graphics::matplot(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1,
      ...
    )
  } else {
    Ygr <- cbind(Y, rep(NA, n))
    if (is.null(X)) {
      X <- cbind(rep(1, n))
    }
    if (method == "L") {
      tmp <- sapply(1:ng, function(s) {
        piik(theta, i, s, ng, X) *
          gkZIP_cpp(beta, nu, 0, s - 1, nbeta, nnu, A, Y, TCOV, delta, nw)
      })
    } else {
      tmp <- sapply(1:ng, function(s) {
        theta[s] *
          gkZIP_cpp(beta, nu, 0, s - 1, nbeta, nnu, A, Y, TCOV, delta, nw)
      })
    }
    ncol <- which.max(tmp / sum(tmp))
    a <- dec * stats::runif(n * period, 0, 1) * 2 * 3.14159
    r <- dec * stats::runif(n * period, 0, 0.01)
    decx <- r * cos(a)
    decy <- r * sin(a)
    d <- max(decx)
    dy <- max(decy)
    plot(
      A[1, ] + decx[1:period],
      Y[1, ] + decy[1:period],
      type = "b",
      ylim = c(min(Y, na.rm = T), max(Y, na.rm = T)),
      xlim = c(Time[1], Time[2]),
      pch = 16,
      col = cols1[ncol],
      ...
    )
    Ygr[1, period + 1] <- ncol
    for (i in 2:n) {
      if (method == "L") {
        tmp <- sapply(1:ng, function(s) {
          piik(theta, i, s, ng, X) *
            gkZIP_cpp(beta, nu, i - 1, s - 1, nbeta, nnu, A, Y, TCOV, delta, nw)
        })
      } else {
        tmp <- sapply(1:ng, function(s) {
          theta[s] *
            gkZIP_cpp(beta, nu, i - 1, s - 1, nbeta, nnu, A, Y, TCOV, delta, nw)
        })
      }
      ncol <- which.max(tmp / sum(tmp))
      graphics::lines(
        A[i, ] + decx[((i - 1) * period + 1):(i * period)],
        Y[i, ] + decy[i],
        type = "b",
        pch = 16,
        col = grDevices::adjustcolor(cols1[ncol], alpha.f = alpha)
      )
      graphics::points(
        A[i, ] + decx[((i - 1) * period + 1):(i * period)],
        Y[i, ] + decy[i],
        pch = 16,
        col = cols1[ncol]
      )
      Ygr[i, period + 1] <- ncol
    }
    vec <- c()
    for (k in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        exp(sum(beta[[k]] * pas[s]**(0:(nbeta[k] - 1))))
      })
      tmp2 <- sapply(1:length(pas), function(s) {
        exp(sum(nu[[k]] * pas[s]**(0:(nnu[k] - 1)))) /
          (1 + exp(sum(nu[[k]] * pas[s]**(0:(nnu[k] - 1)))))
      })
      vec <- cbind(vec, (1 - tmp2) * tmp)
    }
    graphics::matlines(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1
    )
    for (k in 1:ng) {
      if (mean == TRUE) {
        graphics::lines(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k],
          lty = 2
        )
        graphics::points(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k]
        )
      }
    }
  }
  if (is.matrix(plotcov) == FALSE) {
    plotcov <- matrix(plotcov, ncol = period, byrow = TRUE)
  }
  if (any(plotcov != 0)) {
    pas <- 1:period
    for (k in 1:ng) {
      for (row in 1:nrow(plotcov)) {
        tmp <- sapply(1:length(pas), function(s) {
          exp(
            sum(beta[[k]] * pas[s]**(0:(nbeta[k] - 1))) +
              sum(delta[[k]] * plotcov[row, s])
          )
        })
        tmp2 <- sapply(1:length(pas), function(s) {
          exp(sum(nu[[k]] * pas[s]**(0:(nnu[k] - 1)))) /
            (1 + exp(sum(nu[[k]] * pas[s]**(0:(nnu[k] - 1)))))
        })
        vec <- (1 - tmp2) * tmp
        graphics::lines(
          pas,
          vec,
          type = "l",
          pch = 16,
          col = cols2[k],
          lwd = 2,
          lty = 2
        )
      }
    }
  }
}
#################################################################################
plotTrajNL <- function(
  beta,
  sigma,
  theta,
  plotcov,
  Y,
  A,
  X,
  mean,
  alpha,
  ng,
  n,
  Time,
  col,
  degre,
  ymin,
  ymax,
  method,
  fct,
  main,
  TCOV,
  ...
) {
  period <- length(Time)
  if (is.null(plotcov)) {
    plotcov <- matrix(rep(0, period), nrow = 1)
  }
  nbeta <- degre + 1
  ntheta <- length(theta) / ng
  betatmp <- beta
  j <- 1
  beta <- list()
  nw <- 0
  for (i in 1:length(nbeta)) {
    beta[[i]] <- betatmp[j:sum(nbeta[1:i])]
    j <- sum(nbeta[1:i]) + 1
  }
  if (length(col) == 2 * ng) {
    cols1 <- col[1:ng]
    cols2 <- col[(ng + 1):(2 * ng)]
  } else {
    cols1 <- grDevices::gray.colors(
      ng,
      start = 0.3,
      end = 0.7,
      gamma = 2.2,
      alpha = 0.5
    )
    cols2 <- rep("black", ng)
  }
  pas <- seq(Time[1], Time[period], (Time[period] - Time[1]) / 100)
  if (is.null(Y) | is.null(A)) {
    vec <- c()
    for (i in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        fct(pas[s], beta[[i]], TCOV)
      })
      vec <- cbind(vec, tmp)
    }
    graphics::matplot(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1,
      ...
    )
  } else {
    Ygr <- cbind(Y, rep(NA, n))
    if (is.null(X)) {
      X <- cbind(rep(1, n))
    }
    if (method == "L") {
      tmp <- sapply(1:ng, function(s) {
        piik(theta, i, s, ng, X) * gkNL(beta, sigma, i, s, TCOV, A, Y, fct)
      })
    } else {
      tmp <- sapply(1:ng, function(s) {
        theta[s] * gkNL(beta, sigma, i, s, TCOV, A, Y, fct)
      })
    }
    ncol <- which.max(tmp / sum(tmp))
    plot(
      A[1, ],
      Y[1, ],
      type = "b",
      ylim = c(min(Y), max(Y)),
      pch = 16,
      col = cols1[ncol],
      ...
    )
    Ygr[1, period + 1] <- ncol
    for (i in 2:n) {
      if (method == "L") {
        tmp <- sapply(1:ng, function(s) {
          piik(theta, i, s, ng, X) * gkNL(beta, sigma, i, s, TCOV, A, Y, fct)
        })
      } else {
        tmp <- sapply(1:ng, function(s) {
          theta[s] * gkNL(beta, sigma, i, s, TCOV, A, Y, fct)
        })
      }
      ncol <- which.max(tmp / sum(tmp))
      graphics::lines(A[i, ], Y[i, ], type = "b", pch = 16, col = cols1[ncol])
      Ygr[i, period + 1] <- ncol
    }
    vec <- c()
    for (k in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        fct(pas[s], beta[[k]], TCOV)
      })
      vec <- cbind(vec, tmp)
    }
    graphics::matlines(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1
    )
    for (k in 1:ng) {
      if (mean == TRUE) {
        graphics::lines(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k],
          lty = 2
        )
        graphics::points(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k]
        )
      }
    }
  }
}
##########################################################################################
# plot BETA function
##########################################################################################
plotTrajBETA <- function(
  beta,
  phi,
  theta,
  delta,
  Y,
  A,
  X,
  TCOV,
  ng,
  n,
  Time,
  dec,
  col,
  degre,
  plotcov,
  mean,
  alpha,
  method,
  degre.phi,
  ...
) {
  period <- length(Time)
  if (is.null(plotcov)) {
    plotcov <- matrix(rep(0, period), nrow = 1)
  }
  nbeta <- degre + 1
  nphi <- degre.phi + 1
  ntheta <- length(theta) / ng
  j <- 1
  betatmp <- beta
  beta <- list()
  for (i in 1:ng) {
    beta[[i]] <- betatmp[j:sum(nbeta[1:i])]
    j <- sum(nbeta[1:i]) + 1
  }
  j <- 1
  phitmp <- phi
  phi <- list()
  for (i in 1:ng) {
    phi[[i]] <- phitmp[j:sum(nphi[1:i])]
    j <- sum(nphi[1:i]) + 1
  }
  if (any(is.na(delta))) {
    nw <- 0
    delta <- rep(list(0), ng)
  } else {
    nw <- length(delta) / ng
    deltatmp <- delta
    ndeltacum <- cumsum(c(0, rep(nw, ng)))
    delta <- list()
    for (i in 1:ng) {
      delta[[i]] <- deltatmp[(ndeltacum[i] + 1):(ndeltacum[i + 1])]
    }
  }
  if (length(col) == 2 * ng) {
    cols1 <- col[1:ng]
    cols2 <- col[(ng + 1):(2 * ng)]
  } else {
    cols1 <- grDevices::gray.colors(
      ng,
      start = 0.3,
      end = 0.7,
      gamma = 2.2,
      alpha = 0.5
    )
    cols2 <- rep("black", ng)
  }
  pas <- seq(Time[1], Time[period], (Time[period] - Time[1]) / 100)
  if (is.null(Y) | is.null(A)) {
    vec <- c()
    for (i in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        exp(sum(beta[[i]] * pas[s]**(0:(nbeta[i] - 1))))
      })
      vec <- cbind(vec, tmp / (1 + tmp))
    }
    graphics::matplot(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1,
      ...
    )
  } else {
    Ygr <- cbind(Y, rep(NA, n))
    if (is.null(X)) {
      X <- cbind(rep(1, n))
    }
    if (method == "L") {
      tmp <- sapply(1:ng, function(s) {
        piik(theta, i, s, ng, X) *
          gkBETA_cpp(beta, phi, 0, s - 1, nbeta, nphi, A, Y, TCOV, delta, nw)
      })
    } else {
      tmp <- sapply(1:ng, function(s) {
        theta[s] *
          gkBETA_cpp(beta, phi, 0, s - 1, nbeta, nphi, A, Y, TCOV, delta, nw)
      })
    }
    ncol <- which.max(tmp / sum(tmp))
    xlim <- c(min(A), max(A))
    ylim <- c(min(Y, na.rm = TRUE), max(Y, na.rm = TRUE))
    plot(
      A[1, ],
      Y[1, ],
      type = "b",
      ylim = c(min(Y, na.rm = T), max(Y, na.rm = T)),
      pch = 16,
      col = cols1[ncol],
      ...
    )
    Ygr[1, period + 1] <- ncol
    for (i in 2:n) {
      if (method == "L") {
        tmp <- sapply(1:ng, function(s) {
          piik(theta, i, s, ng, X) *
            gkBETA_cpp(
              beta,
              phi,
              i - 1,
              s - 1,
              nbeta,
              nphi,
              A,
              Y,
              TCOV,
              delta,
              nw
            )
        })
      } else {
        tmp <- sapply(1:ng, function(s) {
          theta[s] *
            gkBETA_cpp(
              beta,
              phi,
              i - 1,
              s - 1,
              nbeta,
              nphi,
              A,
              Y,
              TCOV,
              delta,
              nw
            )
        })
      }
      ncol <- which.max(tmp / sum(tmp))
      graphics::lines(
        A[i, ],
        Y[i, ],
        type = "b",
        pch = 16,
        col = grDevices::adjustcolor(cols1[ncol], alpha.f = alpha)
      )
      graphics::points(A[i, ], Y[i, ], pch = 16, col = cols1[ncol])
      Ygr[i, period + 1] <- ncol
    }
    vec <- c()
    for (k in 1:ng) {
      tmp <- sapply(1:length(pas), function(s) {
        exp(sum(beta[[k]] * pas[s]**(0:(nbeta[k] - 1))))
      })
      vec <- cbind(vec, tmp / (1 + tmp))
    }
    graphics::matlines(
      pas,
      vec,
      type = "l",
      pch = 16,
      col = cols2,
      lwd = 4,
      lty = 1
    )
    for (k in 1:ng) {
      if (mean == TRUE) {
        graphics::lines(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k],
          lty = 2
        )
        graphics::points(
          A[1, ],
          colMeans(Ygr[Ygr[, period + 1] == k, -(period + 1)]),
          col = cols2[k]
        )
      }
    }
  }
  if (is.matrix(plotcov) == FALSE) {
    plotcov <- matrix(plotcov, ncol = period, byrow = TRUE)
  }
  if (any(plotcov != 0)) {
    pas <- 1:period
    for (k in 1:ng) {
      for (row in 1:nrow(plotcov)) {
        tmp <- sapply(1:length(pas), function(s) {
          exp(
            sum(beta[[k]] * pas[s]**(0:(nbeta[k] - 1))) +
              sum(delta[[k]] * plotcov[row, s])
          )
        })
        vec <- tmp / (1 + tmp)
        graphics::lines(
          pas,
          vec,
          type = "l",
          pch = 16,
          col = cols2[k],
          lwd = 2,
          lty = 2
        )
      }
    }
  }
}
###################################################################################
# Generic function to plot the trajectory
####################################################################################
#' @title Plot Group Trajectories from a trajeR Model
#'
#' @description A generic plotting function for objects of class `Trajectory.*`. It visualizes both individual and average group trajectories.
#'
#' @param Obj A trajectory object returned by the \code{trajeR} function.
#' @param ... Other parameters passed to the specific plot methods.
#'
#' @return This function is called for its side effect of creating a plot. It does not return any value.
#' @export
#'
#' @examples
#' # See examples for specific methods like plotrajeR.Trajectory.CNORM
#' data <- read.csv(system.file("extdata", "CNORM2gr.csv", package = "trajeR"))
#' sol <- trajeR(Y = data[, 2:6], A = data[, 7:11], Risk = data[, 12, drop = FALSE],
#'               degre = c(2, 2), Model = "CNORM", Method = "EM")
#' plotrajeR(sol)
#'
plotrajeR <- function(Obj, ...) {
  UseMethod("plotrajeR", Obj)
}

###################################################################################
# modification of plot's method for class trajectory.CNORM
####################################################################################
#' @title Plot Method for CNORM Trajectory Objects
#'
#' @description Plots the trajectories for a Censored Normal model fitted by `trajeR`.
#'
#' @param Obj An object of class `Trajectory.CNORM`.
#' @param Y The original response variable matrix used to fit the model. If provided, individual trajectories are plotted.
#' @param A The original time variable matrix used to fit the model.
#' @param Risk Optional. The original risk covariate matrix used to fit the model.
#' @param plotcov An optional vector or matrix of time-dependent covariates to visualize their effect on trajectories.
#' @param col An optional vector of colors. If its length is `2 * ng`, the first `ng` colors are for individual trajectories and the next `ng` are for the mean trajectories.
#' @param mean Logical. If `TRUE`, the mean trajectory for each group (based on assigned individuals) is plotted as a dashed line.
#' @param alpha A numeric value between 0 and 1 for the transparency of individual trajectory lines.
#' @param ... Other graphical parameters passed to `plot`.
#'
#' @return a graphic.
#' @export
#'
plotrajeR.Trajectory.CNORM <- function(
  Obj,
  plotcov = NULL,
  col = "black",
  Y = NULL,
  A = NULL,
  Risk = NULL,
  mean = FALSE,
  alpha = 1,
  ...
) {
  if (!is.null(Y)) {
    Y <- data.matrix(Y)
  }
  if (!is.null(A)) {
    A <- data.matrix(A)
  }
  plotTrajCensoredNormalLikelihood(
    beta = Obj$beta,
    sigma = Obj$sigma,
    theta = Obj$theta,
    delta = Obj$delta,
    plotcov = plotcov,
    Y = Y,
    A = A,
    X = Risk,
    mean = mean,
    alpha = alpha,
    ng = Obj$groups,
    n = Obj$Size,
    Time = Obj$Time,
    col = col,
    degre = Obj$degre,
    ymin = Obj$min,
    ymax = Obj$max,
    method = Obj$Method,
    ...
  )
}
###################################################################################
# modification of plot's method for class trajectory.LOGIT
####################################################################################
#' @title Plot Method for LOGIT Trajectory Objects
#'
#' @description Plots the trajectories for a Logistic model fitted by `trajeR`.
#'
#' @param Obj An object of class `Trajectory.LOGIT`.
#' @param Y The original response variable matrix. If provided, individual data points are jittered and plotted.
#' @param A The original time variable matrix.
#' @param Risk Optional. The original risk covariate matrix.
#' @param plotcov An optional vector or matrix of time-dependent covariates.
#' @param dec A numeric value controlling the amount of jitter applied to data points for better visualization.
#' @param col An optional vector of colors for individual points and mean trajectories.
#' @param mean Logical. If `TRUE`, plots the mean of the observed data for each group.
#' @param alpha A numeric value for the transparency of individual trajectory lines.
#' @param ... Other graphical parameters passed to `plot`.
#'
#' @return a graphic.
#' @export
#'

plotrajeR.Trajectory.LOGIT <- function(
  Obj,
  plotcov = NULL,
  dec = 0,
  col = "black",
  Y = NULL,
  A = NULL,
  Risk = NULL,
  mean = FALSE,
  alpha = 1,
  ...
) {
  if (!is.null(Y)) {
    Y <- data.matrix(Y)
  }
  if (!is.null(A)) {
    A <- data.matrix(A)
  }
  plotTrajLOGIT(
    beta = Obj$beta,
    theta = Obj$theta,
    delta = Obj$delta,
    plotcov = plotcov,
    Y = Y,
    A = A,
    X = Risk,
    mean = mean,
    alpha = alpha,
    ng = Obj$groups,
    n = Obj$Size,
    Time = Obj$Time,
    dec = dec,
    col = col,
    degre = Obj$degre,
    ...
  )
}
###################################################################################
# modification of plot's method for class trajectory.POIS
####################################################################################
#' @title Plot Method for POIS Trajectory Objects
#'
#' @description Plots the trajectories for a Poisson model fitted by `trajeR`.
#'
#' @param Obj An object of class `Trajectory.POIS`.
#' @param Y The original response variable matrix. If provided, individual trajectories are plotted.
#' @param A The original time variable matrix.
#' @param Risk Optional. The original risk covariate matrix.
#' @param plotcov An optional vector or matrix of time-dependent covariates.
#' @param dec A numeric value controlling the amount of jitter applied to data points.
#' @param col An optional vector of colors.
#' @inheritParams trajeR
#' @param mean Logical. If `TRUE`, plots the mean trajectory for each group.
#' @param alpha A numeric value for the transparency of individual trajectory lines.
#' @param ... Other graphical parameters passed to `plot`.
#'
#' @return a graphic.
#' @export
#'

plotrajeR.Trajectory.POIS <- function(
  Obj,
  plotcov = NULL,
  dec = 0,
  col = "black",
  Y = NULL,
  A = NULL,
  Risk = NULL,
  TCOV = NULL,
  mean = FALSE,
  alpha = 1,
  ...
) {
  if (!is.null(Y)) {
    Y <- data.matrix(Y)
  }
  if (!is.null(A)) {
    A <- data.matrix(A)
  }
  plotTrajPOIS(
    beta = Obj$beta,
    theta = Obj$theta,
    delta = Obj$delta,
    plotcov = plotcov,
    Y = Y,
    A = A,
    X = Risk,
    TCOV = TCOV,
    mean = mean,
    alpha = alpha,
    method = Obj$Method,
    ng = Obj$groups,
    n = Obj$Size,
    Time = Obj$Time,
    dec = dec,
    col = col,
    degre = Obj$degre,
    ...
  )
}
###################################################################################
# modification of plot's method for class trajectory.ZIP
####################################################################################
#' @title Plot Method for ZIP Trajectory Objects
#'
#' @description Plots the trajectories for a Zero-Inflated Poisson (ZIP) model fitted by `trajeR`.
#'
#' @param Obj An object of class `Trajectory.ZIP`.
#' @param Y The original response variable matrix.
#' @param A The original time variable matrix.
#' @param Risk Optional. The original risk covariate matrix.
#' @param plotcov An optional vector or matrix of time-dependent covariates.
#' @param dec A numeric value controlling the amount of jitter applied to data points.
#' @param col An optional vector of colors.
#' @inheritParams trajeR
#' @param mean Logical. If `TRUE`, plots the mean trajectory for each group.
#' @param alpha A numeric value for the transparency of individual trajectory lines.
#' @param ... Other graphical parameters passed to `plot`.
#'
#' @return a graphic.
#' @export
#'
plotrajeR.Trajectory.ZIP <- function(
  Obj,
  plotcov = NULL,
  dec = 0,
  col = "black",
  Y = NULL,
  A = NULL,
  Risk = NULL,
  TCOV = NULL,
  mean = FALSE,
  alpha = 1,
  ...
) {
  if (!is.null(Y)) {
    Y <- data.matrix(Y)
  }
  if (!is.null(A)) {
    A <- data.matrix(A)
  }
  plotTrajZIP(
    beta = Obj$beta,
    nu = Obj$nu,
    theta = Obj$theta,
    delta = Obj$delta,
    plotcov = plotcov,
    Y = Y,
    A = A,
    X = Risk,
    TCOV = TCOV,
    mean = mean,
    alpha = alpha,
    ng = Obj$groups,
    n = Obj$Size,
    Time = Obj$Time,
    dec = dec,
    col = col,
    degre = Obj$degre,
    method = Obj$Method,
    degre.nu = Obj$degre.nu,
    period = Obj$period,
    ...
  )
}
###################################################################################
# modification of plot's method for class trajectory.NL
####################################################################################
#' @title Plot Method for Non-Linear Trajectory Objects
#'
#' @description Plots the trajectories for a non-linear model fitted by `trajeR`.
#'
#' @param Obj An object of class `Trajectory.NL`.
#' @param Y The original response variable matrix.
#' @param A The original time variable matrix.
#' @param Risk Optional. The original risk covariate matrix.
#' @param plotcov An optional vector or matrix of time-dependent covariates.
#' @param col An optional vector of colors.
#' @inheritParams trajeR
#' @param mean Logical. If `TRUE`, plots the mean trajectory for each group.
#' @param alpha A numeric value for the transparency of individual trajectory lines.
#' @param ... Other graphical parameters passed to `plot`.
#'
#' @return a graphic.
#' @export
#'
plotrajeR.Trajectory.NL <- function(
  Obj,
  plotcov = NULL,
  col = "black",
  Y = NULL,
  A = NULL,
  Risk = NULL,
  mean = FALSE,
  alpha = 1,
  TCOV = NULL,
  ...
) {
  if (!is.null(Y)) {
    Y <- data.matrix(Y)
  }
  if (!is.null(A)) {
    A <- data.matrix(A)
  }
  plotTrajNL(
    beta = Obj$beta,
    sigma = Obj$sigma,
    theta = Obj$theta,
    plotcov = plotcov,
    Y = Y,
    A = A,
    X = Risk,
    TCOV = TCOV,
    mean = mean,
    alpha = alpha,
    ng = Obj$groups,
    n = Obj$Size,
    Time = Obj$Time,
    col = col,
    degre = Obj$degre,
    ymin = Obj$min,
    ymax = Obj$max,
    method = Obj$Method,
    fct = Obj$fct,
    ...
  )
}
###################################################################################
# modification of plot's method for class trajectory.BETA
####################################################################################
#' @title Plot Method for BETA Trajectory Objects
#'
#' @description Plots the trajectories for a Beta model fitted by `trajeR`.
#'
#' @param Obj An object of class `Trajectory.BETA`.
#' @param Y The original response variable matrix.
#' @param A The original time variable matrix.
#' @param Risk Optional. The original risk covariate matrix.
#' @param plotcov An optional vector or matrix of time-dependent covariates.
#' @param col An optional vector of colors.
#' @inheritParams trajeR
#' @param mean Logical. If `TRUE`, plots the mean trajectory for each group.
#' @param alpha A numeric value for the transparency of individual trajectory lines.
#' @param ... Other graphical parameters passed to `plot`.
#'
#' @return a graphic.
#' @export
#'
plotrajeR.Trajectory.BETA <- function(
  Obj,
  plotcov = NULL,
  col = "black",
  Y = NULL,
  A = NULL,
  Risk = NULL,
  TCOV = NULL,
  mean = FALSE,
  alpha = 1,
  ...
) {
  if (!is.null(Y)) {
    Y <- data.matrix(Y)
  }
  if (!is.null(A)) {
    A <- data.matrix(A)
  }
  plotTrajBETA(
    beta = Obj$beta,
    phi = Obj$phi,
    theta = Obj$theta,
    delta = Obj$delta,
    plotcov = plotcov,
    Y = Y,
    A = A,
    X = Risk,
    TCOV = TCOV,
    mean = mean,
    alpha = alpha,
    ng = Obj$groups,
    n = Obj$Size,
    Time = Obj$Time,
    dec = 1,
    col = col,
    degre = Obj$degre,
    method = Obj$Method,
    degre.phi = Obj$degre.phi,
    ...
  )
}