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