R/plotfit.fd.R

In fda: Functional Data Analysis

#plotfit <- function (x, ...){
#  UseMethod("plotfit")
#}

plotfit.fdSmooth <- function(y, argvals, fdSm, rng = NULL,
                       index = NULL, nfine = 101, residual = FALSE,
                       sortwrd = FALSE, titles=NULL,  ylim=NULL,
                       ask=TRUE, type=c("p", "l")[1+residual],
                       xlab=NULL, ylab=NULL, sub=NULL, col=1:9,
                       lty=1:9, lwd=1, cex.pch=1, axes=NULL, ...) {
  plotfit.fd(y, argvals, fdSm$fd, rng = rng, index = index,
             nfine = nfine, residual = residual,
             sortwrd = sortwrd, titles=titles,  ylim=ylim,
             ask=ask, type=c("p", "l")[1+residual],
             xlab=xlab, ylab=ylab, sub=sub, col=1:9, lty=1:9,
             lwd=1, cex.pch=1, axes=axes, ...)
}

plotfit.fd <- function(y, argvals, fdobj, rng = NULL,
                       index = NULL, nfine = 101, residual = FALSE,
                       sortwrd = FALSE, titles=NULL,  ylim=NULL,
                       ask=TRUE, type=c("p", "l")[1+residual],
                       xlab=NULL, ylab=NULL, sub=NULL, col=1:9, lty=1:9,
                       lwd=1, cex.pch=1, axes=NULL, ...)
{
  #PLOTFIT plots discrete data along with a functional data object for
  #  fitting the data.  It is designed to be used after DATA2FD,
  #  SMOOTH.FD or SMOOTH.BASIS to check the fit of the data offered by
  #  the FD object.
  
  #  Arguments:
  #  Y        ... the data used to generate the fit
  #  ARGVALS  ... discrete argument values associated with data
  #  fdobj       ... a functional data object for fitting the data
  #  RNG      ... a range of argument values to be plotted
  #  INDEX    ... an index for plotting subsets of the curves
  #       (either sorted or not)
  #  NFINE    ... number of points to use for plotting curves
  #  RESIDUAL ... if TRUE, the residuals are plotted instead of
  #         the data plus curve
  #  SORTWRD  ... sort plots by mean square error
  #  TITLES   ... vector of title strings for curves
  
  # Last modified 18 January 2020 by Jim Ramsay
  
  ##
  ## 1.  Basic checks
  ##
  #  check third argument, FDOBJ
  
  if (!(is.fd(fdobj) || is.fdPar(fdobj)))  stop(
    "Third argument is neither a functional data nor a functional parameter object.")
  if (is.fdPar(fdobj)) fdobj <- fdobj$fd
  
  #  extract basis object from FDOBJ
  basisobj <- fdobj$basis
  #  set range of arguments
  if (is.null(rng)) rng <- basisobj$rangeval
  #  extract FDNAMES from FDOBJ
  fdnames <- fdobj$fdnames
  yName   <- substring(deparse(substitute(y)), 1, 33)
  fdName  <- paste(substring(deparse(substitute(fdobj)), 1, 22),
                   "$coef", sep='')
  ##
  ## 2.  Use 'checkDims' to reconcile y and fdoj$coef
  ##
  #  The default fdnames may not work well
  defaultNms <- c(fdnames[2], fdnames[3], x='x')
  if((length(defaultNms[[2]])<2) && !is.null(names(defaultNms))
     && !is.na(names(defaultNms)[2]))
    defaultNms[[2]] <- names(defaultNms)[2]
  subset <- checkDims3(y, fdobj$coef, defaultNames = defaultNms,
                       xName=yName, yName=fdName)
  y <- subset$x
  fdobj$coef <- subset$y
  #  number of argument values
  n <- dim(y)[1]
  #  number of replications
  nrep <- dim(y)[2]
  #  number  of variables
  nvar <- dim(y)[3]
  curveno <- 1:nrep
  #  names for argument values
  argname <- names(fdnames)[[1]]
  if (is.null(argname)) {
    if (is.null(fdnames[[1]])) argname <- "Argument Value"
    else                       argname <- fdnames[[1]]
  }
  if (is.null(xlab)) xlab <- argname
  #  names for replicates
  casenames <- dimnames(y)[[2]]
  #  names for variables
  varnames  <- dimnames(y)[[3]]
  if (is.null(axes)) {
    if (is.null(fdobj$basis$axes)) {
      Axes  <- TRUE
      axFun <- FALSE
    } else {
      if (!inherits(fdobj$basis$axes, "list"))
        stop("fdobj$basis$axes must be a list;  ",
             "class(fdobj$basis$axes) = ", class(fdobj$basis$axes))
      if (!(inherits(fdobj$basis$axes[[1]], "character") ||
            inherits(fdobj$basis$axes[[1]], "function")))
        stop("fdobj$basis$axes[[1]] must be either a function or the ",
             "name of a function;  class(fdobj$basis$axes[[1]]) = ",
             class(fdobj$basis$axes[[1]]))
      Axes   <- FALSE
      axFun  <- TRUE
      axList <- c(fdobj$basis$axes, ...)
    }
  } else {
    if (is.logical(axes)) {
      Axes  <- axes
      axFun <- FALSE
    }
    else {
      if (!inherits(axes, "list"))
        stop("axes must be a logical or a list;  class(axes) = ",
             class(axes))
      if (!(inherits(axes[[1]], "character") || inherits(axes[[1]],
                                                         "function")))
        stop("axes[[1]] must be either a function or the ",
             "name of a function;  class(axes[[1]]) = ",
             class(axes[[1]]))
      Axes <- FALSE
      axFun <- TRUE
      axList <- c(axes, ...)
    }
  }
  dots <- list(...)
  if (is.null(titles) && ("main" %in% names(dots)))
    titles <- dots$main
  ##
  ## 3.  Computed fitted values for argvals and a fine mesh
  ##
  yhat.  <- eval.fd(argvals, fdobj)
  yhat   <- as.array3(yhat.)
  res    <- y - yhat
  MSE    <- apply(res^2,c(2,3),mean)
  dimnames(MSE) <- list(casenames, varnames)
  MSEsum <- apply(MSE,1,sum)
  #  compute fitted values for fine mesh of values
  xfine  <- seq(rng[1], rng[2], len=nfine)
  yfine  <- array(eval.fd(xfine, fdobj),c(nfine,nrep,nvar))
  #  sort cases by MSE if desired?????
  if (sortwrd && nrep > 1) {
    MSEind <- order(MSEsum)
    y      <- y    [,MSEind,, drop=FALSE]
    yhat   <- yhat [,MSEind,, drop=FALSE]
    yfine  <- yfine[,MSEind,, drop=FALSE]
    res    <- res  [,MSEind,, drop=FALSE]
    MSE    <- MSE  [ MSEind,, drop=FALSE]
    casenames  <- casenames[MSEind]
    titles     <- titles[MSEind]
  }
  
  #  set up fit and data as 3D arrays, selecting curves in INDEX
  
  if (is.null(index)) index <- 1:nrep
  #
  nrepi <- length(index)
  y     <- y    [,index,, drop=FALSE]
  yhat  <- yhat [,index,, drop=FALSE]
  res   <- res  [,index,, drop=FALSE]
  yfine <- yfine[,index,, drop=FALSE]
  MSE   <- MSE  [ index,, drop=FALSE]
  casenames <- casenames[index]
  titles    <- titles[index]
  
  # How many plots on a page?
  nOnOne <- 1   #  only one plot is default
  # types of plots
  if (ask & ((nvar*nrepi/nOnOne) > 1))
    cat('Multiple plots:  Click in the plot to advance to the next plot\n')
  col <- rep(col, length=nOnOne)
  lty <- rep(lty, length=nOnOne)
  lwd <- rep(lwd, length=nOnOne)
  cex.pch <- rep(cex.pch, length=nOnOne)
  
  #  select values in ARGVALS, Y, and YHAT within RNG
  
  argind    <- ((argvals >= rng[1]) & (argvals <= rng[2]))
  argvals   <- argvals[argind]
  casenames <- casenames[argind]
  y         <- y   [argind,,, drop=FALSE]
  yhat      <- yhat[argind,,, drop=FALSE]
  res       <- res [argind,,, drop=FALSE]
  n         <- length(argvals)
  
  xfiind <- ((xfine >= rng[1]) & (xfine <= rng[2]))
  xfine  <- xfine[xfiind]
  yfine  <- yfine[xfiind,,, drop=FALSE]
  nfine  <- length(xfine)
  ##
  ## 4.  Plot the results
  ##
  ndigit = abs(floor(log10(min(c(MSE)))) - 1)
  if (is.null(sub))
    sub <- paste("  RMS residual =", round(sqrt(MSE),ndigit))
  if (length(sub) != length(MSE)){
    warning('length(sub) = ', length(sub), ' != ',
            length(MSE), ' = ', length(MSE), '; forcing equality')
    sub <- rep(sub, length=length(MSE))
  }
  if (is.null(dim(sub))){
    dim(sub) <- dim(MSE)
  }
  if (!all(dim(sub)==dim(MSE))) {
    warning('dim(sub) = ', dim(sub), " != dim(MSE) = ",
            paste(dim(MSE), collapse=', '), ';  forcing equality.')
    dim(sub) <- dim(MSE)
  }
  # 'ask' is controlled by 'nOnOne' ...
  op <- par(ask=FALSE)
  # Don't ask for the first plot,
  # but if ask==TRUE, do ask for succeeding plots
  on.exit(par(op))
  # Reset 'ask' on.exit (normal or otherwise)
  iOnOne <- 0
  if (residual) {
    #  plot the residuals
    if(is.null(ylim))ylim <- range(res)
    if(is.null(ylab))ylab=paste('Residuals for', varnames)
    for (j in 1:nvar) {
      for (i in 1:nrepi){
        if(iOnOne %in% c(0, nOnOne)[1:(1+ask)]){
          if(is.null(xlab))xlab <- argname
          plot(rng, ylim, type="n", xlab=xlab,
               ylab=ylab[j], axes=Axes, ...)
          if(axFun)
            do.call(axList[[1]], axList[-1])
          par(ask=ask)
          abline(h=0, lty=4, lwd=2)
          if(nOnOne==1){
            {
              if (is.null(titles))
                title(main=casenames[i],sub=sub[i, j])
              else title(main=titles[i], sub=sub[i, j])
            }
          }
          iOnOne <- 0
        }
        iOnOne <- iOnOne+1
        lines(argvals, res[,i,j], type=type,
              col=col[iOnOne], lty=lty[iOnOne],
              lwd=lwd[iOnOne], cex=cex.pch[iOnOne])
      }
    }
  } else {
    #  plot the data and fit
    if(is.null(ylim))ylim <- range(c(c(y),c(yfine)))
    varname <- names(fdnames)[[3]]
    if (is.null(varname)) {
      if (is.null(fdnames[[3]])) rep("Function Value", nvar)
      else                       varname <- fdnames[[3]]
    }
    if(is.null(ylab))ylab <- varname
    if (nvar == 1) {
      for (i in 1:nrepi) {
        if(iOnOne %in% c(0, nOnOne)[1:(1+ask)]){
          plot(rng, ylim, type="n", xlab=xlab,
               ylab=ylab, axes=Axes, ...)
          if(axFun)
            do.call(axList[[1]], axList[-1])
          par(ask=ask)
          if(nOnOne==1){
            if(is.null(titles)) title(main=casenames[i],
                                      sub=sub[i, 1])
            else title(main=titles[i], sub=sub[i, 1])
          }
          iOnOne <- 0
        }
        iOnOne <- iOnOne+1
        points(argvals, y[,i,1], type=type,
               xlab=xlab, ylab=varnames, col=col[iOnOne],
               lty=lty[iOnOne], lwd=lwd[iOnOne],
               cex=cex.pch[iOnOne])
        lines(xfine, yfine[,i,1], col=col[iOnOne],
              lty=lty[iOnOne], lwd=lwd[iOnOne])
      }
    } else {
      if (ask) {
        aski = FALSE
        for (i in 1:nrepi) {
          par(mfrow=c(nvar,1),ask=aski)
          for (j in 1:nvar) {
            plot(rng, ylim, type="n", xlab=xlab,ylab=varnames[j], axes=Axes)
            if (axFun) do.call(axList[[1]], axList[-1])
            if (j == 1) {
              if (is.null(titles)) {
                title(paste("Record",i))
              } else {
                title(main=titles[i], sub=sub[i, j])
              }
            }
            points(argvals, y[,i,j], xlab=xlab, ylab=varnames[j])
            lines(xfine, yfine[,i,j])
          }
          aski = TRUE
        }
      } else {
        askj <- FALSE
        for (j in 1:nvar) {
          par(mfrow=c(1,1),ask=askj)
          matplot(argvals, y[,,j], type="p", pch="o", ylim=ylim, xlab=xlab,
                  ylab=varnames[j])
          if (axFun) do.call(axList[[1]], axList[-1])
          matlines(xfine, yfine[,,j])
          askj <- TRUE
        }
      }
    }
  }
  invisible(NULL)
}