R/1TM_Qsun_fit.R

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

#' Frequentist inference for the one thermal mass model (1TM) including
#' incoming solar radiation
#'
#' Brief description
#'
#' @param data A dataset in the format of either \code{\link{cwall_east}}
#'   or \code{\link{cwall_north}}.
#' @param air A logical scalar.  Should we use air (\code{TRUE}) or
#'   surface (\code{FALSE}) temperatures?
#' @param init A numeric vector of length 4.  Optional initial estimates of
#'   the parameters R1, R2, R3 and C1 of the 1TM + Qsun.
#' @param tau A numeric scalar. Time interval between successive measurements,
#'   in hours.
#' @param ep A numeric scalar.  Lower bound imposed on the (positive)
#'   parameters R1, R2 and C1.  Passed to the \code{lower} argument of
#'   \code{\link[stats]{nls}}.
#' @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 use_gnls A logical scalar.  In the \code{hetero = FALSE}) case
#'   should we use \code{\link[nlme]{gnls}} to fit the model
#'   (\code{use_gnls == TRUE}) or \code{\link[stats]{nls}}
#'   (\code{use_gnls == FALSE})?  This is provided as a means to check that
#'   these two functions give equivalent results.
#' @param ... Further arguments to be pased to \code{\link[stats]{nls}}.
#' @details Add details
#' @return An object of class \code{"nls"} returned from
#'   \code{\link[stats]{nls}}.
#' @references Gori, V. (2017) A novel method for the estimation of
#'   thermophysical properties of walls from short and seasonally
#'   independent in-situ surveys. Doctoral thesis. UCL (University
#'   College London).
#'   \url{http://discovery.ucl.ac.uk/1568418/}
#' @examples
#' ### cwall_east --------
#'
#' # Constant error variance
#' fit3 <- one_tm_sun(data = cwall_east, hetero = FALSE)
#' fit3
#' logLik(fit3)
#' plot(fit3, form = resid(., type = "p") ~ fitted(.) | in_out)
#'
#' # Different error variances for Qin and Qout
#' fit4 <- one_tm_sun(data = cwall_east)
#' fit4
#' plot(fit4, form = resid(., type = "p") ~ fitted(.) | in_out)
#' logLik(fit4)
#' @export
one_tm_sun <- function(data, air = TRUE, init = NULL, tau = 1 / 12,
                       hetero = TRUE, ep = 1e-6, nls_control = NULL,
                       gnls_control = NULL, use_gnls = FALSE,...) {
  # Check that ep is positive
  if (ep <= 0) {
    stop("ep must be positive")
  }
  # Create a data frame: Qin data (top half) and Qout data (bottom half)
  data_df <- one_tm_sun_create_data_df(data, air, tau)
  f <- Qp ~ one_tm_sun_grad_fn(x1, x2, x3, x4, x5, Qp1, R1, R2, R3, C1,
                               in_out, tau)
  if (is.null(init)) {
    init <- c(R1 = 0.04891, R2 =  0.20447, R3 = 0.04286, C1 = 4.83928)
  }
  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))
  }
  if (hetero) {
    # Fit the model with different error variances but error correlation = 0
    sigma_ratio <- 3
    fit <- nlme::gnls(f, data = data_df, start = init,
                      weights = nlme::varIdent(c(inside = sigma_ratio ^ 2),
                                               form = ~ 1 | in_out),
                      control = gnls_control, ...)
  } else {
    if (use_gnls) {
      fit <- nlme::gnls(f, data = data_df, start = init,
                        control = gnls_control, ...)
    } else {
      fit <- stats::nls(f, data = data_df, start = init,
                        control = nls_control, ...)
    }
  }
  # Make "walls" the first class
  class(fit) <- c("walls", class(fit))
  # Manually add the data to the object returned
  fit$data <- data_df
  attr(fit$data, "which_data") <- "Qin and Qout"
  # If we used nlme::gnls() the reverse the order of the classes of the
  # object returned
  if (class(fit)[2] == "gnls") {
    class(fit) <- c("walls", "gls", "gnls")
  }
  # I have found that unless I do these things plot(object) doesn't work
  # in cases where it is necessary to access the original data frame used
  # in the call to nlme::gnls()
  # Save the name of the function, for use in summary.walls().
  fit$walls_fn <- "one_tm_sun"
  return(fit)
}