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