Nothing
## FR->FR accessors : http://glmm.wikidot.com/pkg-comparison
# lists of parameter names estimated by the two methods
.get_methods_disp <- function(object, family=object$family, families=object$families, phi.object=object$phi.object, mv_it=NULL) {
iterativeEsts <- character(0)
optimEsts <- character(0)
if ( ! is.null(families)) {
for (mv_it in seq_along(families)) {
methods_disp <- .get_methods_disp(object, family=families[[mv_it]], families=NULL, phi.object=phi.object[[mv_it]], mv_it=mv_it)
iterativeEsts <- c(iterativeEsts,methods_disp$iterativeEsts)
optimEsts <- c(optimEsts,methods_disp$optimEsts)
}
return(list(iterativeEsts=unique(iterativeEsts),optimEsts=unique(optimEsts)))
}
## What the HLfit call says:
if ( object$models[["eta"]] != "etaGLM") {
if ( any(object$lambda.object$type=="outer")) optimEsts <- c(optimEsts,"lambda")
if ( any(object$lambda.object$type=="inner")) iterativeEsts <- c(iterativeEsts,"lambda")
}
if ( (family$family %in% c("gaussian","Gamma"))) { # exponential families with a "phi" dispersion parameter
if ( ! is.null(phi.object$fixef) ) {
iterativeEsts <- c(iterativeEsts,"phi")
} else if ( identical(attr(phi.object$phi_outer,"type"),"var")) {
if (is.null(mv_it)) {
optimEsts <- c(optimEsts,"phi")
} else optimEsts <- c(optimEsts,paste("phi",mv_it,sep="_"))
}
# confint (fitme HLfitoide -> 3e cas)
}
optimNames <- names(attr(object,"optimInfo")$LUarglist$canon.init) ## will contain eg NB_shape
corrPars <- get_ranPars(object,which="corrPars")
optimEsts <- c(optimEsts,intersect(optimNames,c(names(corrPars),"COMP_nu","NB_shape","beta_prec")))
iterativeEsts <- c(iterativeEsts,setdiff(optimNames,optimEsts))
return(list(iterativeEsts=iterativeEsts,optimEsts=optimEsts))
}
.make_beta_table <- function(object, p_value="") {
namesOri <- attr(model.matrix(object),"namesOri")
nc <- length(namesOri)
betaOri_cov <- matrix(NA,ncol=nc,nrow=nc,dimnames=list(rownames=namesOri,colnames=namesOri))
beta_cov <- .get_beta_cov_any_version(object)
betaOri_cov[colnames(beta_cov),colnames(beta_cov)] <- beta_cov
beta_se <- sqrt(diag(betaOri_cov))
fixef_z <- object$fixef/beta_se
beta_table <- cbind(Estimate=object$fixef,"Cond. SE"=beta_se,"t-value"=fixef_z)
if (p_value=="Wald") {
beta_table <- cbind(beta_table,"p-value"=1-pchisq(fixef_z^2,df=1))
} else if (p_value) { # t-test for each coefficient
dft <- df.residual(object)
beta_table <- cbind(beta_table,"p-value"=2*(1-pt(abs(fixef_z), df=dft)))
}
rownames(beta_table) <- names(object$fixef)
return(beta_table)
}
.legend_lambda <- function(object, in_table, type="") {
## (1) analyse distributions
randfams <- object$rand.families
rff <- sapply(randfams, getElement, name="family")
# info about links:
rfl <- sapply(rff, switch,
# "Beta" = "log(lambda)", ## u ~ Beta(1/(2lambda),1/(2lambda))
# "inverse.Gamma" = "log(lambda)", ## u ~ I-Gamma(1+1/lambda,1/lambda)
# "log(lambda)" ## gaussian, or gamma with u ~ Gamma(sh=1/lambda,sc=lambda)
"log(lambda)" ## currently constant except as modified below for adjd
)
check_adjd <- any(unlist(lapply(object$lambda.object$coefficients_lambdaS, names))=="adjd")
if (check_adjd) {
if (object$spaMM.version < "2.2.116") {
whichadj <- attr(attr(object$ZAlist,"ranefs"),"type") %in% c("adjacency")
} else whichadj <- attr(object$ZAlist,"exp_ranef_types") %in% c("adjacency")
rfl[whichadj] <- "inverse[ lambda_i =var(V'u) ]"
rff <- rff[!whichadj] # and remove such terms from the family info
}
#
if (type=="family") {
# print family info except for adjd
urff <- unique(rff)
abyss <- lapply(urff, switch,
"Beta" = cat("lambda = 4 var(u)/(1 - 4 var(u)) for u ~ Beta[1/(2*lambda),1/(2*lambda)]; \n"),
"inverse.Gamma" = cat("lambda = var(u)/(1 + var(u)) for u ~ inverse-Gamma(sh=1+1/lambda, rate=1/lambda); \n"),
"Gamma" = cat("lambda = var(u) for u ~ Gamma(sh=1/lambda, sc=lambda); \n"),
"gaussian" = cat("lambda = var(u) for u ~ Gaussian; \n")
)
} else if (type=="link") {
# print link info
urfl <- unique(rfl[in_table])
cat(paste0(" --- Coefficients for ",paste(urfl,collapse=" or "),":"))
}
invisible(NULL)
}
# .degree <- rawToChar(as.raw(c(194L,176L))) # cf HelpersMG::char().
# with integers returned by HelpersMG::asc(), which is strtoi(charToRaw(x), 16L) also using only primitive fns.
## string compar by swithch) appears to fail:
# .ML_with_Lap_obs <- paste0("by ML (p",.degree,"_v approximation of logL).")
# .REML_with_Lap_obs <- paste0("by REML (p",.degree,"_bv approximation of restricted logL).")
.get_obsInfo <- function(object) { # back-compatibility utility
obsInfo <- object$how$obsInfo
if (is.null(obsInfo)) obsInfo <- FALSE
obsInfo
}
.MLmess <- function(object,ranef=FALSE) {
if (object$models[["eta"]]=="etaGLM") {
return("by ML.")
} else if (identical(attr(object$models,"LMMbool"),TRUE)) { # v3.5.9+
return("by ML.")
} else if (object$spaMM.version < "3.5.9" && ## back compat
object$family$family=="gaussian" &&
object$family$link=="identity" &&
all(attr(object$rand.families,"lcrandfamfam")=="gaussian")) {
return("by ML.")
} else {
if (object$HL[1]=='SEM') {
return("by stochastic EM.")
} else if (object$HL[1L]==1L) {
if (.get_obsInfo(object)) {
return("by ML (P_v approximation of logL).")
} else return("by ML (p_v approximation of logL).")
} else if (object$HL[1L]==0L) {
if (ranef) { # info for inner-estimated ranPars objective fn
if (.get_obsInfo(object)) {
return("by ML (P_v approximation of logL).")
} else return("by ML (p_v approximation of logL).")
} else return("by h-likelihood approximation.")
}
}
}
## FR->FR il faudrait distinguer EQL approx of REML ?
.REMLmess <- function(object,return_message=TRUE) {
## 'object' has no 'processed' element but its processed$REMLformula was copied to 'REMLformula' element.
## ./. It is by default NULL if REML was used, but may be an explicit non-default formula
## ./. It is an explicit formula if ML was used
if (return_message) {
if (is.null(object$REMLformula)) { ## default REML case
if (object$HL[1]=='SEM') {
resu <- ("by stochastic EM.")
} else if (object$models[["eta"]]=="etaGLM" && object$family$family =="gaussian") {
resu <- ("by REML.")
} else if ( identical(attr(object$models,"LMMbool"),TRUE)) {
resu <- ("by REML.")
} else if (object$spaMM.version < "3.5.9" &&
object$family$family=="gaussian" &&
object$family$link=="identity" &&
all(attr(object$rand.families,"lcrandfamfam")=="gaussian")) {
resu <- ("by REML.")
} else if (.get_obsInfo(object)) {
resu <- "by REML (P_bv approximation of restricted logL)."
} else resu <- "by REML (p_bv approximation of restricted logL)."
} else {
if (identical(attr(object$REMLformula,"isML"),TRUE)) { ## FALSE also if object created by spaMM <1.9.15
resu <- (.MLmess(object, ranef=TRUE))
} else { ## if nontrivial REML formula was used...
resu <- ("by non-standard REML")
attr(resu,"fixeformFromREMLform") <- .stripRanefs(object$REMLformula)
}
}
} else return( is.null(object$REMLformula) && object$HL[1]!='SEM')
return(resu)
}
summary.HLfitlist <- function(object, ...) {
if (object[[1]]$spaMM.version<"1.11.57") {
warning("This fit object was created with spaMM version<1.11.57, and is no longer supported.\n Please recompute it.")
} ## warning added in v2.2.43 2017/10/28
cumlvls <- character(0)
for (it in seq_along(object)) {
lvl <- attr(object[[it]]$data,"identifier")
cat(crayon::underline(paste0("Category ",lvl," vs. other one(s)")))
if (it>1L) cat(crayon::underline(paste0(" except ",paste(cumlvls, collapse=", "))))
cumlvls <- c(cumlvls,lvl)
cat("\n")
summary.HLfit(object[[it]])
}
cat(" ======== Global likelihood values ========\n")
zut <- attr(object,"APHLs")
cat(paste0(names(zut),": ",signif(unlist(zut),6), collapse="; "),"\n")
invisible(object)
}
.prettify_family <- function(family,linkstring="") {
famfam <- family$family
if ( ! is.null(withArgs <- attr(famfam,"withArgs"))) {
withArgs <- eval(withArgs,envir=environment(family$aic))
legend <- paste0(withArgs, "( ",linkstring, family$link," )")
} else legend <- paste0(famfam, "( ",linkstring, family$link," )")
if (identical(family$zero_truncated, TRUE)) legend <- paste("0-truncated", legend)
return(legend)
}
.get_compact_cov_mats <- function(strucList) {
isRandomSlope <- attr(strucList,"isRandomSlope")
if (any(isRandomSlope)) {
cov.mats <- vector("list",length(strucList))
for (it in seq_along(strucList)) {
if ( isRandomSlope[it] ) {
cov.mats[[it]] <- attr(strucList[[it]],"latentL_blob")$compactcovmat
}
## keep NULL slots for other elements as expected by summary.HLfit
}
return(cov.mats)
} else return(NULL)
}
# conversion of columns of data frames:
.aschar_adhoc <- function(data, digits=4) {
for (colit in seq_len(ncol(data))) if (is.numeric(colval <- data[,colit])) {
charvec <- as.character(signif(colval,digits))
charvec[is.na(charvec)] <- "" ## the reason why we create an ad hoc fn.
data[,colit] <- charvec
}
return(data)
}
.display_raw_lambdas <- function(in_pointLambda, row_map, lambda.object, hy_ranges) {
# if (details$ranCoefs) {
# displaypos <- innerlambda_pos
# } else displaypos <- setdiff(innerlambda_pos, random_slope_pos)
# displayrows <- unlist(row_map[displaypos])
#if ( ! details$ranCoefs) displaypos <- setdiff(displaypos, random_slope_pos)
displaypos <- which(in_pointLambda)
displayrows <- unlist(row_map[displaypos])
nicertypes <- lambda.object$type[displaypos]
posf <- ( nicertypes=="fixed")
posfh <- ( nicertypes=="fix_hyper") # see paste0(hy_lam_type,"_hyper") in .calc_initial_init_lambda()
posoh <- ( nicertypes=="outer_hyper") # again, see paste0(hy_lam_type,"_hyper") in .calc_initial_init_lambda()
nicertypes[posf] <- "[fixed]"
nicertypes[posfh] <- "[fixed]" # "[fix_hyper]" # need to distinguish in summary not obvious
for (it in seq_along(hy_ranges)) {
hy_range <- hy_ranges[[it]]
if (nicertypes[hy_range[1L]] == "outer_hyper") {
nicertypes[hy_range[1L]] <- "" ## outer-optim
nicertypes[hy_range[-1L]] <- ("[set by hyper]")
}
}
nicertypes[ ! (posf | posfh | posoh )] <- ""
nicertypes <- rep(nicertypes, unlist(lapply(row_map[displaypos], length)))
if ( ! is.null(displayrows)) {
lambda_list <- lambda.object$lambda_list # this includes a vector of 1's for ranCoefs with Xi_ncol>1L, and the actual lambda for Xi_ncol=1
if (is.null(lambda_list)) lambda_list <- lambda.object$lambda # back-compatibility fix (from test on v2.5.0 object stored in package vullioud2018)
print_lambda <- unlist(lambda_list)
cat(paste(" ",
names(displayrows)," : ",
signif(print_lambda[displayrows],4),
nicertypes,
collapse="\n"),"\n")
}
}
.lambda_table_fn <- function(namesTerms, object, lambda.object, linklam_coeff_list=lambda.object$coefficients_lambdaS) {
nrand <- length(namesTerms)
nrows <- unlist(lapply(namesTerms,length))
cum_nrows <- cumsum(c(0,nrows))
names(cum_nrows) <- NULL
row_map <- lapply(nrows,seq)
for (it in seq_len(length(row_map))) row_map[[it]] <- row_map[[it]]+cum_nrows[it]
# first construct a table including NA's for some coefficients (not "inner" estimated), then remove these rows
term_names <- names(namesTerms)
repGroupNames <- rep(term_names, sapply(namesTerms,length))
## i.e. for namesTerms = list(Subject=c("(Intercept)", "Days")), repGroupNames[[1]] is c("Subject", "Subject")
lambda_table <- data.frame(Group=repGroupNames,Term=unlist(namesTerms))
in_table <- rep(FALSE,nrand)
in_pointLambda <- rep(TRUE,nrand)
maxnrow <- cum_nrows[nrand+1] ## maxnrow should = nrow(lambda_table)
#
############### better initialization of the table (03/2021)
template <- lambda.object$lambda_list # RHS has elements for each of the rows of the table, except the "adjd" coeff. Hence:
if (is.null(template)) template <- lambda.object$lambda # back compat fix (pointed by fits in vullioud2018 by version 2.5.0)
for (it in seq_len(length(namesTerms))) {
if ("adjd" %in% namesTerms[[it]]) {
template[[it]][2] <- NA
}
}
summ_variances <- data.frame("Var."=matrix(.unlist(template),ncol=1L,nrow=maxnrow))
# So now this df has a col with:
# lambda values for 'simple lambdas', and ranCoefs without corr matrix (Xi_ncol=1L)
# '1's for the lambdas of ranCoefs with corr matrix (Xi_ncol>1)
# a variance and an NA for adjacency terms
# Previously the column was initiated with NA's and lambdas for ranCoefs without corr matrix (Xi_ncol=1L) were never filled in => defective VarCorr().
################
cov.mats <- .get_compact_cov_mats(object$strucList)
if ( length(cov.mats)) { ## fixme ? rename cov.mats to refer to ranCoefs ?
#.varcorr <- function(nrows, maxnrow, cov.mats, in_table, in_pointLambda, cum_nrows) {
summ_corr_cols <- data.frame(matrix(NA,ncol=max(nrows-1L),nrow=maxnrow))
for (mt in seq_len(length(cov.mats))) {
m <- cov.mats[[mt]]
if ( ! is.null(m)) { # Xi_ncol>1L...
# if ( # kappa(m)>1e14 ||
# min(eigen(m,only.values = TRUE)$values)<1e-6) { # lme4:::rePCA.merMod uses svd() [less assuming but beware of sign]
# # Motivations for the check discussed in Bates et al 2015 on arXiv, but other views are possible...
# message(paste0("Random-coefficient term '",term_names[[mt]],"' has a practically singular covariance matrix"))
# }
in_table[mt] <- TRUE
in_pointLambda[mt] <- FALSE
inrows <- cum_nrows[mt]+(1:nrow(m))
summ_variances[inrows,1] <- diag(m)
m <- cov2cor(m)
for (it in (2:nrow(m))) {
for (jt in (1:(it-1))) {
summ_corr_cols[cum_nrows[mt]+it,jt] <- m[it,jt]
}
}
}
}
colnames(summ_corr_cols) <- rep("Corr.",ncol(summ_corr_cols))
# }
random_slope_ncol_geq_1_pos <- which(attr(object$strucList,"isRandomSlope"))
random_slope_ncol_geq_1_rows <- unlist(row_map[ random_slope_ncol_geq_1_pos ])
} else {
random_slope_ncol_geq_1_rows <- random_slope_ncol_geq_1_pos <- integer(0)
summ_corr_cols <- NULL
}
if ( ! is.null(linklam_coeff_list)) {
in_table[which( ! unlist(lapply(linklam_coeff_list,is.null)))] <- TRUE
lambda_table <- cbind(lambda_table, Estimate=unlist(linklam_coeff_list),
"Cond.SE"=lambda.object$lambda_se)
info_rows <- which(! is.na(lambda_table$Estimate)) ## must be evaluated before the next line sets more NAs
} else info_rows <- NULL
for (it in seq_len(length(namesTerms))) {
if ("adjd" %in% namesTerms[[it]]) {
in_pointLambda[it] <- FALSE
in_table[it] <- TRUE
}
}
lambda_table <- cbind(lambda_table, summ_variances)
if ( length(cov.mats)) lambda_table <- cbind(lambda_table, summ_corr_cols)
lambda_table <- structure(lambda_table,
class=c("lambda_table",class(lambda_table)), info_rows=info_rows,
random_slope_ncol_geq_1_rows=random_slope_ncol_geq_1_rows, # rows of terms with Xi_ncol>1
random_slope_ncol_geq_1_pos=random_slope_ncol_geq_1_pos,
in_pointLambda=in_pointLambda, # defines which ranefs in the "Variance parameters ('lambda')" table in summary
# (not adjd nor Xi_ncol>1-ranCoefs)
row_map=row_map,
in_table=in_table)
return(lambda_table)
}
.print_adj_ests <- function(object, namesTerms, linklam) {
ncoeffs <- attr(object$ZAlist,"Xi_cols") ## Xi_cols elements may be not 1 for ranCoefs terms
any_adjd <- FALSE
for (it in seq_len(length(namesTerms))) if ("adjd" %in% namesTerms[[it]]) {
ncoeffs[it] <- 2L
any_adjd <- TRUE
}
if (any_adjd) {
cat(" --- Variance parameters ('lambda'):\n") ## NOT the (Gamma-GLM) coefficients and SEs
}
cum_ncoeffs <- c(0,cumsum(ncoeffs))
for (it in seq_len(length(namesTerms))) {
if ("adjd" %in% namesTerms[[it]]) {
namenames <- names(namesTerms[it])
pos <- cum_ncoeffs[it]+1L
cat(paste("Estimate of rho (",namenames,"CAR): ",
signif( - linklam[pos+1L]/linklam[pos],4),"\n"))
cat(paste("Estimate of lambda factor (",namenames,"CAR): ",
with(object$lambda.object,signif(linkinvS[[rand_to_glm_map[it]]](linklam[pos]),4)),"\n"))
}
}
return(any_adjd)
}
.print_lambda_table <- function(object,lambda_table, details, namesTerms, linklam_coeff_list, digits=4) {
attribs <- attributes(lambda_table)
# strings for screen output:
lambda_table <- .aschar_adhoc(lambda_table, digits=digits)
if (length(attribs$random_slope_ncol_geq_1_rows)) {
keep <- which( ! colnames(lambda_table) %in% c("Estimate","Cond.SE")) ## may be null if only outer... + fixed
cov_table <- lambda_table[attribs$random_slope_ncol_geq_1_rows,keep] ## EXCLUDES the "Estimate","Cond.SE" cols
cat(" --- Random-coefficients Cov matrices:\n")
print(cov_table, digits = 4, row.names = FALSE)
}
#
keep <- which( colnames(lambda_table) %in% c("Group","Term","Estimate","Cond.SE")) ## may be null if only outer... + fixed
lambda_table <- lambda_table[,keep]
#
# point lambda's (+adjd rho) estimates, NOT the (Gamma-GLM) coefficients and SEs
linklam <- unlist(linklam_coeff_list)
any_adjd <- .print_adj_ests(object, namesTerms, linklam)
if (any(attribs$in_pointLambda)) {
if ( ! any_adjd) cat(" --- Variance parameters ('lambda'):\n")
.legend_lambda(object, type="family") ## print info about the meaning of lambda according to the rand family (not rand link)
.display_raw_lambdas(in_pointLambda=attribs$in_pointLambda,
row_map=attribs$row_map,
object$lambda.object,
object$ranef_info$hyper_info$ranges
) ## not the table with SEs / covariances
}
#
# NOW the (Gamma-GLM) coefficients and SEs
info_rows <- attribs$info_rows
if ( ! details$ranCoefs) info_rows <- setdiff(info_rows, attribs$random_slope_ncol_geq_1_rows) ## removes random-coef info except if detaisl are requested
lambda_table <- lambda_table[ info_rows ,]
in_table <- attribs$in_table
if (nrow(lambda_table)) {
if (any(in_table) && ! is.null(linklam_coeff_list[in_table])) { ## (mixed cov.mat and) coefficients_lambdaS output
.legend_lambda(object, in_table, type="link")
cat("\n")
# } else { ## only ranCoefs
# cat("Variances and correlation for random-coefficient terms:\n")
} ## else only other terms with outer lambda (eg Matern term)
print(lambda_table,digits=4,row.names=FALSE)
}
invisible(lambda_table) ## not currently correct as it is truncated
}
# # Nice names for the raw lambda vector in the HLfit object:
# .rename_res_lambda <- function (object) {
# lambda.object <- object$lambda.object
# namesTerms <- lambda.object$print_namesTerms ## list of vectors of variable length
# linklam_coeff_list <- lambda.object$coefficients_lambdaS
# lambda_table <- .lambda_table_fn(namesTerms, object, lambda.object,linklam_coeff_list) ## uses object$strucList
# attribs <- attributes(lambda_table)
# displayrows <- unlist(attribs$row_map)
# lambda <- unlist(lambda.object$lambda_list)
# structure(lambda, names=names(displayrows))
# }
# # => this function is not used. Later I wrote:
# .print_lambda <- function(lambda_list) { # to have the right names for a 'lambda' vector by unlisting à 'lambda_list'.
# for (rd in seq_along(lambda_list)) {
# if (length(lambda_list[[rd]])==1L && names(lambda_list[[rd]])=="(Intercept)") names(lambda_list[[rd]]) <- NULL
# # in which case unlisting will paste the group label (say "Subject2") and NULL => "Subject2.
# # In all other cases unlisting lambda_list will paste the unique group label
# # (say "Subject1") and the coefficient names (say "(Intercept)" and "age") => "Subject1.(Intercept)" and "Subject1.age"
# }
# unlist(lambda_list)
# }
.table_glm_phi <- function(glm_phi, loc_p_phi, phi.object, summ, object, family, # 'family' is response fam (of mv_it submodel if mv)
phiform=.get_phiform(object), resid.family=eval(.get_phifam(object))
) {
phi_se <- summary(glm_phi,dispersion=1)$coefficients[(loc_p_phi+1L):(2L*loc_p_phi)]
## note dispersion set to 1 to match SmythHV's 'V_1' method, which for a log link has steps:
#SmythHVsigd <- as.vector(sqrt(2)*phi_est);SmythHVG <- as.vector(phi_est); tmp <- SmythHVG / SmythHVsigd
## tmp is here sqrt(2) !
#if (length(tmp)>1) {SmythHVZstar <- diag(tmp) %*% X_disp} else SmythHVZstar <- tmp * X_disp
#SmythHVcovmat <- solve(.ZtWZ(SmythHVZstar,(1-lev_phi))); phi_se <- sqrt(diag(SmythHVcovmat)) print(phi_se)
#
phi_table <- cbind(phi.object$fixef,phi_se)
colnames(phi_table) <- c("Estimate", "Cond. SE")
rownames(phi_table) <- namesX_disp <- names(phi.object$fixef)
summ$phi_table <- phi_table
phiinfo <- resid.family$link
if (phiinfo=="identity") {phiinfo="phi "} else {phiinfo <- paste0(phiinfo,"(phi) ")}
phiinfo <- paste("Coefficients for",phiinfo,paste0(phiform,collapse=" ")," :\n")
cat(phiinfo)
print(phi_table,4)
dispOffset <- attr(phiform,"off")
if (!is.null(dispOffset)) dispOffset <- unique(dispOffset)
if (length(namesX_disp)==1 && namesX_disp[1]=="(Intercept)" && length(dispOffset)<2L) {
# constant phi: we can display it
phi_est <- (phi.object$fixef)
if (length(dispOffset)==1L) phi_est <- phi_est+dispOffset
phi_est <- resid.family$linkinv(phi_est)
if (family$family=="Gamma") {
cat(paste("Estimate of phi: ",signif(phi_est,4),"\n"))
## la var c'est phi mu^2...
} else cat(paste("Estimate of phi=residual var: ",signif(phi_est,4),"\n"))
} ## else phi not constant; We don't try to display it
wa <- glm_phi$warnmess
if ( ! is.null(wa)) {
if (wa=="glm.fit: algorithm did not converge") {
cat("glm.fit for estimation of phi SE did not converge; this suggests\n")
cat(" non-identifiability of some phi (and possibly also lambda) coefficients.\n")
} else {
cat("warning in glm.fit for estimation of phi SE: \n")
cat(wa,"\n")
}
}
return(summ)
}
.is_unit <- function(pw) {
if (is.null(is_unit <- attr(pw, "is_unit"))) {
is_unit <- (identical(attr(pw,"unique"),TRUE) && pw[1]==1L)
}
is_unit
}
.summary_phi_object <- function(object, phi.object=object$phi.object, family=object$family,
pw=object$prior.weights, summ, phimodel, mv_it=NULL) { # 'mv_it' needed for non-trivial .get_glm_phi(object, it=it)
if (family$family %in% c("gaussian","Gamma")) {
if (! is.null(mv_it)) {
cat(crayon::underline("* response", mv_it))
if (family$family=="Gamma") {
cat(" (Gamma) residual var = phi * mu^2:\n")
} else cat(" (gaussian) residual variance: \n")
} else {
cat(" -")
if (family$family=="Gamma") {
cat("-- Residual variation ( var = phi * mu^2 ) --\n")
} else cat("------------- Residual variance ------------\n")
}
if ( ! .is_unit(pw)) cat(paste("Prior weights:",
paste(signif(pw[1:min(5,length(pw))],6),collapse=" "),
"...\n"))
if ( ! is.null(phi_outer <- phi.object$phi_outer)) {
if ( identical(attr(phi_outer,"type"),"fix") ) {
if (length(phi_outer)==1L) {
cat(paste("phi was fixed to",signif(phi_outer,6),"\n"))
} else if (phimodel=="phiGLM") { # fixed by an offset in a resid.formula (see .preprocess_phi_model)
phiform <- .get_phiform(object, mv_it)
if (length(phiform)==2) phiform <- as.formula(paste('"phi"',paste(phiform,collapse=" ")))
cat(paste("phi was fixed [through ",deparse(phiform),"] to",
paste(signif(phi_outer[1:min(5,length(phi_outer))],6),collapse=" "),"...\n"))
} else cat(paste("phi was fixed to",paste(signif(phi_outer[1:min(5L,length(phi_outer))],6),collapse=" "),"...\n"))
} else {
if (length(phi_outer)==1L) {
cat(paste("phi estimate was",signif(phi_outer,6),"\n"))
} else cat(paste("phi was estimated.\n"))
}
#summ$phi_outer <- phi.object$phi_outer ## mv: {actually from phi.object[[mv_it]]; overwrites previous mv_it}
# spaMM does nothing of it anyway so I remove
} else {
if (! is.null(mv_it)) {
is_hlfit <- ! is.null(object$resid_fits[[mv_it]])
} else is_hlfit <- ! is.null(object$resid_fit)
if (is_hlfit) {
if (phimodel=="phiHGLM") {
if (! is.null(mv_it)) {
cat(paste0("Residual dispersion model includes random effects:\n use summary(<fit object>$resid_fits[[",mv_it,"]]) to display results.\n"))
} else cat("Residual dispersion model includes random effects:\n use summary(<fit object>$resid_fit) to display results.\n")
} else if (! is.null(mv_it)) {
cat(paste0("Residual dispersion model is HLfit object:\n use summary(<fit object>$resid_fits[[",mv_it,"]]) to display results.\n"))
} else cat("Residual dispersion model is HLfit object:\n use summary(<fit object>$resid_fit) to display results.\n") # _F I X M E__ that is reasonable output but try extracting info from the fits?
} else if ((loc_p_phi <- length(phi.object$fixef))) { # there are phi params (possibly outer estimated),
glm_phi <- .get_glm_phi(object, mv_it=mv_it)
summ <- .table_glm_phi(glm_phi, loc_p_phi, phi.object, summ, object, family=family,
phiform=.get_phiform(object,mv_it), resid.family=eval(.get_phifam(object,mv_it)))
} else {
# Other cases should have been handled by identical(attr(phi_outer,"type"),"fix") && phimodel=="phiGLM" above
warning("Unexpected case in .summary_phi_object(): maybe harmless, but please contact the maintainer.") # warning set on 3.7.24 03/2021
}
}
} else if (family$family %in% c("beta_resp","betabin", "COMPoisson", "negbin1", "negbin2")) {
not_pw_1 <- ! .is_unit(pw) # But relevant only for beta_resp and betabin (which have pw)? (__F I X M E___?? Don't confuse with resid.model feature which works for negbin[] models too)
has_dispenv_beta <- ( ! is.null(beta <- (disp_env <- family$resid.model)$beta))
if (not_pw_1 || has_dispenv_beta) {
if (! is.null(mv_it)) {
cat(crayon::underline("* response", mv_it))
info <- switch(family$family,
"beta_resp" = " dispersion model for beta_resp\n",
"betabin" = " dispersion model for betabin\n",
"COMPoisson"= " dispersion model for COMPoisson\n",
"negbin1" = " dispersion model for negbin1\n",
"negbin2" = " dispersion model for negbin2\n",
"family dispersion parameter"
)
} else info <- switch(family$family,
"beta_resp" = "------- Dispersion model for beta_resp --------\n",
"betabin" = "------- Dispersion model for betabin ----------\n",
"COMPoisson"= "------- Dispersion model for COMPoisson -------\n",
"negbin1" = "------- Dispersion model for negbin1 ----------\n",
"negbin2" = "------- Dispersion model for negbin2 ----------\n",
"family dispersion parameter"
)
cat(info)
}
if (not_pw_1) cat(paste("Prior weights:",
paste(signif(pw[1:min(5,length(pw))],6),collapse=" "),
"...\n"))
if (has_dispenv_beta) {
info <- switch(family$family,
"beta_resp" = "prec",
"betabin" = "prec",
"COMPoisson"= "nu",
"negbin1"="shape",
"negbin2"="shape",
"family dispersion parameter"
)
info <- paste0("Coefficients for log(",info,")", deparse(disp_env$resid.formula)," :\n")
dispcoef_table <- cbind(disp_env$beta, NA)
colnames(dispcoef_table) <- c("Estimate", "Cond. SE")
summ$dispcoef_table <- dispcoef_table
cat(info)
print(dispcoef_table,6)
}
} ## else binomial or poisson, no dispersion param
return(summ)
}
.check_inverse_pred <- function(object, family=object$family) {
if (family$link=="inverse") {
if (family$family=="Gamma" &&
identical(attr(object$eta,"any_neg_eta"),TRUE)
) {
"Non-positive predictions implied by final fitted coefficients in Gamma(inverse)-response model.\n"
} else if (family$family=="gaussian" &&
identical(attr(object$eta,"any_tiny_eta"),TRUE)
) {
"Practically zero predictions implied by final fitted coefficients in gaussian(inverse)-response model.\n"
}
}
} # returns NULL or message
`summary.HLfit` <- function(object, details=FALSE, max.print=100L, verbose=TRUE, ...) {
parent_from_there <- parent.frame()
if ( ! verbose) {
mc <- match.call(expand.dots = TRUE)
mc$verbose <- TRUE
capture.output({silent <- eval(mc, envir=parent_from_there)})
return(silent)
}
oldopt <- options(max.print=max.print)
if (is.null(names(details))) details <- structure(rep(details,2),names=c("ranCoefs","p_value")) ## handle FALSE or TRUE input
details <- as.list(details)
if (is.null(details[["ranCoefs"]])) details["ranCoefs"] <- FALSE
if (is.null(details[["p_value"]])) details["p_value"] <- "" ## a string such as "Wald"
if (is.null(details[["digits"]])) details["digits"] <- 4 # _F I X M E__ document this? Only used in .print_lambda_table() so rather cryptic.
models <- object$models
lcrandfamfam <- attr(object$rand.families,"lcrandfamfam")
randfamfamlinks <- unlist(lapply(object$rand.families, .prettify_family))
randfamlinks <- unlist(lapply(object$rand.families, getElement, name="link"))
summ <- list()
form <- formula.HLfit(object, which="hyper")
if (is.list(form)) {
for (mv_it in seq_along(form)) {cat("formula_",mv_it,": ",sep=""); print(form[[mv_it]],showEnv=FALSE)}
} else {
cat("formula: ")
print(form,showEnv=FALSE)
}
#
# HLchar <- paste(as.character(object$HL),collapse="")
# cat(paste0("[code: ",HLchar,"]"," method: "))
messlist <- list()
if (length(object$fixef)) messlist[["fixed"]] <- .MLmess(object)
messlist[["ranef"]] <- .REMLmess(object)
summ$formula <- object$formula
summ$REMLformula <- object$REMLformula
## Distinguishing iterative algo within HLfit and numerical maximization outside HLfit
locblob <- .get_methods_disp(object)
iterativeEsts <- locblob$iterativeEsts
optimEsts <- locblob$optimEsts
RE_Ests <- c(iterativeEsts,optimEsts)
RE_Ests <- gsub("phi_.*", "phi", RE_Ests) # for mv
RE_Ests <- unique(RE_Ests)
len <- length(RE_Ests)
lenIt <- length(iterativeEsts) ## p_v maxim, or p_bv(X.Re)
lenOpt <- length(optimEsts) ## p_v maxim, or p_bv(X.Re)
if (len > 1) RE_Ests <- paste(c(paste(RE_Ests[-len],collapse=", "),RE_Ests[len]),collapse=" and ")
#if (lenIt > 1) iterativeEsts <- paste(c(paste(iterativeEsts[-lenIt],collapse=", "),iterativeEsts[lenIt]),collapse=" and ") # nothing is done of that
## The next block outputs classes of estimated parameters and method, and further code prints apart those that are outer estimated
## However, if there is no iterative est, only the latter 'outer' message (without the words 'Laplace... approximation') is output
## So if no fixed effect is estimated, no message with the words 'Laplace... approximation' is output (subliminal imperfection)
if (lenIt) { # if there are inner estimated parameter,s we output a message distinct from the 'outer' pars message, listing all estimated disp params
if (messlist[["ranef"]]=="by REML.") {## REMLmess has checked that this is a LMM
cat("REML: Estimation of ")
cat(RE_Ests);
tab <- " "
} else if (messlist[["ranef"]]=="by ML.") {
cat("ML: Estimation of ")
cat(RE_Ests);
tab <- " "
} else {
cat("Estimation of ")
cat(RE_Ests);
tab <-""
}
## the aim of the tab is to align "Estimation of..." vertically
cat(" ")
cat(messlist[["ranef"]])
if ( messlist[["ranef"]]=="by non-standard REML") {
cat("\n");cat(tab);cat(" based on fixed-effects model: ")
print(attr(messlist[["ranef"]],"fixeformFromREMLform"),showEnv=FALSE)
} else cat("\n") ## normal output for standard REML formula
} else tab <- ""
if (length(object$fixef)) {
cat(tab)
cat("Estimation of fixed effects ")
cat(messlist[["fixed"]]);
cat("\n")
}
if (lenOpt > 1) optimEsts <- paste(c(paste(optimEsts[-lenOpt],collapse=", "),optimEsts[lenOpt]),collapse=" and ")
if (lenOpt > 0) {
objective <- attr(object,"optimInfo")$objective
if(is.null(objective)) {
#stop("attr(object,'optimInfo')$objective is missing: malformed object.")
## may happen when one refits an HLCorcall =>
## optimEsts ultimately deduced by its $processed$ranPars as in the original fit,
## but no optimInfo in the refit.
} else {
objString <- switch(objective,
p_bv= "'outer' REML, maximizing restricted logL",
p_v= "'outer' ML, maximizing logL",
cAIC= "'outer' minimization of cAIC",
paste("'outer' maximization of",objective)
)
outst <- paste0("Estimation of ",optimEsts," by ",objString,".\n")
cat(outst)
}
}
if (is.null(object$family)) {
cat("Families: ")
nfam <- length(object$families)
famst <- character(nfam)
inverse_pred_mess <- vector("list", nfam)
for (mv_it in seq_along(object$families)) {
famst[mv_it] <- paste0(mv_it,": ", .prettify_family(object$families[[mv_it]], linkstring = "") , sep="")
inverse_pred_mess[[mv_it]] <- .check_inverse_pred(object=object, family=object$families[[mv_it]])
# : not elegant: assumes that .muetafn() passes the attrs checked by .check_inverse_pred (OK so far...),
# and that they correspond to a single family-link combination (OK so far...).
}
inverse_pred_mess <- .unlist(inverse_pred_mess) # NULL of "character" vector
cat(paste0(famst, collapse="; "))
cat("\n")
} else {
inverse_pred_mess <- .check_inverse_pred(object)
cat("family:", .prettify_family(object$family, linkstring = "link = ") , "\n")
}
summ$family <- object$family
if (length(object$fixef)==0L) {
cat("No fixed effect\n")
} else {
cat(" ------------ Fixed effects (beta) ------------\n")
beta_table <- .make_beta_table(object, p_value=details$p_value)
print(beta_table,4)
if (prod(dim(beta_table))>max.print && missing(max.print) ) {
cat(' [ Use the summary.HLfit() "max.print" argument\n to control getOption("max.print") at this point ].\n')
#: problem dealt here is that the print() message referring to getOption("max.print") suggests
# that the user can control this by setting the option globally, but this will fail 'locally'.
}
summ$beta_table <- beta_table
}
if (models[["eta"]]=="etaHGLM") {
cat(" --------------- Random effects ---------------\n")
urff <- unique(lcrandfamfam)
urffl <- unique(randfamfamlinks)
if (length(urffl)==1L) {
cat("Family:", .prettify_family(object$rand.families[[1]], linkstring = "link = ") , "\n")
} else {
cat("Families(links):", paste(randfamfamlinks,collapse=", "), "\n")
}
corrPars <- get_ranPars(object,which="corrPars")
corrPars <- corrPars[sort(names(corrPars))]
cP <- unlist(corrPars)
if ( ! is.null(cP) ) {
moreargs <- .get_moreargs(object)
control_dists <- lapply(moreargs,getElement,name="control.dist")
dist_methods <- lapply(control_dists,getElement,name="dist.method")
if (any((udm <- unlist(dist_methods))!="Euclidean")) {
ocd <- rep("Euclidean",length(dist_methods))
names(ocd) <- names(dist_methods)
ocd[names(udm)] <- udm
cat("Distance(s):",paste(ocd, collapse=","),"\n")
}
cat(" --- Correlation parameters:")
corrFixNames <- names(unlist(corrPars[which(attr(corrPars,"type")=="fix")]))
if (length(corrFixNames)>1) {
cat(" [",paste0(corrFixNames,collapse=",")," were fixed]")
} else if (length(corrFixNames)==1L) cat(" [",corrFixNames," was fixed]")
cat("\n")
print(cP)
}
lambda.object <- object$lambda.object
namesTerms <- lambda.object$print_namesTerms ## list of vectors of variable length
if (any(object$models[["lambda"]] == "lamHGLM")) {
stop("voir ici dans summary.HLfit")
} else {
linklam_coeff_list <- lambda.object$coefficients_lambdaS ## used beyond the next line
summ$lambda_table <- lambda_table <- .lambda_table_fn(namesTerms, object, lambda.object,linklam_coeff_list)
.print_lambda_table(object,lambda_table, details=details, namesTerms, linklam_coeff_list,
digits=details$digits)
}
groups_n <- unlist(lapply(object$ZAlist,ncol))/attr(object$ZAlist,"Xi_cols")
cat(paste0("# of obs: ",nrow(object$data),"; # of groups: ",
paste0(names(namesTerms),", ",groups_n, collapse="; ")
), "\n")
}
##
if (length(vec_nobs <- object$vec_nobs)) { # fitmv case; then object$phi.object must be a list of phi objects
if (any(.unlist(lapply(object$families,`[[`, "family")) %in% c("gaussian","Gamma")) ||
! all(sapply(object$prior.weights, attr, "is_unit")) ||
length(setdiff(object$models$rdispar, c("","rdiOff")))) {
cat("-------------- Residual variation -------------\n")
for (mv_it in seq_along(form)) {
summ <- .summary_phi_object(object, phi.object=object$phi.object[[mv_it]],
family=object$families[[mv_it]],
pw=object$prior.weights[[mv_it]],
summ=summ, phimodel=models[["phi"]][[mv_it]], mv_it=mv_it)
}
}
} else summ <- .summary_phi_object(object, summ=summ, phimodel=models[["phi"]])
if (object$HL[1]==0L) {
validnames <- intersect(names(object$APHLs),c("hlik","p_v","p_bv"))
} else validnames <- intersect(names(object$APHLs),c("p_v","p_bv"))
if (length(validnames)) { ## may be 0 in SEM...
likelihoods <- unlist(object$APHLs[validnames]) # NULL if no validnames
if ( models[["eta"]]=="etaHGLM"){
APHLlegend <- c(hlik=" h-likelihood:",
p_v= if (.get_obsInfo(object)) {
paste0("logL (P_v(h)):")
} else "logL (p_v(h)):",
p_bv= if (.get_obsInfo(object)) {
paste0("Re.logL (P_b,v(h)):")
} else "Re.logL (p_b,v(h)):")
} else APHLlegend <- c(p_v="logL :",
p_bv="log restricted-lik :")
names(likelihoods) <- APHLlegend[validnames]
if ( is.null(object$distinctX.Re)) {
## standard REML
} else {
whichp_bv <- which(validnames=="p_bv")
if (ncol(object$distinctX.Re)==0L) {
likelihoods <- likelihoods[-whichp_bv] ## ML
} else names(likelihoods)[whichp_bv] <- "(Non-standard?) ReL:"
}
} else likelihoods <- numeric(0)
# messlist[["ranef"]]
logLapp <- object$APHLs$logLapp
if (!is.null(logLapp)) {
locli <- list(logLapp[1]) ## [1] removes attribute
names(locli)[1] <- attr(logLapp,"method")
likelihoods <- c(likelihoods,locli)
}
cat(" ------------- Likelihood values -------------\n")
astable <- as.matrix(likelihoods);colnames(astable)<-"logLik";
print(astable)
summ$likelihoods <- likelihoods
if (length(object$warnings) ) silent <- sapply(object$warnings, cat, "\n")
cat(paste(inverse_pred_mess, collapse="\n"))
options(oldopt)
invisible(summ)
}
print.HLfit <- function(x,...) {
summary(x,...)
invisible(x)
}
print.HLfitlist <- function(x,...) {
summary(x,...)
invisible(x)
}
div_info <- function(object, ...) {
ranges <- NULL # avoid bug when no divinfo
if ( ! is.null(ranFixes <- object$divinfo$high_kappa$ranFixes)) { # remains NULL if problem too large for fast condnum computation
corr_info <- .get_from_ranef_info(object)
ranFixes <- lapply(ranFixes, .canonizeRanPars, corr_info=corr_info,checkComplete = FALSE, rC_transf = .spaMM.data$options$rC_transf)
if (length(ranFixes)>1L) {
ranges <- t(apply(do.call(rbind,lapply(ranFixes,unlist)),2L,range))
if (is.null(ranges)) {
message("Near-singular Hessian matrix, but could not be related to specific parameter values.")
} else {
message("Near-singular Hessian matrix for the following ranges of parameters:")
colnames(ranges) <- c("lower","upper")
print(ranges)
message("It may be worth refitting in a constrained parameter space excluding some of these ranges.")
}
} else if (length(ranFixes)) {
message("Near-singular Hessian matrix for the following parameters:")
print(unlist(ranFixes))
message("It may be worth refitting in a constrained parameter space excluding these parameters.")
}
}
if ( ! is.null(ranFixes <- object$divinfo$unknown$ranFixes)) { # if processed$envir$ranFix was NULL there are explicit list(NULL) here
corr_info <- .get_from_ranef_info(object)
ranFixes <- lapply(ranFixes, .canonizeRanPars, corr_info=corr_info,checkComplete = FALSE, rC_transf = .spaMM.data$options$rC_transf)
if (length(ranFixes)>1L) {
ranges <- t(apply(do.call(rbind,lapply(ranFixes,unlist)),2L,range))
if (is.null(ranges)) {
message("A non-convergence issue was detected but could not be diagnosed\n nor related to specific parameter values.")
} else {
message("Uncharacterized numerical issue for the following ranges of parameters:")
colnames(ranges) <- c("lower","upper")
print(ranges)
message("It may be worth refitting in a constrained parameter space excluding some of these ranges.")
}
} else if (length(ranFixes)) {
message("Uncharacterized numerical issue for the following parameters:")
print(unlist(ranFixes))
message("It may be worth refitting in a constrained parameter space excluding these parameters.")
}
}
invisible(ranges)
}
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.