R/sample_diagnostics.R

In paulnorthrop/walls: Statistical Inference for the Thermophysical Properties of Walls from In-Situ Surveys

#' Summarise a walls object
#'
#' \code{summary} method for class "walls".
#'
#' @param object an object of class "walls", a result of a call to
#'   \code{\link{one_tm}}, \code{\link{one_tm_sun}} or \code{\link{two_tm}}.
#' @param param A character scalar indicating which parameterization should
#'   be used in the results for the thermal resistance parameters:
#'   \itemize{
#'     \item{"original": }{(R1, R2, R3)}
#'     \item{"R": }{(sumR = R1 + R2 + R3, fracR1 = R1 / sumR,
#'       fracR2 = R2 / sumR)}
#'     \item{"U": }{(Uvalue = 1 / sumR, fracR1 = R1 / sumR,
#'       fracR2 = R2 / sumR)}
#'  }
#' @param nls_control An optional list of control settings to be passed to
#'   \code{\link[stats]{nls}}.
#' @param gnls_control An optional list of control settings to be passed to
#'   \code{\link[nlme]{gnls}}.
#' @param ... Additional arguments passed on to the \code{\link{print}}.
#' @details Note: if \code{param = "R"} or \code{param = "U"} then
#'   the model is refitted under the new parameterization and therefore may
#'   take a few seconds.
#' @examples
#' fit2 <- two_tm(data = cwall_east)
#' summary(fit2)
#' summary(fit2, param = "R")
#' summary(fit2, param = "U")
#'
#' fit1s <- one_tm_sun(data = cwall_east)
#' summary(fit1s)
#' summary(fit1s, param = "R")
#' summary(fit1s, param = "U")
#' @seealso \code{\link{one_tm}}, \code{\link{one_tm_sun}} and
#'   \code{\link{two_tm}}
#' @export
summary.walls <- function(object, param = c("original", "R", "U"),
                          nls_control = NULL, gnls_control = NULL, ...) {
  if (!inherits(object, "walls")) {
    stop("use only with \"walls\" objects")
  }
  param <- match.arg(param)
  # Is the (secondary) class "nls" or "gls"
  class_object <- class(object)[2]
  # Get the S3 print method
  summary_fn <- getS3method("summary", class_object)
  # If param = "original" then call the nls() or gnls() method, as appropriate
  # Otherwise, refit the model under the required parameterization
  if (param == "original") {
    object <- summary_fn(object)
    class(object) <- paste("summary.", class_object, sep = "")
    return(object)
  }
  # Use the same default control parameters as the fitting functions
  if (is.null(gnls_control$nlsTol)) {
    gnls_control <- c(gnls_control, list(nlsTol = 100))
  }
  if (is.null(nls_control$nlsTol)) {
    nls_control <- c(nls_control, list(maxiter = 200))
  }
  # Extract the (original) parameter estimates
  if (class_object == "nls") {
    ests <- object$m$getAllPars()
  } else {
    ests <- object$coefficients
  }
  if (param == "R") {
    if (object$walls_fn == "two_tm") {
      sumR <- sum(ests[c("R1", "R2", "R3")])
      fracR1 <- ests["R1"] / sumR
      fracR2 <- ests["R2"] / sumR
      init <- c(sumR, fracR1, fracR2, ests[c("C1", "C2")])
      names(init) <- c("sumR", "fracR1", "fracR2", "C1", "C2")
      f <- Qp ~ R_two_tm_grad_fn(x1, x2, x3, Qp1, sumR, fracR1, fracR2, C1,
                                 C2, in_out, tau)
      if (class_object == "nls") {
        fit <- stats::nls(f, data = object$data, start = init, ...)
      } else {
        sigma_ratio <- coef(object$modelStruct$varStruct, uncons = FALSE,
                            allCoef = FALSE)
        fit <- nlme::gnls(f, data = object$data, start = init,
                          weights = nlme::varIdent(c(inside = sigma_ratio ^ 2),
                                                   form = ~ 1 | in_out),
                          control = gnls_control, ...)
      }
      object <- summary_fn(fit)
    } else if (object$walls_fn == "one_tm_sun") {
      sumR <- sum(ests[c("R1", "R2", "R3")])
      fracR1 <- ests["R1"] / sumR
      fracR2 <- ests["R2"] / sumR
      init <- c(sumR, fracR1, fracR2, ests["C1"])
      names(init) <- c("sumR", "fracR1", "fracR2", "C1")
      f <- Qp ~ R_one_tm_sun_grad_fn(x1, x2, x3, x4, x5, Qp1, sumR, fracR1,
                                     fracR2, C1, in_out, tau)
      if (class_object == "nls") {
        fit <- stats::nls(f, data = object$data, start = init, ...)
      } else {
        sigma_ratio <- coef(object$modelStruct$varStruct, uncons = FALSE,
                            allCoef = FALSE)
        fit <- nlme::gnls(f, data = object$data, start = init,
                          weights = nlme::varIdent(c(inside = sigma_ratio ^ 2),
                                                   form = ~ 1 | in_out),
                          control = gnls_control, ...)
      }
      object <- summary_fn(fit)
    }
  } else if (param == "U") {
    if (object$walls_fn == "two_tm") {
      sumR <- sum(ests[c("R1", "R2", "R3")])
      Uvalue <- 1 / sumR
      fracR1 <- ests["R1"] / sumR
      fracR2 <- ests["R2"] / sumR
      init <- c(Uvalue, fracR1, fracR2, ests[c("C1", "C2")])
      names(init) <- c("Uvalue", "fracR1", "fracR2", "C1", "C2")
      f <- Qp ~ U_two_tm_grad_fn(x1, x2, x3, Qp1, Uvalue, fracR1, fracR2, C1,
                                 C2, in_out, tau)
      if (class_object == "nls") {
        fit <- stats::nls(f, data = object$data, start = init, ...)
      } else {
        sigma_ratio <- coef(object$modelStruct$varStruct, uncons = FALSE,
                            allCoef = FALSE)
        fit <- nlme::gnls(f, data = object$data, start = init,
                          weights = nlme::varIdent(c(inside = sigma_ratio ^ 2),
                                                   form = ~ 1 | in_out),
                          control = gnls_control, ...)
      }
      object <- summary_fn(fit)
    } else if (object$walls_fn == "one_tm_sun") {
      sumR <- sum(ests[c("R1", "R2", "R3")])
      Uvalue <- 1 / sumR
      fracR1 <- ests["R1"] / sumR
      fracR2 <- ests["R2"] / sumR
      init <- c(Uvalue, fracR1, fracR2, ests["C1"])
      names(init) <- c("Uvalue", "fracR1", "fracR2", "C1")
      f <- Qp ~ U_one_tm_sun_grad_fn(x1, x2, x3, x4, x5, Qp1, Uvalue, fracR1,
                                     fracR2, C1, in_out, tau)
      if (class_object == "nls") {
        fit <- stats::nls(f, data = object$data, start = init, ...)
      } else {
        sigma_ratio <- coef(object$modelStruct$varStruct, uncons = FALSE,
                            allCoef = FALSE)
        fit <- nlme::gnls(f, data = object$data, start = init,
                          weights = nlme::varIdent(c(inside = sigma_ratio ^ 2),
                                                   form = ~ 1 | in_out),
                          control = gnls_control, ...)
      }
      object <- summary_fn(fit)
    }
  }
  class(object) <- paste("summary.", class_object, sep = "")
  return(object)
}

#' Print a walls object
#'
#' \code{print} method for class "walls".
#'
#' @param x an object of class "walls", a result of a call to
#'   \code{\link{one_tm}}, \code{\link{one_tm_sun}} or \code{\link{two_tm}}.
#' @param ... Additional arguments.
#' @examples
#' in_air <- one_tm(data = cwall_east)
#' summary(in_air)
#' @seealso \code{\link{one_tm}}, \code{\link{one_tm_sun}} and
#'   \code{\link{two_tm}}
#' @export
print.walls <- function(x, ...) {
  if (!inherits(x, "walls")) {
    stop("use only with \"walls\" objects")
  }
  # Replace the data with a description of which data were used, i.e.
  # Qin only, Qout only or Qin and Qout.
  # Otherwise, print.nls() will print all the data, which we don't want
  # and print.gls will just give the name of the dataframe in which the
  # data sat, which is not very informative
  #
  # Is the (secondary) class "nls" or "gls"
  class_x <- class(x)[2]
  # Get the S3 print method
  print_fn <- getS3method("print", class_x)
  # Extract the description of the response data
  describe_data <- as.symbol(attr(x$data, "which_data"))
  # Save the data to be replaced and replace this with the description
  if (class_x == "gls") {
    save_data <- x$call$data
    x$call$data <- describe_data
  } else if (class_x == "nls") {
    save_data <- x$data
    x$data <- describe_data
  }
  # Produce the print
  print_fn(x)
  # Restore the data in the object to it's original state
  if (class_x == "gls") {
    x$call$data <- save_data
  } else if (class_x == "nls")  {
    x$data <- save_data
  }
  invisible(x)
}