Nothing
fitmv_body <- function(processed,
init=list(),
fixedS,
lower=list(),upper=list(),
# control.dist=list(), ## info in processed
control=list(), ## optimizer, <optimizer controls>, precision
nb_cores=NULL,
... ## cf dotnames processing below
) {
dotlist <- list(...) ## forces evaluations, which makes programming easier...
# if (is.list(processed)) { # actually FALSE in mv fits
# proc1 <- processed[[1L]]
# } else
proc1 <- processed
verbose <- proc1$verbose
HLnames <- (c(names(formals(HLCor)),names(formals(HLfit)),
names(formals(mat_sqrt)),names(formals(make_scaled_dist)))) ## cf parallel code in HLCor.obj
## fill HLCor.args
good_dotnames <- intersect(names(dotlist),HLnames) ## those specifically for the called fns as def'd by HLnames
if (length(good_dotnames)) {
HLCor.args <- dotlist[good_dotnames]
} else HLCor.args <- list()
## replace some HLCor.args members
if ( is.list(processed) ) {
pnames <- names(processed[[1]])
} else pnames <- names(processed)
for (st in pnames) HLCor.args[st] <- NULL
# 'processed' may be modified below, then will be copied in HLCor.args (and then removed from this envir for safety)
optim.scale <- control[["optim.scale"]]
if (is.null(optim.scale)) optim.scale="transformed" ## currently no public alternative
sparse_precision <- proc1$is_spprec
#
init$lambda <- .reformat_lambda(init$lambda, nrand=length(processed$ZAlist),
namesTerms=attr(processed$ZAlist,"namesTerms"), full_lambda=FALSE)
init$phi <- .reformat_phi(init$phi, n_models=length(processed$unmerged), full_phi=FALSE) # input list -> provide names to list -> output named list
# Automatically reformatting lower,upper, init NB_shape and COMP_nu would be less trivial as it should have the following format:
if (length(na.omit(init$NB_shape))!=length(names(init$NB_shape)) ||
length(na.omit(lower$NB_shape))!=length(names(lower$NB_shape)) ||
length(na.omit(upper$NB_shape))!=length(names(upper$NB_shape))
) stop("Any non-NULL [init|lower|upper]$NB_shape should be a vector of numeric values (no NA),\n named by index of submodels (e.g. c('2'=1)).")
if (length(na.omit(init$beta_prec))!=length(names(init$beta_prec)) ||
length(na.omit(lower$beta_prec))!=length(names(lower$beta_prec)) ||
length(na.omit(upper$beta_prec))!=length(names(upper$beta_prec))
) stop("Any non-NULL [init|lower|upper]$beta_prec should be a vector of numeric values (no NA),\n named by index of submodels (e.g. c('2'=1)).")
if (length(na.omit(init$COMP_nu))!=length(names(init$COMP_nu)) ||
length(na.omit(lower$COMP_nu))!=length(names(lower$COMP_nu)) ||
length(na.omit(upper$COMP_nu))!=length(names(upper$COMP_nu))
) stop("Any non-NULL [init|lower|upper]$NB_shape should be a vector of numeric values (no NA),\n named by index of submodels (e.g. c('2'=1)).")
#
user_init_optim <- init # more explicit name; will serve as template for canon.init which will serve as template for names of lewer, upper, initvec
optim_blob <- .calc_optim_args_mv(processed, map_rd_mv=attr(processed$ZAlist, "map_rd_mv"),
user_init_optim=user_init_optim, fixedS=fixedS, user.lower=lower, user.upper=upper, verbose=verbose, optim.scale=optim.scale)
# .calc_optim_args() for a familiy with outer disp param (negbin...) may tend to outer optimize lambda when it will tend to inner optimize it for other distribs
# so we have a bit incorent choices across submodels, that could be improved __FIXME__
# modify HLCor.args and <>bounds; ## distMatrix or uniqueGeo potentially added to HLCor.args:
# init <- optim_blob$inits$`init` ## list; keeps all init values, all in untransformed scale
init.optim <- optim_blob$inits$`init.optim` ## list; subset of all estimands, as name implies, and in transformed scale
init.HLfit <- optim_blob$inits$`init.HLfit` ## list; subset as name implies
fixed <- optim_blob$fixed
corr_types <- optim_blob$corr_types
LUarglist <- optim_blob$LUarglist
moreargs <- LUarglist$moreargs
LowUp <- optim_blob$LowUp
lower <- LowUp$lower ## list ! which elements may have length >1 !
upper <- LowUp$upper ## list !
#
if ( FALSE && # following is a block for fitme_body which is not suited for mv fits (with proc1$residProcesseds, not proc1$residProcessed !)
# This means that outer_optim_resid is far from ready in mv fits
! is.null(residproc1 <- proc1$residProcessed) &&
identical(spaMM.getOption("outer_optim_resid"),TRUE)) {
#
# spaMM.getOption("outer_optim_resid") being a private option, NULL By default, THIS IS NOT RUN
#
## Problem is that outer optim at the mean model level is useful if we can avoid computation of the leverages
## But here anyway we need the leverages of the 'mean' model to define the resid model response.
resid_optim_blob <- .calc_optim_args(proc_it=residproc1, processed=proc1,
user_init_optim=proc1$residModel[["init"]], fixed=proc1$residModel$fixed, ## all user input must be in proc1$residModel
lower=proc1$residModel$lower, upper=proc1$residModel$upper, ## all user input must be in proc1$residModel
verbose=c(SEM=FALSE), optim.scale=optim.scale, For="fitme")
resid_init.optim <- resid_optim_blob$inits$`init.optim` ## list; subset of all estimands, as name implies, and in transformed scale
proc1$residModel$`init.HLfit` <- resid_optim_blob$inits$`init.HLfit` ## list; subset as name implies
proc1$residModel$fixed <- resid_optim_blob$fixed
##### residproc1$corr_types <- resid_optim_blob$corr_types
# resid_user.lower <- resid_optim_blob$user.lower
# resid_user.upper <- resid_optim_blob$user.upper
init.optim <- c(init.optim,list(resid=resid_init.optim))
lower <- c(lower,list(resid=resid_optim_blob$LowUp$lower))
upper <- c(upper,list(resid=resid_optim_blob$LowUp$upper))
LowUp <- list(lower=lower,upper=upper)
moreargs <- c(moreargs,list(resid=resid_optim_blob$LUarglist$moreargs))
}
#
processedHL1 <- proc1$HL[1]
needHLCor_specific_args <- (length(unlist(lower$corrPars, use.names = FALSE)) ||
length(intersect(corr_types,c("Matern","Cauchy","adjacency","AR1","corrMatrix", "IMRF","corrFamily"))))
if (needHLCor_specific_args) {
HLcallfn.obj <- "HLCor.obj"
HLcallfn <- "HLCor"
control.dist <- vector("list",length(moreargs))
for (nam in names(moreargs)) control.dist[[nam]] <- moreargs[[nam]]$control.dist
HLCor.args[["control.dist"]] <- control.dist ## always reconstructed locally, not in the fitme_body call
} else {
HLcallfn.obj <- "HLfit.obj"
HLcallfn <- "HLfit"
}
HLCor.args$fixed <- fixed ## to be modified by objective function (at least in HLCor case)
HLCor.args$init.HLfit <- init.HLfit
HLCor.args$processed <- processed ## for the <...>.obj and <...>_body functions
##
anyHLCor_obj_args <- HLCor.args
## HLCor.obj uses a vector + skeleton
initvec <- unlist(init.optim)
if (needHLCor_specific_args) {
anyHLCor_obj_args$skeleton <- structure(init.optim, moreargs=moreargs, ## moreargs is a list over ranefs
type=relist(rep("fix",length(initvec)),init.optim) )
} else {
anyHLCor_obj_args$skeleton <- structure(init.optim,
type=relist(rep("fix",length(initvec)),init.optim))
}
.assignWrapper(anyHLCor_obj_args$processed,
paste0("return_only <- \"",proc1$objective,"APHLs\""))
if (length(initvec)) {
augZXy_phi_est <- NULL ## the value assumed by later code when augZXy was not used
if (.safe_true(verbose["getCall"][[1L]])) { ## to get an optim call with its initial value. Then HLcallfn is called and its call returned.
## confint -> get_HLCorcall needs an HLCor call with the following ranFix
ranPars_in_refit <- structure(.modify_list(fixed,init.optim),
# I label this "F I X" as a TAG for this modif type attribute:
type=.modify_list(.relist_rep("fix",fixed), #attr(fixed,"type"),
relist(rep("outer",length(initvec)),init.optim)) )
} else {
use_SEM <- (!is.null(processedHL1) && processedHL1=="SEM")
time2 <- Sys.time()
if (use_SEM) {
if (is.null(proc1$SEMargs$control_pmvnorm$maxpts)) {
.assignWrapper(processed,"SEMargs$control_pmvnorm$maxpts <- quote(250L*nobs)")
} ## else default visible in SEMbetalambda
## its names should match the colnames of the data in Krigobj = the parameters of the likelihood surface. Current code maybe not general.
loclist <- list(anyHLCor_obj_args=anyHLCor_obj_args, ## contains $processed
LowUp=LowUp,init.corrHLfit=init.optim, ## F I X M E usage of user_init_optim probably not definitive
control.corrHLfit=control,
verbose=verbose[["iterateSEM"]],
nb_cores=nb_cores)
optr <- .probitgemWrap("iterateSEMSmooth",arglist=loclist, pack="probitgem") # do.call("iterateSEMSmooth",loclist)
optPars <- relist(optr$par,init.optim)
if (!is.null(optPars)) attr(optPars,"method") <-"optimthroughSmooth"
refit_info <- control[["refit"]] ## may be a NULL/ NA/ boolean or a list of booleans
if ( is.null(refit_info) || is.na(refit_info)) refit_info <- FALSE ## alternatives are TRUE or an explicit list or NULL
} else {
optPars <- .new_locoptim(init.optim, ## try to use gradient? But neither minqa nor _LN_BOBYQA use gradients. optim() can
LowUp,
control, objfn_locoptim=.objfn_locoptim,
HLcallfn.obj=HLcallfn.obj, anyHLCor_obj_args=anyHLCor_obj_args,
user_init_optim=user_init_optim,
grad_locoptim=NULL, verbose=verbose[["TRACE"]])
refit_info <- attr(optPars,"refit_info") ## 'derives' from control[["refit"]] with modification ## may be NULL but not NA
}
optim_time <- .timerraw(time2)
augZXy_phi_est <- proc1$augZXy_env$phi_est ## may be NULL: if phi was not estimated by augZXy
refit_args <- .get_refit_args(fixed, optPars, processed, moreargs, proc1, refit_info, HLCor.args, augZXy_phi_est)
HLCor.args <- refit_args$HLCor.args
ranPars_in_refit <- refit_args$ranPars_in_refit
} ## end if ...getCall... else
#
# refit_info is list if so provided by user, else typically boolean. An input NA should have been converted to something else (not documented).
if (needHLCor_specific_args) attr(ranPars_in_refit,"moreargs") <- moreargs
HLCor.args$fixed <- ranPars_in_refit
# initiate *phi* and resid models for refit, up to random effect
residProcesseds <- processed$residProcesseds
for (mv_it in seq_along(residProcesseds)) {
if (length(residinfo_it <- residProcesseds[[mv_it]]$port_env$port_fit_values)) { # implies that there was a phi-resid.model
# next line typically creates a list, from possibly nonfirst element... untidy, but this seems OK.
HLCor.args$init.HLfit$phi[[mv_it]] <- residinfo_it$fv # The phi of the mean-response model is initiated by fitted values of the resid.model
processed$residModels[[mv_it]]$init.HLfit <- residinfo_it$init_HLfit # may have fixef, v_h, and more generally the elements of get_inits_from_fit(phifit,inner_lambdas = TRUE)$init.HLfit
}
}
# At this point there may be null elements in the init.HLfit$phi list, but .denullify() should handle this.
# 'cover_residM_reinit' fit in test-mv-extra covers this code.
} else if (len_ranPars <- length(unlist(HLCor.args$fixed, use.names = FALSE))){ ## Set attribute
HLCor.args$fixed <- structure(HLCor.args$fixed,
type = relist(rep("fix", len_ranPars), HLCor.args$fixed),
moreargs=moreargs) ## moreargs needed if user handles fixed(<transformed params>) ('hyper' tests)
}
#
# not local to anyHLCor_obj_args$processed: change processed globally
.assignWrapper(HLCor.args$processed,"return_only <- NULL")
.assignWrapper(HLCor.args$processed,"verbose['warn'] <- TRUE") ## important!
hlcor <- do.call(HLcallfn,HLCor.args) ## recomputation post optimization, or only computation if length(initvec)=0, or the HLCorcall
if (is.call(hlcor)) {
if (length(initvec)) {
attr(hlcor,"optimInfo") <- list(LUarglist=LUarglist, init.optim=init.optim,
objective=proc1$objective,
augZXy_phi_est=augZXy_phi_est,
rC_transf=.spaMM.data$options$rC_transf)
}
return(hlcor) ## HLCorcall
} else {
if ( processed$fitenv$prevmsglength) { # there was output for a phi-resid.model. The fit object may then be printed...
cat("\n")
processed$fitenv$prevmsglength <- 0L
}
}
# hlcor<- .update_ranef_info(hlcor, moreargs=moreargs)
if (length(initvec)) {
attr(hlcor,"optimInfo") <- list(LUarglist=optim_blob$LUarglist, optim.pars=optPars,
objective=proc1$objective,
augZXy_phi_est=augZXy_phi_est, ## gives info to interpret optim.pars in confint.HLfit()
optim_time=optim_time,
rC_transf=.spaMM.data$options$rC_transf) ## processed was erased for safety
locoptr <- attr(optPars,"optr")
if (attr(optPars,"method")=="nloptr") {
if (locoptr$status<0L) hlcor$warnings$optimMessage <- paste0("nloptr() message: ",
locoptr$message," (status=",locoptr$status,")")
} else if ( attr(optPars,"method")=="optim" ) {
if (locoptr$convergence) hlcor$warnings$optimMessage <- paste0("optim() message: ",locoptr$message,
" (convergence=",locoptr$convergence,")")
} else if ( attr(optPars,"method")== "optimthroughSmooth") {
# provide logL estimate from the smoothing, to be used rather than the hlcor logL :
logLapp <- optr$value
attr(logLapp,"method") <- " logL (smoothed)"
hlcor$APHLs$logLapp <- logLapp
}
hlcor$warnings$suspectRho <- .check_suspect_rho(corr_types, ranPars_in_refit, LowUp)
if ( ! is.null(PQLdivinfo <- processed$envir$PQLdivinfo)) {
hlcor$divinfo <- PQLdivinfo
hlcor$warnings$divinfo <- "Numerical issue detected; see div_info(<fit object>) for more information."
warning(hlcor$warnings$divinfo)
}
}
## substantial effect on object size! :
lsv <- c("lsv",ls())
rm(list=setdiff(lsv,"hlcor"))
##
return(hlcor)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.