Nothing
.is_2_in_1 <- function(X1,X2, tol=1e-10) {
qrX <- qr(X1)
proj2in1 <- X1 %*% solve(qrX,X2) # is.qr.solve()
return(all(abs(proj2in1 -X2)<tol))
}
# deprecated:
.process_ranef_case <- function(object, object2, nest="") {
l1 <- logLik(object)
l2 <- logLik(object2)
testlik <- unique(c(names(l1),names(l2)))
if (length(testlik)!=1L) stop(paste("Unable to determine a unique objective function from logLik() names.\n",
"Check that both fits are fitted by the same 'method'."))
if (nest=="2in1" || (nest!="1in2" && l1>l2)) {
nullfit <- object2
fullfit <- object
} else { # nest=="1in2" || l2>l1
nullfit <- object
fullfit <- object2
}
return(list(fullfit=fullfit,nullfit=nullfit,test_obj=testlik,df=NA))
}
.guess_Rnest <- function(object, object2, Xnest) {
dfR1 <- sum(unlist(object$dfs, use.names=FALSE)) - object$dfs$pforpv
dfR2 <- sum(unlist(object2$dfs, use.names=FALSE)) - object2$dfs$pforpv
if (dfR1>dfR2) {
Rnest <- "2in1"
if (is.null(Xnest)) { # fixed effects are identical
message(paste("spaMM did not ascertain whether the models are nested in their random-effect specifications.\n",
" You are responsible for that."))
# comparison of random-coefficients, (LHS|RHS), vs (1|RHS) leads here and it's a bit sad that spaMM does not diagnose that, but... not now.
# (same in 1in2 case)
} else warning(paste("Fixed-effects specifications differ between the models, and spaMM did not ascertain\n",
" whether the models are similarly nested in their random-effect specifications.\n",
" You are responsible for that. spaMM will guess nestedness from number of parameters."), immediate. = TRUE)
} else if (dfR1<dfR2) {
Rnest <- "1in2"
if (is.null(Xnest)) { # fixed effects are identical
message(paste("spaMM did not ascertain whether the models are nested in their random-effect specifications.\n",
" You are responsible for that."))
} else warning(paste("Fixed-effects specifications differ between the models, and spaMM did not ascertain\n",
" whether the models are similarly nested in their random-effect specifications.\n",
" You are responsible for that. spaMM will guess nestedness from number of parameters."), immediate. = TRUE)
} else if (is.null(Xnest)) { # fixed effects are identical, dfs are identical
stop("The models compared appear to have the same number of parameters, hence may be identical or non-nested.",
call.=FALSE)
} else {
Rnest <- NULL
warning(paste("The models have the same number of parameters except for fixed effects, but spaMM did not ascertain\n",
" whether the random-effect specifications are identical. You are responsible for that."), immediate. = TRUE)
}
list(Rnest=Rnest, dfs=c(dfR1, dfR2))
}
.compare_phimodels_structures <- function(phifit1, phifit2, df1, df2, df=abs(df1-df2)) {
if (inherits(phifit1,"HLfit")) {
if (inherits(phifit2,"HLfit")) {
cmp <- .compare_model_structures(phifit1,phifit2, args_are_phimodels=TRUE)
if (cmp$anyranef) return(list(warning="Mixed-effect residual dispersion model found. No chi2 LRT performed."))
# df <- cmp$df
# df1 <- cmp$df1
# df2 <- cmp$df2
## boostrap tests still possible
} else {
# # phifit2 is a phiScal... the input df2 may be 0 or 1
# the input df1 may be object$dfs$p_fixef_phi[[mv_it]] which is then NA
if (phifit1$models$eta=="etaHGLM") {
return(list(warning="Mixed-effect residual dispersion model found. No chi2 LRT performed."))
} else if (df2==1L && "(Intercept)" %in% names(fixef(phifit1))) {
df1 <- sum(.unlist(phifit1$dfs))
} else df1 <- NA
if (length(attr(phifit1$ZAlist,"exp_ranef_strings"))) df <- NA
}
} else { # phifit1 is a phiScal...
if (inherits(phifit2,"HLfit")) {
if (phifit2$models$eta=="etaHGLM") {
return(list(warning="Mixed-effect residual dispersion model found. No chi2 LRT performed."))
} else if (df1==1L && "(Intercept)" %in% names(fixef(phifit2))) {
df2 <- sum(.unlist(phifit2$dfs))
} else df2 <- NA
if (length(attr(phifit2$ZAlist,"exp_ranef_strings"))) df <- NA
}
}
return(c(df1,df2,df))
}
.is_nested_family <- function(famfam, famfam2, object2, mv_it=NULL) {
(
(famfam=="poisson" &&
famfam2 %in% c("negbin1","negbin2","COMPoisson") &&
length(intersect(unlist(.get_methods_disp(object2, mv_it=mv_it)), c("NB_shape","COMP_nu","beta_prec"))) # detection fitted par as in .calc_dfs
)
) || (
famfam=="binomial" &&
famfam2=="betabin" &&
length(intersect(unlist(.get_methods_disp(object2, mv_it=mv_it)), c("beta_prec")))
)
}
# Handles mv fits, their resid models included.
.compare_model_structures <- function(object,object2, args_are_phimodels, boot.repl) {
if (inherits(object,"HLfitlist") || inherits(object2,"HLfitlist")) {
warning("This does not yet work on HLfitlist objects", immediate.=TRUE, call.=FALSE)
}
X.pv1 <- model.matrix(object)
X.pv2 <- model.matrix(object2)
X1 <- attr(X.pv1,"namesOri") ## need to track NA beta's
X2 <- attr(X.pv2,"namesOri")
if (length(X1)==0L) {
REML1 <- NULL ## compatible with both ML or REML tests
} else REML1 <- (object$APHLs$p_v != object$APHLs$p_bv)
if (length(X2)==0L) {
REML2 <- NULL ## idem
} else REML2 <- (object2$APHLs$p_v != object2$APHLs$p_bv)
REML <- unique(c(REML1,REML2))
Fnest <- NULL
if (! args_are_phimodels) {
meth1 <- object$HL
meth2 <- object2$HL
if (! identical(meth1,meth2) || length(REML)>1L ) {
stop("object fitted by different methods cannot be compared.", call.=FALSE)
}
if (mv_model <- is.null(object$family)) {
nestedfam <- all(sapply(object$families, `[[`, x="link")==sapply(object2$families, `[[`, x="link"))
} else nestedfam <- (object$family$link==object2$family$link)
if ( ! nestedfam) stop("Models not nested (distinct family links).", call.=FALSE)
if (mv_model) {
famfams1 <- sapply(object$families, `[[`, x="family")
famfams2 <- sapply(object2$families, `[[`, x="family")
if ( any(famfams1!=famfams2)) {
for (mv_it in seq_along(famfams1)) {
if (famfams1[mv_it]!= famfams2[mv_it]) {
if (.is_nested_family(famfams1[mv_it], famfams2[mv_it], object2, mv_it=mv_it)) {
Fnest <- unique(c(Fnest, "1in2"))
} else if (.is_nested_family(famfam2, famfam, object, mv_it=mv_it)) {
Fnest <- unique(c(Fnest, "2in1"))
}
if (length(Fnest)>1L) stop("Models may not be nested (distinct families).", call.=FALSE)
}
}
}
} else if (object$family$family!=object2$family$family) {
famfam <- object$family$family
famfam2 <- object2$family$family
if (.is_nested_family(famfam, famfam2, object2)) {
Fnest <- "1in2"
} else if (.is_nested_family(famfam2, famfam, object)) {
Fnest <- "2in1"
} else stop("Models may not be nested (distinct families).", call.=FALSE)
}
}
if ( ! is.null(X1)) X1 <- sapply(strsplit(X1,':'), function(x) paste(sort(x),collapse=':')) ## JBF 2015/02/23: sort variables in interaction terms before comparison
if ( ! is.null(X2)) X2 <- sapply(strsplit(X2,':'), function(x) paste(sort(x),collapse=':'))
dX12 <- setdiff(X1,X2)
dX21 <- setdiff(X2,X1)
if (length(dX12) && length(dX21)) {
if (.is_2_in_1(X1=X.pv1, X2=X.pv2)) {
Xnest <- "2in1"
} else if (.is_2_in_1(X1=X.pv2, X2=X.pv1)) {
Xnest <- "1in2"
} else if (args_are_phimodels) {
stop("Fixed effects seem non-nested for residual-dispersion model.", call.=FALSE)
} else stop("Fixed effects seem non-nested.", call.=FALSE)
} else if (length(dX12)) {
Xnest <- "2in1"
} else if (length(dX21)) {
Xnest <- "1in2"
} else {
message("Fixed effects terms appear identical for both models.")
Xnest <- NULL
}
#
if (object$spaMM.version < "2.2.116") {
ranterms1 <- attr(object$ZAlist,"ranefs")
} else ranterms1 <- attr(object$ZAlist,"exp_ranef_strings")
if (object2$spaMM.version < "2.2.116") {
ranterms2 <- attr(object2$ZAlist,"ranefs")
} else ranterms2 <- attr(object2$ZAlist,"exp_ranef_strings")
randist1 <- lapply(object$rand.families, function(v) paste(paste(v)[1:2],collapse="")) ## makes a string from each $family and $link
randist2 <- lapply(object2$rand.families, function(v) paste(paste(v)[1:2],collapse="")) ## makes a string from each $family and $link
if (args_are_phimodels && # nested call compare_model_str -> compare_phimodel_str -> compare_model_str
(length(randist1) || length(randist2))) {# either of the two residual-disp models is a mixed-effect model
# => no objective function maximized => no LRT, not even a bootstrap one
# df <- NA # inhibits asymptotic LRT
return(list(anyranef=TRUE))
}
#
ranterms1 <- paste(ranterms1,randist1) ## joins each term and its distrib
ranterms2 <- paste(ranterms2,randist2) ## joins each term and its distrib
dR12 <- setdiff(ranterms1,ranterms2)
dR21 <- setdiff(ranterms2,ranterms1)
if (length(dR12) && length(dR21)) { # no obvious nested structure on ranefs. Trying to guess 'Rnest'
Rnest_blob <- .guess_Rnest(object, object2, Xnest)
Rnest <- Rnest_blob$Rnest
#.process_ranef_case(object, object2, Xnest=Xnest)
} else if (length(dR12)) {
Rnest <- "2in1"
} else if (length(dR21)) {
Rnest <- "1in2"
} else Rnest <- NULL
XRnest <- c(Xnest,Rnest)
uXRnest <- unique(XRnest)
if (length(unique(c(Fnest,uXRnest)))>1L) {
# if (args_are_phimodels) {
# warning("Models not nested (opposite nestings for fixed and random terms in residual-dispersion model). No test performed.", call.=FALSE)
# } else
warning("Models not nested (opposite nestings for fixed and random terms\n or response families). No test performed.", call.=FALSE)
return(list(fullfit=NULL,nullfit=NULL,test_obj=NULL,df=NA, df1=NA, df2=NA))
}
# ELSE
if ( ! args_are_phimodels) { #check their residual-dispersion models ##################################### ( length(unest)==0L) { # no difference between models found yet. Check residual dispersion models
phimodel1 <- object$models[["phi"]]
phimodel2 <- object2$models[["phi"]]
n_mv <- length(phimodel1)
df_mat <- matrix(NA, nrow=3L, ncol=n_mv)
if (n_mv>1L) {
for (mv_it in seq_len(n_mv)) {
df_col <- .compare_phimodels_structures(phifit1 = residVar(object, which="fit", submodel = mv_it),
phifit2 = residVar(object2, which="fit", submodel = mv_it),
df1=object$dfs$p_fixef_phi[[mv_it]],
df2=object2$dfs$p_fixef_phi[[mv_it]])
if (is.list(df_col)) {
warning(df_col$warning, call.=FALSE)
return(list(fullfit=NULL,nullfit=NULL,test_obj=NULL,df=NA,
warning=df_col$warning))
} else if (anyNA(df_col[1:2])) warning(paste("apparently non-nested residual-dispersion models for submodel",mv_it), call.=FALSE)
df_mat[,mv_it] <- df_col
}
} else {
df_col <- .compare_phimodels_structures(phifit1 = residVar(object, which="fit", submodel = NULL),
phifit2 = residVar(object2, which="fit", submodel = NULL),
df1=object$dfs$p_fixef_phi,
df2=object2$dfs$p_fixef_phi)
if (is.list(df_col)) {
warning(df_col$warning, call.=FALSE)
return(list(fullfit=NULL,nullfit=NULL,test_obj=NULL,df=NA))
} else df_mat[,1L] <- df_col
}
if (anyNA(df_mat[1:2,])) {
warning("apparently non-nested residual-dispersion models", call.=FALSE)
return(list(fullfit=NULL,nullfit=NULL,test_obj=NULL,df=NA))
} else if (all(df_mat[1,]==df_mat[2,])) {
Pnest <- NULL
} else if (all(df_mat[1,]>=df_mat[2,])) {
Pnest <- "2in1"
} else if (all(df_mat[1,]<=df_mat[2,])) {
Pnest <- "1in2"
} else {
warning(paste("The two models appear non-nested (opposite nesting of residual dispersion models?).\n No test performed."), call.=FALSE)
return(list(fullfit=NULL,nullfit=NULL,test_obj=testlik,df=NA))
}
XRPnest <- c(XRnest,Pnest)
unest <- unique(XRPnest)
if (length(unest)==2L) {
warning("Models not nested (opposite nestings for mean-response and residual-dispersion models).\n No test performed.", call.=FALSE)
return(list(fullfit=NULL,nullfit=NULL,test_obj=NULL,df=NA))
} else if (length(unest)==0L) {
warning("The two models appear equivalent (except perhaps for offsets). No test performed.", call.=FALSE)
return(list(fullfit=NULL,nullfit=NULL,test_obj=NULL,df=0))
} else if (length(uXRnest)==0L) {
message(paste("The two models appear equivalent except for residual dispersion models."))
}
# overwrite the p_fixef_phi (with ma further be a list for mv fits) with a single scalar for all params of phi models
dfs1 <- object$dfs
dfs1$p_fixef_phi <- NULL
dfs1$p_phi <- .calc_p_rdisp(object)
dfs2 <- object2$dfs
dfs2$p_fixef_phi <- NULL
dfs2$p_phi <- .calc_p_rdisp(object2)
df1 <- sum(unlist(dfs1, use.names=FALSE))
df2 <- sum(unlist(dfs2, use.names=FALSE))
df <- abs(df1-df2)
if (unest=="1in2") {
fullm <- object2
nullm <- object
} else {
fullm <- object
nullm <- object2
}
if (length(XRnest)==2L) {
message("Models differing both by in their fixed and in their random terms. ")
if (REML) warning("LRT comparing REML fits with different fixed-effect conditions is highly suspect",
immediate.=TRUE, call.=FALSE)
testlik <- "p_v"
} else if (length(uXRnest)==0L) { # that means that the difference is in Pnest
if (REML) {
testlik <- "p_bv"
} else {
message("Note: ML fits used to compare different residual-dispersion models.")
testlik <- "p_v"
}
} else {
if (is.null(Rnest) && is.null(Pnest)) { ## fixed effect test
if (REML) {
## checking the comparability of REML fits (must be an hack for non-standard REML)
if ( ! is.null(fullm$distinctX.Re) ) {
df.f.Re <- ncol(fullm$distinctX.Re)
} else df.f.Re <- ncol(fullm$`X.pv`)
if ( ! is.null(nullm$distinctX.Re) ) {
df.n.Re <- ncol(nullm$distinctX.Re)
} else df.n.Re <- ncol(nullm$`X.pv`)
if ( df.f.Re != df.n.Re ) {
warning("LRT comparing REML fits with different fixed-effect conditions is highly suspect",
immediate.=TRUE, call.=FALSE)
}
}
testlik <- "p_v"
} else {
if ( ! is.null(Rnest)) { ## random effect test (case of models differing only through resid.model was captured above )
if (length(dR12) && length(dR21) &&
# comparing two random-effect models, possibly differing in other components but not in a clearly conflicting way
min(Rnest_blob$dfs)>0L) {
warnmess <- "Asymptotic chi-square tests comparing different\n random-effect specifications may be poorly behaved."
if (boot.repl) {
message(paste("Note:",warnmess))
} else warning(warnmess, immediate. = TRUE, call.=FALSE)
df <- abs(df1-df2)
RLRbool <- FALSE
} else {
# detect comparison of LM to LMM with single random effect, and no heteroscedasticity:
RLRbool <- (unique(fullm$models$phi)=="phiScal" &&
fullm$family$family =="gaussian" &&
(fullm$models[["eta"]]=="etaHGLM" &&
all(attr(fullm$rand.families,"lcrandfamfam")=="gaussian")) &&
nullm$models[["eta"]]=="etaGLM")
df <- NA # inhibits asymptotic LRT
}
} else { # no longer sure it can be reached
df <- NA # inhibits asymptotic LRT
RLRbool <- FALSE
}
if (REML) {
testlik <- "p_bv"
if (is.null(Xnest) && RLRbool) message("Note: procedures from package 'RLRsim' may be useful for this test. See help('get_RLRsim_args').")
} else {
testlik <- "p_v"
if (is.null(Xnest)) message("Note: ML fits used to compare different random-effects or residual-dispersion models.")
if (RLRbool) message("Note: procedures from package 'RLRsim' may be useful for this test. See help('get_RLRsim_args').")
}
}
}
return(list(fullfit=fullm,nullfit=nullm,test_obj=testlik,df=df,
df1=df1,df2=df2))
} else { # args_are_phimodels: (fixed-effects only, mixed-effect models caught early above)
df1 <- sum(unlist(object$dfs, use.names=FALSE))
df2 <- sum(unlist(object2$dfs, use.names=FALSE))
df <- abs(df1-df2)
return(list(df1=df1,df2=df2, df=df))
}
}
.add_famPars_outer <- function(parlist, fitobject, names_u_c_inits=NULL, type_attr) {
if (is.null(names_u_c_inits)) {
canon.init <- attr(fitobject,"optimInfo")$LUarglist$canon.init ## includes user init
names_u_c_inits <- names(unlist(canon.init))
}
# For each of the following family param there may already be attr(outer_ests,"type")$tr<fam par> but let's not assume that messy thing
if ( ! is.null(families <- fitobject$families)) {
for (mv_it in seq_along(families)) {
fam_it <- families[[mv_it]]
char_mv_it <- as.character(mv_it)
if ((parname <- paste("NB_shape",mv_it,sep=".")) %in% names_u_c_inits) {
parlist[["NB_shape"]][char_mv_it] <- environment(fam_it$aic)$shape # <vector element> <-
if (type_attr) attr(parlist,"type")[["NB_shape"]][char_mv_it] <- "outer"
} else if ((parname <- paste("COMP_nu",mv_it,sep=".")) %in% names_u_c_inits) {
parlist[["COMP_nu"]][char_mv_it] <- environment(fam_it$aic)$nu # <vector element> <-
if (type_attr) attr(parlist,"type")[["COMP_nu"]][char_mv_it] <- "outer"
} else if ((parname <- paste("beta_prec",mv_it,sep=".")) %in% names_u_c_inits) {
parlist[["beta_prec"]][char_mv_it] <- environment(fam_it$aic)$prec # <vector element> <-
if (type_attr) attr(parlist,"type")[["beta_prec"]][char_mv_it] <- "outer"
} else if (! is.null(vec <- canon.init$rdisPars[[char_mv_it]])) {
fitted_beta <- fam_it$resid.model$beta[names(vec)]
parlist[["rdisPars"]][char_mv_it] <- list(fitted_beta)
if (type_attr) {
typ <- rep("outer", length(fitted_beta))
names(typ) <- names(vec)
attr(parlist,"type")[["rdisPars"]][char_mv_it] <- typ
}
}
}
return(parlist)
}
##
if ("NB_shape" %in% names_u_c_inits) {
parlist$NB_shape <- environment(fitobject$family$aic)$shape
if (type_attr) attr(parlist,"type")$NB_shape <- "outer"
}
if ("COMP_nu" %in% names_u_c_inits) {
parlist$COMP_nu <- environment(fitobject$family$aic)$nu
if (type_attr) attr(parlist,"type")$COMP_nu <- "outer"
}
if ("beta_prec" %in% names_u_c_inits) {
parlist$beta_prec <- environment(fitobject$family$aic)$prec
if (type_attr) attr(parlist,"type")$beta_prec <- "outer"
}
if (! is.null(vec <- canon.init[["rdisPars"]])) {
fitted_beta <- fitobject$family$resid.model$beta[names(vec)]
parlist["rdisPars"] <- list(fitted_beta)
if (type_attr) {
typ <- rep("outer", length(fitted_beta))
names(typ) <- names(vec)
attr(parlist,"type")$rdisPars <- typ
}
}
parlist
}
# Construct new outer inits from outer fitted values and/or initial outer values of input fit.
.get_outer_inits_from_fit <- function(fitobject, keep_canon_user_inits) {
canon.init <- attr(fitobject,"optimInfo")$LUarglist$canon.init ## includes sanitized user init
if (FALSE) {
outer_ests <- get_ranPars(fitobject,lambda_names = "") ## "CorrEst_and_RanFix only"
# which means that only these params were conceived to be controlled in the new outer inits.
# But get_ranPars(, which=NULL) is not formally defined, so has been changing... the original comment is no longer true.
## If the alternative not valid in the long run, this get_ranPars(.) calls should be modified.
} else {
optimInfo <- attr(fitobject,"optimInfo")
outer_ests <- optimInfo$optim.pars # transparent (but potentially more comprehensive set of params)
if ( ! is.null(outer_ests)) {
attr(outer_ests,"moreargs") <- optimInfo$LUarglist$moreargs # necess to canonize Matern params...
outer_ests <- .canonizeRanPars(outer_ests, corr_info=fitobject$ranef_info$sub_corr_info,
checkComplete=FALSE, rC_transf=.spaMM.data$options$rC_transf)
}
couter_est <- c(outer_ests) # drops fancy attributes
}
if (keep_canon_user_inits &&
length(user_inits <- .reformat_corrPars(getCall(fitobject)$init,corr_families=fitobject$corr_info$corr_families))
) { # keep them (as interpreted in canon.init: minimum phi is 1e-4, etc) in return value
# => remove the fitted values from the nullranPars used to modify_list
# => keep them in 'removand' list of pars to remove from nullranPars
# => exclude them from 'not_user_inits_names' to remove from 'removand' !
in_init_not_user_inits <- setdiff(names(unlist(canon.init)), names(unlist(user_inits))) # names, excluding those of parameters with user inits
in_user_inits <- setdiff(names(unlist(outer_ests)), in_init_not_user_inits) ## it rem
## removand: user_inits, fixed, or inner optimized corrPars
# locinit will retain parameters that were outer optimized without an explicit user init
if ( is.null(in_user_inits)) {
locinit <- outer_ests
} else { # replace initial value by [fitted values, except those that had a user_init]
# => the locinit retains the sanitized user init
not_in_user_inits <- .remove_from_cP(outer_ests, u_names=in_user_inits)
locinit <- .modify_list(canon.init, not_in_user_inits) ## loses attributes
}
} else locinit <- outer_ests
return(locinit)
}
.get_rC_inits_from_hlfit <- function(hlfit,
type # set type=NULL to retain all types
) {
reinit <- .get_compact_cov_mats(hlfit$strucList)
seq_rd <- seq_along(reinit)
if (length(seq_rd)) {
for (rd in seq_rd) {
if ( ! is.null(ini_mix <- reinit[[rd]])) {
ini_corr <- cov2cor(ini_mix)
lowerblocF <- lower.tri(ini_corr,diag=FALSE)
ini_mix[lowerblocF] <- ini_corr[lowerblocF] # replaces covby corr
reinit[[rd]] <- ini_mix[lower.tri(ini_mix,diag=TRUE)] ## mix cov/corr in vector form
} else reinit[rd] <- NA # the syntax understood by fitting functions
}
names(reinit) <- paste(seq_rd)
if ( ! is.null(type)) {
reinit[ ! hlfit$lambda.object$type %in% type] <- NA
reinit <- reinit[ ! is.na(reinit)] # otherwise for "outer" type the NA ends in the optimizer's init...
}
}
return(reinit)
}
# for inner CAR it returns the single lambda factor which is OK
.get_lambdas_notrC_from_hlfit <- function(hlfit, type, keep_names=(type=="adhoc"),
isRandomSlope=attr(hlfit$strucList,"isRandomSlope")) {
lambdas <- hlfit$lambda.object$lambda_list
seq_rd <- seq_along(lambdas)
if (length(seq_rd)) {
any_ranCoef <- FALSE
for (rd in seq_rd) {
if (isRandomSlope[[rd]]) {
lambdas[[rd]] <- NA # the syntax understood by fitting functions
any_ranCoef <- TRUE
} else {
# pretty renaming of simple lambda:
if (type=="adhoc") if (length(lambdas[[rd]])==1L && names(lambdas[[rd]])=="(Intercept)") names(lambdas[[rd]]) <- NULL
}
}
lambdas <- unlist(lambdas)
if ( ! keep_names) names(lambdas) <- paste(seq_rd)
if (type !="adhoc") lambdas[hlfit$lambda.object$type!=type] <- NA
lambdas <- lambdas[ ! is.na(lambdas)] # otherwise for "outer" type the NA ends in the optimizer's init...
# for "inner" NA's are harmless but not necessary when names are paste(seq_rd)
if (any_ranCoef && type=="adhoc") lambdas <- structure(lambdas,
message="Random-coefficient variances removed. Use e.g. VarCorr() to get them.")
}
return(lambdas) # vector wwith NAs for everything not wanted
}
# to initiate a *f*ullfit # , inner_lambdas= TRUE
get_inits_from_fit <- function(from, template=NULL, to_fn=NULL, inner_lambdas=FALSE) { # 'to_fn' may differ from that of 'from' and 'to'
new_outer_inits <- .get_outer_inits_from_fit(fitobject=from, keep_canon_user_inits = FALSE)
# check fromfn and to_fn
if (is.null(to_fn)) {
# recent objects should have how$fnname. Otherwise, .get_bare_fnname() may not return a valid name.
fromfn <- .get_bare_fnname.HLfit(from)
if (is.null(template)) {
to_fn <- fromfn
} else if (inherits(template,"HLfit")) {
to_fn <- .get_bare_fnname.HLfit(template)
} else stop("Invalid 'template' argument.")
} else if (to_fn=="fitme_body") { ## using to_fn to modify fromfn...
## ad hoc fix for residModel: fitme_body is called directly so the final object's call is to HLCor of HLfit
fromfn <- "fitme"
} else fnname <- .get_bare_fnname.HLfit(from)
# Inner-estimated lambda and ranCoefs (_F I X M E__ could add phi: amusing has this was never done... inner estimated mv phi exist, incidentally)
init.HLfit <- NULL
rC_inner_inits <- .get_rC_inits_from_hlfit(from, type="inner") # (yes, inner, not inner_ranCoefs)
if (length(rC_inner_inits) ) init.HLfit <- list(ranCoefs=rC_inner_inits)
if (inner_lambdas) {
lambda_inner_inits <- .get_lambdas_notrC_from_hlfit(from, type="inner")
if (length(lambda_inner_inits)) init.HLfit <- c(init.HLfit, list(lambda=lambda_inner_inits))
}
#
if (to_fn %in% c("fitme", "fitmv", "fitme_body")) {
new_inits <- list(init=new_outer_inits,init.HLfit=init.HLfit)
} else if (to_fn=="corrHLfit") {
new_inits <- list(init.corrHLfit=new_outer_inits,init.HLfit=init.HLfit)
} else new_inits <- list(init.HLfit=init.HLfit)
# Add initial value for fixed effects
if (length(fixef_from <- na.omit(fixef(from)))) {
# fixef() returns a vector with NA's ifor non-estimable parameters;
# but these NA should not reach the code initializing eta as X.beta in HLfit_body, using beta from the inits... => na.omit
if (inherits(template,"HLfit")) { # This was motivated by the Leucadendron_hard.R bootstrap replicates.
new_HLfit_inits <- na.omit(fixef(template))
new_HLfit_inits[names(new_HLfit_inits)] <- 0
shared_names <- intersect(names(new_HLfit_inits), names(fixef_from))
new_HLfit_inits[shared_names] <- fixef_from[shared_names]
new_HLfit_inits <- list(fixef=new_HLfit_inits)
new_inits[["init.HLfit"]] <- c(new_inits[["init.HLfit"]], new_HLfit_inits)
} else new_inits[["init.HLfit"]] <- c(new_inits[["init.HLfit"]], list(fixef=fixef_from))
}
return(new_inits)
}
.update_control <- function(fit_call, optim_boot, from_fn=NULL) {
if (is.null(from_fn)) from_fn <- paste(fit_call[[1]])
if (from_fn=="fitme") {
ctrl_opt <- fit_call[["control"]]
if (is.null(ctrl_opt)) {
ctrl_opt <- list(optimizer=optim_boot)
} else ctrl_opt[["optimizer"]] <- optim_boot
} else if (from_fn=="corrHLfit") {
ctrl_opt <- fit_call[["control.corrHLfit"]]
if (is.null(ctrl_opt)) {
ctrl_opt <- list(optimizer=optim_boot)
} else ctrl_opt[["optimizer"]] <- optim_boot
} else ctrl_opt <- NULL
return(ctrl_opt)
}
eval_replicate <- function(y) { # no additional arguments, to ease parallel programming => next lines instead
# the function will be called within e.g. pbapply so it's useless to refer to parent.frame() here
enclosing_env <- parent.env(environment()) ## this is not necess for the code to run, but for the CRAN checks not to complain
nullfit <- get("nullfit", enclosing_env)
fullfit <- get("fullfit", enclosing_env)
#dotargs <- get("dotargs", enclosing_env)
test_obj <- get("test_obj", enclosing_env)
debug. <- get("debug.", enclosing_env)
# .condition <- get(".condition", enclosing_env)
null_call <- getCall(nullfit)
verbose <- null_call$verbose
verbose["phifit"] <- FALSE
null_fit_fn <- .get_bare_fnname.HLfit(nullfit, call.=null_call)
full_fit_fn <- .get_bare_fnname.HLfit(fullfit)
newinits <- .get_outer_inits_from_fit(fitobject=nullfit, keep_canon_user_inits = FALSE) # for next new_nullfit
ctrl_opt <- .update_control(fit_call=null_call, optim_boot=.spaMM.data$options$optim_boot, from_fn=null_fit_fn) # need .safe_opt when newinits are at bound.
if (null_fit_fn=="fitme") {
new_args <- list(init=newinits, control=ctrl_opt, verbose=verbose)
} else if (null_fit_fn=="corrHLfit") {
new_args <- list(init.corrHLfit=newinits, control.corrHLfit=ctrl_opt, verbose=verbose)
} else new_args <- list(verbose=verbose)
# pbbly never good not to use try(). It's the handling of try-error that may vary
# debug(pbapply) (or pblapply ?) may be useful to
if (debug.==2L) {# Shuld be prevented by spaMM's calling fns in a parallel session
re_nullfit <- do.call(update_resp, c(list(object=nullfit, newresp = y),new_args))
} else {
re_nullfit <- try(do.call(update_resp, c(list(object=nullfit, newresp = y),new_args)))
if (inherits(re_nullfit,"try-error")) { ## (debug.= TRUE or 1L) to return error info in parallel mode: return the try-error object
if (debug.) {
utils::dump.frames(dumpto="dump_on_re_nullfit", to.file=TRUE) # but doesn't stop
# .check_bootreps() checks the first attr of the first element of the returned list,
# looking for a try-error => 'null' should be thefirst element here.
return(list(null=structure(NA,re_nullfit=re_nullfit),
full=structure(NA,info="no fullfit")))
} else return(c(full=NA, null=NA)) # attributes would be lost at the level of the pbapply() closure.
# cf apply -> array -> as.vector -> loses all attributes as doc'ed in apply and as.vector
} ## ELSE continue
}
logL_re_null <- logLik(re_nullfit,which=test_obj)
# Allows the user to control the starting values of the re_fullfit:
# edotargs <- dotargs
# for (st in setdiff(names(edotargs),"prior.weights")) {
# edotargs[[st]] <- eval(edotargs[[st]],env=environment()) ## evaluate the promises in the current execution envir
# }
# #
newinits <- get_inits_from_fit(from=re_nullfit, template=fullfit, to_fn=full_fit_fn )
lens <- rep(NA, length(newinits))
for (it in seq_along(newinits)) lens[it] <- length(newinits[[it]])
newinits <- newinits[lens>0L]
subsets_inits <- .all.subsets(newinits) ## limited scope: could effectively do this is a nested way (with a skeleton+relist technique)
# subsets_inits is always a list of lists, with at least one element (minimally: list(list()))
for (it in seq_along(subsets_inits)) {
inits <- subsets_inits[[it]]
if (full_fit_fn=="fitme") {
new_args <- list(init=inits[["init"]], init.HLfit=inits[["init.HLfit"]], control=ctrl_opt, verbose=verbose)
} else if (full_fit_fn=="corrHLfit") {
new_args <- list(init.corrHLfit=inits[["init.corrHLfit"]], init.HLfit=inits[["init.HLfit"]],
control.corrHLfit=ctrl_opt, verbose=verbose)
} else new_args <- list(init.HLfit=inits[["init.HLfit"]], verbose=verbose)
if (debug.==2) {
re_fullfit <- do.call(update_resp, c(list(object=fullfit, newresp = y),new_args)) # may stop on error
} else {
re_fullfit <- try(do.call(update_resp, c(list(object=fullfit, newresp = y),new_args)))
if (inherits(re_fullfit,"try-error")) { ## (debug.= TRUE or 1L) to return error info in parallel mode: return the try-error object
if (debug.) {
utils::dump.frames(dumpto="dump_on_re_fullfit", to.file=TRUE) # but doesn't stop
return(list(full=structure(NA,re_fullfit=re_fullfit), # but the attributes seem lost at the level of the pbapply() closure.
null=structure(logL_re_null,re_nullfit=re_nullfit)))
} else return(c(full=NA, null=logL_re_null)) # attributes would be lost at the level of the pbapply() closure.
} ## ELSE continue
}
logL_re_full <- logLik(re_fullfit,which=test_obj)
if (logL_re_full<logL_re_null && ## re_fullfit poor but ...
(logL_re_full>logL_re_null-5e-02 && .check_outer_at_bound(re_fullfit)) # ...'not too poor' if at bound => automatically correct
) logL_re_full <- logL_re_null
if (logL_re_full>logL_re_null-1e-04) break # break the for loop.
}
# if the fullfit was even poorer than the 5e-02 threshold the LR stat is still negative. Correct it if outer_at_bound (___F I X M E____ warning?):
if (logL_re_full<logL_re_null &&
.check_outer_at_bound(re_fullfit)) logL_re_full <- logL_re_null
resu <- c(full=logL_re_full,null=logL_re_null)
# if ( ! is.null(.condition)) {
# condition <- with(list(nullfit=re_nullfit, fullfit=re_fullfit), eval(.condition))
# resu <- c(resu, condition=condition)
# }
return(resu)
}
.eval_replicate2 <- function(y) {
enclosing_env <- parent.env(environment())
nullfit <- get("nullfit",enclosing_env)
fullfit <- get("fullfit",enclosing_env)
# dotargs <- get("dotargs",enclosing_env)
test_obj <- get("test_obj",enclosing_env)
debug. <- get("debug.", enclosing_env)
# .condition <- get(".condition", enclosing_env)
null_call <- getCall(nullfit)
null_fit_fn <- .get_bare_fnname.HLfit(nullfit, call.=null_call)
full_fit_fn <- .get_bare_fnname.HLfit(fullfit)
conv_full <- conv_null <- FALSE
best_logL_full <- best_logL_null <- prev_logL_full <- prev_logL_null <- -Inf
# Allows the user to control the starting values of the initial new_nullfit
# edotargs <- dotargs
# for (st in setdiff(names(edotargs),"prior.weights")) {
# edotargs[[st]] <- eval(edotargs[[st]],env=environment()) ## evaluate the promises in the current execution envir
# }
ctrl_opt <- .update_control(fit_call=null_call, optim_boot=.spaMM.data$options$optim_boot, from_fn=null_fit_fn) # need .safe_opt when newinits are at bound.
newinits <- .get_outer_inits_from_fit(fitobject=nullfit, keep_canon_user_inits = FALSE) # for next new_nullfit
if (null_fit_fn=="fitme") {
new_args <- list(init=newinits, control=ctrl_opt)
} else if (null_fit_fn=="corrHLfit") {
new_args <- list(init.corrHLfit=newinits, control.corrHLfit=ctrl_opt)
} else new_args <- NULL
while ( TRUE ) {
if (debug.==2) {# Shuld be prevented by spaMM's calling fns in a parallel session
new_nullfit <- do.call(update_resp, c(list(object=nullfit, newresp = y),new_args))
} else {
new_nullfit <- try(do.call(update_resp, c(list(object=nullfit, newresp = y),new_args)))
if (inherits(new_nullfit,"try-error")) { ## (debug.= TRUE or 1L) to return error info in parallel mode: return the try-error object
if (debug.) {
utils::dump.frames(dumpto="dump_on_new_nullfit", to.file=TRUE) # but doesn't stop
return(list(full=structure(NA,new_nullfit=new_nullfit), # but the attributes seem lost at the level of the pbapply() closure.
null=structure(NA,info="no nullfit")))
} else return(c(full=NA, null=NA)) # attributes would be lost at the level of the pbapply() closure.
} ## ELSE continue
}
logL_new_null <- logLik(new_nullfit,which=test_obj)
#cat(logL_new_null)
conv_null <- (abs(logL_new_null - prev_logL_null)<1e-4)
if (logL_new_null>best_logL_null) { # always true the first time
best_logL_null <- logL_new_null
best_nullfit <- new_nullfit
}
if (conv_null) break # no point in refitting the full model if the new inits from null fit don't change
# ELSE
prev_logL_null <- logL_new_null
newinits <- get_inits_from_fit(from=best_nullfit, template=fullfit, to_fn=full_fit_fn) # using default 'inner_lambdas' argument
lens <- rep(NA, length(newinits))
for (it in seq_along(newinits)) lens[it] <- length(newinits[[it]])
newinits <- newinits[lens>0L]
subsets_inits <- .all.subsets(newinits) ## limited scope: could effectively do this is a nested way (with a skeleton+relist technique)
# subsets_inits is always a list of lists, with at least one element (minimally: list(list()))
for (it in seq_along(subsets_inits)) {
inits <- subsets_inits[[it]]
if (full_fit_fn=="fitme") {
new_args <- list(init=inits[["init"]], init.HLfit=inits[["init.HLfit"]], control=ctrl_opt)
} else if (full_fit_fn=="corrHLfit") {
new_args <- list(init.corrHLfit=inits[["init.corrHLfit"]], init.HLfit=inits[["init.HLfit"]], control.corrHLfit=ctrl_opt)
} else new_args <- list(init.HLfit=inits[["init.HLfit"]])
if (debug.==2) {
new_fullfit <- do.call(update_resp, c(list(object=fullfit, newresp = y),new_args)) # may stop on error
} else {
new_fullfit <- try(do.call(update_resp, c(list(object=fullfit, newresp = y),new_args)))
if (inherits(new_fullfit,"try-error")) { ## (debug.= TRUE or 1L) to return error info in parallel mode: return the try-error object
if (debug.) {
utils::dump.frames(dumpto="dump_on_new_fullfit", to.file=TRUE) # but doesn't stop
return(list(full=structure(NA,new_fullfit=new_fullfit), # but the attributes seem lost at the level of the pbapply() closure.
null=structure(logL_new_null,new_nullfit=new_nullfit)))
} else return(c(full=NA, null=logL_new_null)) # attributes would be lost at the level of the pbapply() closure.
} ## ELSE continue
}
logL_new_full <- logLik(new_fullfit,which=test_obj)
#cat(" ",logL_new_full,"\n")
conv_full <- (abs(logL_new_full - prev_logL_full)<1e-4)
if (logL_new_full>best_logL_full) { # always true the first time
best_fullfit <- new_fullfit
best_logL_full <- logL_new_full
if (logL_new_full>best_logL_null) break # the for loop
}
}
if (conv_full) break # no point in refitting the null model if the new inits from full fit don't change
# ELSE
prev_logL_full <- logL_new_full
newinits <- .get_outer_inits_from_fit(fitobject=best_fullfit, keep_canon_user_inits = FALSE) # for next new_nullfit
if (null_fit_fn=="fitme") {
new_args <- list(init=newinits, control=ctrl_opt)
} else if (null_fit_fn=="corrHLfit") {
new_args <- list(init.corrHLfit=newinits, control.corrHLfit=ctrl_opt)
} else new_args <- NULL
} # end while()
# print(logLik(new_fullfit,which=test_obj) - logLik(new_nullfit,which=test_obj))
resu <- c(full=best_logL_full,null=best_logL_null)
return(resu)
}
.check_bootreps <- function(bootreps) {
if (is.list(bootrep <- bootreps[[1]])) { # debug.=TRUE:
# If error in LRT() with debug.=TRUE and not doSNOW (at least), the replicate result is a list.
if (anyNA(bootrep)) { # ...and the list contains NA
first_failed <- which(is.na(bootrep))[1] # first of pair of fits
failure <- attributes(bootrep[[first_failed]])[[1]]
if (inherits(failure,"try-error") ) { # ...and ideally the first attr is a try error
condmess <- conditionMessage(attr(failure,"condition"))
if (length(grep("could not find",condmess))) {
firstpb <- strsplit(condmess,"\"")[[1]][2]
cat(crayon::bold(paste0(
"Hmmm. It looks like some variables were not passed to the parallel processes.\n",
"Maybe add ",firstpb," = ",firstpb," to spaMM_boot()'s 'fit_env' argument?\n"
)))
} else {
message(crayon::bold("If debug.=TRUE was used, a dump file may have been saved."))
}
stop(condmess,call. = FALSE)
} else { # ideally this does not happen, but just in case...
cat(crayon::bold(
"Hmmm. An unanticipated issue occurred. Here is the struture of the first replicate:\n",
))
str(bootrep)
}
}
} else if (anyNA(bootreps)) { # debug.=FALSE:
bootreps <- na.omit(bootreps)
if ( ! nrow(bootreps)) {
warnmess <- paste0("All bootstrap replicate(s) apparently failed.\n",
"To diagnose issues specific to parallel execution, argument\n",
"'debug.=TRUE' in LRT() or spaMM_boot() call may be useful\n",
"(it may generate a dump file from a child process for inspection).")
warning(warnmess, immediate.=TRUE, call.=FALSE)
} else {
n_omitted <- length(attr(bootreps,"na.action"))
warnmess <- paste0(n_omitted," bootstrap replicate(s) apparently failed and are omitted for p-value computations.")
}
warnmess
}
}
.warn_ests_at_bound <- function(fitobject) {
pars_at_bound <- NULL
outer_at_bound <- .check_outer_at_bound(fitobject, bool=FALSE)
if (length(outer_at_bound)) pars_at_bound <- paste(names(outer_at_bound),'=',signif(outer_at_bound), collapse=", ")
inner_lambdas <-.get_lambdas_notrC_from_hlfit(fitobject, type="inner") # __F I X M E___ not all inner params checked.
# inner ranCoefs not very interesting. Detecting inner phi at bound seems contrived (moreover with mv fits).
inner_at_bound <- (inner_lambdas<1e-6)
if (any(inner_at_bound)) {
inner_at_bound <- inner_lambdas[which(inner_at_bound)]
pars_at_bound <- paste(
pars_at_bound,
paste0("lambda.",names(inner_at_bound),' = ',signif(inner_at_bound), collapse=", "),
collapse=", "
)
warnmess <- paste("Full model fitted at bound of parameter ranges, for:\n",
pars_at_bound,"\n",
"=> Asymptotic test is unreliable (even with Bartlett correction).")
warning(warnmess, immediate.=TRUE, call.=FALSE)
}
}
LRT <- function(object,object2,boot.repl=0L,# nb_cores=NULL,
resp_testfn=NULL, simuland=eval_replicate,
# .condition = NULL, ## bc expected by simuland, but not operational,
...) {
if (nrow(object$data)!=nrow(object2$data)) {
stop("models were not both fitted to the same size of dataset.")
}
#if (length(list(...))) warning("...' arguments are currently ignored in LRT()", immediate. = TRUE)
# which is a bit unfortunate ( say ...control=list(optimizer="bobyqa")) but makes parallelisation so much more straightforward...
info <- .compare_model_structures(object,object2,
args_are_phimodels=FALSE,
boot.repl=boot.repl)
fullfit <- info$fullfit
fix_neg_LRT <- TRUE
if (is.null(fullfit)) { # could not determine which is nullfit and which is fullfit
# => nested resid models involving resid HGLM(s); or possibly non nested models.
if (boot.repl) {
fix_neg_LRT <- FALSE
fullfit <- object
nullfit <- object2
# (__F I X M E___ afficher la promise pour object2? mais alors in faut resoudre les ..n potentiels)
message(crayon::bold(paste0("Tentatively using 'object2' as \"null\" model\n",
" of the test, and thus as sample-generating model for bootstrap:")))
} else return(NULL)
} else nullfit <- info$nullfit
df <- info$df
if ( ! is.na(df) ) .warn_ests_at_bound(fitobject=fullfit)
test_obj <- info$test_obj
LRTori <- 2*(logLik(fullfit,which=test_obj)-logLik(nullfit,which=test_obj))
if (is.na(df)) {
resu <- list(nullfit=nullfit,fullfit=fullfit,basicLRT = data.frame(chi2_LR=LRTori,df=NA,p_value=NA))
} else {
pvalue <- 1-pchisq(LRTori,df=df) ## but not valid for testing null components of variance
resu <- list(nullfit=nullfit,fullfit=fullfit,basicLRT = data.frame(chi2_LR=LRTori,df=df,p_value=pvalue)) ## format appropriate for more tests
}
if (boot.repl) {
if (boot.repl<99L && ! is.na(df)) message("Note: It is recommended to set boot.repl>=99 for Bartlett correction")
# dotargs <- match.call(expand.dots = FALSE)$... ## produce a pairlist of (essentially) promises. No quote() needed
isdebugd <- isdebugged(simuland) # bc the assignment of environment drops this status
environment(simuland) <- environment() # enclosing env(simuland) <- evaluation env(LRT)
if (isdebugd) debug(simuland)
bootblob <- spaMM_boot(object=nullfit,nsim = boot.repl,
simuland=simuland,
#nb_cores = nb_cores,
resp_testfn = resp_testfn,
#aslistfull=aslistfull, aslistnull=aslistnull#, simbData=simbData,
# debug.=debug., now in the ...
type="marginal", # mandatory arg of spaMM_boot()
#, control.foreach=list(.errorhandling="pass")
...
)
bootblob$warnings$n_omitted <- .check_bootreps(bootblob$bootreps)
resu <- .add_boot_results(bootblob, resu, LRTori, df, test_obj, fix_neg_LRT=fix_neg_LRT)
}
class(resu) <- c("fixedLRT",class(resu))
return(resu)
}
.anova.lm <- function (object, ...) {
if (!inherits(object, "HLfit"))
warning("calling anova.lm(<fake-HLfit-object>) ...")
w <- weights(object, type="prior")
ssr <- sum(if (is.null(w)) residuals.HLfit(object,type="response")^2 else w * residuals.HLfit(object,type="response")^2)
mss <- sum(if (is.null(w)) fitted(object)^2 else w *
fitted(object)^2)
if (ssr < 1e-10 * mss)
warning("ANOVA F-tests on an essentially perfect fit are unreliable")
dfr <- df.residual(object)
p <- object$dfs$pforpv
if (p > 0L) {
qr_sXaug <- object$envir$qr_X
comp <- crossprod(qr.Q(qr_sXaug),object$y) # object$effects[p1] # the better tibco doc says it is t(Q) %*% y
# but given it is the QR for scaled X varaibles,
if (is.null(qr_sXaug)) stop("HLfit object does not have a 'qr' factorization of the model matrix.")
# sub_pivot_ori_X <- attr(object$X.pv,"rankinfo")$whichcols
asgn <- attr(model.matrix(object),"assign")# [sub_pivot_ori_X]
nmeffects <- c("(Intercept)", attr(terms(object), "term.labels"))
tlabels <- nmeffects[1 + unique(asgn)]
ss <- c(vapply(split(comp^2, asgn), sum, 1), ssr)
df <- c(lengths(split(asgn, asgn)), dfr)
}
else {
ss <- ssr
df <- dfr
tlabels <- character()
}
ms <- ss/df
f <- ms/(ssr/dfr)
P <- pf(f, df, dfr, lower.tail = FALSE)
table <- data.frame(df, ss, ms, f, P)
table[length(P), 4:5] <- NA ## row of residual SS
dimnames(table) <- list(c(tlabels, "Residuals"), c("Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)"))
#if (attr(object$terms, "intercept")) table <- table[-1, ]
table <- table[ ! rownames(table) == "(Intercept)", ]
structure(table, heading = c("Analysis of Variance Table\n",
paste("Response:", deparse((formula(object))[[2L]]))), ## formula ## FIXME use formula.HLfit()
class = c("anova", "data.frame"))
}
.null.deviance <- function(object, intercept=attr(terms(object), "intercept")) { # see glm.fit
if (length(offset <- model.offset(model.frame(object))) && intercept > 0L) {
# in that case it is necessary to fit the model with intercept and offset to obtain the sensible null deviance.
# This is what glm() does, calling [if (length(offset) && attr(mt, "intercept") > 0L) { fit2 <- ...]
termsv <- terms(object) # sufficient for fixed-effect model
offset_term <- (as.character(attr(termsv, "variables"))[-1L])[[attr(termsv, "offset")]]
newform <- as.formula(paste(" . ~ 1 +",offset_term))
null_dev_fit <- update(object, formula.= newform )
deviance(null_dev_fit)
} else {
# that's what is returned as $null.deviance by glm.fit,
pw <- weights(object, type="prior")
wtdmu <- if (intercept) {
sum(pw * object$y)/sum(pw)
} else object$family$linkinv(model.offset(model.frame(object)))
sum(object$family$dev.resids(object$y, wtdmu, pw))
}
}
.df.null <- function(object, intercept=attr(terms(object), "intercept")) {
n.ok <- length(object$y) - sum(weights(object, type="prior") == 0)
nulldf <- n.ok - as.integer(intercept)
}
.drop.terms <- function (termobj, dropx = NULL, keep.response = FALSE) {
if (is.null(dropx)) {
termobj
} else {
if (!inherits(termobj, "terms"))
stop(gettextf("'termobj' must be a object of class %s",
dQuote("terms")), domain = NA)
itcp <- attr(termobj, "intercept")
if ( ( ! length(dropped <- attr(termobj, "term.labels")[-dropx])) &&
itcp) dropped <- '1'
newformula <- reformulate(dropped, response = if (keep.response) termobj[[2L]],
intercept = itcp, env = environment(termobj))
result <- terms(newformula, specials = names(attr(termobj, "specials")))
response <- attr(termobj, "response")
dropOpt <- if (response && !keep.response)
c(response, dropx + length(response))
else dropx + max(response)
if (!is.null(predvars <- attr(termobj, "predvars"))) {
attr(result, "predvars") <- predvars[-(dropOpt +
1)]
}
if (!is.null(dataClasses <- attr(termobj, "dataClasses"))) {
attr(result, "dataClasses") <- dataClasses[-dropOpt]
}
result
}
}
.anova.glm <- function(object, ..., dispersion = NULL, test = NULL) {
doscore <- !is.null(test) && test == "Rao"
x <- model.matrix(object)
if (is.null(asgn <- attr(x,"assignOri"))) asgn <- attr(x,"assign")
nvars <- max(0, asgn)
resdev <- resdf <- NULL
termsv <- terms(object)
if (doscore) {
score <- numeric(nvars) # misnomer: regressors, not variables
subx <- x[, asgn == 0L, drop = FALSE]
# E.g., the design matrix 'x' may have a col for intercept, two cols for a first factor, two cols for a second factor
# asgn is 0 1 1 2 2 and there are 3 terms in the euation. dropterns removes the terms from the equation
# and x[, asgn <= i, drop = FALSE] removes the corresponding blocs of cols
# 0 always stands for the Intercept, not for the first term.
# This first doscore fit uses the original GLM family ect to generate residuals used in the next doscore fits which are LMs
# x = x[, asgn == 0, drop = FALSE], y = y, weights = object$prior.weights,
# start = object$start, offset = object$offset, family = object$family,
# drop.terms fails to remove all terms because because attr(,"term.labels") does not include the intercept
# hence calling the patch fn .drop.terms()
newform <- .drop.terms(termsv, seq_len(nvars), keep.response = TRUE)
refit <- update(object, formula.= newform)
r <- residuals(refit, type="working")
w <- weights(refit, type="working")
icpt <- attr(termsv, "intercept")
}
if (nvars > 1 || doscore) {
method <- object$method
for (i in seq_len(max(nvars - 1L, 0))) {
newform <- drop.terms(termsv, (i+1L):nvars, keep.response = TRUE)
refit <- update(object, formula.= newform)
if (doscore) {
subx <- x[, asgn <= i, drop = FALSE]
zz <- glm.fit(x=subx, y = r, weights = w, intercept = icpt)
score[i] <- zz$null.deviance - zz$deviance
r <- residuals(refit, type="working")
w <- weights(refit, type="working")
}
resdev <- c(resdev, deviance(refit))
resdf <- c(resdf, df.residual(refit))
}
if (doscore) {
zz <- glm.fit(x=x, y = r, weights = w, intercept = icpt)
score[nvars] <- zz$null.deviance - zz$deviance
}
}
resdf <- c(.df.null(object), resdf, df.residual(object))
resdev <- c(.null.deviance(object), resdev, deviance(object))
table <- data.frame(c(NA, -diff(resdf)), c(NA, pmax(0, -diff(resdev))),
resdf, resdev)
tl <- attr(termsv, "term.labels")
if (length(tl) == 0L)
table <- table[1, , drop = FALSE]
dimnames(table) <- list(c("NULL", tl), c("Df", "Deviance",
"Resid. Df", "Resid. Dev"))
if (doscore)
table <- cbind(table, Rao = c(NA, score))
varlist <- attr(termsv, "variables")
title <- paste0("Analysis of Deviance Table", "\n\nModel: ",
object$family$family, ", link: ", object$family$link,
"\n\nResponse: ", as.character(varlist[-1L])[1L], "\n\nTerms added sequentially (first to last)\n\n")
df.dispersion <- Inf
if (is.null(test)) {
# Trying to integrate new default in stats:::anova.glm
fam <- object$family
if (is.null(dispersion)) {
dispersion <- residVar(object,"fit")
# for Gamma GLM anova.glm uses the MME estimate based on Pearson residuals, dispersion <- summary(object, dispersion = dispersion)$dispersion
if (object$dfs$p_fixef_phi==0L) { # __F I X M E___ too narrow conditions for test <- FALSE ?
df.dispersion <- Inf
test <- "Chisq"
# We actually tested phi model = phiScal in the calling function
} else if (is.null(fam$dispersion)) { # cautious when not a stats:: family
test <- FALSE
} else { #
test <- "F"
df.dispersion <- df.residual(object)
}
} else test <- "Chisq" # if dispersion explicitly fixed by user
} else {
if (is.null(dispersion)) {
dispersion <- residVar(object,"fit")
# for Gamma GLM anova.glm uses the MME estimate based on Pearson residuals, dispersion <- summary(object, dispersion = dispersion)$dispersion
df.dispersion <- if (object$dfs$p_fixef_phi==0L) {
# we actually tested phi model = phiScal in the calling function
Inf
} else df.residual(object)
}
}
if (!isFALSE(test)) {
if (test == "F" && df.dispersion == Inf) {
famfam <- object$family$family
if (famfam == "binomial" || fam == "poisson")
warning(gettextf("using F test with a '%s' family is inappropriate",
famfam), domain = NA)
else warning("using F test with a fixed dispersion is inappropriate")
}
table <- stat.anova(table = table, test = test, scale = dispersion,
df.scale = df.residual(object), n = NROW(x))
}
structure(table, heading = title, class = c("anova", "data.frame"))
}
# LU facto with 1 on diag of L (Doolittle's scheme)
.lu_doo <- function(x, eps = 1e-06) {
stopifnot(ncol(x) > 1L)
n <- nrow(x)
L <- U <- matrix(0, nrow = n, ncol = n)
diag(L) <- rep(1, n)
for (i in 1:n) {
ip1 <- i + 1
im1 <- i - 1
for (j in 1:n) {
U[i, j] <- x[i, j]
if (im1 > 0) {
for (k in 1:im1) {
U[i, j] <- U[i, j] - L[i, k] * U[k, j]
}
}
}
if (ip1 <= n) {
for (j in ip1:n) {
L[j, i] <- x[j, i]
if (im1 > 0) {
for (k in 1:im1) {
L[j, i] <- L[j, i] - L[j, k] * U[k, i]
}
}
L[j, i] <- if (abs(U[i, i]) < eps)
0
else L[j, i]/U[i, i]
}
}
}
L[abs(L) < eps] <- 0
U[abs(U) < eps] <- 0
list(L = L, U = U)
}
.term2colX <- function (term_names, itcp, asgn) {
col_terms <- if (itcp)
c("(Intercept)", term_names)[asgn + 1L]
else term_names[asgn[asgn > 0L]]
nm <- union(unique(col_terms), term_names)
res <- lapply(setNames(as.list(nm), nm), function(x) numeric(0L))
map <- split(seq_along(col_terms), col_terms)
res[names(map)] <- map
res[nm]
}
.get_type1_contrasts <- function (model,
termsv=terms(model),
X=model.matrix(model),
asgn=attr(X,"assign"), ## "for each column in the matrix ... the term in the formula which gave rise to the column" (for mv, see below)
rankinfo=attr(X,"rankinfo"),
colrange=seq_len(ncol(X))) {
if (inherits(termsv,"list")) {
resu <- vector("list",length(termsv))
col_ranges <- attr(X,"col_ranges")
for (mv_it in seq_along(termsv)) {
colrange <- col_ranges[[mv_it]]
# here asgn must be a list with elements "for each column in the [sub]matrix ... the term in the [sub]formula which gave rise to the column"
resu[[mv_it]] <- .get_type1_contrasts(model, termsv=termsv[[mv_it]], X=X, asgn=asgn[[mv_it]],
rankinfo=rankinfo[[mv_it]], colrange=colrange)
names(resu[[mv_it]]) <- paste0(names(resu[[mv_it]]),"_",mv_it)
}
return(unlist(resu,use.names = TRUE, recursive = FALSE))
}
p <- length(colrange)
if (p == 0L) return(list(matrix(numeric(0L), nrow = 0L)))
itcp <- attr(termsv, "intercept")
if (p == 1L && itcp) return(NULL) # return(list(matrix(numeric(0L), ncol = ncol(X)))) # if only an intercept term
L <- if (p == 1L) { matrix(1L,, ncol = ncol(X)) } else t(.lu_doo(crossprod(X))$L)
dimnames(L) <- list(colnames(X), colnames(X))
term_names <- attr(termsv, "term.labels")
term_names <- term_names[unique(asgn[rankinfo$whichcols])] # line handling singular design matrix (the present X already being reduced)
ind.list <- .term2colX(term_names=term_names, itcp=itcp, asgn= asgn)[term_names]
lapply(ind.list, function(rows) L[colrange[rows], , drop = FALSE])
}
.term_contain <- function (term, factors, dataClasses, term_names)
{
get_vars <- function(term) rownames(factors)[factors[, term] == 1L]
contain <- function(F1, F2) {
all(vars[[F1]] %in% vars[[F2]]) &&
length(setdiff(vars[[F2]], vars[[F1]])) > 0L && setequal(numerics[[F1]], numerics[[F2]])
}
vars <- lapply(setNames(term_names, term_names), get_vars)
numerics <- lapply(vars, function(varnms) varnms[which(dataClasses[varnms] == "numeric")])
sapply(term_names, function(term_nm) contain(term, term_nm))
}
.containment <- function (termsv) {
data_classes <- attr(termsv, "dataClasses")
term_names <- attr(termsv, "term.labels")
factor_mat <- attr(termsv, "factors")
lapply(setNames(term_names, term_names), function(term) {
term_names[.term_contain(term, factor_mat, data_classes, term_names)]
})
}
.get_type2_contrasts <- function (model,
termsv=terms(model),
X=model.matrix(model),
asgn=attr(X,"assign"),
rankinfo=attr(X,"rankinfo"),
colrange=seq_len(ncol(X))) {
if (inherits(termsv,"list")) {
resu <- vector("list",length(termsv))
col_ranges <- attr(X,"col_ranges")
for (mv_it in seq_along(termsv)) {
colrange <- col_ranges[[mv_it]]
resu[[mv_it]] <- .get_type2_contrasts(model, termsv=termsv[[mv_it]], X=X, asgn=asgn[[mv_it]],
rankinfo=rankinfo[[mv_it]],colrange=colrange)
names(resu[[mv_it]]) <- paste0(names(resu[[mv_it]]),"_",mv_it)
}
return(unlist(resu,use.names = TRUE, recursive = FALSE))
}
data_classes <- attr(termsv, "dataClasses")
term_names <- attr(termsv, "term.labels")
#
which <- term_names
if (length(colrange) <= 1L || length(term_names) <= 1L)
return(.get_type1_contrasts(model, termsv=termsv, X=X, asgn=asgn, rankinfo=rankinfo, colrange=colrange)) # possibly the arguments for the mv_it submodel
#
else stopifnot(is.character(which), all(which %in% term_names))
which <- setNames(as.list(which), which)
is_contained <- .containment(model)
itcp <- attr(termsv, "intercept") > 0
asgn <- unique(setdiff(asgn,0L)) # unique(asgn[rankinfo$whichcols])
term_names <- term_names[asgn]
which <- if (itcp) {c("(Intercept)", term_names)} else term_names
term2colX <- split(seq_along(which), which)[unique(which)]
Llist <- lapply(which, function(term) {
subX_termcols <- unlist(term2colX[c(term, is_contained[[term]])])
X_termcols <- colrange[subX_termcols]
colperm <- c(setdiff(seq_len(ncol(X)),X_termcols), X_termcols)
Xnew <- X[,colperm]
Lc <- t(.lu_doo(crossprod(Xnew))$L)
dimnames(Lc) <- list(colnames(Xnew), colnames(Xnew))
Lc[colperm %in% X_termcols, colnames(X), drop = FALSE]
})
names(Llist) <- which
Llist
}
# Not well tested
.get_type3_contrasts <- function(model,
termsv=terms(model),
Xorig=model.matrix(model),
asgn=attr(Xorig,"assign"),
rankinfo=attr(Xorig,"rankinfo"),
colrange=seq_len(ncol(Xorig))) {
if (inherits(termsv,"list")) {
resu <- vector("list",length(termsv))
col_ranges <- attr(Xorig,"col_ranges")
for (mv_it in seq_along(termsv)) {
colrange <- col_ranges[[mv_it]]
resu[[mv_it]] <- .get_type3_contrasts(model, termsv=termsv[[mv_it]], Xorig=Xorig, asgn=asgn[[mv_it]],
rankinfo=rankinfo[[mv_it]],colrange=colrange)
names(resu[[mv_it]]) <- paste0(names(resu[[mv_it]]),"_",mv_it)
}
return(unlist(resu,use.names = TRUE, recursive = FALSE))
}
data_classes <- attr(termsv, "dataClasses")
term_names <- attr(termsv, "term.labels")
which <- term_names
if (length(colrange) <= 1L || length(term_names) <= 1L)
return(.get_type1_contrasts(model, termsv=termsv, X=Xorig, asgn=asgn, rankinfo=rankinfo, colrange=colrange)) # possibly the arguments for the mv_it submodel
#
else stopifnot(is.character(which), all(which %in% term_names))
codings <- unlist(attr(Xorig, "contrast"))
if (length(codings) > 0 && all(is.character(codings)) &&
all(codings %in% c("contr.treatment"))) {
.extract_contrasts_type3 <- get("extract_contrasts_type3", asNamespace("lmerTest"))
return(.extract_contrasts_type3(model, X = Xorig))
}
.get_contrast_coding <- get("get_contrast_coding", asNamespace("lmerTest"))
Contrasts <- .get_contrast_coding(model, contrasts = "contr.treatment")
X <- model.matrix(terms(model), data = model.frame(model), contrasts.arg = Contrasts)
.ensure_full_rank <- get("ensure_full_rank", asNamespace("lmerTest"))
X <- .ensure_full_rank(X, silent = TRUE, test.ans = FALSE)
.extract_contrasts_type3 <- get("extract_contrasts_type3", asNamespace("lmerTest"))
type3ctr <- .extract_contrasts_type3(model, X = X)
map <- zapsmall(ginv(X) %*% Xorig)
rownames(map) <- colnames(X)
lapply(type3ctr[which], function(L) L %*% map)
}
.anova_fallback <- function(fitobject, type="2", rhs=NULL, test="Chisq.", ...) {
beta <- na.omit(fixef(fitobject))
beta_cov <- vcov(fitobject)
# Note that 'type' is formally documented for LMMs only, which do not use this code.
# So the following switches are not part of API, and not formally tested.
# (The original idea behing the .get_type[]_contrasts() was to perform the LMM tests using spaMM code rather than lmerTest)
Llist <- switch(type,
"I" = .get_type1_contrasts(fitobject),
"1" = .get_type1_contrasts(fitobject),
"II" = .get_type2_contrasts(fitobject),
"2" = .get_type2_contrasts(fitobject),
"III" = stop("type-3 contrasts not implemented (and may never be)"),
# .get_type3_contrasts(fitobject), #
"3" = stop("type-3 contrasts not implemented (and may never be)"),
# .get_type3_contrasts(fitobject), #
stop("Unknown ANOVA 'type' speification")
)
# no line for Intercept: as in lmerTest output
numtable <- matrix(NA,nrow=length(Llist), ncol = 3)
colnames(numtable) <- c("Df", test, paste("Pr(>", test, ")", sep = ""))
for (it in seq_along(Llist)) {
L <- Llist[[it]]
if (is.null(dim(L))) L <- t(L)
if (is.null(rhs)) rhs <- rep(0, nrow(L))
q <- NROW(L)
Lbeta <- L %*% beta - rhs
vcov_Lbeta <- L %*% tcrossprod(beta_cov,L)
waldX2 <- as.vector(t(Lbeta) %*% solve(vcov_Lbeta,Lbeta))
p <- pchisq(waldX2, q, lower.tail = FALSE)
numtable[it, ] <- c(q, waldX2, p)
}
resu <- as.data.frame(numtable)
rownames(resu) <- names(Llist)
roman.type <- if (type %in% c("1","2","3")) {utils::as.roman(as.integer(type))} else type
attr(resu, "heading") <- paste(test, "tests for each term (type", roman.type, "contrasts)")
attr(resu, "hypotheses") <- Llist
class(resu) <- c("anova", "data.frame")
resu
}
anova.HLfit <- function(object, object2=NULL, type="2", method="", ...) {
if (is.null(object2)) {
if (method == "") { # defaults
models <- object$models[c("eta","phi")]
if (length(models$phi)==1L && # exclude mv fits => all in .anova_fallback()
models$phi %in% c("phiScal","")) {
if (models$eta=="etaGLM") {
if ( ! type %in% c("2","II")) {
warning("anova.HLfit() for (G)LMs ignores 'type' consistently with anova.[g]lm().")
}
if (object$family$family=="gaussian" && object$family$link=="identity") {
return(.anova.lm(object, ...))
} else return(.anova.glm(object, ...))
} else if (object$family$family=="gaussian" && object$family$link=="identity") { # LMM
if (requireNamespace("lmerTest",quietly=TRUE)) { # if the package is available
# if (length(attr(object$X.pv,"namesOri")) == ncol(object$X.pv)) {
lmlt <- as_LMLT(object, ...)
# next line calls spaMM:::anova.LMLT ->
# as(as(object,"LMLTslots"),"LMLTinternal")
# and lmerTest:::anova.lmerModLmerTest as LMLTinternal contains = c("LMLTslots","lmerModLmerTest")
return(anova(lmlt, type=type)) # other possible arguments not currently meaningful
# } else {
# message(paste("Original model has a singular design matrix: lmerTest's F tests cannot be computed\n",
# "(you can fix this by correcting the model formula)"))
# }
} else if ( ! identical(spaMM.getOption("lmerTest_warned"),TRUE)) {
message("If the lmerTest package were installed, a traditional anova table could be computed.")
.spaMM.data$options$lmerTest_warned <- TRUE
}
}
}
}
# Fallback for non-LM, GLM or LMM fits: (only case where spaMM's contrast code is used).
return(.anova_fallback(fitobject=object, type=type, ...)) # dots may contain 'rhs' and (later) 'test'
} else {
## anova treated as alias for LRT()
LRT(object,object2, ...)
}
}
# anova is an S3 class, LMLT an S4 class
## Aim of this method is to ensure that the lmerTest package is loaded when anova() is called
# bc it may be called when object has class LMLT, but no def may yet be accessible for such objects.
## See the 'ZZZ' version for detailed explanations of the code.
anova.LMLT <- function(object, ...) {
if (is.null(getClassDef("LMLT", where = .spaMM.data$class_cache, inherits = FALSE))) {
# is as_LMLT has not been called since restarting R session-> lmerTest presumably not loaded
if (requireNamespace("lmerTest",quietly=TRUE)) {
# Hack to define object not of class LMLT (-> infinite recursion) but with similar "contains", using only default coerce methods
setClass(Class="LMLT", contains = c("LMLTslots","lmerModLmerTest"), where=.spaMM.data$class_cache)
setClass(Class="LMLTinternal", contains = c("LMLTslots","lmerModLmerTest"), where=.spaMM.data$class_cache)
} else message("If the lmerTest package were installed, a traditional anova table could be computed.")
}
object <- as(as(object,"LMLTslots"),"LMLTinternal")
anova(object, ...)
}
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.