Nothing
R/postfit_internals.R
In spaMM: Mixed-Effect Models, with or without Spatial Random Effects
Defines functions
.u_h_v_h_from_v_h
.make_CorrEst_and_RanFix
.calc_newZACvar
.get_beta_cov_info
.calc_beta_cov_info_from_sXaug
.calc_beta_cov_info_others
.calc_invV_factors
.get_logdispObject
.calc_logdispObject
.calc_d2hdv2_info
..calc_d2hdv2_info
.calc_dvdloglamMat_new
.calc_dvdlogphiMat_new
.get_info_crits
DoF
.calc_p_rdisp
.calc_cAIC_pd_others
.calc_cAIC_pd_from_sXaug
.calc_cAIC_pd_spprec
.calc_pd_product
Documented in
DoF
.calc_pd_product <- function(tcrossfac_v_beta_cov, Md2clikdvb2, blockSize=5000L) {
if ((nc <- ncol(tcrossfac_v_beta_cov))>(blockSize)) {
## We reached this point by sparse matrix computations. We need to save memory in the following dense computations,
message(paste0("Conditional-AIC computation requires operations on a large matrix (square with dimension ",nc,"),\n",
"which may take a bit of time. Use 'also_cAIC=FALSE' to avoid it."))
slices <- unique(c(seq(0L,nc,blockSize),nc))
nslices <- length(slices)-1L
it <- 0L ## 'global definition' for Rcmd check
foreach_args <- list(it = seq_len(nslices), .combine = "sum")
foreach_blob <- do.call(foreach::foreach,foreach_args)
abyss <- foreach::`%do%`(foreach_blob, Sys.setenv(LANG = "en"))
barstyle <- eval(spaMM.getOption("barstyle"))
progrbar_setup <- .set_progrbar(max = nslices, style = barstyle, char="s")
pd <- foreach::`%do%`(foreach_blob, {
slice <- (slices[it]+1L):slices[it+1L]
tmp <- t(.crossprod(tcrossfac_v_beta_cov, Md2clikdvb2[,slice]))
if (barstyle) progrbar_setup$progress(it)
return(sum(tcrossfac_v_beta_cov[slice,] * tmp))
})
if (barstyle) close(progrbar_setup$pb)
# We could parallelize using %dopar% (twice) but is that worth the overhead?
} else {
# logic of following code is
# pd = sum(diag(solve(Md2hdbv2,Md2clikdbv2[c(113:114,1:112),])))
# = sum(diag((tcrossprod(R_invMd2hdvb2)[c(113:114,1:112),] %*% Md2clikdbv2[sort.list(c(113:114,1:112)),])))
# = sum(diag(R_invMd2hdvb2[c(113:114,1:112),] %*% (crossprod(R_invMd2hdvb2, Md2clikdbv2[sort.list(c(113:114,1:112)),]))))
# = sum((R_invMd2hdvb2[c(113:114,1:112),] * t(crossprod(R_invMd2hdvb2, Md2clikdbv2[sort.list(c(113:114,1:112)),]))))
# but we directly use v,b order rather than b,v
pd <- t(.crossprod(tcrossfac_v_beta_cov,Md2clikdvb2))
pd <- sum(tcrossfac_v_beta_cov * pd)
# that is actually the logic of .traceAB:
# pd = .traceAB(tcrossfac_v_beta_cov, t(tcrossfac_v_beta_cov),Md2clikdvb2,diag(nrow =ncol(Md2clikdvb2)))
# but the matrices are square so no use here
}
return(pd)
}
.calc_cAIC_pd_spprec <- function(object) {
H_w.resid <- .get_H_w.resid(object)
#ZAL <- .compute_ZAL(XMatrix=object$strucList, ZAlist=object$ZAlist,as_matrix=.eval_as_mat_arg.HLfit(object))
ZAL <- get_ZALMatrix(object, force_bind=FALSE) # allows ZAXlist; but if a non-ZAXlist is already in the $envir, no effect; + we will solve(chol_Q, Diagonal()) so the gain is not obvious
if ( ncol(X_ori <- model.matrix(object)) ) {
M12 <- .crossprod(ZAL, .Dvec_times_m_Matrix(H_w.resid, X_ori), as_mat=TRUE)
Md2clikdvb2 <- rbind2(cbind2(as.matrix(.safe_ZtWZwrapper(ZAL,H_w.resid)), M12), ## this .ZtWZwrapper() takes time
cbind2(t(M12), as.matrix(.safe_ZtWZwrapper(X_ori,H_w.resid))))
# _FIXME_ any way to avoid formation of this matrix ? Or otherwise message() ?
} else {
Md2clikdvb2 <- as.matrix(.safe_ZtWZwrapper(ZAL,H_w.resid))
}
tcrossfac_v_beta_cov <- .calc_Md2hdvb2_info_spprec(X.pv=X_ori, envir=object$envir,
which="tcrossfac_v_beta_cov")
# not triang if we used sparse QR. Following code should not assume triangularity
pd <- .calc_pd_product(tcrossfac_v_beta_cov, Md2clikdvb2)
return(pd)
}
# From Matrix 1.6-2, "Methods for generic functions rbind2, cbind2, %*%, %&%, crossprod, and tcrossprod determine the class of the result using more strict rules, designed to avoid “surprising” coercions where possible. Notably, operations involving RsparseMatrix now return an RsparseMatrix in more cases. TsparseMatrix and diagonalMatrix may be handled as CsparseMatrix or as RsparseMatrix, depending on context. "
# Maybe some old as.matrix() coertions were previously useful.
# But now I can remove some of then with great gains.
.calc_cAIC_pd_from_sXaug <- function(object) {
if (is.matrix(beta_cov_info <- object$envir$beta_cov_info) || ## matrix is old format, should be a list now
is.null(tcrossfac_beta_v_cov <- beta_cov_info$tcrossfac_beta_v_cov)) {
tcrossfac_beta_v_cov <- .get_beta_cov_info(object)$tcrossfac_beta_v_cov
}
X_ori <- model.matrix(object)
pforpv <- ncol(X_ori)
nc <- ncol(tcrossfac_beta_v_cov)
n_u_h <- nc - pforpv
seqp <- seq_len(pforpv)
perm <- c(pforpv+seq_len(n_u_h), seqp)
tcrossfac_v_beta_cov <- tcrossfac_beta_v_cov[perm,,drop=FALSE] # useful to keep it for predVar computations?
H_w.resid <- .get_H_w.resid(object)
#ZAL <- .compute_ZAL(XMatrix=object$strucList, ZAlist=object$ZAlist,as_matrix=.eval_as_mat_arg.HLfit(object))
ZAL <- get_ZALMatrix(object, force_bind=FALSE) # allows ZAXlist; but if a non-ZAXlist is already in the $envir, no effect; + we will solve(chol_Q, Diagonal()) so the gain is not obvious
if ( ncol(X_ori) ) {
M12 <- .crossprod(ZAL, .Dvec_times_m_Matrix(H_w.resid, X_ori), as_mat=TRUE)
if (.spaMM.data$options$Matrix_old) {
Md2clikdvb2 <- rbind2(cbind2(as.matrix(.safe_ZtWZwrapper(ZAL,H_w.resid)), M12),
cbind2(t(M12), as.matrix(.safe_ZtWZwrapper(X_ori,H_w.resid))))
message("If this step was slow, updating the Matrix package to a version >=1.6-2 might be helpful.")
} else { # AIC(bigranefs) showed the gain of removing the as.matrix()s.
Md2clikdvb2 <- rbind2(cbind2(.safe_ZtWZwrapper(ZAL,H_w.resid), M12),
cbind2(t(M12), .safe_ZtWZwrapper(X_ori,H_w.resid)))
}
# Old comments:
# .ZtWZwrapper(ZAL...) may take time
} else {
Md2clikdvb2 <- as.matrix(.safe_ZtWZwrapper(ZAL,H_w.resid))
}
# not triang if we used sparse QR. Following code should not assume triangularity
pd <- .calc_pd_product(tcrossfac_v_beta_cov, Md2clikdvb2)
return(pd)
}
# Here, using qr.Md2hdbv2$pivot implies that qr() was run on a base matrix.
# So as.matrix() may remain necessary unless further changes are made.
.calc_cAIC_pd_others <- function(X.pv, ZAL, w.resid, d2hdv2, blockSize=1000L) {
if ( ncol(X.pv) ) { ## the projection matrix for the response always includes X even for REML!
hessnondiag <- .crossprod(ZAL, .Dvec_times_m_Matrix(w.resid, X.pv))
Md2hdbv2 <- as.matrix(rbind2(cbind2(.safe_ZtWZwrapper(X.pv,w.resid), t(hessnondiag)),
cbind2(hessnondiag, - d2hdv2)))
Md2clikdbv2 <- as.matrix(rbind2(cbind2(.safe_ZtWZwrapper(X.pv,w.resid), t(hessnondiag)),
cbind2(hessnondiag, .safe_ZtWZwrapper(ZAL,w.resid))))
} else {
Md2hdbv2 <- - d2hdv2
Md2clikdbv2 <- as.matrix(.safe_ZtWZwrapper(ZAL,w.resid))
}
if (inherits(Md2hdbv2,"diagonalMatrix")) {
pd <- sum(diag(Md2clikdbv2)/diag(Md2hdbv2)) ## is sum(diag(solve(Md2hdbv2,Md2clikdbv2)))
} else {
## dans un LMM avec estimation ML, pd = sum(lev_phi), mais pas de simplif plus generale
if ((nc <- ncol(Md2hdbv2))>(blockSize)) {
message(paste0("Conditional-AIC computation requires operations on a large dense matrix (square with dimension ",nc,"),\n",
"which may take a lot of time. Use 'also_cAIC=FALSE' to avoid it."))
}
## if we reached this point with a huge dense matrix, then we have huge memory, hence saving memory may not be the issue,
# But we might save computation by computing only the requested diagonal (one qr, one backsolve, one sum(. * .))
# using Md2hdbv2=QRP, tr=sum_i (invP invR Qt(Md2hdbv2))_ii
qr.Md2hdbv2 <- try(qr(Md2hdbv2))
if (inherits(qr.Md2hdbv2,"try-error")) {
warning("Computation of cAIC/GoF df's failed because the information matrix appears singular.")
pd <- NA
} else {
solveR <- try(backsolve(qr.R(qr.Md2hdbv2),diag(nrow =nc)))
if (inherits(solveR,"try-error")) {
warning("Computation of cAIC/GoF df's failed because the information matrix appears singular.")
# determinant(qrR,logarithm=FALSE)$modulus > 1e-14
pd <- NA
} else {
# using inv(RP)= inv(R[,perm <- sort.list(.$pivot)]) = inv(R)[perm,] (as in h9[,perm] %*% solve(h9)[perm,] for any perm)
# though use of pivot may occur only for nearly singular matrices (poorly doc)
pd <- sum(solveR[sort.list(qr.Md2hdbv2$pivot),] * t(qr.qty(qr.Md2hdbv2,Md2clikdbv2)))
}
}
}
return(pd)
}
# Consistent with Saefken et al, only in terms of theta and mu. The linear predictor eta and link do not appear.
# Simulations for that paper use predictions eta only in the poisson(log) case;
# In the exponential (->Gamma(log)) they do use the theta deduced as -1/mu, not the eta. => cAIC4:::conditionalBootstrap is odd...
.calc_boot_AIC_dfs <- function (object, nsim, type="residual", seed=NULL, # (_F I X M E__) does not handle mv fits ? Could be docu'ed at least
fit_env=NULL, ...) {
family <- family(object)
if ( ! family$flags$exp) {
stop("Bootstrap bias correction not implemented for families not from GLM (exponential family) class.") # assuming $exp methods are always available for GLMs (but see negbin2_dvl)
} else if (family$family %in% c("negbin","negbin2")) {
warning("Bootstrap method does not strictly apply for negbin models.")
# Ad hoc version of .theta.mu.canonical that gives correct results in Poisson limit:
# as if the canonical parm was -log(1+shape/mu) + log(shape)
# This is not the commonly used formula, but still seems consistent with general theory
# (theta def'ed up to some forms of constants)
# devel/negbin_link has a version of HLfit_internals with this 'revised' canonical link.
thetafn <- function(x) {
refit <- update_resp(object, newresp=x)
muFREQS <- predict(refit, type="response")
shape <- environment(family(refit)$aic)$shape # allows shape to depend on refitted value
th <- -log(1+shape/muFREQS)+log(shape)
th
}
} else {
thetafn <- function(x) {
refit <- update_resp(object, newresp=x)
muFREQS <- predict(refit, type="response") # predict(refit(object, newresp = x))
.theta.mu.canonical(muFREQS,family(object)) # would return mu for non-GLMs
}
}
bootsims <- simulate(object, nsim = nsim, type = type, verbose=FALSE, seed=seed) # the corresponding lmer code returns a data frame
thetas <- dopar(bootsims, thetafn, fit_env=fit_env, ...)
if (inherits(thetas[[1]],"simpleError")) {
stop(paste("From dopar() call: ",
thetas[[1]]$message))
} else if (inherits(thetas,"try-error")) {
stop('dopar() call failed: use control=list(.errorhandling = "pass")\n to obtain more informative error message')
}
# if (is.factor(bootsims[1])) dataMatrix <- as.numeric(dataMatrix) - 1 # in cAIC4:::conditionalBootstrap
bootsims <- bootsims - rowMeans(bootsims)
phis <- residVar(object, which="phi") # For Gamma(), the phis are still those of he canonical form of the exponential family,
# not the full residual variance (name 'get_residVar' is ambiguous)
bootBC <- sum(colSums(thetas * bootsims/phis))/(nsim-1) # Handles the case Where the phis are heteroscedastic,
return(bootBC)
}
.calc_p_rdisp <- function(object, dfs=object$dfs) {
p_phi <- dfs[["p_fixef_phi"]]
if ( ! is.null(resid_fits <- object$resid_fits)) { # i.e mv fit =>$resid_fit*s*
p_phi <- sum(na.omit(unlist(p_phi)),
sum(unlist(lapply(resid_fits, `[[`, x="dfs"), recursive = TRUE, use.names = FALSE)) )
} else if ( ! is.null(resid_fit <- object$resid_fit)) {
# input p_phi (above) is typically set to NA, and will be ignored
p_phi <- sum(.unlist(resid_fit$dfs)) ## phi_pd is relevant only for measuring quality of prediction by the resid_fit!
} else p_phi <- sum(.unlist(dfs[["p_fixef_phi"]])) # .unlist() for mv
names_est_ranefPars <- unlist(.get_methods_disp(object))
fam_disp_parsnames <- intersect(names_est_ranefPars,c("NB_shape","COMP_nu","beta_prec",
"Tw_index","Tw_link"))
if (length(fam_disp_parsnames)) {
p_GLM_family <- length( # compatible with mv:
.unlist(.get_outer_inits_from_fit(object, keep_canon_user_inits = FALSE)[fam_disp_parsnames]))
p_phi <- p_phi+p_GLM_family ## effectively part of the model for residual error structure
}
p_phi
}
DoF <- function(object) {
dfs <- object$dfs
if ( ! inherits(dfs,"list")) dfs <- as.list(dfs) ## back compatibility
dfs[["p_fixef"]] <- dfs[["pforpv"]]
dfs[["p_rdisp"]] <- .calc_p_rdisp(object, dfs)
dfs["pforpv"] <- dfs["p_phi"] <- dfs["p_fixef_phi"] <- list(NULL) # make sure to erase any previously pre-existing redundant info
if (object$models[[1]]=="etaHGLM") {
if (is.null(p_corrPars <- dfs[["p_corrPars"]])) { ## back compatibility code,
dfs[["p_corrPars"]] <- length(which(unlist(attr(object$CorrEst_and_RanFix,"type")$corrPars, use.names = FALSE) %in% c("outer","var")))
# This is not the full count for up-to date spaMM (hyper param are missing),
# but should be OK for old objects to which this back compat code applies.
} # ELSE up do date objects should already have a correct p_corrPars
}
unlist(dfs)
}
.get_info_crits <- function(object, also_cAIC=TRUE, nsim=0L, ...) {
if (is.null(info_crits <- object$envir$info_crits) || (also_cAIC && is.null(info_crits[["cAIC"]]))) {
dof <- DoF(object)
p_rdisp <- dof[["p_rdisp"]]
APHLs <- object$APHLs
info_crits <- list()
# poisson-Gamma and negbin should have similar similar mAIC => NB_shape as one df or lambda as one df
forAIC <- APHLs
if (object$models[[1]]=="etaHGLM") {
if (object$HL[1]=="SEM") {
forAIC <- list(p_v=APHLs$logLapp,p_bv=APHLs$logLapp,clik=APHLs$clik)
}
# if standard ML: there is an REMLformula ~ 0 (or with ranefs ?); processed$X.Re is 0-col matrix
# if standard REML: REMLformula is NULL: $X.Re is X.pv, processed$X.Re is NULL
# non standard REML: other REMLformula: $X.Re and processed$X.Re identical, and may take essentially any value
X.pv <- model.matrix(object)
alldfs <- sum(dof)
info_crits$mAIC <- -2*forAIC$p_v + 2 *(alldfs)
info_crits$dAIC <- -2*forAIC$p_bv + 2 * (alldfs- dof[["p_fixef"]]) ## HaLM07 (eq 10) focussed for dispersion params
# including the rho param of an AR model
if (also_cAIC) {
if (.is_spprec_fit(object)) {
pd <- .calc_cAIC_pd_spprec(object)
} else if ( ! is.null(object$envir$sXaug)) {
pd <- .calc_cAIC_pd_from_sXaug(object)
} else { # SEM
H_w.resid <- .get_H_w.resid(object)
ZAL <- .compute_ZAL(XMatrix=object$strucList, ZAlist=object$ZAlist,as_matrix=.eval_as_mat_arg.HLfit(object))
d2hdv2 <- .calcD2hDv2(ZAL,H_w.resid,object$w.ranef) ## update d2hdv2= - t(ZAL) %*% diag(w.resid) %*% ZAL - diag(w.ranef)
pd <- .calc_cAIC_pd_others(X.pv, ZAL, H_w.resid, d2hdv2)
}
info_crits$GoFdf <- length(object$y) - pd ## <- nobs minus # df absorbed in inference of ranefs
## eqs 4,7 in HaLM07
info_crits$cAIC <- -2*forAIC$clik + 2*(pd+p_rdisp) ## no p_lambda ! Presumably if resid disp is declared as individual-level ranef,
# p_rdisp is reduced to 0 but pd is increased by as much so the result is still correct...
}
# print(c(pd,p_phi))
# but Yu and Yau then suggest caicc <- -2*clik + ... where ... involves d2h/db d[disp params] and d2h/d[disp params]2
} else { ## fixed effect model
info_crits$cAIC <- info_crits$mAIC <- -2*forAIC$p_v+2*(dof[["p_fixef"]]+p_rdisp)
# => sets cAIC so that (also_cAIC && is.null(info_crits[["cAIC"]])) becomes FALSE although also_cAIC is true by default...
}
object$envir$info_crits <- info_crits
} else p_rdisp <- NULL
if ( nsim>0L ) {
if (is.null(p_rdisp)) p_rdisp <- .calc_p_rdisp(object)
b_pd <- .calc_boot_AIC_dfs(object, nsim=nsim, ...)
if (object$models[[1]]=="etaHGLM") {
object$envir$info_crits$b_cAIC <- -2*object$APHLs$clik + 2*(b_pd+p_rdisp)
} else object$envir$info_crits$b_cAIC <- -2*object$APHLs$p_v + 2*(b_pd+p_rdisp) # as documented
}
return(object$envir$info_crits)
}
# Two ad hoc functions implying large atrix computations post-fit. In-fit .hatvals2std_lev() has comments about avoiding these matrix computations there. Not immediately transposable here.
.calc_dvdlogphiMat_new <- function(dh0deta,ZAL,
sXaug,d2hdv2_info=NULL ## either one
) {
## cf calcul dhdv, but here we want to keep each d/d phi_i distinct hence not sum over observations i
if (inherits(ZAL,"ZAXlist")) { # appeared in first try gaussian("logit")... with ARp()...
neg.d2h0_dv_dlogphi <- .crossprod(ZAL,.sparseDiagonal(x= drop(dh0deta), shape="g"))
} else neg.d2h0_dv_dlogphi <- .m_Matrix_times_Dvec(t(ZAL), drop(dh0deta)) # n_u_h*nobs: each ith column is a vector of derivatives wrt v_k# dh0dv <- t(ZAL) %*% diag(as.vector(dh0deta))
if (is.null(d2hdv2_info)) { # call by HLfit_body
dvdlogphiMat <- get_from_MME(sXaug,"solve_d2hdv2",B=neg.d2h0_dv_dlogphi)
} else if (inherits(d2hdv2_info,"CHMfactor")) { # CHM of ***-*** d2hdv2
dvdlogphiMat <- solve(d2hdv2_info, - neg.d2h0_dv_dlogphi) # efficient without as.matrix()!
} else if (inherits(d2hdv2_info,"qr") || inherits(d2hdv2_info,"sparseQR") ) {
dvdlogphiMat <- solve(d2hdv2_info, as.matrix(neg.d2h0_dv_dlogphi)) # rXn
} else if (is.environment(d2hdv2_info)) {
# dvdlogphiMat <- d2hdv2_info %*% neg.d2h0_dv_dlogphi # rXn
rhs <- d2hdv2_info$chol_Q %*% neg.d2h0_dv_dlogphi
rhs <- solve(d2hdv2_info$G_CHMfactor,rhs)
dvdlogphiMat <- - .crossprod(d2hdv2_info$chol_Q,rhs) # don't forget '-'
} else { ## then d2hdv2_info is ginv(d2hdv2) or a sparse matrix inverse of (d2hdv2) (spprec code will provide a dsCMatrix)
if (inherits(d2hdv2_info,"dsCMatrix")) {
d2hdv2_info <- as(d2hdv2_info,"dgeMatrix") ## more efficient if inv_d2hdv2 is math-dense
# It would be nice to store only the half matrix but then as( - d2hdv2_info, "dpoMatrix") and reversing sign afterwards.
}
dvdlogphiMat <- d2hdv2_info %*% neg.d2h0_dv_dlogphi # rXn
}
return(dvdlogphiMat)
}
.calc_dvdloglamMat_new <- function(neg.d2f_dv_dloglam, # a .sparseDiagonal or the result of .bdiag()
d2hdv2_info) {
if (inherits(d2hdv2_info,"CHMfactor")) {# CHM of ***-*** d2hdv2
dvdloglamMat <- solve(d2hdv2_info, - neg.d2f_dv_dloglam) # rXr !
} else if (inherits(d2hdv2_info,"qr") || inherits(d2hdv2_info,"sparseQR") ) { ## much slower than using CHMfactor
if (nrow(neg.d2f_dv_dloglam)>5000L) message("[one-time solve()ing of large matrix, which may be slow]")
dvdloglamMat <- solve(d2hdv2_info, neg.d2f_dv_dloglam) # rXr !
} else if (is.environment(d2hdv2_info)) {
# dvdloglamMat <- solve(d2hdv2_info, diag( neg.d2f_dv_dloglam )) # rXr !
rhs <- d2hdv2_info$chol_Q %*% neg.d2f_dv_dloglam
rhs <- solve(d2hdv2_info$G_CHMfactor, rhs)
dvdloglamMat <- - .crossprod(d2hdv2_info$chol_Q, rhs) # don't forget '-'
} else { ## then d2hdv2_info is ginv(d2hdv2) or some other form of inverse
# This case happens at least in my ressp devel example (tiny ranCoefs spprec)
if (inherits(d2hdv2_info,"dsCMatrix")) d2hdv2_info <- as(d2hdv2_info, "generalMatrix")
dvdloglamMat <- d2hdv2_info %*% neg.d2f_dv_dloglam
}
# I returned as.matrix(dvdloglamMat) a long time ago, and found it terribly inefficient in bigranefs case;
# but bigranefs example no longer runs calls .calc_dvdloglamMat_new()
return(dvdloglamMat)
} ## square matrix, by the formulation of the algo
..calc_d2hdv2_info <- function(object, ZAL) {
envir <- object$envir
if (is.null(factor_inv_Md2hdv2 <- envir$factor_inv_Md2hdv2)) { ## code not using promises
# we reach here at least in simple dense-correlation random-slope fit
# (where d2hdv2 may be sparse, and resulting d2hdv2_info too)
if ( ! is.null(envir$G_CHMfactor)) {
if (length(object$y) > ncol(envir$chol_Q)) { # We precompute inverse(d2hdv2) so that .calc_dvdlogphiMat_new() has ONE costly (r*r) %*% (r*n).
if ( is.null(envir$invL_G.P)) { # allowing for old objects that did not have this element
envir$factor_inv_Md2hdv2 <- Matrix::solve(envir$G_CHMfactor, as(envir$G_CHMfactor,"pMatrix") %*% envir$chol_Q,system="L") ## dgCMatrix
} else envir$factor_inv_Md2hdv2 <- envir$invL_G.P %*% envir$chol_Q ## dgCMatrix
# envir$factor_inv_Md2hdv2 <- Matrix::drop0(LL,tol=.Machine$double.eps)
d2hdv2_info <- - .crossprod(envir$factor_inv_Md2hdv2)
} else d2hdv2_info <- envir # then .calc_dvdlogphiMat_new() has TWO (r*r) %*% (r*n) (plus an efficient solve())
} else {
d2hdv2 <- .calcD2hDv2(ZAL,.get_H_w.resid(object),object$w.ranef)
if (inherits(d2hdv2,"sparseMatrix")) {
d2hdv2_info <- .silent_W_E(Cholesky( - d2hdv2,LDL=FALSE,perm=TRUE )) # '-' !
# this d2hdv2_info seems used only in solve(, default system="A") so perm is automatically handled.
if (inherits(d2hdv2_info, "CHMfactor")) {
#d2hdv2_info <- structure(d2hdv2_info, BLOB=list2env(list(), parent=environment(.solve_CHM)))
rank <- ncol(d2hdv2)
} else {
# oldMDCopt <- options(Matrix.warnDeprecatedCoerce = 0) # qr(<dsC>) problem in Matrix v1.4.2
d2hdv2_info <- qr(d2hdv2, tol=spaMM.getOption("qrTolerance")) # sparseQR
# options(oldMDCopt)
rank <- sum(abs(diag(qrR(d2hdv2_info,backPermute=FALSE)))>1e-7)
}
} else { # d2hdv2 was a matrix or a Matrix in dense format (dgeMatrix...)
d2hdv2_info <- qr(d2hdv2, tol=spaMM.getOption("qrTolerance")) # qr
rank <- d2hdv2_info$rank
}
if (rank<ncol(d2hdv2)) {
d2hdv2_info <- .force_solve(d2hdv2) # inverse(d2hdv2)
} # else we keep the QR facto.
}
} else {
# we reach here at least in simple ssprec random-slope fit ... so dense result here for spprec
d2hdv2_info <- - .crossprodCpp_d(as.matrix(factor_inv_Md2hdv2), yy=NULL) # inverse(d2hdv2)
# .crossprodCpp_d much faster than Matrix::crossprod(dsCMatrix...);
# still so for get_predVar(fit.Frailty) 2025/01/12;
# and still faster than any of the as_mat variants, eg .Rcpp_crossprod(factor_inv_Md2hdv2, BB = NULL, as_mat=TRUE)
# bc .crossprodCpp_d(dense) is fast relative to .Rcpp_crossprod(sparse), irrespective of as_mat
}
return(d2hdv2_info)
}
.calc_d2hdv2_info <- function(object, ZAL) {
if (TRUE) return(..calc_d2hdv2_info(object, ZAL)) # rediscovered long afterwards... so ./.
## The code below first checks for a "factor_inv_Md2hdv2" promise. But
## there is no such promise in the current object bc HLfit_body() has run
## .init_promises_spprec(sXaug, non_X_ones=FALSE, nullify_X_ones =TRUE)
## and non_X_ones=FALSE implies that there is no "factor_inv_Md2hdv2" promise.
#
# ELSE older version
envir <- object$envir
if ( ! is.null(envir$G_CHMfactor)) { ## spprec code
if (! .is_evaluated("factor_inv_Md2hdv2", envir)) {
if (length(object$y) > ncol(envir$chol_Q)) { # We precompute inverse(d2hdv2) so that .calc_dvdlogphiMat_new() has ONE costly (r*r) %*% (r*n).
d2hdv2_info <- - .crossprod(envir$factor_inv_Md2hdv2)
} else d2hdv2_info <- envir # then .calc_dvdlogphiMat_new() has TWO (r*r) %*% (r*n) (plus an efficient solve())
} else d2hdv2_info <- - .crossprodCpp_d(as.matrix(envir$factor_inv_Md2hdv2), yy=NULL)
} else {
if (is.null(factor_inv_Md2hdv2 <- envir$factor_inv_Md2hdv2)) {
d2hdv2 <- .calcD2hDv2(ZAL,object$w.resid,object$w.ranef)
if (inherits(d2hdv2,"sparseMatrix")) {
d2hdv2_info <- .silent_W_E(Cholesky( - d2hdv2,LDL=FALSE,perm=TRUE )) # '-' !
if (inherits(d2hdv2_info, "CHMfactor")) {
#d2hdv2_info <- structure(d2hdv2_info, BLOB=list2env(list(), parent=environment(.solve_CHM)))
rank <- ncol(d2hdv2)
} else {
d2hdv2_info <- qr(d2hdv2, tol=spaMM.getOption("qrTolerance")) # sparseQR
rank <- sum(abs(diag(qrR(d2hdv2_info,backPermute=FALSE)))>1e-7)
}
} else { # d2hdv2 was a matrix or a Matrix in dense format (dgeMatrix...)
d2hdv2_info <- qr(d2hdv2, tol=spaMM.getOption("qrTolerance")) # qr
rank <- d2hdv2_info$rank
}
if (rank<ncol(d2hdv2)) {
d2hdv2_info <- .force_solve(d2hdv2) # inverse(d2hdv2)
} # else we keep the QR facto.
} else d2hdv2_info <- - .crossprodCpp_d(as.matrix(envir$factor_inv_Md2hdv2), yy=NULL) # inverse(d2hdv2) # much faster than Matrix::crossprod(dsCMatrix...)
# and still faster than any of the as_mat variants, eg .Rcpp_crossprod(factor_inv_Md2hdv2, BB = NULL, as_mat=TRUE)
#
}
return(d2hdv2_info)
}
.calc_logdispObject <- function(object, envir=object$envir, force_fixed,
# promises:
strucList=object$strucList,
isRandomSlope= attr(strucList,"isRandomSlope")) {
dvdloglamMat <- envir$dvdloglamMat
dvdloglamMat_needed <- ( is.null(dvdloglamMat) &&
# The fact that next line is commented allows computation of predVar effect for random slope lambdas):
# all(unlist(attr(object$ZAlist,"namesTerms"))=="(Intercept)") && ## (1|.) or CAR or Matern
(force_fixed ||
any( ! object$lambda.object$type %in% c("fixed","fix_ranCoefs","fix_hyper"))) ) ## some lambda params were estimated
dvdlogphiMat <- envir$dvdlogphiMat
dvdlogphiMat_needed <- (is.null(dvdlogphiMat) &&
( any((phimodel <- object$models[["phi"]])=="phiScal") ||
force_fixed))
dvdlogphiMat_needed <- dvdlogphiMat_needed || identical(envir$forcePhiComponent,TRUE) ## hack for code testing !
if (dvdloglamMat_needed || dvdlogphiMat_needed) {
ZAL <- get_ZALMatrix(object, force_bind = ! (.is_spprec_fit(object)) )
d2hdv2_info <- .calc_d2hdv2_info(object, ZAL) # may be a qr object, or not (SPPREC). F I X M E a gentle message for long computations ?
}
if (dvdloglamMat_needed) { # $\partial_\tau \bv$ in latex doc.
cum_n_u_h <- attr(.get_u_h(object),"cum_n_u_h")
psi_M <- rep(attr(object$rand.families,"unique.psi_M"),diff(cum_n_u_h))
u_h <- .get_u_h(object)
dlogfthdth <- (psi_M - u_h)/object$lambda.object$lambda_est ## the d log density of th(u)
neg.d2f_dv_dloglam <- .calc_neg_d2f_dv_dloglam(dlogfthdth, cum_n_u_h,
lcrandfamfam=attr(object$rand.families,"lcrandfamfam"),
rand.families=object$rand.families, u_h=u_h)
## I tried using the private .fix_RHS_dvdloglamMat() here instead of the following simple conversion.
## See further comments there.
neg.d2f_dv_dloglam <- .sparseDiagonal(x= neg.d2f_dv_dloglam, shape="g")
dvdloglamMat <- .calc_dvdloglamMat_new(neg.d2f_dv_dloglam,
d2hdv2_info=d2hdv2_info) ## d2hdv2_info is either a qr factor or the inverse as a matrix or an environment
# envir$dvdloglamMat <- dvdloglamMat # needed only for debugging
}
if (dvdlogphiMat_needed) {
muetablob <- object$muetablob
# .get_H_w.resid() rather than .get_w.resid here. See section mentioning "dvdloglamMat" in the long doc.
if ( ! is.null(envir$G_CHMfactor)) { # spprec; possibly generalisable code not using math-dense ZAL
# rhs <- .Matrix_times_Dvec(t(envir$sXaug$AUGI0_ZX$ZAfix), - dh0deta) # efficient
# rhs <- solve(envir$G_CHMfactor,rhs,system="A") # efficient
unW_dh0deta <- (object$y-muetablob$mu)/muetablob$dmudeta ## (soit Bin -> phi fixe=1, soit BinomialDen=1)
if (is.null(envir$invG_ZtW)) {
if (is.null(envir$ZtW)) {
# cf comments in .old_calc_Md2hdvb2_info_spprec_by_r22()
envir$ZtW <- t(.Dvec_times_m_Matrix(.get_H_w.resid(envir=envir), envir$sXaug$AUGI0_ZX$ZAfix))
}
envir$invG_ZtW <- solve(envir$G_CHMfactor,
envir$ZtW, system="A") # hardly avoidable has there is no comparable operation elsewhere (check "A")
}
rhs <- .Matrix_times_Dvec(envir$invG_ZtW, -unW_dh0deta) # efficient
dvdlogphiMat <- .crossprod(envir$chol_Q, rhs) # _FIXME_ bottleneck in large spprec but .crossprodCpp_d not useful here
} else {
dh0deta <- ( .get_H_w.resid(object) *(object$y-muetablob$mu)/muetablob$dmudeta ) ## (soit Bin -> phi fixe=1, soit BinomialDen=1)
dvdlogphiMat <- .calc_dvdlogphiMat_new(dh0deta=dh0deta, ZAL=ZAL,
d2hdv2_info=d2hdv2_info ## either a qr factor or a matrix inverse or envir
)
}
# envir$dvdlogphiMat <- dvdlogphiMat # needed only for debugging
}
invV_factors <- .calc_invV_factors(object, r_x_r_needed = any(isRandomSlope)) ## n_x_r and r_x_n in repres of invV as diag(w.resid)- [n_x_r %*% r_x_n]
if (length(phimodel)>1L) {
.calc_logdisp_cov_mv(object, dvdloglamMat=dvdloglamMat, ## square matrix, by the formulation of the algo
dvdlogphiMat=dvdlogphiMat, invV_factors=invV_factors,
force_fixed=force_fixed)
} else
.calc_logdisp_cov(object, dvdloglamMat=dvdloglamMat, ## square matrix, by the formulation of the algo
dvdlogphiMat=dvdlogphiMat, invV_factors=invV_factors,
force_fixed=force_fixed)
}
.get_logdispObject <- function(object) { ##
envir <- object$envir
if (is.null(envir$logdispObject) && object$models[["eta"]]=="etaHGLM" ) {
envir$logdispObject <- .calc_logdispObject(object, envir=envir, force_fixed=FALSE)
}
return(envir$logdispObject)
} # if dvdloglamMat or dvdlogphiMat were computed ex-tempo, they are NOT saved.
## This provides factor n_x_r and r_x_n of the representation of invV as diag(w.resid)- [n_x_r %*% r_x_n = t(Ztw) %*% invG.ZtW] (nXr %*% rxn)
## slow computation the one time .get_logdispObject() is called, for variances$disp (no need to store the result in an $envir)
.calc_invV_factors <- function(object,
# $r_x_r returned for spprec or if 'r_x_r_needed' (set up for ranCoefs):
r_x_r_needed=FALSE) { ## used by .get_logdispObject
## Store inv( G=[{precmat=inv(L invWranef Lt)} +ZtWZ] ) as two matrix nXr and rXn rather than their nXn product
# F I X M E yet there will be cases where n<r and then it's better to store the n x n product !
if (.is_spprec_fit(object)) {
envir <- object$envir
## code clearly related to .Sigsolve_sparsePrecision() algorithm:
if (is.null(envir$invG_ZtW)) {
if (is.null(envir$ZtW)) {
# cf comments in .old_calc_Md2hdvb2_info_spprec_by_r22()
H_w.resid <- .get_H_w.resid(envir=envir)
envir$ZtW <- t(.Dvec_times_m_Matrix(H_w.resid, envir$sXaug$AUGI0_ZX$ZAfix))
}
object$envir$invG_ZtW <- solve(envir$G_CHMfactor, envir$ZtW, system="A") # hardly avoidable has there is no comparable operation elsewhere (check "A")
}
RES <- list(n_x_r=t(envir$ZtW), r_x_n=as.matrix(envir$invG_ZtW)) # requires both being kept in the envir
ZAfix <- .get_ZAfix(object, as_matrix=FALSE)
RES$r_x_r <- object$envir$invG_ZtW %*% ZAfix
return(RES)
} else {
ZAfix <- .get_ZAfix(object, as_matrix=FALSE)
H_w.resid <- .get_H_w.resid(object)
if (.is_identity(ZAfix)) {
# FIXME inelegant ZAL computation only to test if it is diagonal
ZAL <- .compute_ZAL(XMatrix=object$strucList, ZAlist=object$ZAlist,as_matrix=.eval_as_mat_arg.HLfit(object))
if (inherits(ZAL,"diagonalMatrix")) {
RES <- list(n_x_r=Diagonal(x = H_w.resid),
r_x_n=Diagonal(x = H_w.resid/(H_w.resid + object$w.ranef/diag(ZAL)^2)))
return(RES)
} ## ELSE
} ## ELSE
#
invL <- .get_invL_HLfit(object) # t(tcrossfac(precision_matrix))
wrZ <- .Dvec_times_m_Matrix(H_w.resid, ZAfix) # suppressMessages(sweep(t(ZA), 2L, w.resid,`*`))
if (TRUE) {
ZtwrZ <- .crossprod(ZAfix, wrZ, as_sym=TRUE) ## seems more precise and we must compute wrZ anyway
} else ZtwrZ <- .safe_ZtWZwrapper(ZAfix,.get_H_w.resid(object)) ## -> calls .crossprod(., y=NULL) affects numerical precision of twolambda test in test-predVar.R
## is general ZtwrZ should be sparse, while invL may not. Large dense invL will away be a problem
## for small Z, ZtwrZ may bedsy, but this an un-intersting subcase that does not call for a special treatment
if (inherits(ZtwrZ,"dsCMatrix") || inherits(ZtwrZ,"dsyMatrix")) {
# (NB: dsy+dsy faster than dsC+dsy but maybe just because of conversion from dsC to dsy ?) => no obvious improvement
if (.is_identity(invL)) {
precmat <- .symDiagonal(x=object$w.ranef) # dsC+dsC
} else precmat <- .safe_ZtWZwrapper(invL,object$w.ranef) # dsC+whatever
} else if (inherits(ZtwrZ,"ddiMatrix")) {
if (.is_identity(invL)) {
precmat <- Diagonal(x=object$w.ranef) # ddi+ddi...
} else precmat <- .safe_ZtWZwrapper(invL,object$w.ranef) # ddi+whatever
} else { # ("matrix") may never occur
if (.is_identity(invL)) {
precmat <- diag(x=object$w.ranef)
} else precmat <- .safe_ZtWZwrapper(invL,object$w.ranef)
message("Possibly inefficient code in .calc_invV_factors().")
}
## try to sum dsC or to sum dense matrix but not to mix types..., and to avoid formation of a large nxn matrix:
if (inherits(precmat,"dsCMatrix") && inherits(ZtwrZ,"dsCMatrix")) {
Gmat <- .dsCsum(precmat,ZtwrZ)
} else Gmat <- precmat + ZtwrZ
invG_ZtW <- tryCatch(solve(Gmat, t(wrZ)),error=function(e) e)
if (inherits(invG_ZtW,"simpleError")) { ## but that should be well behaved when precmat is.
invG <- .ginv(as.matrix(Gmat)) ## FIXME quick patch at least
invG_ZtW <- .tcrossprod(invG, wrZ)
}
RES <- list(n_x_r=wrZ, r_x_n=invG_ZtW)
if (r_x_r_needed) RES$r_x_r <- invG_ZtW %*% ZAfix
return(RES)
}
}
.calc_beta_cov_info_others <- function(wAugX=NULL, AUGI0_ZX, ZAL, ww) { ## post-fit fn
if (is.null(wAugX)) {
if (is.null(ZAL)) { # GLM... or LLM
if (any(ww<0)) { # ... LLM
negHess <- crossprod(AUGI0_ZX$X.pv, .Dvec_times_m_Matrix( ww, AUGI0_ZX$X.pv))
cholH <- tryCatch(chol(negHess),error=function(e) e)
if (inherits(cholH, "simpleError")) { # => brute regularization
warning("logLik presumably not maximized (information matrix is not positive definite at attained estimates).", immediate. = TRUE)
beta_cov <- tcrossfac_beta_v_cov <- matrix(NA, ncol=ncol(negHess), nrow=ncol(negHess))
} else {
beta_cov <- chol2inv(cholH)
tcrossfac_beta_v_cov <- t(cholH)
}
colnames(beta_cov) <- rownames(beta_cov) <- colnames(negHess) # sigh
return(list(beta_cov=beta_cov, tcrossfac_beta_v_cov=tcrossfac_beta_v_cov))
} else wAugX <- .calc_wAugX(XZ_0I=AUGI0_ZX$X.pv,sqrt.ww=sqrt(ww))
# => ... may have zero cols... X.pv being a zero-col *m*atrix, in which case .spaMM.data$options$matrix_method is selected below
# and then get_from_MME(<matrix_method>, "beta_cov_info_from_wAugX") must handle this case.
} else {
XZ_0I <- .calc_XZ_0I(AUGI0_ZX=AUGI0_ZX,ZAL) # ZAL is not ZAXlist since .calc_beta_cov_info_others not called in spprec
wAugX <- .calc_wAugX(XZ_0I=XZ_0I,sqrt.ww=sqrt(ww))
}
} ## wAugX is in XZ_OI order
if (inherits(wAugX,"Matrix")) {
sXaug_method <- .spaMM.data$options$Matrix_method
} else sXaug_method <- .spaMM.data$options$matrix_method
.silent_M_E(untrace(sXaug_method, where=asNamespace("spaMM"))) # try() bc this fails when called by Infusion
# test: Infusion tests -> ... -> .predict_body -> ;.. -> .calc_beta_cov_info_others -> untrace -> def_sXaug_EigenDense_QRP_Chol_scaled not found
# Note that the function is in exportPattern, explicitly exporting/impporting it does not help; attaching spaMM seems required to avoid untrace's error..
# hack to recycle sXaug code; all weights are 1 or unit vectors as the order is not that assumed by sXaug_method.
wAugX <- do.call(sXaug_method,list(Xaug=wAugX, weight_X=rep(1,nrow(AUGI0_ZX$X.pv)),
w.ranef=rep(1,ncol(AUGI0_ZX$I)), ## we need at least its length for get_from Matrix methods
H_global_scale=1))
attr(wAugX,"AUGI0_ZX") <- AUGI0_ZX # cf use of trDiag in .sXaug_Matrix_QRP_CHM_scaled -> solve_R_scaled
beta_cov_info <- get_from_MME(wAugX,"beta_cov_info_from_wAugX")
return(beta_cov_info)
}
.calc_beta_cov_info_from_sXaug <- function(BLOB, sXaug, tcrossfac) {
pforpv <- attr(sXaug,"pforpv")
X_scaling <- sqrt(rep(attr(sXaug,"H_global_scale"),pforpv))
X_scale <- attr(sXaug,"scaled:scale") # this is (! spprec) code and the X.pv attribute has been copied here
if ( ! is.null(X_scale)) X_scaling <- X_scaling/X_scale
diagscalings <- c(1/sqrt(attr(sXaug,"w.ranef")), X_scaling)
if ( ! is.null(BLOB$sortPerm)) { # depending on method used for QR facto
sctPmat <- sparseMatrix(seq_along(BLOB$sortPerm), BLOB$sortPerm, x=1)
sctPmat <- .Dvec_times_Matrix(diagscalings, sctPmat) # Pmat <- .Matrix_times_Dvec(Pmat,diagscalings)
tcrossfac_v_beta_cov <- sctPmat %*% tcrossfac
} else tcrossfac_v_beta_cov <- .Dvec_times_m_Matrix(diagscalings, tcrossfac) ## loses colnames...
dgC_good <- (inherits(tcrossfac_v_beta_cov, "dgCMatrix") &&
.calc_denseness(tcrossfac_v_beta_cov, relative=TRUE)<0.35 )
if ( ! dgC_good) tcrossfac_v_beta_cov <- as.matrix(tcrossfac_v_beta_cov) # bigranefs.R shows that conversion is not always good.
rownames(tcrossfac_v_beta_cov) <- colnames(sXaug) ## necessary for summary.HLfit, already lost in BLOB$R_scaled
seqp <- seq_len(pforpv)
beta_pos <- attr(sXaug,"n_u_h")+seqp
beta_v_order <- c(beta_pos,seq(attr(sXaug,"n_u_h")))
tcrossfac_beta_v_cov <- tcrossfac_v_beta_cov[beta_v_order,,drop=FALSE]
if (inherits(sXaug,"dtCMatrix")) tcrossfac_beta_v_cov <- as(tcrossfac_beta_v_cov, "sparseMatrix")
beta_cov <- as.matrix(.tcrossprod(tcrossfac_beta_v_cov[seqp,,drop=FALSE])) ## assignment in .make_beta_table() assumes a dense matrix
return( list(beta_cov=beta_cov,
#beta_v_cov=beta_v_cov,
tcrossfac_beta_v_cov=tcrossfac_beta_v_cov) )
}
.get_beta_cov_info <- function(res) {
# Provide list(beta_cov=., tcrossfac_beta_v_cov=.)
if (.is_spprec_fit(res)) {
if (is.null(dcdb_p4m <- res$envir$sXaug$AUGI0_ZX$dcdb_p4m)) dcdb_p4m <- res$X.pv # p4m spprec has $dcdb_p4m
return(.calc_beta_cov_info_spprec(X.pv=dcdb_p4m, envir=res$envir))
} else if (! is.null(res$envir$sXaug)) { # excludes SEM *and* fixed-effect models
if (prod(dim(res$envir$sXaug))>1e7) message("[one-time computation of covariance matrix, which may be slow]")
if (res$spaMM.version<"2.7.34") attr(res$envir$sXaug,"scaled:scale") <- attr(res$X.pv,"scaled:scale")
res$envir$beta_cov_info <- get_from_MME(res$envir$sXaug,which="beta_cov_info_from_sXaug")
} else { ## older code, generic for non-spprec.
# What this [block -> .calc_beta_cov_info_others()] provides is a beta_cov + full beta_v_cov
if (is.matrix(beta_cov_info <- res$envir$beta_cov_info) || ## old format (old fit object), should be a list now
is.null(beta_cov_info$tcrossfac_beta_v_cov)) { ## Use .calc_beta_cov_info_others -> get_from_MME()
# F I X M E test useless if call to .get_beta_cov_info() already dependent on the test
# Further, can $tcrossfac_beta_v_cov be NULL if $beta_cov_info is list ?
ZAL <- get_ZALMatrix(res, force_bind=FALSE) # note that .calc_beta_cov_info_others -> ... -> rbind2() but ZAL should not be a ZAXlist here
n_u_h <- length(res$w.ranef) # may be 0: this code may be run for LMs.
pforpv <- ncol(res$X.pv)
if (inherits(ZAL,"Matrix")) {
AUGI0_ZX <- list(I=.sparseDiagonal(n=n_u_h, shape="g"),
trDiag=..trDiagonal(n=n_u_h+pforpv), # trDiag needed by solve_R_scaled promise (from summary(<probitgem fit>)).
ZeroBlock=Matrix(0,nrow=n_u_h,ncol=pforpv),X.pv=res$X.pv)
} else { # ZAL is dense or is NULL (what if itis dense and X.pv is sparse?)
AUGI0_ZX <- list(I=diag(nrow=n_u_h),ZeroBlock=matrix(0,nrow=n_u_h,ncol=pforpv),X.pv=res$X.pv)
if (inherits(res$X.pv,"Matrix")) AUGI0_ZX$trDiag <- ..trDiagonal(n=n_u_h+pforpv) # useful for LM with sparse X,
# cf use of trDiag in .sXaug_Matrix_QRP_CHM_scaled -> solve_R_scaled
}
res$envir$beta_cov_info <- .calc_beta_cov_info_others(AUGI0_ZX=AUGI0_ZX,ZAL=ZAL,ww=c(.get_H_w.resid(res),res$w.ranef)) ## with beta_v_cov attribute
}
}
return(res$envir$beta_cov_info)
}
.calc_newZACvar <- function(newZAlist,cov_newLv_oldv_list) {
newZACvarlist <- vector("list",length(newZAlist))
for (new_rd in seq_along(newZAlist)) newZACvarlist[[new_rd]] <- newZAlist[[new_rd]] %ZA*gI% cov_newLv_oldv_list[[new_rd]]
newZAC <- .ad_hoc_cbind(newZACvarlist, as_matrix=FALSE)
if (inherits(newZAC,"sparseMatrix") && .calc_denseness(newZAC, relative=TRUE)>0.35) newZAC <- as.matrix(newZAC)
newZAC
}
## Aggregate info on corrpars, inner-estimated and inner-fixed.
## $corrPars is only for info in messages() and return value, (?!)
.make_CorrEst_and_RanFix <- function(ranFix, ## has "fix", "outer", and also "var" values ! code corrently assumes "var" <=> corr_est
corr_est
) {
# When fitting function was HLCor, the type attribute has not been added to the corrPars
if (is.null(attr(ranFix,"type"))) {
rFtype <- .relist_rep("fix", ranFix)
# HLCor may have been called through fitme/fitmv with inner estimation of lambda...
# Test case is first fitmv_body fit for BbyP_sp: initvec is NULL, inner estim of lambda
# get_fittedPars -> ... -> .get_ranPars_notPhi() finally gets the correct info from lambda.object$type.
# but get_ranPars gets info from $corrEst_and_ranFix which gets info from here
# This means we need to distinguish inner from fix here.
rFtype["lambda"] <- .modify_list(rFtype["lambda"],
relist("inner",ranFix["lambda"])[is.na(ranFix["lambda"])])
attr(ranFix,"type") <- rFtype
}
if ( ! is.null(corr_est)) {
ranFix <- structure(.modify_list(ranFix,corr_est),
type=.modify_list(attr(ranFix,"type"),.relist_rep("var",corr_est)))
}
return(ranFix) ## correlation params + whatever else was in ranFix
}
# Has become a postfit fn; to recycle processed$u_h_v_h_from_v_h one would have to remove processed from its defining envir.
.u_h_v_h_from_v_h <- function(v_h, rand.families, cum_n_u_h, lower.v_h, upper.v_h,
u_list=vector("list", length(rand.families))) {
if(!is.null(lower.v_h)) {v_h[v_h<lower.v_h] <- lower.v_h}
if(!is.null(upper.v_h)) {v_h[v_h>upper.v_h] <- upper.v_h}
nrand <- length(rand.families)
for (it in seq_len(nrand)) {
u.range <- (cum_n_u_h[it]+1L):(cum_n_u_h[it+1L])
u_list[[it]] <- rand.families[[it]]$linkinv(v_h[u.range])
if (any(is.infinite(u_list[[it]]))) {
warning("infinite random values ('u_h') were constrained to finite range.")
u_list[[it]] <- pmin(.Machine$double.xmax, pmax(-.Machine$double.xmax,u_list[[it]]) )
}
}
u_h <- .unlist(u_list)
## if there were box constr, v_h may have been modified, we put it in return value
if ( ! (is.null(lower.v_h) && is.null(upper.v_h))) attr(u_h,"v_h") <- v_h
return(u_h)
}
Try the spaMM package in your browser
Any scripts or data that you put into this service are public.
spaMM documentation built on Sept. 10, 2026, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .u_h_v_h_from_v_h .make_CorrEst_and_RanFix .calc_newZACvar .get_beta_cov_info .calc_beta_cov_info_from_sXaug .calc_beta_cov_info_others .calc_invV_factors .get_logdispObject .calc_logdispObject .calc_d2hdv2_info ..calc_d2hdv2_info .calc_dvdloglamMat_new .calc_dvdlogphiMat_new .get_info_crits DoF .calc_p_rdisp .calc_cAIC_pd_others .calc_cAIC_pd_from_sXaug .calc_cAIC_pd_spprec .calc_pd_product
Documented in DoF
.calc_pd_product <- function(tcrossfac_v_beta_cov, Md2clikdvb2, blockSize=5000L) {
if ((nc <- ncol(tcrossfac_v_beta_cov))>(blockSize)) {
## We reached this point by sparse matrix computations. We need to save memory in the following dense computations,
message(paste0("Conditional-AIC computation requires operations on a large matrix (square with dimension ",nc,"),\n",
"which may take a bit of time. Use 'also_cAIC=FALSE' to avoid it."))
slices <- unique(c(seq(0L,nc,blockSize),nc))
nslices <- length(slices)-1L
it <- 0L ## 'global definition' for Rcmd check
foreach_args <- list(it = seq_len(nslices), .combine = "sum")
foreach_blob <- do.call(foreach::foreach,foreach_args)
abyss <- foreach::`%do%`(foreach_blob, Sys.setenv(LANG = "en"))
barstyle <- eval(spaMM.getOption("barstyle"))
progrbar_setup <- .set_progrbar(max = nslices, style = barstyle, char="s")
pd <- foreach::`%do%`(foreach_blob, {
slice <- (slices[it]+1L):slices[it+1L]
tmp <- t(.crossprod(tcrossfac_v_beta_cov, Md2clikdvb2[,slice]))
if (barstyle) progrbar_setup$progress(it)
return(sum(tcrossfac_v_beta_cov[slice,] * tmp))
})
if (barstyle) close(progrbar_setup$pb)
# We could parallelize using %dopar% (twice) but is that worth the overhead?
} else {
# logic of following code is
# pd = sum(diag(solve(Md2hdbv2,Md2clikdbv2[c(113:114,1:112),])))
# = sum(diag((tcrossprod(R_invMd2hdvb2)[c(113:114,1:112),] %*% Md2clikdbv2[sort.list(c(113:114,1:112)),])))
# = sum(diag(R_invMd2hdvb2[c(113:114,1:112),] %*% (crossprod(R_invMd2hdvb2, Md2clikdbv2[sort.list(c(113:114,1:112)),]))))
# = sum((R_invMd2hdvb2[c(113:114,1:112),] * t(crossprod(R_invMd2hdvb2, Md2clikdbv2[sort.list(c(113:114,1:112)),]))))
# but we directly use v,b order rather than b,v
pd <- t(.crossprod(tcrossfac_v_beta_cov,Md2clikdvb2))
pd <- sum(tcrossfac_v_beta_cov * pd)
# that is actually the logic of .traceAB:
# pd = .traceAB(tcrossfac_v_beta_cov, t(tcrossfac_v_beta_cov),Md2clikdvb2,diag(nrow =ncol(Md2clikdvb2)))
# but the matrices are square so no use here
}
return(pd)
}
.calc_cAIC_pd_spprec <- function(object) {
H_w.resid <- .get_H_w.resid(object)
#ZAL <- .compute_ZAL(XMatrix=object$strucList, ZAlist=object$ZAlist,as_matrix=.eval_as_mat_arg.HLfit(object))
ZAL <- get_ZALMatrix(object, force_bind=FALSE) # allows ZAXlist; but if a non-ZAXlist is already in the $envir, no effect; + we will solve(chol_Q, Diagonal()) so the gain is not obvious
if ( ncol(X_ori <- model.matrix(object)) ) {
M12 <- .crossprod(ZAL, .Dvec_times_m_Matrix(H_w.resid, X_ori), as_mat=TRUE)
Md2clikdvb2 <- rbind2(cbind2(as.matrix(.safe_ZtWZwrapper(ZAL,H_w.resid)), M12), ## this .ZtWZwrapper() takes time
cbind2(t(M12), as.matrix(.safe_ZtWZwrapper(X_ori,H_w.resid))))
# _FIXME_ any way to avoid formation of this matrix ? Or otherwise message() ?
} else {
Md2clikdvb2 <- as.matrix(.safe_ZtWZwrapper(ZAL,H_w.resid))
}
tcrossfac_v_beta_cov <- .calc_Md2hdvb2_info_spprec(X.pv=X_ori, envir=object$envir,
which="tcrossfac_v_beta_cov")
# not triang if we used sparse QR. Following code should not assume triangularity
pd <- .calc_pd_product(tcrossfac_v_beta_cov, Md2clikdvb2)
return(pd)
}
# From Matrix 1.6-2, "Methods for generic functions rbind2, cbind2, %*%, %&%, crossprod, and tcrossprod determine the class of the result using more strict rules, designed to avoid “surprising” coercions where possible. Notably, operations involving RsparseMatrix now return an RsparseMatrix in more cases. TsparseMatrix and diagonalMatrix may be handled as CsparseMatrix or as RsparseMatrix, depending on context. "
# Maybe some old as.matrix() coertions were previously useful.
# But now I can remove some of then with great gains.
.calc_cAIC_pd_from_sXaug <- function(object) {
if (is.matrix(beta_cov_info <- object$envir$beta_cov_info) || ## matrix is old format, should be a list now
is.null(tcrossfac_beta_v_cov <- beta_cov_info$tcrossfac_beta_v_cov)) {
tcrossfac_beta_v_cov <- .get_beta_cov_info(object)$tcrossfac_beta_v_cov
}
X_ori <- model.matrix(object)
pforpv <- ncol(X_ori)
nc <- ncol(tcrossfac_beta_v_cov)
n_u_h <- nc - pforpv
seqp <- seq_len(pforpv)
perm <- c(pforpv+seq_len(n_u_h), seqp)
tcrossfac_v_beta_cov <- tcrossfac_beta_v_cov[perm,,drop=FALSE] # useful to keep it for predVar computations?
H_w.resid <- .get_H_w.resid(object)
#ZAL <- .compute_ZAL(XMatrix=object$strucList, ZAlist=object$ZAlist,as_matrix=.eval_as_mat_arg.HLfit(object))
ZAL <- get_ZALMatrix(object, force_bind=FALSE) # allows ZAXlist; but if a non-ZAXlist is already in the $envir, no effect; + we will solve(chol_Q, Diagonal()) so the gain is not obvious
if ( ncol(X_ori) ) {
M12 <- .crossprod(ZAL, .Dvec_times_m_Matrix(H_w.resid, X_ori), as_mat=TRUE)
if (.spaMM.data$options$Matrix_old) {
Md2clikdvb2 <- rbind2(cbind2(as.matrix(.safe_ZtWZwrapper(ZAL,H_w.resid)), M12),
cbind2(t(M12), as.matrix(.safe_ZtWZwrapper(X_ori,H_w.resid))))
message("If this step was slow, updating the Matrix package to a version >=1.6-2 might be helpful.")
} else { # AIC(bigranefs) showed the gain of removing the as.matrix()s.
Md2clikdvb2 <- rbind2(cbind2(.safe_ZtWZwrapper(ZAL,H_w.resid), M12),
cbind2(t(M12), .safe_ZtWZwrapper(X_ori,H_w.resid)))
}
# Old comments:
# .ZtWZwrapper(ZAL...) may take time
} else {
Md2clikdvb2 <- as.matrix(.safe_ZtWZwrapper(ZAL,H_w.resid))
}
# not triang if we used sparse QR. Following code should not assume triangularity
pd <- .calc_pd_product(tcrossfac_v_beta_cov, Md2clikdvb2)
return(pd)
}
# Here, using qr.Md2hdbv2$pivot implies that qr() was run on a base matrix.
# So as.matrix() may remain necessary unless further changes are made.
.calc_cAIC_pd_others <- function(X.pv, ZAL, w.resid, d2hdv2, blockSize=1000L) {
if ( ncol(X.pv) ) { ## the projection matrix for the response always includes X even for REML!
hessnondiag <- .crossprod(ZAL, .Dvec_times_m_Matrix(w.resid, X.pv))
Md2hdbv2 <- as.matrix(rbind2(cbind2(.safe_ZtWZwrapper(X.pv,w.resid), t(hessnondiag)),
cbind2(hessnondiag, - d2hdv2)))
Md2clikdbv2 <- as.matrix(rbind2(cbind2(.safe_ZtWZwrapper(X.pv,w.resid), t(hessnondiag)),
cbind2(hessnondiag, .safe_ZtWZwrapper(ZAL,w.resid))))
} else {
Md2hdbv2 <- - d2hdv2
Md2clikdbv2 <- as.matrix(.safe_ZtWZwrapper(ZAL,w.resid))
}
if (inherits(Md2hdbv2,"diagonalMatrix")) {
pd <- sum(diag(Md2clikdbv2)/diag(Md2hdbv2)) ## is sum(diag(solve(Md2hdbv2,Md2clikdbv2)))
} else {
## dans un LMM avec estimation ML, pd = sum(lev_phi), mais pas de simplif plus generale
if ((nc <- ncol(Md2hdbv2))>(blockSize)) {
message(paste0("Conditional-AIC computation requires operations on a large dense matrix (square with dimension ",nc,"),\n",
"which may take a lot of time. Use 'also_cAIC=FALSE' to avoid it."))
}
## if we reached this point with a huge dense matrix, then we have huge memory, hence saving memory may not be the issue,
# But we might save computation by computing only the requested diagonal (one qr, one backsolve, one sum(. * .))
# using Md2hdbv2=QRP, tr=sum_i (invP invR Qt(Md2hdbv2))_ii
qr.Md2hdbv2 <- try(qr(Md2hdbv2))
if (inherits(qr.Md2hdbv2,"try-error")) {
warning("Computation of cAIC/GoF df's failed because the information matrix appears singular.")
pd <- NA
} else {
solveR <- try(backsolve(qr.R(qr.Md2hdbv2),diag(nrow =nc)))
if (inherits(solveR,"try-error")) {
warning("Computation of cAIC/GoF df's failed because the information matrix appears singular.")
# determinant(qrR,logarithm=FALSE)$modulus > 1e-14
pd <- NA
} else {
# using inv(RP)= inv(R[,perm <- sort.list(.$pivot)]) = inv(R)[perm,] (as in h9[,perm] %*% solve(h9)[perm,] for any perm)
# though use of pivot may occur only for nearly singular matrices (poorly doc)
pd <- sum(solveR[sort.list(qr.Md2hdbv2$pivot),] * t(qr.qty(qr.Md2hdbv2,Md2clikdbv2)))
}
}
}
return(pd)
}
# Consistent with Saefken et al, only in terms of theta and mu. The linear predictor eta and link do not appear.
# Simulations for that paper use predictions eta only in the poisson(log) case;
# In the exponential (->Gamma(log)) they do use the theta deduced as -1/mu, not the eta. => cAIC4:::conditionalBootstrap is odd...
.calc_boot_AIC_dfs <- function (object, nsim, type="residual", seed=NULL, # (_F I X M E__) does not handle mv fits ? Could be docu'ed at least
fit_env=NULL, ...) {
family <- family(object)
if ( ! family$flags$exp) {
stop("Bootstrap bias correction not implemented for families not from GLM (exponential family) class.") # assuming $exp methods are always available for GLMs (but see negbin2_dvl)
} else if (family$family %in% c("negbin","negbin2")) {
warning("Bootstrap method does not strictly apply for negbin models.")
# Ad hoc version of .theta.mu.canonical that gives correct results in Poisson limit:
# as if the canonical parm was -log(1+shape/mu) + log(shape)
# This is not the commonly used formula, but still seems consistent with general theory
# (theta def'ed up to some forms of constants)
# devel/negbin_link has a version of HLfit_internals with this 'revised' canonical link.
thetafn <- function(x) {
refit <- update_resp(object, newresp=x)
muFREQS <- predict(refit, type="response")
shape <- environment(family(refit)$aic)$shape # allows shape to depend on refitted value
th <- -log(1+shape/muFREQS)+log(shape)
th
}
} else {
thetafn <- function(x) {
refit <- update_resp(object, newresp=x)
muFREQS <- predict(refit, type="response") # predict(refit(object, newresp = x))
.theta.mu.canonical(muFREQS,family(object)) # would return mu for non-GLMs
}
}
bootsims <- simulate(object, nsim = nsim, type = type, verbose=FALSE, seed=seed) # the corresponding lmer code returns a data frame
thetas <- dopar(bootsims, thetafn, fit_env=fit_env, ...)
if (inherits(thetas[[1]],"simpleError")) {
stop(paste("From dopar() call: ",
thetas[[1]]$message))
} else if (inherits(thetas,"try-error")) {
stop('dopar() call failed: use control=list(.errorhandling = "pass")\n to obtain more informative error message')
}
# if (is.factor(bootsims[1])) dataMatrix <- as.numeric(dataMatrix) - 1 # in cAIC4:::conditionalBootstrap
bootsims <- bootsims - rowMeans(bootsims)
phis <- residVar(object, which="phi") # For Gamma(), the phis are still those of he canonical form of the exponential family,
# not the full residual variance (name 'get_residVar' is ambiguous)
bootBC <- sum(colSums(thetas * bootsims/phis))/(nsim-1) # Handles the case Where the phis are heteroscedastic,
return(bootBC)
}
.calc_p_rdisp <- function(object, dfs=object$dfs) {
p_phi <- dfs[["p_fixef_phi"]]
if ( ! is.null(resid_fits <- object$resid_fits)) { # i.e mv fit =>$resid_fit*s*
p_phi <- sum(na.omit(unlist(p_phi)),
sum(unlist(lapply(resid_fits, `[[`, x="dfs"), recursive = TRUE, use.names = FALSE)) )
} else if ( ! is.null(resid_fit <- object$resid_fit)) {
# input p_phi (above) is typically set to NA, and will be ignored
p_phi <- sum(.unlist(resid_fit$dfs)) ## phi_pd is relevant only for measuring quality of prediction by the resid_fit!
} else p_phi <- sum(.unlist(dfs[["p_fixef_phi"]])) # .unlist() for mv
names_est_ranefPars <- unlist(.get_methods_disp(object))
fam_disp_parsnames <- intersect(names_est_ranefPars,c("NB_shape","COMP_nu","beta_prec",
"Tw_index","Tw_link"))
if (length(fam_disp_parsnames)) {
p_GLM_family <- length( # compatible with mv:
.unlist(.get_outer_inits_from_fit(object, keep_canon_user_inits = FALSE)[fam_disp_parsnames]))
p_phi <- p_phi+p_GLM_family ## effectively part of the model for residual error structure
}
p_phi
}
DoF <- function(object) {
dfs <- object$dfs
if ( ! inherits(dfs,"list")) dfs <- as.list(dfs) ## back compatibility
dfs[["p_fixef"]] <- dfs[["pforpv"]]
dfs[["p_rdisp"]] <- .calc_p_rdisp(object, dfs)
dfs["pforpv"] <- dfs["p_phi"] <- dfs["p_fixef_phi"] <- list(NULL) # make sure to erase any previously pre-existing redundant info
if (object$models[[1]]=="etaHGLM") {
if (is.null(p_corrPars <- dfs[["p_corrPars"]])) { ## back compatibility code,
dfs[["p_corrPars"]] <- length(which(unlist(attr(object$CorrEst_and_RanFix,"type")$corrPars, use.names = FALSE) %in% c("outer","var")))
# This is not the full count for up-to date spaMM (hyper param are missing),
# but should be OK for old objects to which this back compat code applies.
} # ELSE up do date objects should already have a correct p_corrPars
}
unlist(dfs)
}
.get_info_crits <- function(object, also_cAIC=TRUE, nsim=0L, ...) {
if (is.null(info_crits <- object$envir$info_crits) || (also_cAIC && is.null(info_crits[["cAIC"]]))) {
dof <- DoF(object)
p_rdisp <- dof[["p_rdisp"]]
APHLs <- object$APHLs
info_crits <- list()
# poisson-Gamma and negbin should have similar similar mAIC => NB_shape as one df or lambda as one df
forAIC <- APHLs
if (object$models[[1]]=="etaHGLM") {
if (object$HL[1]=="SEM") {
forAIC <- list(p_v=APHLs$logLapp,p_bv=APHLs$logLapp,clik=APHLs$clik)
}
# if standard ML: there is an REMLformula ~ 0 (or with ranefs ?); processed$X.Re is 0-col matrix
# if standard REML: REMLformula is NULL: $X.Re is X.pv, processed$X.Re is NULL
# non standard REML: other REMLformula: $X.Re and processed$X.Re identical, and may take essentially any value
X.pv <- model.matrix(object)
alldfs <- sum(dof)
info_crits$mAIC <- -2*forAIC$p_v + 2 *(alldfs)
info_crits$dAIC <- -2*forAIC$p_bv + 2 * (alldfs- dof[["p_fixef"]]) ## HaLM07 (eq 10) focussed for dispersion params
# including the rho param of an AR model
if (also_cAIC) {
if (.is_spprec_fit(object)) {
pd <- .calc_cAIC_pd_spprec(object)
} else if ( ! is.null(object$envir$sXaug)) {
pd <- .calc_cAIC_pd_from_sXaug(object)
} else { # SEM
H_w.resid <- .get_H_w.resid(object)
ZAL <- .compute_ZAL(XMatrix=object$strucList, ZAlist=object$ZAlist,as_matrix=.eval_as_mat_arg.HLfit(object))
d2hdv2 <- .calcD2hDv2(ZAL,H_w.resid,object$w.ranef) ## update d2hdv2= - t(ZAL) %*% diag(w.resid) %*% ZAL - diag(w.ranef)
pd <- .calc_cAIC_pd_others(X.pv, ZAL, H_w.resid, d2hdv2)
}
info_crits$GoFdf <- length(object$y) - pd ## <- nobs minus # df absorbed in inference of ranefs
## eqs 4,7 in HaLM07
info_crits$cAIC <- -2*forAIC$clik + 2*(pd+p_rdisp) ## no p_lambda ! Presumably if resid disp is declared as individual-level ranef,
# p_rdisp is reduced to 0 but pd is increased by as much so the result is still correct...
}
# print(c(pd,p_phi))
# but Yu and Yau then suggest caicc <- -2*clik + ... where ... involves d2h/db d[disp params] and d2h/d[disp params]2
} else { ## fixed effect model
info_crits$cAIC <- info_crits$mAIC <- -2*forAIC$p_v+2*(dof[["p_fixef"]]+p_rdisp)
# => sets cAIC so that (also_cAIC && is.null(info_crits[["cAIC"]])) becomes FALSE although also_cAIC is true by default...
}
object$envir$info_crits <- info_crits
} else p_rdisp <- NULL
if ( nsim>0L ) {
if (is.null(p_rdisp)) p_rdisp <- .calc_p_rdisp(object)
b_pd <- .calc_boot_AIC_dfs(object, nsim=nsim, ...)
if (object$models[[1]]=="etaHGLM") {
object$envir$info_crits$b_cAIC <- -2*object$APHLs$clik + 2*(b_pd+p_rdisp)
} else object$envir$info_crits$b_cAIC <- -2*object$APHLs$p_v + 2*(b_pd+p_rdisp) # as documented
}
return(object$envir$info_crits)
}
# Two ad hoc functions implying large atrix computations post-fit. In-fit .hatvals2std_lev() has comments about avoiding these matrix computations there. Not immediately transposable here.
.calc_dvdlogphiMat_new <- function(dh0deta,ZAL,
sXaug,d2hdv2_info=NULL ## either one
) {
## cf calcul dhdv, but here we want to keep each d/d phi_i distinct hence not sum over observations i
if (inherits(ZAL,"ZAXlist")) { # appeared in first try gaussian("logit")... with ARp()...
neg.d2h0_dv_dlogphi <- .crossprod(ZAL,.sparseDiagonal(x= drop(dh0deta), shape="g"))
} else neg.d2h0_dv_dlogphi <- .m_Matrix_times_Dvec(t(ZAL), drop(dh0deta)) # n_u_h*nobs: each ith column is a vector of derivatives wrt v_k# dh0dv <- t(ZAL) %*% diag(as.vector(dh0deta))
if (is.null(d2hdv2_info)) { # call by HLfit_body
dvdlogphiMat <- get_from_MME(sXaug,"solve_d2hdv2",B=neg.d2h0_dv_dlogphi)
} else if (inherits(d2hdv2_info,"CHMfactor")) { # CHM of ***-*** d2hdv2
dvdlogphiMat <- solve(d2hdv2_info, - neg.d2h0_dv_dlogphi) # efficient without as.matrix()!
} else if (inherits(d2hdv2_info,"qr") || inherits(d2hdv2_info,"sparseQR") ) {
dvdlogphiMat <- solve(d2hdv2_info, as.matrix(neg.d2h0_dv_dlogphi)) # rXn
} else if (is.environment(d2hdv2_info)) {
# dvdlogphiMat <- d2hdv2_info %*% neg.d2h0_dv_dlogphi # rXn
rhs <- d2hdv2_info$chol_Q %*% neg.d2h0_dv_dlogphi
rhs <- solve(d2hdv2_info$G_CHMfactor,rhs)
dvdlogphiMat <- - .crossprod(d2hdv2_info$chol_Q,rhs) # don't forget '-'
} else { ## then d2hdv2_info is ginv(d2hdv2) or a sparse matrix inverse of (d2hdv2) (spprec code will provide a dsCMatrix)
if (inherits(d2hdv2_info,"dsCMatrix")) {
d2hdv2_info <- as(d2hdv2_info,"dgeMatrix") ## more efficient if inv_d2hdv2 is math-dense
# It would be nice to store only the half matrix but then as( - d2hdv2_info, "dpoMatrix") and reversing sign afterwards.
}
dvdlogphiMat <- d2hdv2_info %*% neg.d2h0_dv_dlogphi # rXn
}
return(dvdlogphiMat)
}
.calc_dvdloglamMat_new <- function(neg.d2f_dv_dloglam, # a .sparseDiagonal or the result of .bdiag()
d2hdv2_info) {
if (inherits(d2hdv2_info,"CHMfactor")) {# CHM of ***-*** d2hdv2
dvdloglamMat <- solve(d2hdv2_info, - neg.d2f_dv_dloglam) # rXr !
} else if (inherits(d2hdv2_info,"qr") || inherits(d2hdv2_info,"sparseQR") ) { ## much slower than using CHMfactor
if (nrow(neg.d2f_dv_dloglam)>5000L) message("[one-time solve()ing of large matrix, which may be slow]")
dvdloglamMat <- solve(d2hdv2_info, neg.d2f_dv_dloglam) # rXr !
} else if (is.environment(d2hdv2_info)) {
# dvdloglamMat <- solve(d2hdv2_info, diag( neg.d2f_dv_dloglam )) # rXr !
rhs <- d2hdv2_info$chol_Q %*% neg.d2f_dv_dloglam
rhs <- solve(d2hdv2_info$G_CHMfactor, rhs)
dvdloglamMat <- - .crossprod(d2hdv2_info$chol_Q, rhs) # don't forget '-'
} else { ## then d2hdv2_info is ginv(d2hdv2) or some other form of inverse
# This case happens at least in my ressp devel example (tiny ranCoefs spprec)
if (inherits(d2hdv2_info,"dsCMatrix")) d2hdv2_info <- as(d2hdv2_info, "generalMatrix")
dvdloglamMat <- d2hdv2_info %*% neg.d2f_dv_dloglam
}
# I returned as.matrix(dvdloglamMat) a long time ago, and found it terribly inefficient in bigranefs case;
# but bigranefs example no longer runs calls .calc_dvdloglamMat_new()
return(dvdloglamMat)
} ## square matrix, by the formulation of the algo
..calc_d2hdv2_info <- function(object, ZAL) {
envir <- object$envir
if (is.null(factor_inv_Md2hdv2 <- envir$factor_inv_Md2hdv2)) { ## code not using promises
# we reach here at least in simple dense-correlation random-slope fit
# (where d2hdv2 may be sparse, and resulting d2hdv2_info too)
if ( ! is.null(envir$G_CHMfactor)) {
if (length(object$y) > ncol(envir$chol_Q)) { # We precompute inverse(d2hdv2) so that .calc_dvdlogphiMat_new() has ONE costly (r*r) %*% (r*n).
if ( is.null(envir$invL_G.P)) { # allowing for old objects that did not have this element
envir$factor_inv_Md2hdv2 <- Matrix::solve(envir$G_CHMfactor, as(envir$G_CHMfactor,"pMatrix") %*% envir$chol_Q,system="L") ## dgCMatrix
} else envir$factor_inv_Md2hdv2 <- envir$invL_G.P %*% envir$chol_Q ## dgCMatrix
# envir$factor_inv_Md2hdv2 <- Matrix::drop0(LL,tol=.Machine$double.eps)
d2hdv2_info <- - .crossprod(envir$factor_inv_Md2hdv2)
} else d2hdv2_info <- envir # then .calc_dvdlogphiMat_new() has TWO (r*r) %*% (r*n) (plus an efficient solve())
} else {
d2hdv2 <- .calcD2hDv2(ZAL,.get_H_w.resid(object),object$w.ranef)
if (inherits(d2hdv2,"sparseMatrix")) {
d2hdv2_info <- .silent_W_E(Cholesky( - d2hdv2,LDL=FALSE,perm=TRUE )) # '-' !
# this d2hdv2_info seems used only in solve(, default system="A") so perm is automatically handled.
if (inherits(d2hdv2_info, "CHMfactor")) {
#d2hdv2_info <- structure(d2hdv2_info, BLOB=list2env(list(), parent=environment(.solve_CHM)))
rank <- ncol(d2hdv2)
} else {
# oldMDCopt <- options(Matrix.warnDeprecatedCoerce = 0) # qr(<dsC>) problem in Matrix v1.4.2
d2hdv2_info <- qr(d2hdv2, tol=spaMM.getOption("qrTolerance")) # sparseQR
# options(oldMDCopt)
rank <- sum(abs(diag(qrR(d2hdv2_info,backPermute=FALSE)))>1e-7)
}
} else { # d2hdv2 was a matrix or a Matrix in dense format (dgeMatrix...)
d2hdv2_info <- qr(d2hdv2, tol=spaMM.getOption("qrTolerance")) # qr
rank <- d2hdv2_info$rank
}
if (rank<ncol(d2hdv2)) {
d2hdv2_info <- .force_solve(d2hdv2) # inverse(d2hdv2)
} # else we keep the QR facto.
}
} else {
# we reach here at least in simple ssprec random-slope fit ... so dense result here for spprec
d2hdv2_info <- - .crossprodCpp_d(as.matrix(factor_inv_Md2hdv2), yy=NULL) # inverse(d2hdv2)
# .crossprodCpp_d much faster than Matrix::crossprod(dsCMatrix...);
# still so for get_predVar(fit.Frailty) 2025/01/12;
# and still faster than any of the as_mat variants, eg .Rcpp_crossprod(factor_inv_Md2hdv2, BB = NULL, as_mat=TRUE)
# bc .crossprodCpp_d(dense) is fast relative to .Rcpp_crossprod(sparse), irrespective of as_mat
}
return(d2hdv2_info)
}
.calc_d2hdv2_info <- function(object, ZAL) {
if (TRUE) return(..calc_d2hdv2_info(object, ZAL)) # rediscovered long afterwards... so ./.
## The code below first checks for a "factor_inv_Md2hdv2" promise. But
## there is no such promise in the current object bc HLfit_body() has run
## .init_promises_spprec(sXaug, non_X_ones=FALSE, nullify_X_ones =TRUE)
## and non_X_ones=FALSE implies that there is no "factor_inv_Md2hdv2" promise.
#
# ELSE older version
envir <- object$envir
if ( ! is.null(envir$G_CHMfactor)) { ## spprec code
if (! .is_evaluated("factor_inv_Md2hdv2", envir)) {
if (length(object$y) > ncol(envir$chol_Q)) { # We precompute inverse(d2hdv2) so that .calc_dvdlogphiMat_new() has ONE costly (r*r) %*% (r*n).
d2hdv2_info <- - .crossprod(envir$factor_inv_Md2hdv2)
} else d2hdv2_info <- envir # then .calc_dvdlogphiMat_new() has TWO (r*r) %*% (r*n) (plus an efficient solve())
} else d2hdv2_info <- - .crossprodCpp_d(as.matrix(envir$factor_inv_Md2hdv2), yy=NULL)
} else {
if (is.null(factor_inv_Md2hdv2 <- envir$factor_inv_Md2hdv2)) {
d2hdv2 <- .calcD2hDv2(ZAL,object$w.resid,object$w.ranef)
if (inherits(d2hdv2,"sparseMatrix")) {
d2hdv2_info <- .silent_W_E(Cholesky( - d2hdv2,LDL=FALSE,perm=TRUE )) # '-' !
if (inherits(d2hdv2_info, "CHMfactor")) {
#d2hdv2_info <- structure(d2hdv2_info, BLOB=list2env(list(), parent=environment(.solve_CHM)))
rank <- ncol(d2hdv2)
} else {
d2hdv2_info <- qr(d2hdv2, tol=spaMM.getOption("qrTolerance")) # sparseQR
rank <- sum(abs(diag(qrR(d2hdv2_info,backPermute=FALSE)))>1e-7)
}
} else { # d2hdv2 was a matrix or a Matrix in dense format (dgeMatrix...)
d2hdv2_info <- qr(d2hdv2, tol=spaMM.getOption("qrTolerance")) # qr
rank <- d2hdv2_info$rank
}
if (rank<ncol(d2hdv2)) {
d2hdv2_info <- .force_solve(d2hdv2) # inverse(d2hdv2)
} # else we keep the QR facto.
} else d2hdv2_info <- - .crossprodCpp_d(as.matrix(envir$factor_inv_Md2hdv2), yy=NULL) # inverse(d2hdv2) # much faster than Matrix::crossprod(dsCMatrix...)
# and still faster than any of the as_mat variants, eg .Rcpp_crossprod(factor_inv_Md2hdv2, BB = NULL, as_mat=TRUE)
#
}
return(d2hdv2_info)
}
.calc_logdispObject <- function(object, envir=object$envir, force_fixed,
# promises:
strucList=object$strucList,
isRandomSlope= attr(strucList,"isRandomSlope")) {
dvdloglamMat <- envir$dvdloglamMat
dvdloglamMat_needed <- ( is.null(dvdloglamMat) &&
# The fact that next line is commented allows computation of predVar effect for random slope lambdas):
# all(unlist(attr(object$ZAlist,"namesTerms"))=="(Intercept)") && ## (1|.) or CAR or Matern
(force_fixed ||
any( ! object$lambda.object$type %in% c("fixed","fix_ranCoefs","fix_hyper"))) ) ## some lambda params were estimated
dvdlogphiMat <- envir$dvdlogphiMat
dvdlogphiMat_needed <- (is.null(dvdlogphiMat) &&
( any((phimodel <- object$models[["phi"]])=="phiScal") ||
force_fixed))
dvdlogphiMat_needed <- dvdlogphiMat_needed || identical(envir$forcePhiComponent,TRUE) ## hack for code testing !
if (dvdloglamMat_needed || dvdlogphiMat_needed) {
ZAL <- get_ZALMatrix(object, force_bind = ! (.is_spprec_fit(object)) )
d2hdv2_info <- .calc_d2hdv2_info(object, ZAL) # may be a qr object, or not (SPPREC). F I X M E a gentle message for long computations ?
}
if (dvdloglamMat_needed) { # $\partial_\tau \bv$ in latex doc.
cum_n_u_h <- attr(.get_u_h(object),"cum_n_u_h")
psi_M <- rep(attr(object$rand.families,"unique.psi_M"),diff(cum_n_u_h))
u_h <- .get_u_h(object)
dlogfthdth <- (psi_M - u_h)/object$lambda.object$lambda_est ## the d log density of th(u)
neg.d2f_dv_dloglam <- .calc_neg_d2f_dv_dloglam(dlogfthdth, cum_n_u_h,
lcrandfamfam=attr(object$rand.families,"lcrandfamfam"),
rand.families=object$rand.families, u_h=u_h)
## I tried using the private .fix_RHS_dvdloglamMat() here instead of the following simple conversion.
## See further comments there.
neg.d2f_dv_dloglam <- .sparseDiagonal(x= neg.d2f_dv_dloglam, shape="g")
dvdloglamMat <- .calc_dvdloglamMat_new(neg.d2f_dv_dloglam,
d2hdv2_info=d2hdv2_info) ## d2hdv2_info is either a qr factor or the inverse as a matrix or an environment
# envir$dvdloglamMat <- dvdloglamMat # needed only for debugging
}
if (dvdlogphiMat_needed) {
muetablob <- object$muetablob
# .get_H_w.resid() rather than .get_w.resid here. See section mentioning "dvdloglamMat" in the long doc.
if ( ! is.null(envir$G_CHMfactor)) { # spprec; possibly generalisable code not using math-dense ZAL
# rhs <- .Matrix_times_Dvec(t(envir$sXaug$AUGI0_ZX$ZAfix), - dh0deta) # efficient
# rhs <- solve(envir$G_CHMfactor,rhs,system="A") # efficient
unW_dh0deta <- (object$y-muetablob$mu)/muetablob$dmudeta ## (soit Bin -> phi fixe=1, soit BinomialDen=1)
if (is.null(envir$invG_ZtW)) {
if (is.null(envir$ZtW)) {
# cf comments in .old_calc_Md2hdvb2_info_spprec_by_r22()
envir$ZtW <- t(.Dvec_times_m_Matrix(.get_H_w.resid(envir=envir), envir$sXaug$AUGI0_ZX$ZAfix))
}
envir$invG_ZtW <- solve(envir$G_CHMfactor,
envir$ZtW, system="A") # hardly avoidable has there is no comparable operation elsewhere (check "A")
}
rhs <- .Matrix_times_Dvec(envir$invG_ZtW, -unW_dh0deta) # efficient
dvdlogphiMat <- .crossprod(envir$chol_Q, rhs) # _FIXME_ bottleneck in large spprec but .crossprodCpp_d not useful here
} else {
dh0deta <- ( .get_H_w.resid(object) *(object$y-muetablob$mu)/muetablob$dmudeta ) ## (soit Bin -> phi fixe=1, soit BinomialDen=1)
dvdlogphiMat <- .calc_dvdlogphiMat_new(dh0deta=dh0deta, ZAL=ZAL,
d2hdv2_info=d2hdv2_info ## either a qr factor or a matrix inverse or envir
)
}
# envir$dvdlogphiMat <- dvdlogphiMat # needed only for debugging
}
invV_factors <- .calc_invV_factors(object, r_x_r_needed = any(isRandomSlope)) ## n_x_r and r_x_n in repres of invV as diag(w.resid)- [n_x_r %*% r_x_n]
if (length(phimodel)>1L) {
.calc_logdisp_cov_mv(object, dvdloglamMat=dvdloglamMat, ## square matrix, by the formulation of the algo
dvdlogphiMat=dvdlogphiMat, invV_factors=invV_factors,
force_fixed=force_fixed)
} else
.calc_logdisp_cov(object, dvdloglamMat=dvdloglamMat, ## square matrix, by the formulation of the algo
dvdlogphiMat=dvdlogphiMat, invV_factors=invV_factors,
force_fixed=force_fixed)
}
.get_logdispObject <- function(object) { ##
envir <- object$envir
if (is.null(envir$logdispObject) && object$models[["eta"]]=="etaHGLM" ) {
envir$logdispObject <- .calc_logdispObject(object, envir=envir, force_fixed=FALSE)
}
return(envir$logdispObject)
} # if dvdloglamMat or dvdlogphiMat were computed ex-tempo, they are NOT saved.
## This provides factor n_x_r and r_x_n of the representation of invV as diag(w.resid)- [n_x_r %*% r_x_n = t(Ztw) %*% invG.ZtW] (nXr %*% rxn)
## slow computation the one time .get_logdispObject() is called, for variances$disp (no need to store the result in an $envir)
.calc_invV_factors <- function(object,
# $r_x_r returned for spprec or if 'r_x_r_needed' (set up for ranCoefs):
r_x_r_needed=FALSE) { ## used by .get_logdispObject
## Store inv( G=[{precmat=inv(L invWranef Lt)} +ZtWZ] ) as two matrix nXr and rXn rather than their nXn product
# F I X M E yet there will be cases where n<r and then it's better to store the n x n product !
if (.is_spprec_fit(object)) {
envir <- object$envir
## code clearly related to .Sigsolve_sparsePrecision() algorithm:
if (is.null(envir$invG_ZtW)) {
if (is.null(envir$ZtW)) {
# cf comments in .old_calc_Md2hdvb2_info_spprec_by_r22()
H_w.resid <- .get_H_w.resid(envir=envir)
envir$ZtW <- t(.Dvec_times_m_Matrix(H_w.resid, envir$sXaug$AUGI0_ZX$ZAfix))
}
object$envir$invG_ZtW <- solve(envir$G_CHMfactor, envir$ZtW, system="A") # hardly avoidable has there is no comparable operation elsewhere (check "A")
}
RES <- list(n_x_r=t(envir$ZtW), r_x_n=as.matrix(envir$invG_ZtW)) # requires both being kept in the envir
ZAfix <- .get_ZAfix(object, as_matrix=FALSE)
RES$r_x_r <- object$envir$invG_ZtW %*% ZAfix
return(RES)
} else {
ZAfix <- .get_ZAfix(object, as_matrix=FALSE)
H_w.resid <- .get_H_w.resid(object)
if (.is_identity(ZAfix)) {
# FIXME inelegant ZAL computation only to test if it is diagonal
ZAL <- .compute_ZAL(XMatrix=object$strucList, ZAlist=object$ZAlist,as_matrix=.eval_as_mat_arg.HLfit(object))
if (inherits(ZAL,"diagonalMatrix")) {
RES <- list(n_x_r=Diagonal(x = H_w.resid),
r_x_n=Diagonal(x = H_w.resid/(H_w.resid + object$w.ranef/diag(ZAL)^2)))
return(RES)
} ## ELSE
} ## ELSE
#
invL <- .get_invL_HLfit(object) # t(tcrossfac(precision_matrix))
wrZ <- .Dvec_times_m_Matrix(H_w.resid, ZAfix) # suppressMessages(sweep(t(ZA), 2L, w.resid,`*`))
if (TRUE) {
ZtwrZ <- .crossprod(ZAfix, wrZ, as_sym=TRUE) ## seems more precise and we must compute wrZ anyway
} else ZtwrZ <- .safe_ZtWZwrapper(ZAfix,.get_H_w.resid(object)) ## -> calls .crossprod(., y=NULL) affects numerical precision of twolambda test in test-predVar.R
## is general ZtwrZ should be sparse, while invL may not. Large dense invL will away be a problem
## for small Z, ZtwrZ may bedsy, but this an un-intersting subcase that does not call for a special treatment
if (inherits(ZtwrZ,"dsCMatrix") || inherits(ZtwrZ,"dsyMatrix")) {
# (NB: dsy+dsy faster than dsC+dsy but maybe just because of conversion from dsC to dsy ?) => no obvious improvement
if (.is_identity(invL)) {
precmat <- .symDiagonal(x=object$w.ranef) # dsC+dsC
} else precmat <- .safe_ZtWZwrapper(invL,object$w.ranef) # dsC+whatever
} else if (inherits(ZtwrZ,"ddiMatrix")) {
if (.is_identity(invL)) {
precmat <- Diagonal(x=object$w.ranef) # ddi+ddi...
} else precmat <- .safe_ZtWZwrapper(invL,object$w.ranef) # ddi+whatever
} else { # ("matrix") may never occur
if (.is_identity(invL)) {
precmat <- diag(x=object$w.ranef)
} else precmat <- .safe_ZtWZwrapper(invL,object$w.ranef)
message("Possibly inefficient code in .calc_invV_factors().")
}
## try to sum dsC or to sum dense matrix but not to mix types..., and to avoid formation of a large nxn matrix:
if (inherits(precmat,"dsCMatrix") && inherits(ZtwrZ,"dsCMatrix")) {
Gmat <- .dsCsum(precmat,ZtwrZ)
} else Gmat <- precmat + ZtwrZ
invG_ZtW <- tryCatch(solve(Gmat, t(wrZ)),error=function(e) e)
if (inherits(invG_ZtW,"simpleError")) { ## but that should be well behaved when precmat is.
invG <- .ginv(as.matrix(Gmat)) ## FIXME quick patch at least
invG_ZtW <- .tcrossprod(invG, wrZ)
}
RES <- list(n_x_r=wrZ, r_x_n=invG_ZtW)
if (r_x_r_needed) RES$r_x_r <- invG_ZtW %*% ZAfix
return(RES)
}
}
.calc_beta_cov_info_others <- function(wAugX=NULL, AUGI0_ZX, ZAL, ww) { ## post-fit fn
if (is.null(wAugX)) {
if (is.null(ZAL)) { # GLM... or LLM
if (any(ww<0)) { # ... LLM
negHess <- crossprod(AUGI0_ZX$X.pv, .Dvec_times_m_Matrix( ww, AUGI0_ZX$X.pv))
cholH <- tryCatch(chol(negHess),error=function(e) e)
if (inherits(cholH, "simpleError")) { # => brute regularization
warning("logLik presumably not maximized (information matrix is not positive definite at attained estimates).", immediate. = TRUE)
beta_cov <- tcrossfac_beta_v_cov <- matrix(NA, ncol=ncol(negHess), nrow=ncol(negHess))
} else {
beta_cov <- chol2inv(cholH)
tcrossfac_beta_v_cov <- t(cholH)
}
colnames(beta_cov) <- rownames(beta_cov) <- colnames(negHess) # sigh
return(list(beta_cov=beta_cov, tcrossfac_beta_v_cov=tcrossfac_beta_v_cov))
} else wAugX <- .calc_wAugX(XZ_0I=AUGI0_ZX$X.pv,sqrt.ww=sqrt(ww))
# => ... may have zero cols... X.pv being a zero-col *m*atrix, in which case .spaMM.data$options$matrix_method is selected below
# and then get_from_MME(<matrix_method>, "beta_cov_info_from_wAugX") must handle this case.
} else {
XZ_0I <- .calc_XZ_0I(AUGI0_ZX=AUGI0_ZX,ZAL) # ZAL is not ZAXlist since .calc_beta_cov_info_others not called in spprec
wAugX <- .calc_wAugX(XZ_0I=XZ_0I,sqrt.ww=sqrt(ww))
}
} ## wAugX is in XZ_OI order
if (inherits(wAugX,"Matrix")) {
sXaug_method <- .spaMM.data$options$Matrix_method
} else sXaug_method <- .spaMM.data$options$matrix_method
.silent_M_E(untrace(sXaug_method, where=asNamespace("spaMM"))) # try() bc this fails when called by Infusion
# test: Infusion tests -> ... -> .predict_body -> ;.. -> .calc_beta_cov_info_others -> untrace -> def_sXaug_EigenDense_QRP_Chol_scaled not found
# Note that the function is in exportPattern, explicitly exporting/impporting it does not help; attaching spaMM seems required to avoid untrace's error..
# hack to recycle sXaug code; all weights are 1 or unit vectors as the order is not that assumed by sXaug_method.
wAugX <- do.call(sXaug_method,list(Xaug=wAugX, weight_X=rep(1,nrow(AUGI0_ZX$X.pv)),
w.ranef=rep(1,ncol(AUGI0_ZX$I)), ## we need at least its length for get_from Matrix methods
H_global_scale=1))
attr(wAugX,"AUGI0_ZX") <- AUGI0_ZX # cf use of trDiag in .sXaug_Matrix_QRP_CHM_scaled -> solve_R_scaled
beta_cov_info <- get_from_MME(wAugX,"beta_cov_info_from_wAugX")
return(beta_cov_info)
}
.calc_beta_cov_info_from_sXaug <- function(BLOB, sXaug, tcrossfac) {
pforpv <- attr(sXaug,"pforpv")
X_scaling <- sqrt(rep(attr(sXaug,"H_global_scale"),pforpv))
X_scale <- attr(sXaug,"scaled:scale") # this is (! spprec) code and the X.pv attribute has been copied here
if ( ! is.null(X_scale)) X_scaling <- X_scaling/X_scale
diagscalings <- c(1/sqrt(attr(sXaug,"w.ranef")), X_scaling)
if ( ! is.null(BLOB$sortPerm)) { # depending on method used for QR facto
sctPmat <- sparseMatrix(seq_along(BLOB$sortPerm), BLOB$sortPerm, x=1)
sctPmat <- .Dvec_times_Matrix(diagscalings, sctPmat) # Pmat <- .Matrix_times_Dvec(Pmat,diagscalings)
tcrossfac_v_beta_cov <- sctPmat %*% tcrossfac
} else tcrossfac_v_beta_cov <- .Dvec_times_m_Matrix(diagscalings, tcrossfac) ## loses colnames...
dgC_good <- (inherits(tcrossfac_v_beta_cov, "dgCMatrix") &&
.calc_denseness(tcrossfac_v_beta_cov, relative=TRUE)<0.35 )
if ( ! dgC_good) tcrossfac_v_beta_cov <- as.matrix(tcrossfac_v_beta_cov) # bigranefs.R shows that conversion is not always good.
rownames(tcrossfac_v_beta_cov) <- colnames(sXaug) ## necessary for summary.HLfit, already lost in BLOB$R_scaled
seqp <- seq_len(pforpv)
beta_pos <- attr(sXaug,"n_u_h")+seqp
beta_v_order <- c(beta_pos,seq(attr(sXaug,"n_u_h")))
tcrossfac_beta_v_cov <- tcrossfac_v_beta_cov[beta_v_order,,drop=FALSE]
if (inherits(sXaug,"dtCMatrix")) tcrossfac_beta_v_cov <- as(tcrossfac_beta_v_cov, "sparseMatrix")
beta_cov <- as.matrix(.tcrossprod(tcrossfac_beta_v_cov[seqp,,drop=FALSE])) ## assignment in .make_beta_table() assumes a dense matrix
return( list(beta_cov=beta_cov,
#beta_v_cov=beta_v_cov,
tcrossfac_beta_v_cov=tcrossfac_beta_v_cov) )
}
.get_beta_cov_info <- function(res) {
# Provide list(beta_cov=., tcrossfac_beta_v_cov=.)
if (.is_spprec_fit(res)) {
if (is.null(dcdb_p4m <- res$envir$sXaug$AUGI0_ZX$dcdb_p4m)) dcdb_p4m <- res$X.pv # p4m spprec has $dcdb_p4m
return(.calc_beta_cov_info_spprec(X.pv=dcdb_p4m, envir=res$envir))
} else if (! is.null(res$envir$sXaug)) { # excludes SEM *and* fixed-effect models
if (prod(dim(res$envir$sXaug))>1e7) message("[one-time computation of covariance matrix, which may be slow]")
if (res$spaMM.version<"2.7.34") attr(res$envir$sXaug,"scaled:scale") <- attr(res$X.pv,"scaled:scale")
res$envir$beta_cov_info <- get_from_MME(res$envir$sXaug,which="beta_cov_info_from_sXaug")
} else { ## older code, generic for non-spprec.
# What this [block -> .calc_beta_cov_info_others()] provides is a beta_cov + full beta_v_cov
if (is.matrix(beta_cov_info <- res$envir$beta_cov_info) || ## old format (old fit object), should be a list now
is.null(beta_cov_info$tcrossfac_beta_v_cov)) { ## Use .calc_beta_cov_info_others -> get_from_MME()
# F I X M E test useless if call to .get_beta_cov_info() already dependent on the test
# Further, can $tcrossfac_beta_v_cov be NULL if $beta_cov_info is list ?
ZAL <- get_ZALMatrix(res, force_bind=FALSE) # note that .calc_beta_cov_info_others -> ... -> rbind2() but ZAL should not be a ZAXlist here
n_u_h <- length(res$w.ranef) # may be 0: this code may be run for LMs.
pforpv <- ncol(res$X.pv)
if (inherits(ZAL,"Matrix")) {
AUGI0_ZX <- list(I=.sparseDiagonal(n=n_u_h, shape="g"),
trDiag=..trDiagonal(n=n_u_h+pforpv), # trDiag needed by solve_R_scaled promise (from summary(<probitgem fit>)).
ZeroBlock=Matrix(0,nrow=n_u_h,ncol=pforpv),X.pv=res$X.pv)
} else { # ZAL is dense or is NULL (what if itis dense and X.pv is sparse?)
AUGI0_ZX <- list(I=diag(nrow=n_u_h),ZeroBlock=matrix(0,nrow=n_u_h,ncol=pforpv),X.pv=res$X.pv)
if (inherits(res$X.pv,"Matrix")) AUGI0_ZX$trDiag <- ..trDiagonal(n=n_u_h+pforpv) # useful for LM with sparse X,
# cf use of trDiag in .sXaug_Matrix_QRP_CHM_scaled -> solve_R_scaled
}
res$envir$beta_cov_info <- .calc_beta_cov_info_others(AUGI0_ZX=AUGI0_ZX,ZAL=ZAL,ww=c(.get_H_w.resid(res),res$w.ranef)) ## with beta_v_cov attribute
}
}
return(res$envir$beta_cov_info)
}
.calc_newZACvar <- function(newZAlist,cov_newLv_oldv_list) {
newZACvarlist <- vector("list",length(newZAlist))
for (new_rd in seq_along(newZAlist)) newZACvarlist[[new_rd]] <- newZAlist[[new_rd]] %ZA*gI% cov_newLv_oldv_list[[new_rd]]
newZAC <- .ad_hoc_cbind(newZACvarlist, as_matrix=FALSE)
if (inherits(newZAC,"sparseMatrix") && .calc_denseness(newZAC, relative=TRUE)>0.35) newZAC <- as.matrix(newZAC)
newZAC
}
## Aggregate info on corrpars, inner-estimated and inner-fixed.
## $corrPars is only for info in messages() and return value, (?!)
.make_CorrEst_and_RanFix <- function(ranFix, ## has "fix", "outer", and also "var" values ! code corrently assumes "var" <=> corr_est
corr_est
) {
# When fitting function was HLCor, the type attribute has not been added to the corrPars
if (is.null(attr(ranFix,"type"))) {
rFtype <- .relist_rep("fix", ranFix)
# HLCor may have been called through fitme/fitmv with inner estimation of lambda...
# Test case is first fitmv_body fit for BbyP_sp: initvec is NULL, inner estim of lambda
# get_fittedPars -> ... -> .get_ranPars_notPhi() finally gets the correct info from lambda.object$type.
# but get_ranPars gets info from $corrEst_and_ranFix which gets info from here
# This means we need to distinguish inner from fix here.
rFtype["lambda"] <- .modify_list(rFtype["lambda"],
relist("inner",ranFix["lambda"])[is.na(ranFix["lambda"])])
attr(ranFix,"type") <- rFtype
}
if ( ! is.null(corr_est)) {
ranFix <- structure(.modify_list(ranFix,corr_est),
type=.modify_list(attr(ranFix,"type"),.relist_rep("var",corr_est)))
}
return(ranFix) ## correlation params + whatever else was in ranFix
}
# Has become a postfit fn; to recycle processed$u_h_v_h_from_v_h one would have to remove processed from its defining envir.
.u_h_v_h_from_v_h <- function(v_h, rand.families, cum_n_u_h, lower.v_h, upper.v_h,
u_list=vector("list", length(rand.families))) {
if(!is.null(lower.v_h)) {v_h[v_h<lower.v_h] <- lower.v_h}
if(!is.null(upper.v_h)) {v_h[v_h>upper.v_h] <- upper.v_h}
nrand <- length(rand.families)
for (it in seq_len(nrand)) {
u.range <- (cum_n_u_h[it]+1L):(cum_n_u_h[it+1L])
u_list[[it]] <- rand.families[[it]]$linkinv(v_h[u.range])
if (any(is.infinite(u_list[[it]]))) {
warning("infinite random values ('u_h') were constrained to finite range.")
u_list[[it]] <- pmin(.Machine$double.xmax, pmax(-.Machine$double.xmax,u_list[[it]]) )
}
}
u_h <- .unlist(u_list)
## if there were box constr, v_h may have been modified, we put it in return value
if ( ! (is.null(lower.v_h) && is.null(upper.v_h))) attr(u_h,"v_h") <- v_h
return(u_h)
}
Try the spaMM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.