Nothing
.composite_pred_warn <- local({
composite_pred_warned <- FALSE
function() {
if ( ! environment(.composite_pred_warn)$composite_pred_warned) {
warning(paste("predictions combining 'newdata' and composite random effects\n",
"may fail for most corrFamily and perhaps other terms."),
call.=FALSE)
composite_pred_warned <<- TRUE
}
}
})
.calc_new_X_ZAC_mv <- function(object, newdata=NULL, re.form = NULL,
variances=list(residVar=FALSE, cov=FALSE),invCov_oldLv_oldLv_list,
control=list(), na.action=na.omit) {
locformS <- formula.HLfit(object, which="")
if (inherits(re.form, "formula")) {
re.formS <- vector("list", length(locformS))
for (mv_it in seq_along(locformS)) {
# As the doc says... simulate.HLfit(., type="marginal") -> default re.form = NA -> does not mean that the ranef is absent,
# but that simulation is unconditional on it. re.forms[[mv_it]] = NA should have the same meaning.
re.formS[[mv_it]] <- .update_formula_shared_ranefs(locform=re.form, re.form=locformS[[mv_it]], rm_LHS=TRUE)
}
} else if (length(re.form)==1L && is.na(re.form)) {
re.formS <- rep(NA, length(locformS))
} else re.formS <- re.form # a list or NULL.
# checking variables in the data BEFORE removing "marginalized upon" ranefs.
# These variables are needed in newdata_it <- locdataS[[mv_it]] in simulate.HLfit() -> .calc_ZAlist_newdata_mv()
allvarsS <- vector("list", length(locformS))
for (mv_it in seq_along(locformS)) {
## [-2] important to ignore response variables
allvars_it <- all.vars(.strip_cF_args(locformS[[mv_it]][-2])) ## strip to avoid e.g. 'stuff' being retained as a var from IMRF(..., model=stuff)
if (variances$residVar || control$simulate) {
allvars_it <- unique(c(allvars_it, all.vars(.strip_cF_args(.get_phiform(object, mv_it)))))
allvars_it <- unique(c(allvars_it, all.vars(object$families[[mv_it]]$resid.model$resid.formula)))
}
allvarsS[[mv_it]] <- allvars_it
}
#
## possible change of random effect terms (removing "marginalized upon" ranefs)
for (mv_it in seq_along(locformS)) locformS[[mv_it]] <- .update_formula_shared_ranefs(locformS[[mv_it]], re.formS[[mv_it]], rm_LHS=FALSE)
## matching ranef terms of re.form
if (length(object$ZAlist)) {
if (identical(control$keep_ranef_covs_for_simulate, TRUE) || # : condition for the case
# where only eta_fixed is predicted for marginal simulation, hence re.form is NA ("no prediction for ranef") BUT
# we will need the locdataS with the variables for ranefs, to simulate these ranefs.
# We will need ALSO marginal covariance matrices for the ranefs !! The ZAL in simulate.HLfit() has been correct
# before and after changes in the simulate.HLfit() code 2023/07/23
( re.form_allows_ranefs <- ( is.null(re.form) || # : this condition means that re.form_allows_ranefs may be TRUE despite no ranef in model formula
any( ! sapply(re.formS, .noRanef))) )
) {
map_rd_mv <- attr(object$ZAlist, "map_rd_mv")
ori_exp_ranef_strings <- attr(object$ZAlist,"exp_ranef_strings")
#
newinoldS <- vector("list", length(locformS))
# As the doc says... simulate.HLfit(., type="marginal") -> default re.form = NA -> does not mean that the ranef is absent,
# but that simulation is unconditional on it. re.forms[[mv_it]] = NA should have the same meaning.
# In the present step we retain only the ranefs conditioned upon. The Zlist computed here will be used to put then in
# 'eta_fixed_cond' (in simulate.HLfit).
# The 'marginalized upon' ones are handled by simulate.HLfit -> (get_ZALMatrix or .calc_ZAlist_newdata) + .simulate_ranefs
for (mv_it in seq_along(locformS)) newinoldS[[mv_it]] <- .get_newinold(re.formS[[mv_it]], locformS[[mv_it]],
ori_exp_ranef_strings, rd_in_mv=map_rd_mv[[mv_it]])
newinold <- unique(.unlist(newinoldS))
} else newinold <- NULL
} else newinold <- NULL
#
if (length(newinold)) {
new_exp_ranef_strings <- ori_exp_ranef_strings[newinold]
ranef_form <- as.formula(paste("~",(paste(new_exp_ranef_strings,collapse="+")))) ## effective '.noFixef'
ranefvars <- all.vars(.strip_cF_args(ranef_form))
} else ranefvars <- c()
need_new_design <- ( ( ! is.null(newdata) ) || ! is.null(re.form)) ## newdata or new model
locdata <- .get_locdata(newdata=newdata, locvars=unique(c(.unlist(allvarsS),ranefvars)),
object=object, variances=variances,
na.action= if (need_new_design) {na.pass} else {na.action}) # see comment on evaluation of RESU$newuniqueGeo
#
locdataS <- vector("list", length(allvarsS))
if (need_new_design) {
newX.pv <- eta_fix <- NULL
for (mv_it in seq_along(locformS)) {
allvars_it <- unique(c(allvarsS[[mv_it]], ranefvars)) # all ranefvars + mv_it-specific fixefvars: note difference if ! need_new_design
# presence of ranefvars for all ranefs in each allvars_it allows do.call(rbind, locdataS) later in this fn after selecting them
locdataS[[mv_it]] <- ..get_locdata(locdata, allvars_it, na.action=na.action, mv_it=mv_it)
newX_info <- .get_newX_info(locformS[[mv_it]], locdataS[[mv_it]], object, mv_it=mv_it)
newX.pv <- .merge_Xs(newX.pv, newX_info$newX.pv, mv_it)
eta_fix <- c(eta_fix, newX_info$eta_fix)
}
if ( ! is.null(X2X <- object$X2X)) newX.pv <- newX.pv %*% object$X2X
attr(locdataS,"allvarsS") <- allvarsS # allVarsS (used by map_ranef()) tells which variables were needed for which submodel predictions.
RESU <- list(locdata=locdataS, # locdata in RESU for cbind()ing the predictions in .predict_body();
cum_nobs= cumsum(c(0L,lapply(locdataS, nrow))), # more widely used by .fv_linkinv()
newX.pv=newX.pv, eta_fix=eta_fix)
} else {
for (mv_it in seq_along(allvarsS)) locdataS[[mv_it]] <- ..get_locdata(locdata, allvarsS[[mv_it]], na.action=na.action)
RESU <- list(locdata=locdataS, # locdata in RESU allowing (potential) cbind() with predictions in .predict_body(). (maybe not implemented for mv)
cum_nobs= cumsum(c(0L,lapply(locdataS, nrow))), # more widely used by .fv_linkinv()
newX.pv=model.matrix(object))
}
## so we save 'locdata=locdataS' in RESU, but we will locally modify 'locdataS' by selecting columns in locdataS[[mv_it]]
## and will create a local 'ranefdata' from this locally modified 'locdataS'.
#
## newZAlist and subZAlist appear to have distinct usages since they are created under different conditions.
## subZAlist is a subset of the old ZA, newZAlist contains new ZA
#
if (nrand <- length(newinold)) {
if (length( info_olduniqueGeo <- .get_old_info_uniqueGeo(object) )) {
RESU$newuniqueGeo <- .update_newuniqueGeo(info_olduniqueGeo, newinold, need_new_design, locdata)
# Despite the $newuniqueGeo name, it may be necess without newdata;
# cf preprocess_fix_corr() with NULL 'fixdata' providing info_olduniqueGeo <- fix_info$newuniqueGeo (univariate case).
# There may be a slight suboptimality as it uses 'locdata' produced with na.rm=FALSE.
# But na.rm=TRUE would remove rows used in a submodel, if they have NA's in variables not used in that submodel.
}
RESU$subZAlist <- object$ZAlist[newinold] ## and reordered # a priori for .wrap_calcPredVar
RESU$newinold <- newinold
ori_exp_ranef_types <- attr(object$ZAlist,"exp_ranef_types")
RESU$spatial_old_rd <- which(ori_exp_ranef_types != "(.|.)")
strucList <- object$strucList
if (need_new_design) {
## with newdata we need Evar and then we need nn... if newdata=ori data the Evar (computed with the proper nn) should be 0
# barlist <- .process_bars_mv(predictors=formula.HLfit(object,which = ""),
# map_rd_mv=map_rd_mv,
# as_character=FALSE) ## but default expand =TRUE
# barlist <- structure(barlist[newinold], type=attr(barlist,"type")[newinold])
for (mv_it in seq_along(locdataS)) locdataS[[mv_it]] <- locdataS[[mv_it]][,ranefvars, drop=FALSE]
ranefdata <- do.call(rbind, locdataS)
ori_exp_ranef_terms <- attr(object$ZAlist,"exp_ranef_terms")
# new_exp_ranef_terms <- structure(ori_exp_ranef_terms[newinold], type=attr(ori_exp_ranef_terms,"type")[newinold])
#
# This will first construct a list possibly with excess nonzero elements in the Z's
# This will be corrected in the next loop, cf .Dvec_times_Matrix(newrd_in_obs, newZlist[[new_rd]])
#
raneftypes <- attr(ori_exp_ranef_terms,"type")
if (any(raneftypes[newinold] %in% c("MaternIMRFa","corrFamily") & # tentative
object$ranef_info$is_composite[newinold]) # does not distinguish (0+(mv())) from other composite
) .composite_pred_warn()
newZlist <- .calc_Zlist(exp_ranef_terms=ori_exp_ranef_terms, # subsetting -> new_exp_ranef_terms will be made internally using rd_in_mv arg
data=ranefdata,
rd_in_mv=newinold, # here, must be "conditioned upon" ranefs
rmInt=0L, drop=TRUE,sparse_precision=FALSE,
corr_info=.get_from_ranef_info(object),
#
# ! use levels_type default as is required for simulation !
# In (at least, marginal sim for) univariate case, .calc_new_X_ZAC is not where new Zlist is computed, levels_type= "seq_len" can be used;
# In simulation for mv case, .calc_new_X_ZAC_mv() provides new Zlist, the explicit levels_type arg should not be used.
#
## Old comment:
# levels_type= "seq_len", ## superseded in specific cases: notably,
# ## the same type has to be used by .calc_AMatrix_IMRF() -> .as_factor()
# ## as by .calc_Zmatrix() -> .as_factor() for IMRFs.
# ## This is controlled by package option 'uGeo_levels_type' (default = "data_order" as the most explicit).
# ## The sames functions are called with the same arguments for predict with newdata.
# ##
# ## This means that if there are repeated geo positions in the newdata
# ## we save the time of trying to find them (which perhaps is less justifiable in mv case? __FIXME__)
sub_oldZAlist=object$ZAlist, # subsetting will be made internally
lcrandfamfam=attr(object$rand.families,"lcrandfamfam"))
# As explained above,
# each Z in newZlist then has non-zero rows even for response levels that are not affected by the ranef
# bc it's built from 'ranefdata' built as if all submodels were affected by each (new) ranef
# => need to correct this
newrd_in_obsS <- vector("list", length(newinoldS))
loc_cum_nobs <- RESU$cum_nobs
for(new_rd in seq_along(newinold)) {
for (mv_it in seq_along(newinoldS)) {
nobs_it <- loc_cum_nobs[[mv_it+1L]] - loc_cum_nobs[[mv_it]]
newrd_in_obsS[[mv_it]] <- rep((newinold[new_rd] %in% newinoldS[[mv_it]]), nobs_it)
}
newrd_in_obs <- 1L * .unlist(newrd_in_obsS)
newZlist[[new_rd]] <- structure(.Dvec_times_Matrix(newrd_in_obs, newZlist[[new_rd]]),
is_incid=attr(newZlist[[new_rd]],"is_incid"))
}
amatrices <- .get_new_AMatrices(object, newdata=ranefdata) # .calc_newFrames_ranef(formula=ranef_form,data=ranefdata,fitobject=object)$mf)
## ! complications:
## even if we used perm_Q for *fitting* Matern , the permutation A matrix should not be necessary in building the new correlation matrix, bc ./.
## explicit colnames should handle both cases, so that
## newZAlist <- .calc_normalized_ZAlist( ignoring those A matrices)
## and
## newZAlist <- object$ZAlist
## should be OK.
## But the other Amatrices should be processed before newZACpplist <- .compute_ZAXlist(.) is called
requires_ZCpL <- (attr(newZlist,"exp_ranef_types") %in% c("Matern","Cauchy"))
newZAlist <- .calc_normalized_ZAlist(Zlist=newZlist,
# newZlist has names not necessarily starting at "1"
AMatrices=amatrices[names(newZlist)[ ! requires_ZCpL]],
vec_normIMRF=object$ranef_info$vec_normIMRF,
strucList=strucList)
## must be ordered as parseBars result for the next line to be correct.
attr(newZAlist,"exp_ranef_strings") <- new_exp_ranef_strings ## required pour .compute_ZAXlist to match the ranefs of LMatrix
newZAlist <- .correct_newZAlist_mv_ranCoefs(ZAlist=newZAlist,
cum_nobs=cumsum(c(0L,sapply(locdataS,nrow))))
# : distinct .correct_...() function needed here bc [
# we cannot run an outer loop on mv_it, calling.merge_ZAlists(., ZAlist2 = newZAlist_i...) bc
# expected format of final newZAlist differs from that for the fit: for some ranef types
# we do not try to match the columns of ZA_i's different submodels, instead building larger but diagonal ZA's
# ]
} else {
newZAlist <- object$ZAlist
ranefdata <- object$data
}
RESU$newZAlist <- newZAlist
# We determine which matrices we need for computation of Evar:
need_Cnn <- .calc_need_Cnn(object, newinold, ori_exp_ranef_types, variances, newZAlist)
which_mats <- list(no= need_new_design,
## cov_newLv_newLv_list used in .calc_Evar() whenever newdata, but elements may remain NULL if $cov not requested
## However, for ranCoefs, we need Xi_cols rows for each response's predVar. (FIXME) we store the full matrix.
nn= need_Cnn )
#nrand <- length(newinold)
#
## AT this point both newZAlist and subZAlist may have been reduced to 'newnrand' elements relative to ori object$ZAlist.
## .match_old_new_levels will use it running over newnrand values
## The cov_ lists are reduced too. newinold should be used to construct them
## newZACpplist is reduced.
if (any(unlist(which_mats))
|| any(unlist(variances)) # cov_newLv_oldv_list is always needed for cbind(X.pv,newZAC [which may be ori ZAC]); should ~corr_list when newdata=ori data
) {
blob <- .make_new_corr_lists(object=object,locdata=ranefdata, which_mats=which_mats,
newZAlist=newZAlist, newinold=newinold,
invCov_oldLv_oldLv_list=invCov_oldLv_oldLv_list,
for_mv=TRUE)
RESU <- .update_cov_no_nn(RESU, blob, which_mats, newZAlist)
RESU$newZACpplist <- .compute_ZAXlist(ZAlist=newZAlist, XMatrix=RESU$cov_newLv_oldv_list) ## build from reduced list, returns a reduced list
## This $newZACpplist serves to compute new _point predictions_.
# # this comment may be obsolete : .compute_ZAXlist affects elements of ZAlist that have a ranefs attribute.
# It builds a design matrix to all oldv levels. It does not try to reduce levels.
# But it use newZlist <- .calc_Zlist(...) which may consider a reduced number of levels.
# The function .match_old_new_levels() will perform the match with 'w_h_coeffs' for point prediction.
# it assign values psi_M to new levels of ranefs.
## _Alternatively_ one may construct a newZACvar for _predVar_
# Here we have a slice mechanism (contrary for point pred) hence new ZA with different rows, and the columns of the
# newZACvar constructed there must match those of beta_w_cov for ZWZt_mat_or_diag( <cbind(newX,newZAC)> ,beta_w_cov)
# cov_newLv_oldv_list() provides the info for the expansion from the newZA cols to the oldZA cols.
# In that case one does not need to match levels. .calc_newZACvar() performs a simpler operation than .compute_ZAXlist.
if (need_new_design) RESU$newZACpplist <- .calc_ZAlist(Zlist=RESU$newZACpplist, AMatrices=amatrices[names(newZAlist)[requires_ZCpL]])
}
#
}
return(RESU)
}
.calc_logphiInfo <- function(invV_factors, which, cum_nobs, phi_est, w.resid, spprec) { ## called by .calc_logdisp_cov_mv(), using the result of .calc_invV.dV_info()
n_phi <- length(which)
logphiInfo <- matrix(ncol=n_phi,nrow=n_phi)
lhs_invV.dVdphi_list <- wresid_list <- list()
zerovec <- 0*w.resid
for (colit in seq_along(which)) {
mv_it <- which[colit]
iresp_range <- .subrange(cum_nobs, mv_it)
indic_it <- zerovec
indic_it[iresp_range] <- 1
i_lhs_invV.dVdphi <- .Dvec_times_Matrix(indic_it, invV_factors$n_x_r) # .fill_rhs_invV.dVdlam(template=zerotemplate, urange=uirange[,ilam], invV.dV_info)
lhs_invV.dVdphi_list[[as.character(mv_it)]] <- i_lhs_invV.dVdphi
wresid_list[[as.character(mv_it)]] <- wresid_it <- w.resid*indic_it
if (spprec) {
A <- invV_factors$r_x_n %*% invV_factors$n_x_r
trAB <- sum(A^2)
} else { # .traceAB() differs by checking which is the largest dimension:
trAB <- .traceAB(lA=i_lhs_invV.dVdphi, rA=invV_factors$r_x_n,
B_is_tA=TRUE # (lA is n x r, rA is r x n, so A is square as this cases requires)
)
}
# using the pattern (D-nXr.rXn)^2 = D^2 - 2 D nXr.rXn + (nXr.rXn)^2
logphiInfo[colit,colit] <- sum(wresid_it^2) -2 * .traceDB(wresid_it,i_lhs_invV.dVdphi, invV_factors$r_x_n) + trAB
for (coljt in seq_len(colit-1L)) {
mv_jt <- which[coljt]
wresid_jt <- wresid_list[[as.character(mv_jt)]]
j_lhs_invV.dVdphi <- lhs_invV.dVdphi_list[[as.character(mv_jt)]]
# trace[(D1-l*r) (D2-l*r)], with D1*D2=0:
logphiInfo[coljt,colit] <- - .traceDB(wresid_it,i_lhs_invV.dVdphi, invV_factors$r_x_n) -
.traceDB(wresid_jt,j_lhs_invV.dVdphi, invV_factors$r_x_n) +
.traceAB(i_lhs_invV.dVdphi, invV_factors$r_x_n, t(invV_factors$r_x_n), t(j_lhs_invV.dVdphi))
logphiInfo[colit,coljt] <- logphiInfo[coljt,colit]
}
}
sub_phi_vec <- .unlist(phi_est[which])
logphiInfo <- logphiInfo * (sub_phi_vec %*% t(sub_phi_vec)) # inspired from code for lambda
return(list(logphiInfo=logphiInfo,lhs_invV.dVdphi_list=lhs_invV.dVdphi_list, wresid_list=wresid_list))
}
.calc_logdisp_cov_mv <- function(object, dvdloglamMat=NULL, dvdlogphiMat=NULL, invV_factors=NULL, force_fixed=FALSE) {
lambda.object <- object$lambda.object
strucList <- object$strucList
dispcolinfo <- list()
problems <- list() ## Its elements' names are tested in calcPredVar, and the strings are 'development info'
nrand <- length(strucList)
col_info <- list(nrand=nrand, phi_cols=NULL)
Xi_cols <- attr(object$ZAlist, "Xi_cols")
dwdloglam <- matrix(0,ncol=sum(Xi_cols),nrow=length(object$v_h)) # cols will remain 0 for fixed lambda params
if (force_fixed) {
checklambda <- rep(TRUE, length(lambda.object$type))
} else checklambda <- ( ! (lambda.object$type %in% c("fixed","fix_ranCoefs","fix_hyper")))
if (any(checklambda)) {
exp_ranef_types <- attr(object$ZAlist,"exp_ranef_types")
checkadj <- (exp_ranef_types=="adjacency")
if(any(checkadj)) {
## several blocks of code are "maintained" below for a future dispVar computation for rho
# il me manque dwdrho (et meme dwdloglam pour ce modele ?) donc on inactive les lignes suivantes:
# if (is.null(lambda.object$lambda.fix)) dispnames <- c(dispnames,"loglambda")
# corrFixNames <- names(unlist(object$corrPars[which(attr(corrPars,"type")=="fix")]))
# if (! ("rho" %in% corrFixNames) ) dispnames <- c(dispnames,"rho")
}
if (is.null(dvdloglamMat)) {
## note that .get_logdispObject is computed on request by .get_logdispObject()
problems$stopmiss <- warning("is.null(dvdloglamMat) in a case where it should be available.")
}
dispcolinfo$loglambda <- "loglambda"
#dvdloglam <- matrix(0,nrow=NROW(dvdloglamMat), ncol=sum(Xi_cols))
strucList <- object$strucList
cum_n_u_h <- attr(lambda.object$lambda_list,"cum_n_u_h")
n_u_h <- diff(cum_n_u_h)
cum_Xi_cols <- cumsum(c(0,Xi_cols))
## dwdloglam will include cols of zeros for fixed lambda; matching with reduced logdisp_cov is performed at the end of the function.
for (randit in seq_len(nrand)) { ## ALL ranefs!
range_in_dw <- (cum_n_u_h[randit]+1L):(cum_n_u_h[randit+1L])
if (OLD <- TRUE) { # what about inverting the two operations ? ___F I X M E____
# solve(t(lmatrix) bc to be used as factor to ZAC, not ZA
if ( inherits(strucList[[randit]],"dCHMsimpl")) {
# previous comment: (not expected in default use in mv, for reasons explained in mv, sinc AUG_ZXy presumably FALSE )
# BUT: occurs in univariate case with get_predVar(adjfitsp), where augZXy algo does not appear to be used
for_dw_i <- as(strucList[[randit]], "CsparseMatrix") %*% dvdloglamMat[range_in_dw, ] # i.e L_Q %*% lignes de (t(L_Q) %*% invG %*% L_Q %*% some rhs)
} else if ( ! is.null(lmatrix <- strucList[[randit]])) {
for_dw_i <- solve(t(lmatrix),dvdloglamMat[range_in_dw,]) ## f i x m e for efficiency ? store info about solve(t(lmatrix)) in object ?
} else { ## implicit identity lmatrix
for_dw_i <- dvdloglamMat[range_in_dw,] ## assuming each lambda_i = lambda in each block
}
nblocks_randit <- Xi_cols[randit]
rowranges_in_dw_i <- matrix(seq(n_u_h[randit]),ncol=nblocks_randit) ## this _splits_ seq(n_u_h[randit]) over two columns for a random-slope model
for (row_block in seq_len(nblocks_randit)) { ## half-ranges for random-slope model
rowrange_in_dw_i <- rowranges_in_dw_i[,row_block]
cum_rowrange_in_dw <- rowrange_in_dw_i + cum_n_u_h[randit]
for (randjt in which(checklambda)) { ## NOT all ranefs!
nblocks_randjt <- Xi_cols[randjt]
cum_colrange_in_dw_i <- (cum_n_u_h[randjt]+1L):(cum_n_u_h[randjt+1L])
cum_colranges_in_dw_i <- matrix(cum_colrange_in_dw_i,ncol=nblocks_randjt) ## this _splits_ seq(n_u_h[randit]) over two columns for a random-slope model
for (col_in_colranges_dw_i in nblocks_randjt) { ## half-ranges for random-slope model
cum_col_in_dw <- cum_Xi_cols[randjt]+col_in_colranges_dw_i
cum_cols_in_dw_i <- cum_colranges_in_dw_i[,col_in_colranges_dw_i]
dwdloglam[cum_rowrange_in_dw, cum_col_in_dw] <- rowSums(for_dw_i[rowrange_in_dw_i, cum_cols_in_dw_i,drop=FALSE])
}
}
}
} else {
for_dw_i <- dvdloglamMat[range_in_dw,]
nblocks_randit <- Xi_cols[randit]
rowranges_in_dw_i <- matrix(seq(n_u_h[randit]),ncol=nblocks_randit) ## this _splits_ seq(n_u_h[randit]) over two columns for a random-slope model
for (row_block in seq_len(nblocks_randit)) { ## half-ranges for random-slope model
rowrange_in_dw_i <- rowranges_in_dw_i[,row_block]
cum_rowrange_in_dw <- rowrange_in_dw_i + cum_n_u_h[randit]
for (randjt in which(checklambda)) { ## NOT all ranefs!
nblocks_randjt <- Xi_cols[randjt]
cum_colrange_in_dw_i <- (cum_n_u_h[randjt]+1L):(cum_n_u_h[randjt+1L])
cum_colranges_in_dw_i <- matrix(cum_colrange_in_dw_i,ncol=nblocks_randjt) ## this _splits_ seq(n_u_h[randit]) over two columns for a random-slope model
for (col_in_colranges_dw_i in nblocks_randjt) { ## half-ranges for random-slope model
cum_col_in_dw <- cum_Xi_cols[randjt]+col_in_colranges_dw_i
cum_cols_in_dw_i <- cum_colranges_in_dw_i[,col_in_colranges_dw_i]
dwdloglam[cum_rowrange_in_dw, cum_col_in_dw] <- rowSums(for_dw_i[rowrange_in_dw_i, cum_cols_in_dw_i,drop=FALSE])
}
}
}
if ( inherits(strucList[[randit]],"dCHMsimpl")) {
# previous comment: (not expected in default use in mv, for reasons explained in mv, sinc AUG_ZXy presumably FALSE )
# BUT: occurs in univariate case with get_predVar(adjfitsp), where augZXy algo does not appear to be used
dwdloglam[range_in_dw, ] <- as.matrix(as(strucList[[randit]], "CsparseMatrix") %*%
dwdloglam[range_in_dw, ]) # i.e L_Q %*% lignes de (t(L_Q) %*% invG %*% L_Q %*% some rhs)
} else if ( ! is.null(lmatrix <- strucList[[randit]])) {
dwdloglam[range_in_dw, ] <- as.matrix(solve(t(lmatrix),dwdloglam[range_in_dw,])) ## f i x m e for efficiency ? store info about solve(t(lmatrix)) in object ?
}
}
}
## dwdloglam includes cols of zeros for fixed lambda; matching with reduced logdisp_cov is performed at the end of the function.
ranef_ids <- rep(seq_len(nrand),Xi_cols) ## (repeated for ranCoefs) indices of ranefs, not cols of ranefs
} else ranef_ids <- NULL
###
phimodel <- object$models[["phi"]]
is_phiScalS <- (phimodel=="phiScal" | force_fixed) # vector !
if (any(is_phiScalS) ||
identical(object$envir$forcePhiComponent,TRUE) ## hack for code testing: force dispVar computation as if phi was not fixed.
) { ## semble impliquer pas outer phi.Fix... => no need to test object$phi.object$phi_outer,"type")
phi_est <- object$phi ## no need to get object$phi.object$phi_outer
n_phiScal <- length(which(is_phiScalS))
for (mv_it in which(is_phiScalS)) if (length(phi_est[[mv_it]])!=1L) problems$stopphi <- warning("phimodel=\"phiScal\" but length(phi_est)!=1L.")
if ( ! is.null(dvdlogphiMat)) {
cum_nobs <- attr(object$vec_nobs, "cum_nobs")
dvdlogphiS <- vector("list", length(is_phiScalS))
for (mv_it in seq_along(dvdlogphiS)) { ## there are cols even for fixed phiS
if (is_phiScalS[mv_it]) {
resp_range <- .subrange(cum_nobs, mv_it)
dvdlogphiS[[mv_it]] <- rowSums(dvdlogphiMat[,resp_range, drop=FALSE]) ## using each phi_i = phi # always a vector, even from 0-col matrix
} else dvdlogphiS[[mv_it]] <- NULL # rep(0, nrow(dvdlogphiMat)) ## using each phi_i = phi # always a vector, even from 0-col matrix
}
dvdlogphi <- do.call(cbind, dvdlogphiS)
dwdlogphi <- .calc_invL_coeffs(object,dvdlogphi) # input matrix -> output matrix
col_info$phi_cols <- length(ranef_ids)+seq_len(n_phiScal) ## cols indices for phi
dispcolinfo$logphi <- rep("logphi", n_phiScal)
} else if (object$models[["eta"]]=="etaHGLM") stop("phimodel=='phiScal' but is.null(dvdlogphiMat)")
} else { ## else phimodel="", e.g. binomial
# if binomial or poisson, phimodel=""; warning for other phimodels
if (any(phimodel[!is_phiScalS] !="")) {
problems$structphi <- "phi dispVar component not yet available for phi model != ~1."
if ( ! identical(spaMM.getOption("phi_dispVar_comp_warned"),TRUE)) {
warning(problems$structphi)
.spaMM.data$options$phi_dispVar_comp_warned <- TRUE
}
}
dwdlogphi <- NULL
}
## compute info matrix:
if ((length(dispcolinfo))==0L) {
return(list(problems=problems))
} else {
dwdlogdisp <- cbind(dwdloglam,dwdlogphi) ## typically nobs * 2
attr(dwdlogdisp,"col_info") <- col_info
# cf my documentation, based on McCullochSN08 6.62 and 6.74
# lambda and phi factors enter in dV/dlog(.), computed instead of dV/d(.) to match dwdlog(.) vectors.
#
# use repres of two matrices large A and B, each as (thin) lhs %*% (flat) rhs
ZAL <- get_ZALMatrix(object, force_bind = ! (.is_spprec_fit(object)) )
if ("loglambda" %in% names(dispcolinfo) || "rho" %in% names(dispcolinfo)) {
invV.dV_info <- .calc_invV.dV_info(object, checklambda, invV_factors=invV_factors, ZAL=ZAL) ## $lhs= invV %*% ZALd and $lhs= t(ZALd)
sublambda <- .unlist(invV.dV_info$lambda_list[checklambda])
dispcolinfo$loglambda <- rep("loglambda",length(sublambda))
}
cum_n_disp_pars <- cumsum(c(0,lapply(dispcolinfo,length))) # #ncols for phi, lambda[checklambda], etc.
dispcols <- lapply(seq_along(dispcolinfo), function(varit) {
cum_n_disp_pars[varit]+ seq_along(dispcolinfo[[varit]])
}) ## col ranges for lambda[checklambda], phi, etc (with wrong names)
names(dispcols) <- dispnames <- names(dispcolinfo) ## list names
nrc <- cum_n_disp_pars[length(cum_n_disp_pars)]
#
logdispInfo <- matrix(NA,nrow=nrc,ncol=nrc)
colnames(logdispInfo) <- rownames(logdispInfo) <- .unlist(dispcolinfo)
if ("loglambda" %in% dispnames) {
loglamInfo_blob <- .calc_loglamInfo(invV.dV_info,which=which(checklambda))
logdispInfo[dispcols$loglambda,dispcols$loglambda] <- loglamInfo_blob$loglamInfo
}
if ("rho" %in% dispnames) { ## will occur only when if (any(checkadj)) {...} block above is fixed and active.
# no use of sqrt because adjd can be negative
#invV.dVdrho <- (invV %id*id% ZAL) %*% ( Diagonal(x=lambda*adjd/(denom^2)) %id*id% t(ZAL))
lhs_invV.dVdrho <- .calc_lhs_invV.dVdlam(object, ZAL, invV_factors) # sweep( ZAL,1L,object$w.resid,`*`) - lhs_invV.dVdrho
lambda_adjd <- invV.dV_info$lambda_list[[which(checkadj)]] ## asumes single adjd
rhs_invV.dVdrho <- ( Diagonal(x=lambda_adjd*invV.dV_info$adjd_denom2) %id*id% t(ZAL)) ## FIXME curently meaningful for only one lambda element
#logdispInfo["rho","rho"] <- sum(invV.dVdrho*t(invV.dVdrho))
logdispInfo[dispcols$rho,dispcols$rho] <- .traceAB(lhs_invV.dVdrho,rhs_invV.dVdrho,t(rhs_invV.dVdrho),t(lhs_invV.dVdrho))
if ("loglambda" %in% dispnames) {
sublambda <- .unlist(invV.dV_info$lambda_list)
logdispInfoBlock <- numeric(nrand)
cum_n_u_h <- invV.dV_info$cum_n_u_h
zerotemplate <- rep(0,cum_n_u_h[nrand+1L])
for (randit in which(checklambda)) {
u.range <- (cum_n_u_h[randit]+1L):(cum_n_u_h[randit+1L])
Xi_ncol <- Xi_cols[randit] # say '1 2' for ranCoefs
uirange <- matrix(u.range,ncol=Xi_ncol)
for (ilam in seq_len(Xi_ncol)) {
i_rhs_invV.dVdlam <- loglamInfo_blob$rhs_invV.dVdlam_list[[paste0(randit,"_",ilam)]] #.fill_rhs_invV.dVdlam(template=zerotemplate, urange=uirange[,ilam], invV.dV_info)
colit <- cum_Xi_cols[randit]+ilam
logdispInfoBlock[colit] <- .traceDB(.get_H_w.resid(object), t(rhs_invV.dVdrho),t(lhs_invV.dVdrho)) -
.traceAB(invV.dV_info$lhs_invV.dVdlam, i_rhs_invV.dVdlam,t(rhs_invV.dVdrho),t(lhs_invV.dVdrho))
}
}
logdispInfoBlock <- logdispInfoBlock[which(checklambda)] * sublambda #lambda * sum(invV.dVdlam*t(invV.dVdrho))
logdispInfo[dispcols$loglambda,dispcols$rho] <-
logdispInfo[dispcols$rho,dispcols$loglambda] <- logdispInfoBlock
}
}
## if (! is.null(dwdlogphi)) { ## currently length(phi)==1L && ! is.null(dvdlogphiMat)
if ("logphi" %in% dispnames) { ## more transparent, but error if mismatch of conditions
logphiInfo_blob <- .calc_logphiInfo(invV_factors, which=which(is_phiScalS), cum_nobs=cum_nobs,
phi_est=phi_est, w.resid=.get_H_w.resid(object),
spprec=.is_spprec_fit(object))
logdispInfo[dispcols$logphi,dispcols$logphi] <- logphiInfo_blob$logphiInfo
wresid_list <- logphiInfo_blob$wresid_list
if ("loglambda" %in% dispnames) {
sublambda <- .unlist(invV.dV_info$lambda_list)
# ncol_whichlambda <- sum(Xi_cols[which(checklambda)])
logdispInfoBlock <- matrix( nrow=sum(Xi_cols), ncol=length(wresid_list)) #
cum_n_u_h <- invV.dV_info$cum_n_u_h
zerotemplate <- rep(0,cum_n_u_h[nrand+1L])
col_ids <- rep(FALSE,sum(Xi_cols))
for (randit in which(checklambda)) { # subset of ranefs
u.range <- (cum_n_u_h[randit]+1L):(cum_n_u_h[randit+1L])
Xi_ncol <- Xi_cols[randit] # say '1 2' for ranCoefs
uirange <- matrix(u.range,ncol=Xi_ncol)
for (ilam in seq_len(Xi_ncol)) {
i_rhs_invV.dVdlam <- loglamInfo_blob$rhs_invV.dVdlam_list[[paste0(randit,"_",ilam)]] #.fill_rhs_invV.dVdlam(template=zerotemplate, urange=uirange[,ilam], invV.dV_info)
colit <- cum_Xi_cols[randit]+ilam
col_ids[colit] <- TRUE
# sum(diag(diag(w.resid) %*% lhs_invV.dVdlam)) - sum(diag(lhs_invV.dVdlam %*% invV_factors))
char_mv_jtS <- names(wresid_list)
for (block_jt in seq_along(char_mv_jtS)) {
char_mv_jt <- char_mv_jtS[block_jt]
j_lhs_invV.dVdphi <- logphiInfo_blob$lhs_invV.dVdphi_list[[char_mv_jt]]
logdispInfoBlock[colit,block_jt] <- phi_est[[as.integer(char_mv_jt)]] * (
.traceDB( wresid_list[[char_mv_jt]], invV.dV_info$lhs_invV.dVdlam, i_rhs_invV.dVdlam) -
.traceAB(invV.dV_info$lhs_invV.dVdlam, i_rhs_invV.dVdlam, j_lhs_invV.dVdphi, invV_factors$r_x_n) )
}
}
}
logdispInfoBlock <- .Dvec_times_matrix(sublambda, logdispInfoBlock[which(col_ids),,drop=FALSE])
logdispInfo[dispcols$loglambda,dispcols$logphi] <- logdispInfoBlock
logdispInfo[dispcols$logphi,dispcols$loglambda] <- t(logdispInfoBlock)
}
if ("rho" %in% dispnames) {
char_mv_jtS <- names(wresid_list)
for (jt in seq_along(char_mv_jtS)) {
char_mv_jt <- char_mv_jtS[jt]
j_lhs_invV.dVdphi <- logphiInfo_blob$lhs_invV.dVdphi_list[[char_mv_jt]]
coljt <- dispcols$logphi[jt] # column of full logdispInfo matrix
logdispInfo[dispcols$rho,coljt] <-
logdispInfo[coljt, dispcols$rho] <- phi_est[[as.integer(char_mv_jt)]] * .traceAB(lhs_invV.dVdrho,rhs_invV.dVdrho,
j_lhs_invV.dVdphi, invV_factors$r_x_n)
}
# phi_est * sum(invV.dVdrho * invV)
}
}
logdispInfo <- logdispInfo/2
resu <- .wrap_solve_logdispinfo(logdispInfo, object)
if (any( ! checklambda )) { ## if cols missing from logdisp_cov compared to dwdlogdisp
#ncd <- ncol(dwdlogdisp)
#full_logdisp_cov <- matrix(0,ncd,ncd)
# : alternative ncd below should be equivalent but more self contained
cols_in_logdisp_cov <- rep(checklambda,Xi_cols) ## which cols in dwdloglam match loglambda col in logdisp_cov
if ( ! is.null(dwdlogphi)) cols_in_logdisp_cov <- c(cols_in_logdisp_cov,TRUE) ## col for dwdlogphi
ncd <- length(cols_in_logdisp_cov)
full_logdisp_cov <- matrix(0,ncd,ncd)
full_logdisp_cov[cols_in_logdisp_cov,cols_in_logdisp_cov] <- resu$logdisp_cov
resu$logdisp_cov <- full_logdisp_cov
}
resu$dwdlogdisp <- dwdlogdisp
return(resu)
## more compact than storing ww %*% logdisp_cov %*% t(ww) which is nobs*nobs
}
}
.ugly_na.action_mv <- function(res) {
frames <- attr(res,"frame")
mv_res <- attr(res,"mv")
anyexclude <- FALSE
for (mv_it in seq_along(frames)) {
na.action_it <- attr(frames[[mv_it]],"na.action")
if (inherits(na.action_it,"exclude") ) {
anyexclude <- TRUE
mv_res[[mv_it]] <- .naresid.exclude(na.action_it, mv_res[[mv_it]] ) # INSERT NA's...
}
}
if (anyexclude) {
old_res <- res
res <- as.matrix(unlist(mv_res))
attrs <- attributes(old_res)
attrs <- attrs[setdiff(names(attrs), names(attributes(res)))]
attrs[names(attributes(res))] <- attributes(res)
attrs[["mv"]] <- mv_res
attributes(res) <- attrs
}
res
}
# modification from .make_new_corr_mats_NOT_ranCoef (which cannot be used natively)
.make_new_corr_mats_rhs_composite <- function(
newLv_env,
sub_corr_info,
corr.model=sub_corr_info$corr_types[[old_rd]], # here the one of the RHS !
newZAlist,
old_rd, fix_info,
object, new_rd, which_mats,
ranefs=attr(newZAlist,"exp_ranef_strings"), # "(.|.)" "adjacency" ...,
corrfamily_old_rd=.get_from_ranef_info(object)$corr_families[[old_rd]],
Lnn_not_Cnn, locdata,
for_mv,
spatial_terms=attr(object$ZAlist,"exp_spatial_terms"),
invCov_oldLv_oldLv_list,
latentL_blob) {
# tricky ___F I X M E____
# Previously I used old_kron_Y_Lmat=sub_corr_info$kron_Y_LMatrices[[old_rd]] as the Lmatrix for most of the cases below
# this was OK when newdata=NULL (cf commented code after this function)
# Now I try to construct properly new matrices, but this will fail insome of thecases where it worked.
# Try to identify cases where newdata=NULL (at least when univariate-response...)
namesTerms <- attr(newZAlist,"namesTerms")[[new_rd]]
if ( ! is.null(fix_info)) {
oldlevels <- colnames(fix_info$newZAlist[[old_rd]]) ## old_rd despite the "newZAlist" name: specific to fix_info
} else oldlevels <- colnames(object$ZAlist[[old_rd]])
newlevels <- colnames(newZAlist[[new_rd]]) # may have replicate at this point...
numnamesTerms <- length(namesTerms) ## 2 for random-coef
n_uniq <- length(oldlevels) %/% numnamesTerms
Uoldlevels <- oldlevels[seq_len(n_uniq)]
oldornew <- unique(c(Uoldlevels,newlevels))
repnames <- rep(oldornew, numnamesTerms)
newcols <- repnames %in% newlevels ## handle replicates (don't `[` newoldC using names !)
oldcols <- repnames %in% oldlevels
Unewlevels <- unique(colnames(newZAlist[[new_rd]]))
## rhs_C and compactcovmat for newoldC or newnewC
compactcovmat <- latentL_blob$compactcovmat
if (corr.model =="adjacency") {
if (length(setdiff(Unewlevels,Uoldlevels))) stop(paste0("Found new levels for a '",corr.model,"' random effect."))
### Context:
## We want to compute Z . C . C^{-1/2} v.
## by reproducing the logic of univariate-resp models where
## in permuted case that can be written
## ZP . P' C P . P' C^{-1/2} v [where the P is P_Q]
## = ZP . L_Q^(-T) L_Q^{-1} . L_Q v
## so newZAlist contains permuted ZP, newoldC contains L_Q^(-T) L_Q^{-1} which is permuted,
## and L_Q v is always provided by .calc_invL_coeffs()
### distinct features of this composite adjacency case is that
## (1) difference from composite corrMatrix is that there is no old_kron_Y_Lmat;
## (2) the strucList element represents the Kronecker product, not its LHS =>
## using dimnames is bound to fail since they are repeated.
uuC <- .tcrossprod(object$strucList[[old_rd]], perm=TRUE) ## Might eventually reconstruct permuted (consistent with perm of cols of Z) corrMatrix from its CHM factor
if (which_mats$no) newLv_env$cov_newLv_oldv_list[[new_rd]] <- structure(uuC,
ranefs=ranefs[[new_rd]])
# assign to newLv_env$cov_newLv_oldv_list, cov_newLv_newLv_list, diag_cov_newLv_newLv_list
# although newnew code seems missing, it is not needed (there is a test of predVar cov=TRUE)
} else if (corr.model=="corrMatrix") {
# # "test-predVar-Matern-corrMatrix" shows newdata working with corrMatrix, when newlevels are within oldlevels
# see also test-composite-extra
if (length(setdiff(Unewlevels,Uoldlevels))) stop(paste0("Found new levels for a '",corr.model,"' random effect."))
# all newlevels must be in oldlevels hence the construction of the matrix is a little simplified
uuC <- .tcrossprod(sub_corr_info$kron_Y_LMatrices[[old_rd]], perm=TRUE) ## Can reconstruct permuted (consistent with perm of cols of Z) corrMatrix from its CHM factor
rownames(uuC) <- colnames(uuC) <- Uoldlevels
#if (TRUE) {
.assign_newLv_for_newlevels_corrMatrix(uuC, latentL_blob=latentL_blob,
newlevels=Unewlevels,
newLv_env, new_rd, corr.model, which_mats, ranefs,
Lnn_not_Cnn=Lnn_not_Cnn)
# } else {
# newoldC <- .makelong_kronprod(compactcovmat, kron_Y=uuC[Unewlevels,,drop=FALSE])
# #colnames(newoldC) <- rep(oldlevels, numnamesTerms) ## these will be needed by .match_old_new_levels() or .assign_newLv_for_newlevels_corrMatrix()
# if (which_mats$no) {
# newLv_env$cov_newLv_oldv_list[[new_rd]] <- structure(newoldC, ranefs=ranefs[[new_rd]])
# # if (inherits(newoldC,"bigq")) attr(newLv_env$cov_newLv_oldv_list[[new_rd]],
# # "DIMNAMES") <- list(repnames[newcols],repnames[oldcols]) ## these will be needed by .match_old_new_levels()
# }
#
# if (which_mats$nn[new_rd]) {
# if (Lnn_not_Cnn) {
# compactL <- latentL_blob$design_u
# if (is.null(compactL)) compactL <- # on compAR1fit example
# .wrap_solve_warn(X=t(as.matrix(latentL_blob$compactchol_Q_w)),
# warnmess="Forcing inversion of singular precision matrix.")
# newnewL <- .makelong(compactL,
# longsize=length(Unewlevels)*numnamesTerms,
# kron_Y=old_kron_Y_Lmat[Unewlevels,Unewlevels, drop=FALSE])
# newLv_env$L_newLv_newLv_list[[new_rd]] <- newnewL
# } else {
# newnewC <- .makelong(compactcovmat,
# longsize=length(Unewlevels)*numnamesTerms,
# kron_Y=uuC[Unewlevels,Unewlevels])
# newLv_env$cov_newLv_newLv_list[[new_rd]] <- newnewC #names presumably not used
# #if (inherits(newoldC,"bigq")) attr(diag_cov_newLv_newLv_list[[new_rd]], "DIMNAMES") <- list(repnames[newcols],repnames[newcols])
# }
# } else {
# newnewC <- .makelong(compactcovmat, longsize=length(Unewlevels)*numnamesTerms,
# kron_Y=uuC[Unewlevels,Unewlevels ,drop=FALSE])
# nc <- ncol(newnewC)
# diagPos <- seq.int(1L,nc^2,nc+1L) ## it's not only an efficient syntax, that's the only way to extract the diag from bigq...
# #if (inherits(newoldC,"bigq")) {
# # newLv_env$diag_cov_newLv_newLv_list[[new_rd]] <- newoldC[[]][diagPos[newcols]] ## with the [[]] to access the raw vector
# #} else
# newLv_env$diag_cov_newLv_newLv_list[[new_rd]] <- newnewC[diagPos]
# }
# }
} else if (corr.model == "IMRF") {
# stop("code missing for composite IMRF ranef") # ___F I X M E____ that will indeed fail. Let the actual problem occur
## in this case a new A matrix (by .get_new_AMatrices()) must be computed (elsewhere: it goes in newZAlist)
## but the corr matrix between the nodes is unchanged as node positions do not depend on response position => nn=no
.make_new_corr_lists_IMRFs(newLv_env, which_mats,
object, # shows that direct access to $strucList[[old_rd]] may be useful in generic code
ranefs, new_rd, old_rd)
} else if ((corr.model == "corrFamily" && ! identical(corrfamily_old_rd$levels_type,"time_series")) &&
! is.null(corrfamily_old_rd$make_new_corr_lists)) {
# User must provide corrfamily_old_rd$make_new_corr_lists() code that handles composite ranefs
corrfamily_old_rd$make_new_corr_lists(newLv_env=newLv_env, which_mats=which_mats, ranFix=object$ranFix,
ranefs=ranefs,
newZAlist=newZAlist, new_rd=new_rd, old_rd=old_rd,
corrfamily=corrfamily_old_rd, # easy access to elements added by .preprocess_corrFamily
object=object, # avoid using it... but see comment on IMRFs above
Lnn_not_Cnn=Lnn_not_Cnn # may not be handled in any useful way.
) # Should write in newLv_env
} else if ( ! is.null(corrfamily_old_rd$calc_corr_from_dist) ) {
## all models where a correlation matrix must be computed from a distance matrix => $calc_corr_from_dist() needed
## Notably, Matern...
old_char_rd <- as.character(old_rd)
# olduniqueGeo needed in all cases
if ( ! is.null(fix_info)) {
olduniqueGeo <- fix_info$newuniqueGeo[[old_char_rd]]
} else {
olduniqueGeo <- .get_old_info_uniqueGeo(object, char_rd=old_char_rd)
}
if ( is.data.frame(locdata)) {
geonames <- colnames(olduniqueGeo)
newuniqueGeo <- locdata[,geonames,drop=FALSE] ## includes nesting factor
if (for_mv) { # ...but for mv fits this will be a problem in .compute_ZAXlist() as locnr>locnc
# mv fit
newuniqueGeo <- unique(newuniqueGeo)
rownames(newuniqueGeo) <- .pasteCols(t(newuniqueGeo))
} else if (attr(newZAlist[[new_rd]],"Z_levels_type") != "seq_len") {
# univariate-response model such as Matern(LHS implies ranCoefs|.)
## The location in newuniqueGeo may not be unique despite the name.
## The costs of reducing to unique values may outweight the benefits.
## The 'decision' has been made when creating the new ZA, through its levels type.
## Simple test case where levels_type is seq_len:
## example(bboptim) [Matern with replicate pairs in some locations; locations are not unique in 'newuniqueGeo']
## conversely the predVar computations in block
## 'verif .calc_Evar() with ranCoefs...' in test-devel-predVar-ranCoefs
## stop if unique() is not applied (for the AR1 ranef).
## The Z_levels_type attr should have been used more to secure the code...
newuniqueGeo <- unique(newuniqueGeo)
}
} else { ## locdata is 'preprocessed' list of arrays (tested by get_predCov_var_fix() tests)
newuniqueGeo <- locdata[[as.character(old_rd)]] ## preprocessed, [,geonames,drop=FALSE] not necess ## includes nesting factor
geonames <- colnames(newuniqueGeo)
}
## ... distance matrix and then call to correl fn ...
control_dist_rd <- .get_control_dist(object, old_char_rd)
if (corr.model=="AR1" ||
identical(corrfamily_old_rd$levels_type,"time_series") # identical because Matern, etc. are corrfamily_old_rd without levels_type
# ARp) and ARMA() previously reached here, but no longer as they have a $make_new_corr_lists now (but this does not handle nesting)
) {
# Does not use a kron_Y stored (or not even stored for spprec) in the object: reconstructs matrix from 'distance' matrix.
### older, non nested code:
# if (which$no) resu$uuCnewold <- proxy::dist(newuniqueGeo,olduniqueGeo)
# if (which$nn) resu$uuCnewnew <- proxy::dist(newuniqueGeo)
### new code recycling .get_dist_nested_or_not, but not fastest nor most transparent (fixme?)
onGeo <- rbind(newuniqueGeo,olduniqueGeo) # includes nesting factor
if (object$spaMM.version < "2.2.118") {
blob <- .get_dist_nested_or_not(object$spatial_term, data=onGeo, distMatrix=NULL, uniqueGeo=NULL,
dist.method=control_dist_rd$dist.method, as_matrix=TRUE,
needed=c(distMatrix=TRUE), geo_envir=NULL)
} else {
blob <- .get_dist_nested_or_not(spatial_terms[[old_rd]], data=onGeo, distMatrix=NULL,
uniqueGeo=NULL, ## FIXME provide uniqueGeo 'to save time'(??) ?
dist.method=control_dist_rd$dist.method, as_matrix=TRUE,
needed=c(distMatrix=TRUE), geo_envir=NULL)
}
## we merged old and new so need to get the respective cols (which may overlap)
uli_onGeo <- .ULI(onGeo) # this should give row and columns in the blob ## FIXME how to make sure of that? .get_dist_nested_or_not must use .ULI()
uli_new <- uli_onGeo[seq(nrow(newuniqueGeo))]
uli_old <- uli_onGeo[-seq(nrow(newuniqueGeo))]
if (which_mats$no) uuCnewold <- blob$distMatrix[uli_new,uli_old,drop=FALSE] ## rows match the newZAlist, cols match th u_h
if (which_mats$nn[new_rd]) {
uuCnewnew <- blob$distMatrix[uli_new,uli_new,drop=FALSE]
}
} else if (corr.model %in% c("Cauchy", "Matern")) {
### rho only used to compute scaled distances
rho <- .get_cP_stuff(object$ranFix,"rho", which=old_char_rd)
if ( ! is.null(rho_mapping <- control_dist_rd$rho.mapping)
&& length(rho)>1L ) rho <- .calc_fullrho(rho=rho,coordinates=geonames,rho_mapping=rho_mapping)
# : if control_dist_rd comme from the call (vs from moreargs) rho_mapping may still be NULL
# and then the code assumes that calling .calc_fullrho() is not necessary (that seems OK)
## rows from newuniqueGeo, cols from olduniqueGeo:
txt <- paste(c(spatial_terms[[old_rd]][[2]][[3]])) ## the RHS of the ( . | . ) # c() to handle very long RHS
if (length(grep("%in%",txt))) { # nested geostatistical effect
coord_within <- .extract_check_coords_within(spatial_term=spatial_terms[[old_rd]])
msd.arglist <- list(uniqueGeo=newuniqueGeo[,coord_within,drop=FALSE],
uniqueGeo2=olduniqueGeo[,coord_within,drop=FALSE],
rho=rho,return_matrix=TRUE)
} else msd.arglist <- list(uniqueGeo=newuniqueGeo,uniqueGeo2=olduniqueGeo,
rho=rho,return_matrix=TRUE)
# If control_dist_rd$dist.method is NULL, do not over-write the non-NULL default of make_scaled_dist():
if ( ! is.null(dist.method <- control_dist_rd$dist.method)) msd.arglist$dist.method <- dist.method
if (which_mats$no) {
uuCnewold <- do.call(make_scaled_dist,msd.arglist)
dimnames(uuCnewold) <- list(Unewlevels, Uoldlevels)
} ## ultimately allows products with Matrix ## '*cross*dist' has few methods, not even as.matrix
if (which_mats$nn[new_rd]) {
msd.arglist$uniqueGeo2 <- NULL
if (nrow(msd.arglist$uniqueGeo)==1L) {
uuCnewnew <- matrix(0) ## trivial _distance_ matrix for single point (converted to trivial cov below!)
} else uuCnewnew <- do.call(make_scaled_dist,msd.arglist)
dimnames(uuCnewnew) <- list(Unewlevels, Unewlevels)
}
if (length(grep("%in%",txt))) { # nested geostatistical effect
onGeo <- rbind(newuniqueGeo,olduniqueGeo) # includes nesting factor
isInf <- .get_dist_nested_or_not(spatial_term=spatial_terms[[old_rd]],
data=onGeo, distMatrix=NULL,
uniqueGeo=NULL,
dist.method = dist.method,needed=c(notSameGrp=TRUE),
geo_envir=NULL)$notSameGrp
## we merged old and new so need to get the respective cols (which may overlap)
uli_onGeo <- .ULI(onGeo) # this should give row and columns in the blob ## FIXME how to make sure of that? .get_dist_nested_or_not must use .ULI()
uli_new <- uli_onGeo[seq(nrow(newuniqueGeo))]
uli_old <- uli_onGeo[-seq(nrow(newuniqueGeo))]
if (which_mats$no) {
isInfno <- isInf[uli_new,uli_old,drop=FALSE]
uuCnewold[isInfno] <- Inf ## ultimately allows products with Matrix ## '*cross*dist' has few methods, not even as.matrix
}
if (which_mats$nn[new_rd]) {
isInfnn <- isInf[uli_new,uli_new,drop=FALSE]
uuCnewnew[isInfnn] <- Inf
}
}
} else stop("Unhandled 'corr.model' (make_new_corr_lists() missing from corrFamily descriptor?).")
# ... Finally fill the cov lists ...
if (object$spaMM.version<"2.4.49") {
if (which_mats$no) newLv_env$cov_newLv_oldv_list[[new_rd]] <- structure(.calc_corr_from_dist(uuCnewold, object, corr.model,char_rd=old_char_rd),
corr.model=corr.model,
ranefs=ranefs[[new_rd]])
if (which_mats$nn[new_rd]) newLv_env$cov_newLv_newLv_list[[new_rd]] <- .calc_corr_from_dist(uuCnewnew, object, corr.model,char_rd=old_char_rd)
} else {
if (which_mats$no) {
cov_newLv_oldv <- corrfamily_old_rd$calc_corr_from_dist(
ranFix=object$ranFix, char_rd=old_char_rd, distmat=uuCnewold)
# if (attr(strucList[[old_rd]],"need_gmp") &&
if (inherits(invCov_oldLv_oldLv_list[[old_rd]],"bigq")) cov_newLv_oldv <- structure(gmp::as.bigq(cov_newLv_oldv),
DIMNAMES=dimnames(cov_newLv_oldv))
cov_newLv_oldv <- .makelong_kronprod(compactcovmat, kron_Y=cov_newLv_oldv[Unewlevels,])
newLv_env$cov_newLv_oldv_list[[new_rd]] <- structure(cov_newLv_oldv, corr.model=corr.model, ranefs=ranefs[[new_rd]])
}
if (which_mats$nn[new_rd]) {
uuCnewnew <- corrfamily_old_rd$calc_corr_from_dist(
ranFix=object$ranFix, char_rd=old_char_rd, distmat=uuCnewnew)
if (Lnn_not_Cnn) {
Lmat <- mat_sqrt(m=uuCnewnew)
raw_levels <- .pasteCols(x=t(newuniqueGeo))
rownames(Lmat) <- colnames(Lmat) <- raw_levels # same levels in .as_factor.
compactL <- latentL_blob$design_u
if (is.null(compactL)) compactL <- # on compAR1fit example
.wrap_solve_warn(X=t(as.matrix(latentL_blob$compactchol_Q_w)),
warnmess="Forcing inversion of singular precision matrix.")
newnewL <- .makelong(compactL,
longsize=length(Unewlevels)*numnamesTerms,
kron_Y=Lmat[Unewlevels,Unewlevels])
newLv_env$L_newLv_newLv_list[[new_rd]] <- newnewL
} else {
newnewC <- .makelong(compactcovmat,
longsize=length(Unewlevels)*numnamesTerms,
kron_Y=uuCnewnew)
newLv_env$cov_newLv_newLv_list[[new_rd]] <- newnewC #names presumably not used
}
} else { # (__F I X M E___ ugly inefficient block of code: compute full matrix to get its diag)
uuC <- .tcrossprod(sub_corr_info$kron_Y_LMatrices[[old_rd]], perm=TRUE) ## Can reconstruct permuted (consistent with perm of cols of Z) corrMatrix from its CHM factor
rownames(uuC) <- colnames(uuC) <- Uoldlevels
# uuC <- corrfamily_old_rd$calc_corr_from_dist(
# ranFix=object$ranFix, char_rd=old_char_rd, distmat=uuC)
newnewC <- .makelong(compactcovmat, longsize=length(Unewlevels)*numnamesTerms,
kron_Y=uuC)
nc <- ncol(newnewC)
diagPos <- seq.int(1L,nc^2,nc+1L) ## it's not only an efficient syntax, that's the only way to extract the diag from bigq...
newLv_env$diag_cov_newLv_newLv_list[[new_rd]] <- newnewC[diagPos]
}
}
} # END all models where a correlation matrix must be computed from a distance matrix
}
{ ############## Memo old code:
# namesTerms <- attr(newZAlist,"namesTerms")[[new_rd]]
# # if ( ! is.null(fix_info)) {
# # oldlevels <- colnames(fix_info$newZAlist[[old_rd]]) ## old_rd despite the "newZAlist" name: specific to fix_info
# # } else oldlevels <- colnames(object$ZAlist[[old_rd]])
# # oldlevels <- unique(oldlevels)
# # newlevels <- unique(colnames(newZAlist[[new_rd]]))
# # # "test-predVar-Matern-corrMatrix" shows newdata working with corrMatrix, when newlevels are within oldlevels
# # if (length(setdiff(newlevels,oldlevels))) stop(paste0("Found new levels for a '",corr.model,"' random effect."))
# # kron_Y <- .tcrossprod(kron_Y, perm=TRUE) ## Can reconstruct permuted (consistent with perm of cols of Z) corrMatrix from its CHM factor
# # colnames(kron_Y) <- rownames(kron_Y) <- oldlevels
# kron_Y <- .tcrossprod(kron_Y, perm=TRUE) ## Can reconstruct permuted (consistent with perm of cols of Z) corrMatrix from its CHM factor
# rownames(kron_Y) <- colnames(kron_Y) <- Uoldlevels # oldlevels <- colnames(kron_Y)
# newlevels <- unique(colnames(newZAlist[[new_rd]]))
# # "test-predVar-Matern-corrMatrix" shows newdata working with corrMatrix, when newlevels are within oldlevels
# if (length(setdiff(newlevels,Uoldlevels))) stop(paste0("Found new levels for a '",corr.model,"' random effect."))
# # all newlevels must be in oldlevels hence the construction of the matrix is a little simplified
# compactcovmat <- latentL_blob$compactcovmat
# newoldC <- .makelong_kronprod(compactcovmat, kron_Y=kron_Y[newlevels,])
# #colnames(newoldC) <- rep(oldlevels, numnamesTerms) ## these will be needed by .match_old_new_levels() or .assign_newLv_for_newlevels_corrMatrix()
# #rownames(newoldC) <- rep(newlevels, numnamesTerms)
# if (which_mats$no) {
# newLv_env$cov_newLv_oldv_list[[new_rd]] <- structure(newoldC, ranefs=ranefs[[new_rd]])
# # if (inherits(newoldC,"bigq")) attr(newLv_env$cov_newLv_oldv_list[[new_rd]],
# # "DIMNAMES") <- list(repnames[newcols],repnames[oldcols]) ## these will be needed by .match_old_new_levels()
# }
#
# if (which_mats$nn[new_rd]) {
# if (Lnn_not_Cnn) {
# compactL <- latentL_blob$design_u
# if (is.null(compactL)) compactL <- # on compAR1fit example
# .wrap_solve_warn(X=t(as.matrix(latentL_blob$compactchol_Q_w)),
# warnmess="Forcing inversion of singular precision matrix.")
# newnewL <- .makelong(compactL,
# longsize=length(newlevels)*numnamesTerms,
# kron_Y=kron_Y[newlevels,newlevels])
# newLv_env$L_newLv_newLv_list[[new_rd]] <- newnewL
# } else {
# newnewC <- .makelong(compactcovmat,
# longsize=length(newlevels)*numnamesTerms,
# kron_Y=kron_Y[newlevels,newlevels])
# newLv_env$cov_newLv_newLv_list[[new_rd]] <- newnewC #names presumably not used
# #if (inherits(newoldC,"bigq")) attr(diag_cov_newLv_newLv_list[[new_rd]], "DIMNAMES") <- list(repnames[newcols],repnames[newcols])
# }
# } else {
# newnewC <- .makelong(compactcovmat, longsize=length(newlevels)*numnamesTerms, kron_Y=kron_Y[newlevels,newlevels])
# nc <- ncol(newnewC)
# diagPos <- seq.int(1L,nc^2,nc+1L) ## it's not only an efficient syntax, that's the only way to extract the diag from bigq...
# #if (inherits(newoldC,"bigq")) {
# # newLv_env$diag_cov_newLv_newLv_list[[new_rd]] <- newoldC[[]][diagPos[newcols]] ## with the [[]] to access the raw vector
# #} else
# newLv_env$diag_cov_newLv_newLv_list[[new_rd]] <- newnewC[diagPos]
# }
}
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.