1# These functions are 2# Copyright (C) 1998-2021 T.W. Yee, University of Auckland. 3# All rights reserved. 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18step4vglm <- 19 function (object, scope, # scale = 0, 20 direction = c("both", "backward", "forward"), 21 trace = 1, keep = NULL, steps = 1000, k = 2, 22 ...) { 23 mydeviance <- function(x, ...) { 24 dev <- deviance(x) 25 res <- if (is.null(dev)) extractAIC(x, k = 0)[2L] else dev 26 res 27 } 28 cut.string <- function(string) { 29 if (length(string) > 1L) 30 string[-1L] <- paste0("\n", string[-1L]) 31 string 32 } 33 re.arrange <- function(keep) { 34 namr <- names(k1 <- keep[[1L]]) 35 namc <- names(keep) 36 nc <- length(keep) 37 nr <- length(k1) 38 array(unlist(keep, recursive = FALSE), 39 c(nr, nc), list(namr, namc)) 40 } 41 step.results <- function(models, fit, object) { 42 change <- sapply(models, "[[", "change") 43 rd <- sapply(models, "[[", "deviance") 44 dd <- c(NA, abs(diff(rd))) 45 rdf <- sapply(models, "[[", "df.resid") 46 ddf <- c(NA, diff(rdf)) 47 AIC <- sapply(models, "[[", "AIC") 48 heading <- c("Stepwise Model Path \nAnalysis of Deviance ", 49 "Table", 50 "\nInitial Model:", deparse(formula(object)), 51 "\nFinal Model:", 52 deparse(formula(fit)), "\n") 53 aod <- data.frame(Step = I(change), 54 Df = ddf, 55 Deviance = dd, 56 `Resid. Df` = rdf, 57 `Resid. Dev` = rd, 58 AIC = AIC, 59 check.names = FALSE) 60 attr(aod, "heading") <- heading 61 fit@post$anova <- aod # fit$anova <- aod 62 fit 63 } # step.results 64 65 Terms <- terms(object) 66 object@call$formula <- object@misc$formula <- Terms 67 md <- missing(direction) 68 direction <- match.arg(direction) 69 backward <- is.element(direction, c("both", "backward")) 70 forward <- is.element(direction, c("both", "forward")) 71 if (missing(scope)) { 72 fdrop <- numeric() 73 fadd <- attr(Terms, "factors") 74 if (md) forward <- FALSE 75 } else { 76 if (is.list(scope)) { 77 fdrop <- if (!is.null(fdrop <- scope$lower)) 78 attr(terms(update.formula(object, fdrop)), "factors") else 79 numeric() 80 fadd <- if (!is.null(fadd <- scope$upper)) 81 attr(terms(update.formula(object, fadd)), "factors") 82 } else { 83 fadd <- if (!is.null(fadd <- scope)) 84 attr(terms(update.formula(object, scope)), "factors") 85 fdrop <- numeric() 86 } 87 } 88 models <- vector("list", steps) 89 if (!is.null(keep)) 90 keep.list <- vector("list", steps) 91 n.lm <- nobs(object, type = "lm") 92 n.vlm <- nobs(object, type = "vlm") 93 fit <- object 94 bAIC <- extractAIC(fit, k = k, ...) 95 edf <- bAIC[1L] 96 bAIC <- bAIC[2L] 97 if (is.na(bAIC)) 98 stop("AIC is not defined for this model, so 'step4' ", 99 "cannot proceed") 100 if (bAIC == -Inf) 101 stop("AIC is -infinity for this model, so 'step4' ", 102 "cannot proceed") 103 nm <- 1 104 if (trace) { 105 cat("Start: AIC=", format(round(bAIC, 2)), "\n", 106 cut.string(deparse(formula(fit))), 107 "\n\n", sep = "") 108 flush.console() 109 } 110 models[[nm]] <- list(deviance = mydeviance(fit), 111 df.resid = n.vlm - edf, 112 change = "", AIC = bAIC) 113 if (!is.null(keep)) 114 keep.list[[nm]] <- keep(fit, bAIC) 115 116 while (steps > 0) { 117 steps <- steps - 1 118 AIC <- bAIC 119 ffac <- attr(Terms, "factors") 120 scope <- factor.scope(ffac, list(add = fadd, drop = fdrop)) 121 aod <- NULL 122 change <- NULL 123 124 if (backward && length(scope$drop)) { 125 aod <- drop1(fit, scope$drop, k = k, ...) # trace = trace, 126 rn <- row.names(aod) 127 row.names(aod) <- c(rn[1L], paste("-", rn[-1L])) 128 if (any(aod$Df == 0, na.rm = TRUE)) { 129 zdf <- aod$Df == 0 & !is.na(aod$Df) 130 change <- rev(rownames(aod)[zdf])[1L] 131 } 132 } # if (backward && length(scope$drop)) 133 134 if (is.null(change)) { 135 if (forward && length(scope$add)) { 136 aodf <- add1(fit, scope$add, k = k, ...) # trace = trace, 137 rn <- row.names(aodf) 138 row.names(aodf) <- c(rn[1L], paste("+", rn[-1L])) 139 aod <- if (is.null(aod)) aodf else 140 rbind(aod, aodf[-1, , drop = FALSE]) 141 } 142 attr(aod, "heading") <- NULL 143 nzdf <- if (!is.null(aod$Df)) 144 aod$Df != 0 | is.na(aod$Df) 145 aod <- aod[nzdf, ] 146 if (is.null(aod) || ncol(aod) == 0) 147 break 148 149 nc <- match(c("Cp", "AIC"), names(aod)) 150 nc <- nc[!is.na(nc)][1L] 151 oo <- order(aod[, nc]) 152 if (trace) 153 print(aod[oo, ]) 154 if (oo[1L] == 1) 155 break 156 change <- rownames(aod)[oo[1L]] 157 } # if (is.null(change)) 158 159 160 fit <- update.default(fit, paste("~ .", change), 161 evaluate = FALSE) # update() 162 163 164 fit <- eval.parent(fit) 165 nnew <- nobs(fit, type = "vlm") # use.fallback = TRUE 166 if (all(is.finite(c(n.vlm, nnew))) && nnew != n.vlm) 167 stop("number of rows in use has changed: ", 168 "remove missing values?") 169 170 171 Terms <- terms(fit) 172 bAIC <- extractAIC(fit, k = k, ...) 173 edf <- bAIC[1L] 174 bAIC <- bAIC[2L] 175 if (trace) { 176 cat("\nStep: AIC=", format(round(bAIC, 2)), "\n", 177 cut.string(deparse(formula(fit))), "\n\n", sep = "") 178 flush.console() 179 } 180 if (bAIC >= AIC + 1e-07) 181 break 182 nm <- nm + 1 183 models[[nm]] <- list(deviance = mydeviance(fit), 184 df.resid = n.vlm - edf, 185 change = change, AIC = bAIC) 186 if (!is.null(keep)) 187 keep.list[[nm]] <- keep(fit, bAIC) 188 } # while (steps > 0) 189 190 if (!is.null(keep)) 191 fit$keep <- re.arrange(keep.list[seq(nm)]) 192 step.results(models = models[seq(nm)], fit, object) 193} # step4vglm 194 195 196 197 198 199if (!isGeneric("step4")) 200 setGeneric("step4", function(object, ...) 201 standardGeneric("step4"), 202 package = "VGAM") 203 204 205setMethod("step4", "vglm", 206 function(object, ...) 207 step4vglm(object, ...)) 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222