R/plot.FOBJ.R

In iantsuising/quickfun:

#'@title Visulization of FOBJ
#'@description This function is aimed to plot the functional data from different aspects
#'@param x An FOBJ, generated by \code{FOBJ}
#'@param J The number of PCs to be plotted. Only positive number is allowed
#'@method plot FOBJ
#'@import stats
#'@import graphics
#'@import grDevices
#'@export


plot.FOBJ <- function(x, J = 4){

  if(! 'FOBJ' %in% class(x)){stop("Invalid input class")}

  if(x$K > 4) cat("Only the first 4 groups will be displayed")

  if(J <= 0) {stop("Only positive number allowed")}

  graColors <- c("red", "blue", "black", "darkgreen")
  PointType <- c(3, 17, 5, 0)

  yesorno <- readline(prompt = "To continue hit return/Type n to break")
  if(yesorno == "n"){return("End")}
  else{

      plot(x$t, x$DataGroupby[[1]][1, ],
           main = "Functional Data Curve",
           xlab = " ", ylab = " ",
           ylim = c(min(unlist(x$DataGroupby)) - 0.2 * abs(min(unlist(x$DataGroupby))),
                    max(unlist(x$DataGroupby)) + 0.2 * abs(max(unlist(x$DataGroupby)))),
           type = "l", lty = 1)

      for(i in 1:min(x$K, 4)){
        apply(x$DataGroupby[[i]], 1, FUN = function(y){lines(x$t, y, col = graColors[i], lty = i)})
      }
      legend("topleft", legend = x$y[1:min(x$K, 4)], col = graColors[1:min(x$K, 4)], lty = 1:min(x$K, 4))
  } #Original data

  yesorno <- readline(prompt = "To continue hit return/Type n to break")
  if(yesorno == "n"){return("End")}
  else{
      plot(x$t, apply(x$DataGroupby[[1]], 2, mean),
           main = "Mean curve of each group",
           xlab = " ", ylab = " ",
           ylim = c(min(unlist(x$DataGroupby)) - 0.2 * abs(min(unlist(x$DataGroupby))),
                    max(unlist(x$DataGroupby)) + 0.2 * abs(max(unlist(x$DataGroupby)))),
           type = "l", lty = 1, col = graColors[1])
      for(i in 2:min(x$K, 4)){
        lines(x$t, apply(x$DataGroupby[[i]], 2, mean), col = graColors[i], lty = i)
      }
      legend("topleft", legend = x$y[1:min(x$K, 4)], col = graColors[1:min(x$K, 4)], lty = 1:min(x$K, 4))
   } #Mean curves



  yesorno <- readline(prompt = "To continue hit return/Type n to break")
  if(yesorno == "n"){return("End")}
  else{
          autofit1 <- abs(mean(x$t) / mean(x$poolCov))
          z1 <- autofit1 * x$poolCov
          lim1 <- range(z1)
          len1 <- lim1[2] - lim1[1] + 1
          colorset1 <- terrain.colors(len1, alpha = 0.5)
          color1 <- colorset1[z1 - lim1[1] + 1]
          persp(5 * x$t, 5 * x$t, z1, xlab = "t", ylab =  "s", zlab = " ", col = color1, main = "Pooled Covariance", theta = 45)

  } #Pooled covariance surface
  if(x$K == 2){
    yesorno <- readline(prompt = "To continue hit return/Type n to break")
    if(yesorno == "n"){return("End")}
    else{
    autofit2 <- abs(mean(x$t) / mean(x$singleCov[[1]] - x$singleCov[[2]]))
    z2 <- autofit2 * (x$singleCov[[1]] - x$singleCov[[2]])
    lim2 <- range(z2)
    len2 <- lim2[2] - lim2[1] + 1
    colorset2 <- terrain.colors(len2, alpha = 0.5)
    color2 <- colorset2[z2 - lim2[1] + 1]
    persp(5 * x$t, 5 * x$t, z2, xlab = "t", ylab =  "s", zlab = " ", col = color1, main = "Difference Between Two Covariance Functions", theta = 45)
    }
  }




  for(i in 1:min(J,x$cutPoints)){
    yesorno <- readline(prompt = "Hit return to continue/Type n to break")
    if(yesorno == "n"){return("End")}
    else{
      plx <- numeric()
      ply <- numeric()
      for(ii in 1:min(4, x$K)){
        plx <- append(plx, density(x$scores[[ii]][, i])$x)
        ply <- append(ply, density(x$scores[[ii]][, i])$y)
      }
      plot(density(x$scores[[1]][, i])$x, density(x$scores[[1]][, i])$y,
           type = "l", lty = 1, col = graColors[1], xlab = " ", ylab = " ",
           xlim = range(plx), ylim = range(ply))
      for(j in 2:min(4, x$K)){
        lines(density(x$scores[[j]][, i])$x, density(x$scores[[j]][, i])$y,
              col = graColors[j], lty = j)
      }
      legend("topleft", legend = x$y[1:min(4, x$K)],
             col = graColors[1:min(4, x$K)], lty = 1:min(4, x$K))
      title(main = paste0("Density of scores on FPC", as.character(i), " ", "(", as.character(round(x$proportion[i], 3)) ,")"))
    }
  } #Density of FPCs


  for(i in 1:min(J, x$cutPoints)){
          yesorno <- readline(prompt = "Hit return to continue/Type n to break")
          if(yesorno == "n"){return("End")}
          else{
          theta = seq(0.1, 0.9, 0.1)
          ploty <- numeric()
          for(ii in 1:min(4, x$K)){
            ploty <- append(ploty, x$scores[[ii]][, i])
          }
          plot(theta, quantile(x$scores[[1]][, i], theta), col = graColors[1],
               ylim = range(ploty),
               pch = PointType[1],
               xlab = " ", ylab = " ")
          for(j in 2:min(4, x$K)){
          points(theta, quantile(x$scores[[2]][, i], theta), col = graColors[j], pch = PointType[j])
          }

          legend("topleft", legend = x$y[1:min(4, x$K)],
                 col = graColors[1:min(4, x$K)],
                 pch = PointType[1:min(4, x$K)])
          title(main = paste0("Quantile of scores on FPC", as.character(i), " ",
                              "(", as.character(round(x$proportion[i], 3)),")"))
          }
        } #Quantile of FPCs


  yesorno <- readline(prompt = "To continue hit return/Type n to break")
  if(yesorno == "n"){return("End")}
  else{
            toplotx <- 1:min(J, x$cutPoints)
            plotymean <- numeric()
            for(i in 1:min(4, x$K)){
              plotymean <- append(plotymean, apply(x$scores[[i]][, toplotx], 2, mean))
            }

            plot(toplotx, apply(x$scores[[1]][, toplotx], 2, mean), type = "l", lty = 1,
                 col = graColors[1], xlab = " ", ylab = " ", ylim = range(plotymean),
                 main = "The mean of projection scores from each group",
                 xaxt = "n")
            for(j in 2:min(4, x$K)){
              lines(toplotx, apply(x$scores[[j]][, toplotx], 2, mean), lty = j, col = graColors[j])
            }
            legend("topright", legend = x$y[1:min(4, x$K)], col = graColors[1:min(4, x$K)], pch = rep(1, min(4, x$K)))
            axis(1, 1:min(J, x$cutPoints))
  }  #Mean of FPCs

  yesorno <- readline(prompt = "To continue hit return/Type n to break")
  if(yesorno == "n"){return("End")}
  else{
            plotyvar <- numeric()
            for(i in 1:min(4, x$K)){
              plotyvar <- append(plotyvar, apply(x$scores[[i]][, toplotx], 2, var))
            }

            plot(toplotx, apply(x$scores[[1]][, toplotx], 2, var), type = "l", lty = 1,
                 col = graColors[1], xlab = " ", ylab = " ", ylim = range(plotyvar),
                 main = "The variance of projection scores from each group",
                 xaxt = "n")
            for(j in 2:min(4, x$K)){
              lines(toplotx, apply(x$scores[[j]][, toplotx], 2, var), lty = j, col = graColors[j])
            }
            legend("topright", legend = x$y[1:min(4, x$K)], col = graColors[1:min(4, x$K)], pch = rep(1, min(4, x$K)))
            axis(1, 1:min(J, x$cutPoints))
} #Var of FPCs

  for(i in 1:min(J, x$cutPoints)){
    yesorno <- readline(prompt = "Hit return to continue/Type n to break")
    if(yesorno == "n"){return("End")}
    else{
      plot(x$t, x$phi[, i], type = "l", col = "blue",
           xlab = " ", ylab = " ")
      title(main = paste0("FPC", as.character(i), "(", as.character(round(x$proportion[i], 3)) ,")"))
    }
  } #Plot of FPCs

}