R/2TM_fit.R

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

#' Frequentist inference for the two thermal mass model (2TM)
#'
#' 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 5.  Optional initial estimates of
#'   the parameters R1, R2, R3, C1 and C2 of the 2TM.
#' @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
#' fit1 <- two_tm(data = cwall_east, hetero = FALSE)
#' fit1
#' logLik(fit1)
#' plot(fit1, form = resid(., type = "p") ~ fitted(.) | in_out)
#'
#' # Different error variances for Qin and Qout
#' fit2 <- two_tm(data = cwall_east)
#' fit2
#' plot(fit2, form = resid(., type = "p") ~ fitted(.) | in_out)
#' logLik(fit2)
#' @export
two_tm <- 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 <- two_tm_create_data_df(data, air, tau)
  f <- Qp ~ two_tm_grad_fn(x1, x2, x3, Qp1, R1, R2, R3, C1, C2, in_out, tau)
  # Initial estimates from separate 1TM fits to Qin and Qout data
  in_one_tm <- one_tm(data, flux = "in")
  out_one_tm <- one_tm(data, flux = "out")
  if (is.null(init)) {
    in_ests <- in_one_tm$m$getPars()
    out_ests <- out_one_tm$m$getPars()
    R1in <- in_ests["R1"]
    R2in <- min(in_ests["R2"], out_ests["R1"])
    R3in <- out_ests["R2"]
    C1in <- in_ests["C1"]
    C2in <- out_ests["C1"]
    init <- c(R1in, R2in, R3in, C1in, C2in)
    names(init) <- c("R1", "R2", "R3", "C1", "C2")
  }
  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
    in_sigma <- sigma(in_one_tm)
    out_sigma <- sigma(out_one_tm)
    sigma_ratio <- out_sigma / in_sigma
    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 <- "two_tm"
  return(fit)
}