### Fit a general nonlinear regression model with correlated and/or ### heteroscedastic errors ### ### Copyright 1997-2003 Jose C. Pinheiro, ### Douglas M. Bates ### Copyright 2007-2017 The R Core team # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # A copy of the GNU General Public License is available at # http://www.r-project.org/Licenses/ # gnls <- function(model, data = sys.frame(sys.parent()), params, start, correlation = NULL, weights = NULL, subset, na.action = na.fail, naPattern, control = list(), verbose= FALSE) { finiteDiffGrad <- function(model, data, pars) { dframe <- data.frame(data, pars) base <- eval(model, dframe) nm <- colnames(pars) grad <- array(base, c(length(base), length(nm)), list(NULL, nm)) ssize <- sqrt(.Machine$double.eps) for (i in nm) { diff <- pp <- pars[ , i] diff[pp == 0] <- ssize diff[pp != 0] <- pp[pp != 0] * ssize dframe[[i]] <- pp + diff grad[ , i] <- (base - eval(model, dframe))/diff dframe[[i]] <- pp } grad } ## keeping the call Call <- match.call() ## assigning a new name to the "object" argument form <- model ## control parameters controlvals <- gnlsControl() if (!missing(control)) { controlvals[names(control)] <- control } ## ## checking arguments ## if (!inherits(form, "formula")) stop("'object' must be a formula") if (length(form)!=3) stop("object formula must be of the form \"resp ~ pred\"") ## ## checking if self-starting formula is given ## if (missing(start)) { if (!is.null(attr(eval(form[[3]][[1]]), "initial"))) { nlsCall <- Call[c("","model","data")] nlsCall[[1]] <- quote(stats::nls) names(nlsCall)[2] <- "formula" ## checking if "data" is not equal to sys.frame(sys.parent()) if (is.null(dim(data))) { stop("'data' must be given explicitly to use 'nls' to get initial estimates") } start <- coef(eval(nlsCall)) } else { stop("no initial values for model parameters") } } else { start <- unlist(start) } gnlsModel <- call("-", form[[2]], form[[3]]) ## ## save writing list(...) when only one element ## if (missing(params)) { if (is.null(pNams <- names(start))) { stop("starting estimates must have names when 'params' is missing") } params <- eval(parse(text = paste(paste(pNams, collapse = "+"), "1", sep = "~"))) } if (!is.list(params)) { params <- list(params) } val <- NULL for(i in seq_along(params)) { if (is.name(params[[i]][[2]])) { val <- c(val, list(params[[i]])) } else { ## multiple parameters on left hand side val <- c(val, eval(parse(text = paste("list(", paste(paste(all.vars(params[[i]][[2]]), deparse(params[[i]][[3]]), sep = "~"), collapse=","),")")))) } } params <- as.list(val) pnames <- character(length(params)) for (i in seq_along(params)) { this <- eval(params[[i]]) if (!inherits(this, "formula")) stop ("'params' must be a formula or list of formulae") if (length(this) != 3) stop ("formulae in 'params' must be of the form \"parameter ~ expr\"") if (!is.name(this[[2]])) stop ("formulae in 'params' must be of the form \"parameter ~ expr\"") pnames[i] <- as.character(this[[2]]) } names(params) <- pnames ## ## If data is a pframe, copy the parameters in the frame to frame 1 ## Doesn't exist in R ## if (inherits(data, "pframe")) { ## pp <- parameters(data) ## for (i in names(pp)) { ## assign(i, pp[[i]]) ## } ## attr(data,"parameters") <- NULL ## class(data) <- "data.frame" ## } ## check if correlation is present and has groups groups <- if (!is.null(correlation)) getGroupsFormula(correlation) # else NULL # if (!is.null(correlation)) { # groups <- getGroupsFormula(correlation, asList = TRUE) # if (!is.null(groups)) { # if (length(groups) > 1) { # stop("Only single level of grouping allowed") # } # groups <- groups[[1]] # } else { # if (inherits(data, "groupedData")) { # will use as groups # groups <- getGroupsFormula(data, asList = TRUE) # if (length(groups) > 1) { # ignore it # groups <- NULL # } else { # groups <- groups[[1]] # attr(correlation, "formula") <- # eval(parse(text = paste("~", # deparse(getCovariateFormula(formula(correlation))[[2]]), # "|", deparse(groups[[2]])))) # } # } # } # } else groups <- NULL ## create an gnls structure containing the correlation and weights gnlsSt <- gnlsStruct(corStruct = correlation, varStruct = varFunc(weights)) ## extract a data frame with enough information to evaluate ## form, params, random, groups, correlation, and weights mfArgs <- list(formula = asOneFormula(formula(gnlsSt), form, params, groups, omit = c(pnames, "pi")), data = data, na.action = na.action) if (!missing(subset)) { mfArgs[["subset"]] <- asOneSidedFormula(Call[["subset"]])[[2]] } mfArgs$drop.unused.levels <- TRUE dataMod <- do.call("model.frame", mfArgs) origOrder <- row.names(dataMod) # preserve the original order ## ## Evaluating the groups expression, if needed ## if (!is.null(groups)) { ## sort the model.frame by groups and get the matrices and parameters ## used in the estimation procedures ## always use innermost level of grouping groups <- eval(parse(text = paste("~1", deparse(groups[[2]]), sep = "|"))) grps <- getGroups(dataMod, groups, level = length(getGroupsFormula(groups, asList = TRUE))) ## ordering data by groups ord <- order(grps) grps <- grps[ord] dataMod <- dataMod[ord, ,drop = FALSE] ## revOrder <- match(origOrder, row.names(dataMod)) # putting in orig. order } else grps <- NULL N <- dim(dataMod)[1] # number of observations ## ## evaluating the naPattern expression, if any ## naPat <- if (missing(naPattern)) rep(TRUE, N) else as.logical(eval(asOneSidedFormula(naPattern)[[2]], dataMod)) origOrderShrunk <- origOrder[naPat] dataModShrunk <- dataMod[naPat, , drop=FALSE] yShrunk <- eval(form[[2]], dataModShrunk) if (!is.null(groups)) { ## ordShrunk <- ord[naPat] grpShrunk <- grps[naPat] revOrderShrunk <- match(origOrderShrunk, row.names(dataModShrunk)) } else { grpShrunk <- NULL } ## ## defining list with parameter information ## contr <- list() plist <- vector("list", length(pnames)) names(plist) <- pnames for (nm in pnames) { plist[[nm]] <- TRUE if (deparse(params[[nm]][[3]]) != "1") { plist[[nm]] <- model.matrix(asOneSidedFormula(params[[nm]][[3]]), model.frame(asOneSidedFormula(params[[nm]][[3]]), dataModShrunk)) auxContr <- attr(plist[[nm]], "contrasts") contr <- c(contr, auxContr[is.na(match(names(auxContr), names(contr)))]) } } ## ## Params effects names ## pn <- character(0) currPos <- 0 parAssign <- list() for(nm in pnames) { if (is.logical(p <- plist[[nm]])) { currPos <- currPos + 1 currVal <- list(currPos) pn <- c(pn, nm) names(currVal) <- nm parAssign <- c(parAssign, currVal) } else { currVal <- attr(p, "assign") fTerms <- terms(asOneSidedFormula(params[[nm]][[3]]), data=data) namTerms <- attr(fTerms, "term.labels") if (attr(fTerms, "intercept") > 0) { namTerms <- c("(Intercept)", namTerms) } namTerms <- factor(currVal, labels = namTerms) currVal <- split(order(currVal), namTerms) names(currVal) <- paste(nm, names(currVal), sep = ".") parAssign <- c(parAssign, lapply(currVal, function(el, currPos) { el + currPos }, currPos = currPos)) currPos <- currPos + length(unlist(currVal)) pn <- c(pn, paste(nm, colnames(p), sep = ".")) } } pLen <- length(pn) if (length(start) != pLen) stop ("starting values for parameters are not of the correct length") spar <- start names(spar) <- pn NReal <- sum(naPat) ## ## Creating the params map ## pmap <- list() n1 <- 1 for(nm in pnames) { if (is.logical(p <- plist[[nm]])) { pmap[[nm]] <- n1 n1 <- n1 + 1 } else { pmap[[nm]] <- n1:(n1+ncol(p) - 1) n1 <- n1 + ncol(p) } } ## ## defining the nlFrame, i.e., nlEnv, an environment in R : ## nlEnv <- list2env( list(model = gnlsModel, data = dataMod, plist = plist, beta = as.vector(spar), X = array(0, c(NReal, pLen), list(NULL, pn)), pmap = pmap, N = NReal, naPat = naPat, .parameters = c("beta"), finiteDiffGrad = finiteDiffGrad)) modelExpression <- ~ { pars <- getParsGnls(plist, pmap, beta, N) res <- eval(model, data.frame(data, pars)) if (!length(grad <- attr(res, "gradient"))) { grad <- finiteDiffGrad(model, data, pars)[naPat, , drop = FALSE] } else { grad <- grad[naPat, , drop = FALSE] } res <- res[naPat] for (nm in names(plist)) { gradnm <- grad[, nm] X[, pmap[[nm]]] <- if(is.logical(p <- plist[[nm]])) gradnm else gradnm * p } result <- c(X, res) result[is.na(result)] <- 0 result } modelResid <- ~eval(model, data.frame(data, getParsGnls(plist, pmap, beta, N)))[naPat] w <- eval(modelResid[[2]], envir = nlEnv) ## creating the condensed linear model ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions fixedSigma <- controlvals$sigma > 0 Dims <- list(p = pLen, N = NReal, REML = FALSE) attr(gnlsSt, "conLin") <- list(Xy = array(w, c(NReal, 1), list(row.names(dataModShrunk), deparse(form[[2]]))), dims = Dims, ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions logLik = 0, sigma=controlvals$sigma, auxSigma=0, fixedSigma=fixedSigma) ## additional attributes of gnlsSt attr(gnlsSt, "resp") <- yShrunk attr(gnlsSt, "model") <- modelResid attr(gnlsSt, "local") <- nlEnv attr(gnlsSt, "NReal") <- NReal ## initialization gnlsSt <- Initialize(gnlsSt, dataModShrunk) parMap <- attr(gnlsSt, "pmap") numIter <- 0 # number of iterations nlsSettings <- c(controlvals$nlsMaxIter, controlvals$minScale, controlvals$nlsTol, 0, 0, 0) nlModel <- nonlinModel(modelExpression, nlEnv) repeat { ## alternating algorithm numIter <- numIter + 1 ## GLS step if (needUpdate(gnlsSt)) { # updating varying weights gnlsSt <- update(gnlsSt, dataModShrunk) } if (length(oldPars <- coef(gnlsSt)) > 0) { if (controlvals$opt == "nlminb") { optRes <- nlminb(c(coef(gnlsSt)), function(gnlsPars) -logLik(gnlsSt, gnlsPars), control = list(trace = controlvals$msVerbose, iter.max = controlvals$msMaxIter)) convIter <- optRes$iterations } else { optRes <- optim(c(coef(gnlsSt)), function(gnlsPars) -logLik(gnlsSt, gnlsPars), method = controlvals$optimMethod, control = list(trace = controlvals$msVerbose, maxit = controlvals$msMaxIter, reltol = if(numIter == 0) controlvals$msTol else 100*.Machine$double.eps)) convIter <- optRes$count[2] } aConv <- coef(gnlsSt) <- optRes$par if (verbose) { cat("\n**Iteration", numIter) cat("\n") cat("GLS step: Objective:", format(optRes$value)) print(gnlsSt) } } else { aConv <- oldPars <- NULL } ## NLS step if (is.null(correlation)) { cF <- 1.0 cD <- 1 } else { cF <- corFactor(gnlsSt$corStruct) cD <- Dim(gnlsSt$corStruct) } if (is.null(weights)) { vW <- 1.0 } else { vW <- varWeights(gnlsSt$varStruct) } work <- .C(fit_gnls, thetaNLS = as.double(spar), as.integer(unlist(Dims)), as.double(cF), as.double(vW), as.integer(unlist(cD)), settings = as.double(nlsSettings), additional = double(NReal), as.integer(!is.null(correlation)), as.integer(!is.null(weights)), nlModel, NAOK = TRUE) if (work$settings[4] == 1) { ## convResult <- 2 if (controlvals$returnObject) { warning("step halving factor reduced below minimum in NLS step") } else { stop("step halving factor reduced below minimum in NLS step") } break } oldPars <- c(spar, oldPars) spar[] <- work$thetaNLS if (length(coef(gnlsSt)) == 0 && work$set[5] < controlvals$nlsMaxIter) { break } attr(gnlsSt, "conLin")$Xy[] <- work$additional attr(gnlsSt, "conLin")$logLik <- 0 if (verbose) { cat("\nNLS step: RSS = ", format(work$set[6]), "\n model parameters:") for (i in 1:pLen) cat(format(signif(spar[i]))," ") cat("\n iterations:",work$set[5],"\n") } aConv <- c(spar, aConv) conv <- abs((oldPars - aConv)/ ifelse(abs(aConv) < controlvals$tolerance, 1, aConv)) aConv <- c(max(conv[1:pLen])) names(aConv) <- "params" if (length(conv) > pLen) { conv <- conv[-(1:pLen)] for(i in names(gnlsSt)) { if (any(parMap[,i])) { aConv <- c(aConv, max(conv[parMap[,i]])) names(aConv)[length(aConv)] <- i } } } if (verbose) { cat("\nConvergence:\n") print(aConv) } if ((max(aConv) <= controlvals$tolerance) || (aConv["params"] <= controlvals$tolerance && convIter == 1)) { ## convResult <- 0 break } if (numIter >= controlvals$maxIter) { ## convResult <- 1 if (controlvals$returnObject) { warning("maximum number of iterations reached without convergence") break } else { stop("maximum number of iterations reached without convergence") } } } ## end{ repeat } -------------- ## wraping up ww <- eval(modelExpression[[2]], envir = nlEnv) auxRes <- ww[NReal * pLen + (1:NReal)] attr(gnlsSt, "conLin")$Xy <- array(ww, c(NReal, pLen + 1)) attr(gnlsSt, "conLin") <- c.L <- recalc(gnlsSt) ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions if((sigma <- controlvals$sigma) == 0) { sigma <- sqrt(sum((c.L$Xy[, pLen + 1])^2)/(NReal - pLen)) lsig <- logb(sigma) + 0.5 * logb(1 - pLen/NReal) loglik <- ( - NReal * (1 + logb(2 * pi) + 2 * lsig))/2 + c.L$logLik } else { loglik <- - (NReal * (logb(2 * pi)/2 + logb(sigma)) + sum((c.L$Xy[, pLen + 1])^2) / (2 * sigma^2)) + c.L$logLik lsig <- log(sigma) } ## #### varBeta <- qr(c.L$Xy[ , 1:pLen, drop = FALSE]) if (varBeta$rank < pLen) { ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions print("approximate covariance matrix for parameter estimates not of full rank") return() } attr(parAssign, "varBetaFact") <- varBeta <- sigma * t(backsolve(qr.R(varBeta), diag(pLen))) varBeta <- crossprod(varBeta) dimnames(varBeta) <- list(pn, pn) ## ## fitted.values and residuals (in original order) ## Resid <- resid(gnlsSt) Fitted <- yShrunk - Resid attr(Resid, "std") <- sigma/(varWeights(gnlsSt)) if (!is.null(groups)) { attr(Resid, "std") <- attr(Resid, "std")[revOrderShrunk] Resid[] <- Resid[revOrderShrunk] Fitted[] <- Fitted[revOrderShrunk] grpShrunk[] <- grpShrunk[revOrderShrunk] } names(Resid) <- names(Fitted) <- origOrderShrunk ## getting the approximate var-cov of the parameters ## first making Xy into single column array again attr(gnlsSt, "conLin")$Xy <- array(auxRes, c(NReal, 1)) ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions attr(gnlsSt, "fixedSigma") <- (controlvals$sigma > 0) apVar <- if (controlvals$apVar) gnlsApVar(gnlsSt, lsig, .relStep = controlvals[[".relStep"]], minAbsPar = controlvals[["minAbsParApVar"]]) else "Approximate variance-covariance matrix not available" ## getting rid of condensed linear model and fit oClass <- class(gnlsSt) attributes(gnlsSt) <- attributes(gnlsSt)[!is.na(match(names(attributes(gnlsSt)), ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions c("names","pmap","fixedSigma")))] class(gnlsSt) <- oClass grpDta <- inherits(data, "groupedData") ## ## creating the gnls object ## structure(class = c("gnls", "gls"), list(modelStruct = gnlsSt, dims = Dims, contrasts = contr, coefficients = spar, varBeta = varBeta, ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions sigma = if(controlvals$sigma) controlvals$sigma else sigma, apVar = apVar, logLik = loglik, numIter = numIter, groups = grpShrunk, call = Call, method = "ML", fitted = Fitted, residuals = Resid, plist = plist, pmap = pmap, parAssign = parAssign, na.action = attr(dataMod, "na.action")), ## saving labels and units for plots units = if(grpDta) attr(data, "units"), labels= if(grpDta) attr(data, "labels")) } ### Auxiliary functions used internally in gls and its methods gnlsApVar <- function(gnlsSt, lsigma, conLin = attr(gnlsSt, "conLin"), .relStep = (.Machine$double.eps)^(1/3), minAbsPar = 0, natural = TRUE) { ## calculate approximate variance-covariance matrix of all parameters ## except the coefficients fullGnlsLogLik <- function(Pars, object, conLin, N) { ## logLik as a function of sigma and coef(glsSt) ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions fixedSigma <- attr(object,"fixedSigma") npar <- length(Pars) if (!fixedSigma) { lsigma <- Pars[npar] Pars <- Pars[-npar] } else { lsigma <- log(conLin$sigma) } ####### coef(object) <- Pars conLin <- recalc(object, conLin) conLin[["logLik"]] - N * lsigma - sum(conLin$Xy^2)/(2*exp(2*lsigma)) } fixedSigma <- attr(gnlsSt,"fixedSigma") if (length(gnlsCoef <- coef(gnlsSt)) > 0) { cSt <- gnlsSt[["corStruct"]] if (!is.null(cSt) && inherits(cSt, "corSymm") && natural) { cStNatPar <- coef(cSt, unconstrained = FALSE) class(cSt) <- c("corNatural", "corStruct") coef(cSt) <- log((cStNatPar + 1)/(1 - cStNatPar)) gnlsSt[["corStruct"]] <- cSt gnlsCoef <- coef(gnlsSt) } dims <- conLin$dims N <- dims$N conLin[["logLik"]] <- 0 # making sure ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions Pars <- if(fixedSigma) gnlsCoef else c(gnlsCoef, lSigma = lsigma) # log(sigma) is used as input in contrast to gls val <- fdHess(Pars, fullGnlsLogLik, gnlsSt, conLin, N, .relStep = .relStep, minAbsPar = minAbsPar)[["Hessian"]] if (all(eigen(val, only.values=TRUE)$values < 0)) { ## negative definite - OK val <- solve(-val) ## ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions ## if(fixedSigma && !is.null(dim(val))){ ## Pars <- c(gnlsCoef, lSigma = lsigma) ## npars<-length(Pars) ## val<-cbind(val,rep(0,npars-1)) ## val<-rbind(val,rep(0,npars)) ## } nP <- names(Pars) dimnames(val) <- list(nP, nP) attr(val, "Pars") <- Pars attr(val, "natural") <- natural val } else { ## problem - solution is not maximum "Non-positive definite approximate variance-covariance" } } else { NULL } } ### ### function used to calculate the parameters from ### the params and random effects ### getParsGnls <- function(plist, pmap, beta, N) { pars <- array(0, c(N, length(plist)), list(NULL, names(plist))) for (nm in names(plist)) { pars[, nm] <- if (is.logical(p <- plist[[nm]])) beta[pmap[[nm]]] else p %*% beta[pmap[[nm]]] } pars } ### ### Methods for standard generics ### coef.gnls <- function(object, ...) object$coefficients formula.gnls <- function(x, ...) eval(x$call[["model"]]) getData.gnls <- function(object) { mCall <- object$call data <- eval(mCall$data) if (is.null(data)) return(data) naPat <- eval(mCall$naPattern) if (!is.null(naPat)) { data <- data[eval(naPat[[2]], data), , drop = FALSE] } naAct <- eval(mCall$na.action) if (!is.null(naAct)) { data <- naAct(data) } subset <- mCall$subset if (!is.null(subset)) { subset <- eval(asOneSidedFormula(subset)[[2]], data) data <- data[subset, ] } data } logLik.gnls <- ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions function(object, REML = FALSE, ...) { if (REML) { stop("cannot calculate REML log-likelihood for \"gnls\" objects") } ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions fixSig <- attr(object[["modelStruct"]], "fixedSigma") fixSig <- !is.null(fixSig) && fixSig p <- object$dims$p N <- object$dims$N val <- object[["logLik"]] attr(val, "nobs") <- attr(val, "nall") <- N ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions attr(val, "df") <- p + length(coef(object[["modelStruct"]])) + as.integer(!fixSig) class(val) <- "logLik" val } nobs.gnls <- function(object, ...) object$dims$N predict.gnls <- function(object, newdata, na.action = na.fail, naPattern = NULL, ...) { ## ## method for predict() designed for objects inheriting from class gnls ## if (missing(newdata)) { # will return fitted values return(fitted(object)) } newdata <- data.frame(newdata, check.names = FALSE) mCall <- object$call mfArgs <- list(formula = asOneFormula(formula(object), mCall$params, naPattern, omit = c(names(object$plist), "pi", deparse(getResponseFormula(object)[[2]]))), data = newdata, na.action = na.action, drop.unused.levels = TRUE) dataMod <- do.call("model.frame", mfArgs) ## making sure factor levels are the same as in contrasts contr <- object$contrasts for(i in names(dataMod)) { if (inherits(dataMod[,i], "factor") && !is.null(contr[[i]]) && is.matrix(contr[[i]]) ) { levs <- levels(dataMod[,i]) levsC <- dimnames(contr[[i]])[[1]] if (any(wch <- is.na(match(levs, levsC)))) { stop(sprintf(ngettext(sum(wch), "level %s not allowed for %s", "levels %s not allowed for %s"), paste(levs[wch], collapse = ",")), domain = NA) } attr(dataMod[,i], "contrasts") <- contr[[i]][levs, , drop = FALSE] # if (length(levs) < length(levsC)) { # if (inherits(dataMod[,i], "ordered")) { # dataMod[,i] <- ordered(as.character(dataMod[,i]), levels = levsC) # } else { # dataMod[,i] <- factor(as.character(dataMod[,i]), levels = levsC) # } # } } } N <- nrow(dataMod) ## ## evaluating the naPattern expression, if any ## if (is.null(naPattern)) naPat <- rep(TRUE, N) else naPat <- as.logical(eval(asOneSidedFormula(naPattern)[[2]], dataMod)) ## ## Getting the plist for the new data frame ## ## plist <- object$plist pnames <- names(plist) if (is.null(params <- eval(object$call$params))) { params <- eval(parse(text = paste0(paste(pnames, collapse = "+"), "~ 1"))) } if (!is.list(params)) { params <- list(params) } val <- NULL for(i in seq_along(params)) { if (is.name(params[[i]][[2]])) { val <- c(val, list(params[[i]])) } else { ## multiple parameters on left hand side val <- c(val, eval(parse(text = paste("list(", paste(paste(all.vars(params[[i]][[2]]), deparse(params[[i]][[3]]), sep = "~"), collapse=","),")")))) } } params <- val names(params) <- pnames prs <- coef(object) ## pn <- names(prs) for(nm in pnames) { if (!is.logical(plist[[nm]])) { plist[[nm]] <- model.matrix(asOneSidedFormula(params[[nm]][[3]]), model.frame(asOneSidedFormula(params[[nm]][[3]]), dataMod)) } } modForm <- getCovariateFormula(object)[[2]] val <- eval(modForm, data.frame(dataMod, getParsGnls(plist, object$pmap, prs, N)))[naPat] names(val) <- row.names(newdata) lab <- "Predicted values" if (!is.null(aux <- attr(object, "units")$y)) { lab <- paste(lab, aux) } attr(val, "label") <- lab val } #based on R's update.default update.gnls <- function (object, model., ..., evaluate = TRUE) { call <- object$call if (is.null(call)) stop("need an object with call component") extras <- match.call(expand.dots = FALSE)$... if (!missing(model.)) call$model <- update.formula(formula(object), model.) if(length(extras) > 0) { existing <- !is.na(match(names(extras), names(call))) ## do these individually to allow NULL to remove entries. for (a in names(extras)[existing]) call[[a]] <- extras[[a]] if(any(!existing)) call <- as.call(c(as.list(call), extras[!existing])) } if(evaluate) eval(call, parent.frame()) else call } #update.gnls <- # function(object, model, data = sys.frame(sys.parent()), params, start , # correlation = NULL, weights = NULL, subset, # na.action = na.fail, naPattern, control = list(), # verbose = FALSE, ...) #{ # thisCall <- as.list(match.call())[-(1:2)] # nextCall <- as.list(object$call)[-1] # if (!is.null(thisCall$model)) { # thisCall$model <- update(formula(object), model) # } else { # same model # if (is.null(thisCall$start)) { # thisCall$start <- coef(object) # } # } # if (is.na(match("correlation", names(thisCall))) && # !is.null(thCor <- object$modelStruct$corStruct)) { # thisCall$correlation <- thCor # } # if (is.na(match("weights", names(thisCall))) && # !is.null(thWgt <- object$modelStruct$varStruct)) { # thisCall$weights <- thWgt # } # nextCall[names(thisCall)] <- thisCall # do.call("gnls", nextCall) #} ###*### gnlsStruct - a model structure for gnls fits gnlsStruct <- ## constructor for gnlsStruct objects function(corStruct = NULL, varStruct = NULL) { val <- list(corStruct = corStruct, varStruct = varStruct) val <- val[!sapply(val, is.null)] # removing NULL components # attr(val, "settings") <- attr(val$reStruct, "settings") # attr(val, "resp") <- resp # attr(val, "model") <- model # attr(val, "local") <- local # attr(val, "N") <- N # attr(val, "naPat") <- naPat class(val) <- c("gnlsStruct", "glsStruct", "modelStruct") val } ##*## gnlsStruct methods for standard generics fitted.gnlsStruct <- function(object, ...) attr(object, "resp") - resid(object) Initialize.gnlsStruct <- function(object, data, ...) { if (length(object)) { object[] <- lapply(object, Initialize, data) theta <- lapply(object, coef) len <- lengths(theta) num <- seq_along(len) pmap <- if (sum(len) > 0) outer(rep(num, len), num, "==") else array(FALSE, c(1, length(len))) dimnames(pmap) <- list(NULL, names(object)) attr(object, "pmap") <- pmap if (needUpdate(object)) object <- update(object, data) } object } logLik.gnlsStruct <- function(object, Pars, conLin = attr(object, "conLin"), ...) { coef(object) <- Pars # updating parameter values conLin <- recalc(object, conLin) ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions if(conLin$sigma == 0) { conLin[["logLik"]] - (conLin$dims$N * logb(sum(conLin$Xy^2)))/2 } else { conLin[["logLik"]] - conLin$dims$N * logb(conLin$sigma) - sum(conLin$Xy^2) / (2 * conLin$sigma^2) } } residuals.gnlsStruct <- function(object, ...) { c(eval(attr(object, "model")[[2]], envir = attr(object, "local"))) } gnlsControl <- ## Set control values for iterations within gnls function(maxIter = 50, nlsMaxIter = 7, msMaxIter = 50, minScale = 0.001, tolerance = 1e-6, nlsTol = 0.001, msTol = 1e-7, returnObject = FALSE, msVerbose = FALSE, apVar = TRUE, .relStep = .Machine$double.eps^(1/3), opt = c("nlminb", "optim"), optimMethod = "BFGS", ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions minAbsParApVar = 0.05, sigma=NULL) { ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions if(is.null(sigma)) sigma <- 0 else if(!is.finite(sigma) || length(sigma) != 1 || sigma < 0) stop("Within-group std. dev. must be a positive numeric value") list(maxIter = maxIter, nlsMaxIter = nlsMaxIter, msMaxIter = msMaxIter, minScale = minScale, tolerance = tolerance, nlsTol = nlsTol, msTol = msTol, returnObject = returnObject, msVerbose = msVerbose, apVar = apVar, opt = match.arg(opt), optimMethod = optimMethod, ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions .relStep = .relStep, minAbsParApVar = minAbsParApVar, sigma=sigma) } ## Local Variables: ## ess-indent-offset: 2 ## End: