R/S3_ordgam_plot.R

In ordgam: Additive Model for Ordinal Data using Laplace P-Splines

#' Plot the the additive terms in an <ordgam> object with its credible regions
#'
#' @description Plot the the additive terms in an object generated by \code{ordgam} with its credible regions
#'
#' @usage
#' \method{plot}{ordgam}(x,ngrid=300,ci.level=.95,mfrow=NULL,...)
#'
#' @param x An \link{ordgam.object} generated by \link{ordgam}.
#' @param ngrid An integer indicating the number of gridpoints where the additive terms should be evaluated.
#' @param ci.level Credibility level for the pointwise credible region of the additive terms.
#' @param mfrow (Optional) A vector of the form c(nr, nc). Subsequent figures will be drawn in an nr-by-nc array on the device by rows.
#' @param ... Additional generic plotting arguments.
#'
#' @return In addition to the plots, an invisible object containing information on the estimated additive terms is returned, see the \code{ordgam_additive} function documentation for more details.
#'
#' @author Philippe Lambert \email{p.lambert@uliege.be}
#' @references
#' Lambert, P. and Gressani, 0. (2023) Penalty parameter selection and asymmetry corrections
#' to Laplace approximations in Bayesian P-splines models. Statistical Modelling. <doi:10.1177/1471082X231181173>. Preprint: <arXiv:2210.01668>.
#'
#' @seealso \code{\link{ordgam}}, \code{\link{ordgam.object}}, \code{ordgam.additive}, \code{\link{print.ordregr}}.
#'
#' @examples
#' library(ordgam)
#' data(freehmsData)
#' mod = ordgam(freehms ~ gndr + s(eduyrs) + s(age),
#'              data=freehmsData, descending=TRUE)
#' print(mod)
#' plot(mod)
#'
#' @export
plot.ordgam = function(x, ngrid=300, ci.level=.95, mfrow=NULL,...){
    obj = x
    fhat = ordgam_additive(obj,ngrid=ngrid,ci.level=ci.level) ## Compute additive term + envelope
    if (fhat$J == 0) return(NULL)
    ## Plotting window division
    if (is.null(mfrow)) mfrow=c(1,1)
    maxPlts = prod(mfrow)
    ## Make sure to restore the graphics parameters as they were before calling the function.
    oldpar <- par(no.readonly = TRUE) ; on.exit(par(oldpar))
    dev.new() ; par(mar=c(4,4,1,1),mfrow=mfrow)
    for (j in 1:fhat$J){
        if ((j%/%(maxPlts+1) == 1)) { ## New plotting window if ...
          dev.new() ; par(mar=c(4,4,1,1),mfrow=mfrow)
        }
        lab = names(fhat$f.grid)[j]
        with(fhat$f.grid[[j]], matplot(x, y.mat,type="l",
                                       xlab=lab,ylab=paste("f(",lab,")",sep=""),
                                       lwd=c(2,1,1),lty=c(1,2,2),col=1, ...))
    }
    return(invisible(fhat))
}