R/optimal_parameter.R

Defines functions reml generic_opt optimal_parameter

optimal_parameter <- function(generic_opt,
                              fixed,
                              smp_data,
                              smp_domains,
                              transformation,
                              interval=c(-1,2)) {

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

    # Estimation of optimal lambda parameters
    optimal_parameter <- optimize(generic_opt,
                                  fixed          = fixed,
                                  smp_data       = smp_data,
                                  smp_domains    = smp_domains,
                                  transformation = transformation,
                                  interval       = interval,
                                  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
                        ) {


  #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
             )
        }
      return(optimization)
}



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

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

  sd_transformed_data <- std_data_transformation(fixed          = fixed,
                                                 smp_data       = smp_data,
                                                 transformation = transformation,
                                                 lambda         = lambda
                                                 )
  model_REML <- NULL
  try(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("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)
}
SoerenPannier/emdi documentation built on Oct. 11, 2018, 5:30 p.m.