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
12cao  <-
13  function(formula,
14           family = stop("argument 'family' needs to be assigned"),
15           data = list(),
16           weights = NULL, subset = NULL, na.action = na.fail,
17           etastart = NULL, mustart = NULL, coefstart = NULL,
18           control = cao.control(...),
19           offset = NULL,
20           method = "cao.fit",
21           model = FALSE, x.arg = TRUE, y.arg = TRUE,
22           contrasts = NULL,
23           constraints = NULL,
24           extra = NULL,
25           qr.arg = FALSE, smart = TRUE, ...) {
26  dataname <- as.character(substitute(data))  # "list" if no data=
27  function.name <- "cao"
28
29  ocall <- match.call()
30
31  if (smart)
32    setup.smart("write")
33
34  mt <- terms(formula, data = data)
35  if (missing(data))
36    data <- environment(formula)
37
38  mf <- match.call(expand.dots = FALSE)
39  mf$family <- mf$method <- mf$model <- mf$x.arg <- mf$y.arg <-
40    mf$control <-
41    mf$contrasts <- mf$constraints <- mf$extra <- mf$qr.arg <- NULL
42  mf$coefstart <- mf$etastart <- mf$... <- NULL
43  mf$smart <- NULL
44  mf$drop.unused.levels <- TRUE
45  mf[[1]] <- as.name("model.frame")
46  mf <- eval(mf, parent.frame())
47  if (method == "model.frame")
48    return(mf)
49  na.act <- attr(mf, "na.action")
50
51  xvars <- as.character(attr(mt, "variables"))[-1]
52  if ((yvar <- attr(mt, "response")) > 0)
53    xvars <- xvars[-yvar]
54  xlev <- if (length(xvars) > 0) {
55    xlev <- lapply(mf[xvars], levels)
56    xlev[!sapply(xlev, is.null)]
57  }
58
59  y <- model.response(mf, "numeric")  # model.extract(mf, "response")
60  x <- model.matrix(mt, mf, contrasts)
61  attr(x, "assign") <- attrassigndefault(x, mt)
62  offset <- model.offset(mf)
63  if (is.null(offset))
64    offset <- 0  # yyy ???
65  w <- model.weights(mf)
66  if (!length(w)) {
67    w <- rep_len(1, nrow(mf))
68  } else if (NCOL(w) == 1 && any(w < 0))
69    stop("negative weights not allowed")
70
71  if (is.character(family))
72    family <- get(family)
73  if (is.function(family))
74    family <- family()
75  if (!inherits(family, "vglmff")) {
76    stop("'family = ", family, "' is not a VGAM family function")
77  }
78
79  eval(vcontrol.expression)
80
81  if (!is.null(family@first))
82    eval(family@first)
83
84
85  cao.fitter <- get(method)
86
87
88  deviance.Bestof <- rep_len(NA_real_, control$Bestof)
89  for (tries in 1:control$Bestof) {
90    if (control$trace && (control$Bestof > 1)) {
91      cat(paste("\n========================= Fitting model",
92          tries, "=========================\n"))
93      if (exists("flush.console"))
94        flush.console()
95    }
96    onefit <-
97      cao.fitter(x = x, y = y, w = w, offset = offset,
98                 etastart = etastart, mustart = mustart,
99                 coefstart = coefstart,
100                 family = family,
101                 control = control,
102                 constraints = constraints,
103                 criterion = control$criterion,
104                 extra = extra,
105                 qr.arg = qr.arg,
106                 Terms = mt, function.name = function.name, ...)
107    deviance.Bestof[tries] <- onefit$crit.list$deviance
108    if (tries == 1 ||
109        min(deviance.Bestof[1:(tries-1)]) > deviance.Bestof[tries])
110      fit <- onefit
111  }
112  fit$misc$deviance.Bestof <- deviance.Bestof
113
114  fit$misc$dataname <- dataname
115
116  if (smart) {
117    fit$smart.prediction <- get.smart.prediction()
118    wrapup.smart()
119  }
120
121  answer <-
122  new("rrvgam",
123    "assign"       = attr(x, "assign"),
124    "Bspline"      = fit$Bspline,
125    "call"         = ocall,
126    "coefficients" = fit$coefficients,
127    "criterion"    = fit$crit.list,
128    "family"       = fit$family,
129    "misc"         = fit$misc,
130    "model"        = if (model) mf else data.frame(),
131    "residuals"    = as.matrix(fit$wresiduals),
132    "smart.prediction" = as.list(fit$smart.prediction),
133    "terms"        = list(terms = mt))
134
135  if (!smart)
136    answer@smart.prediction <- list(smart.arg = FALSE)
137
138  if (qr.arg) {
139    class(fit$qr) <- "list"
140    slot(answer, "qr") <- fit$qr
141  }
142  if (length(attr(x, "contrasts")))
143    slot(answer, "contrasts") <- attr(x, "contrasts")
144  if (length(fit$fitted.values))
145    slot(answer, "fitted.values") <- as.matrix(fit$fitted.values)
146  slot(answer, "na.action") <-
147    if (length(na.act)) list(na.act) else list()
148  if (length(offset))
149    slot(answer, "offset") <- as.matrix(offset)
150  if (length(fit$weights))
151    slot(answer, "weights") <- as.matrix(fit$weights)
152  if (x.arg)
153    slot(answer, "x") <- fit$x  # The 'small' design matrix
154  if (length(xlev))
155    slot(answer, "xlevels") <- xlev
156  if (y.arg)
157    slot(answer, "y") <- as.matrix(fit$y)
158
159
160  slot(answer, "control") <- fit$control
161  slot(answer, "extra") <- if (length(fit$extra)) {
162      if (is.list(fit$extra)) fit$extra else {
163          warning("'extra' is not a list, therefore ",
164                  "placing 'extra' into a list")
165          list(fit$extra)
166      }
167  } else list()  # R-1.5.0
168
169  slot(answer, "iter") <- fit$iter
170  fit$predictors <- as.matrix(fit$predictors)  # Must be a matrix
171  dimnames(fit$predictors) <- list(dimnames(fit$predictors)[[1]],
172                                   fit$misc$predictors.names)
173  slot(answer, "predictors") <- fit$predictors
174  if (length(fit$prior.weights))
175    slot(answer, "prior.weights") <- as.matrix(fit$prior.weights)
176
177
178
179
180
181  answer
182}
183attr(cao, "smart") <- TRUE
184
185
186
187