R/optimal_parameter.R

In povmap: Extension to the 'emdi' Package

optimal_parameter <- function(generic_opt,
                              fixed,
                              smp_data,
                              smp_domains,
                              transformation,
                              interval,
                              framework) {

if (transformation != "no" && transformation != "log" &&
      transformation != "ordernorm" && transformation != "arcsin") {
    # no lambda -> no estimation -> no optmimization

    if (transformation == "box.cox" && any(interval == "default")) {
      interval <- c(-1, 2)
    } else if (transformation == "dual" && any(interval == "default")) {
      interval <- c(0, 2)
    } else if (transformation == "log.shift" && any(interval == "default")) {
      # interval = c(min(smp_data[paste(fixed[2])]),
      # max(smp_data[paste(fixed[2])]))
      span <- range(smp_data[paste(fixed[2])])
      if ((span[1] + 1) <= 1) {
        lower <- abs(span[1]) + 1
      } else {
        lower <- 0
      }

      upper <- diff(span) / 2

      interval <- c(lower, upper)
    }

    # Estimation of optimal lambda parameters
    optimal_parameter <- optimize(generic_opt,
      fixed          = fixed,
      smp_data       = smp_data,
      smp_domains    = smp_domains,
      transformation = transformation,
      interval       = interval,
      framework      = framework,
      maximum        = FALSE
    )$minimum
  } else {
    optimal_parameter <- NULL
  }

  return(optimal_parameter)
} # End optimal parameter


# Internal documentation -------------------------------------------------------

# Function generic_opt provides estimation method reml to specifiy
# the optimal parameter lambda. Here its important that lambda is the
# first argument because generic_opt is given to optimize. Otherwise,
# lambda is missing without default.

generic_opt <- function(lambda,
                        fixed,
                        smp_data,
                        smp_domains,
                        transformation,
                        framework) {

  # Definition of optimization function for finding the optimal lambda
  # Preperation to easily implement further methods here
  optimization <- if (TRUE) {
    reml(
      fixed = fixed,
      smp_data = smp_data,
      smp_domains = smp_domains,
      transformation = transformation,
      lambda = lambda,
      framework = framework
    )
  }
  return(optimization)
}



# REML method ------------------------------------------------------------------

reml <- function(fixed = fixed,
                 smp_data = smp_data,
                 smp_domains = smp_domains,
                 transformation = transformation,
                 lambda = lambda,
                 framework = framework) {

  sd_transformed_data <- std_data_transformation(
    fixed = fixed,
    smp_data = smp_data,
    transformation =
      transformation,
    lambda = lambda
  )

  model_REML <- NULL
  try(
    if(!is.null(framework$weights) && framework$weights_type == "nlme_lambda") {

      sd_transformed_data$weights_scaled <- smp_data[,framework$weights] /
        mean(smp_data[,framework$weights], na.rm = TRUE)

      model_REML <- lme(
        fixed = fixed,
        data = sd_transformed_data,
        random = as.formula(paste0("~ 1 | as.factor(", smp_domains, ")")),
        method = "REML",
        keep.data = FALSE,
        weights =
          varComb(varIdent(as.formula(
                      paste0("~ 1 | as.factor(", framework$smp_domains, ")"),)),
                  varFixed(as.formula(paste0("~1/", "weights_scaled")))
        )
      )
    } else {
      model_REML <- lme(
        fixed = fixed,
        data = sd_transformed_data,
        random = as.formula(paste0("~ 1 | as.factor(", smp_domains, ")")),
        method = "REML",
        keep.data = FALSE
      )
    }, silent = TRUE)

  if (is.null(model_REML)) {
    stop(strwrap(prefix = " ", initial = "",
                 "The likelihood does not converge. One reason could be that
                 the interval for the estimation of an optimal transformation
                 parameter is not appropriate. Try another interval. See also
                 help(ebp)."))
  } else {
    model_REML <- model_REML
  }


  log_likelihood <- -logLik(model_REML)

  return(log_likelihood)
}