R/myfitbetaplot.R

In TimeVarConcurrentModel: Concurrent Multivariate Models with Time-Varying Coefficients

###########################
#### Plotting Function ####
###########################

myfitbetaplot = function(betaestlist, betastderrlist=NULL,covnames, argvals=NULL,
                       xlab="", ...)
{
  #  PLOTBETA plots a functional parameter along with confidence
  #  limits
  #  Arguments
  #  BETAESTLIST    ... A list object containing one or more functional
  #                     parameter objects or functional data objects.
  #  BETASTDERRLIST ... A list object containing functional data objects
  #                     for the standard error of the objects in
  #                     BETAESTLIST.
  
  #  Last modified 6 January 2020
  
  #  check BETAESTLIST
  
  ##This function has been modified from the original fda package.
  
  
  if (inherits(betaestlist, "fdPar") || inherits(betaestlist, "fd")) {
    betaestlist = list(betaestlist)
  }
  
  if (!inherits(betaestlist, "list")) {
    stop("BETAESTLIST is not a list, fd, or fdpar object.")
  }
  
  #  check BETASTDERRLIST
  
  if (!is.null(betastderrlist)){
    if (inherits(betastderrlist, "fd")) {
      betastderrlist = list(betastderrlist)
    }
    if (!inherits(betastderrlist, "list")) {
      stop("BETASTDERRLIST is not a list, or fd object.")
    }
  }
  
  #  get range
  
  if (is.fdPar(betaestlist[[1]])) {
    rangeval = betaestlist[[1]]$fd$basis$rangeval
  } else {
    if (is.fd(betaestlist[[1]])) {
      rangeval = betaestlist[[1]]$basis$rangeval
    } else {
      stop(paste("A list does not contain either a functional parameter ",
                 "or a functional data object."))
    }
  }
  
  if (is.null(argvals)) {
    argvals = seq(rangeval[1],rangeval[2],len=51)
  }
  n = length(argvals)
  p = length(betaestlist)
  
  par(ask=T)
  for (j in 1:p) {
    if (is.fdPar(betaestlist[[j]])) {
      betavec = eval.fd(argvals, betaestlist[[j]]$fd)
    } else {
      if (is.fd(betaestlist[[j]])) {
        betavec = eval.fd(argvals, betaestlist[[j]])
      } else {
        stop(
          "BETAESTLIST does not contain a functional parameter or data object.")
      }
    }
    zerofunc = argvals*0
    zeroval  = c(0,0)
    if (is.null(betastderrlist)) {
      plot(argvals, betavec, type="l", xlab=xlab, ylab="",
           xlim=rangeval, ylim=c(min(betavec),max(betavec)), ...)
      lines(rangeval, zeroval, col=1, lty=3)
    }
    else {
      betastderr = eval.fd(argvals, betastderrlist[[j]])
      betavecp   = betavec + 2*betastderr
      betavecm   = betavec - 2*betastderr
      
      if(j==1){
        plot(argvals, betavec, type="l", xlab="t", ylab="Beta_0(t)",
             xlim=rangeval, ylim=c(min(betavecm),max(betavecp)), ...)
        lines(rangeval, zeroval,lty=3, col=2)
        lines(argvals, betavecp, col=1, lwd=1)
        lines(argvals, betavecm, col=1, lwd=1)
      }
      
      else{
        plot(argvals, betavec, type="l", xlab="t", ylab=paste("Beta_",covnames[j-1]," (t)",sep=""),
             xlim=rangeval, ylim=c(min(betavecm),max(betavecp)), ...)
        lines(rangeval, zeroval,lty=3, col=2)
        lines(argvals, betavecp, col=1, lwd=1)
        lines(argvals, betavecm, col=1, lwd=1)
        for(jj in 1:(length(betavecm)-1)){
          if(betavecm[jj] > 0 & betavecm[jj+1] > 0){
            xmin= argvals[jj]
            xmax= argvals[jj+1]
            polygon(c(argvals[argvals>=xmin & argvals<=xmax],rev(argvals[argvals>=xmin & argvals<=xmax])),
                    c(zerofunc[argvals>=xmin & argvals<=xmax],rev(betavecm[argvals>=xmin & argvals<=xmax])),col="#FF3030",density = 10, angle = 45)  
          }
        }  
        for(jj in 1:(length(betavecm)-1)){
          if(betavecp[jj] < 0 & betavecp[jj+1] < 0){
            xmin= argvals[jj]
            xmax= argvals[jj+1]
            polygon(c(argvals[argvals>=xmin & argvals<=xmax],rev(argvals[argvals>=xmin & argvals<=xmax])),                      
                    c(zerofunc[argvals>=xmin & argvals<=xmax],rev(betavecp[argvals>=xmin & argvals<=xmax])),col="#FF3030",density = 10, angle = 45)  
          }
        }
        
      }
      
    }
    if(j==1){title(paste("Intercept"))}
    else{title(paste(covnames[j-1]))}
  }
  
}