### Fit a general nonlinear mixed effects model ### ### Copyright 1997-2003 Jose C. Pinheiro, ### Douglas M. Bates ### Copyright 2006-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/ nlme <- function(model, data = sys.frame(sys.parent()), fixed, random = fixed, groups, start, correlation = NULL, weights = NULL, subset, method = c("ML", "REML"), na.action = na.fail, naPattern, control = list(), verbose= FALSE) { UseMethod("nlme") } nlme.nlsList <- function(model, data = sys.frame(sys.parent()), fixed, random = fixed, groups, start, correlation = NULL, weights = NULL, subset, method = c("ML", "REML"), na.action = na.fail, naPattern, control = list(), verbose= FALSE) { ## control parameters controlvals <- nlmeControl() controlvals[names(control)] <- control thisCall <- as.list(match.call())[-1] ## checking the use of arguments defined within the function if (any(!is.na(match(names(thisCall), c("fixed", "data", "start"))))) { warning("'nlme.nlsList' will redefine 'fixed', 'data', and 'start'") } method <- match.arg(method) REML <- method == "REML" ## add model, data, and optionally groups from the call that created model last.call <- as.list(attr(model, "call"))[-1] last.call$control <- NULL last.call$pool <- NULL thisCall[names(last.call)] <- last.call thisModel <- last.call[["model"]] thisCall[["model"]] <- eval(parse(text=paste( deparse (getResponseFormula (thisModel)[[2]]), c_deparse(getCovariateFormula(thisModel)[[2]]),sep="~"))) ## create "fixed" and "start" cf <- na.omit(coef(model)) start <- list(fixed = unlist(lapply(cf, median, na.rm = TRUE))) pnames <- names(start$fixed) <- names(cf) thisCall[["fixed"]] <- lapply(as.list(pnames), function(el) eval(parse(text = paste(el, 1, sep = "~")))) if (missing(random)) { random <- thisCall[["fixed"]] } reSt <- reStruct(random, data = NULL) if (missing(groups)) { thisCall[["groups"]] <- groups <- getGroupsFormula(model) } if (length(reSt) > 1 || length(groups[[2]]) > 1) { stop("can only fit \"nlsList\" objects with single grouping variable") } ranForm <- formula(reSt)[[1]] if (!is.list(ranForm)) { ranForm <- list(ranForm) } mData <- thisCall[["data"]] if (is.null(mData)) { # will try to construct allV <- unique(unlist(lapply(ranForm, function(el) all.vars(el[[3]])))) if (length(allV) > 0) { alist <- lapply(as.list(allV), as.name) names(alist) <- allV alist <- c(as.list(quote(data.frame)), alist) mode(alist) <- "call" mData <- eval(alist, sys.parent(1)) } } else if (mode(mData) == "name" || mode(mData) == "call") { mData <- eval(mData) } reSt <- reStruct(random, REML = REML, data = mData) names(reSt) <- deparse(groups[[2]]) ## convert list of "name" objects to "character" vector rnames <- sapply(lapply(ranForm, `[[`, 2L), deparse) ## if the random effects are a subset of the coefficients, ## construct initial estimates for their var-cov matrix if (all(match(rnames, pnames, 0))) { madRes <- mad(resid(model), na.rm = TRUE) madRan <- unlist(lapply(cf, mad, na.rm = TRUE)) madRan <- madRan[rnames] if (isInitialized(reSt)) { warning("initial value for 'reStruct' overwritten in 'nlme.nlsList'") } matrix(reSt) <- diag((madRan/madRes)^2, ncol = length(rnames)) } thisCall[["start"]] <- start thisCall[["random"]] <- reSt val <- do.call(nlme.formula, thisCall) val$origCall <- match.call() val } nlme.formula <- function(model, data = sys.frame(sys.parent()), fixed, random, groups, start, correlation = NULL, weights = NULL, subset, method = c("ML", "REML"), na.action = na.fail, naPattern, control = list(), verbose= FALSE) { ## This is the method that actually does the fitting 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() ## control parameters controlvals <- nlmeControl() if (!missing(control)) { controlvals[names(control)] <- control } ## ## checking arguments ## if (!inherits(model, "formula")) stop("'model' must be a formula") if (length(model)!=3) stop("model formula must be of the form \"resp ~ pred\"") method <- match.arg(method) REML <- method == "REML" if (missing(random)) { random <- fixed } reSt <- reStruct(random, REML = REML, data = NULL) if (missing(groups)) { groups <- getGroupsFormula(reSt) } if (is.null(groups)) { if (inherits(data, "groupedData")) { groups <- getGroupsFormula(data) namGrp <- rev(names(getGroupsFormula(data, asList = TRUE))) Q <- length(namGrp) if (length(reSt) != Q) { # may need to repeat reSt if (length(reSt) != 1) stop("incompatible lengths for 'random' and grouping factors") randL <- vector("list", Q) names(randL) <- rev(namGrp) for(i in 1:Q) randL[[i]] <- random reSt <- reStruct(randL, REML = REML, data = NULL) } } else { ## will assume single group groups <- ~ 1 names(reSt) <- namGrp <- "1" } } else { g.exp <- eval(parse(text = paste0("~1 |", deparse(groups[[2]])))) namGrp <- rev(names(getGroupsFormula(g.exp, asList = TRUE))) } names(reSt) <- namGrp ## ## checking if self-starting formula is given ## if (missing(start) && !is.null(attr(eval(model[[3]][[1]]), "initial"))) { nlmeCall <- Call nlsLCall <- nlmeCall[c("","model","data","groups")] nlsLCall[[1]] <- quote(nlme::nlsList) names(nlsLCall)[2] <- "model" nm <- names(nlmeCall) for(i in c("fixed", "data", "groups", "start")) if(i %in% nm) nlmeCall[[i]] <- NULL nlmeCall[[1]] <- quote(nlme::nlme.nlsList) ## checking if "data" is not equal to sys.frame(sys.parent()) if (is.null(dim(data))) { stop("'data' must be given explicitly to use 'nlsList'") } nlsLObj <- eval(nlsLCall) nlmeCall[["model"]] <- quote(nlsLObj) nlmeCall <- as.call(nlmeCall) val <- eval(nlmeCall) val$origCall <- NULL return(val) } if (is.numeric(start)) { # assume it is the fixed effects start <- list(fixed = start) } nlmeModel <- call("-", model[[2]], model[[3]]) ## ## save writing list(...) when only one element ## if (!is.list(fixed)) fixed <- list(fixed) fixed <- do.call(c, lapply(fixed, function(fix.i) { if (is.name(fix.i[[2]])) list(fix.i) else ## multiple parameters on left hand side eval(parse(text = paste0("list(", paste(paste(all.vars(fix.i[[2]]), deparse (fix.i[[3]]), sep = "~"), collapse = ","), ")"))) })) fnames <- lapply(fixed, function(fix.i) { this <- eval(fix.i) if (!inherits(this, "formula")) stop ("'fixed' must be a formula or list of formulae") if (length(this) != 3) stop ("formulae in 'fixed' must be of the form \"parameter ~ expr\"") if (!is.name(this[[2]])) stop ("formulae in 'fixed' must be of the form \"parameter ~ expr\"") as.character(this[[2]]) }) names(fixed) <- fnames ranForm <- formula(reSt) # random effects formula(s) Q <- length(ranForm) # number of groups names(ranForm) <- namGrp rnames <- vector("list", Q) names(rnames) <- namGrp for(i in 1:Q) { rnames[[i]] <- character(length(ranForm[[i]])) for (j in seq_along(ranForm[[i]])) { this <- eval(ranForm[[i]][[j]]) if (!inherits(this, "formula")) stop ("'random' must be a formula or list of formulae") if (length(this) != 3) stop ("formulae in 'random' must be of the form \"parameter ~ expr\"") if (!is.name(this[[2]])) stop ("formulae in 'random' must be of the form \"parameter ~ expr\"") rnames[[i]][j] <- deparse(this[[2]]) } names(ranForm[[i]]) <- rnames[[i]] } ## all parameter names pnames <- unique(c(fnames, unlist(rnames))) ## ## 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 corStruct is present and assign groups to its formula, ## if necessary if (!is.null(correlation)) { if(!is.null(corGrpsForm <- getGroupsFormula(correlation, asList = TRUE))) { corGrpsForm <- unlist(lapply(corGrpsForm, function(el) deparse(el[[2]]))) corQ <- length(corGrpsForm) lmeGrpsForm <- unlist(lapply(splitFormula(groups), function(el) deparse(el[[2]]))) lmeQ <- length(lmeGrpsForm) if (corQ <= lmeQ) { if (any(corGrpsForm != lmeGrpsForm[1:corQ])) { stop("incompatible formulas for groups in \"random\" and \"correlation\"") } if (corQ < lmeQ) { warning("cannot use smaller level of grouping for \"correlation\" than for \"random\". Replacing the former with the latter.") frm <- paste("~", c_deparse(getCovariateFormula(formula(correlation))[[2]]), "|", deparse(groups[[2]])) attr(correlation, "formula") <- eval(parse(text = frm)) } } else { if (any(lmeGrpsForm != corGrpsForm[1:lmeQ])) { stop("incompatible formulas for groups in \"random\" and \"correlation\"") } } } else { ## using the same grouping as in random frm <- paste("~", c_deparse(getCovariateFormula(formula(correlation))[[2]]), "|", deparse(groups[[2]])) attr(correlation, "formula") <- eval(parse(text = frm)) corQ <- lmeQ <- 1 } } else { corQ <- lmeQ <- 1 } ## create an nlme structure containing the random effects model and plug-ins nlmeSt <- nlmeStruct(reStruct = reSt, corStruct = correlation, varStruct = varFunc(weights)) ## extract a data frame with enough information to evaluate ## form, fixed, random, groups, correlation, and weights mfArgs <- list(formula = asOneFormula(formula(nlmeSt), model, fixed, 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 dataMix <- do.call(model.frame, mfArgs) origOrder <- row.names(dataMix) # preserve the original order ## ## Evaluating the groups expression ## grps <- getGroups(dataMix, eval(parse(text = paste("~1", deparse(groups[[2]]), sep = "|")))) N <- dim(dataMix)[1] # number of observations ## ## evaluating the naPattern expression, if any ## if (missing(naPattern)) naPat <- rep(TRUE, N) else naPat <- as.logical(eval(asOneSidedFormula(naPattern)[[2]], dataMix)) origOrderShrunk <- origOrder[naPat] ## ordering data by groups if (inherits(grps, "factor")) { # single level ord <- order(grps) #"order" treats a single named argument peculiarly grps <- data.frame(grps) row.names(grps) <- origOrder names(grps) <- as.character(deparse((groups[[2]]))) } else { ord <- do.call(order, grps) ## making group levels unique for(i in 2:ncol(grps)) { grps[, i] <- as.factor(paste(as.character(grps[, i-1]), as.character(grps[,i]), sep = "/")) } } if (corQ > lmeQ) { ## may have to reorder by the correlation groups ord <- do.call(order, getGroups(dataMix, getGroupsFormula(correlation))) } grps <- grps[ord, , drop = FALSE] dataMix <- dataMix[ord, ,drop = FALSE] ## revOrder <- match(origOrder, row.names(dataMix)) # putting in orig. order naPat <- naPat[ord] # ordered naPat dataMixShrunk <- dataMix[naPat, , drop=FALSE] ## ordShrunk <- ord[naPat] grpShrunk <- grps[naPat,, drop = FALSE] revOrderShrunk <- match(origOrderShrunk, row.names(dataMixShrunk)) yShrunk <- eval(model[[2]], dataMixShrunk) ## ## defining list with parameter information ## contr <- list() plist <- vector("list", length(pnames)) names(plist) <- pnames for (nm in pnames) { this <- list(fixed = !is.null(fixed[[nm]]), random = as.list(lapply(ranForm, function(el, nm) !is.null(el[[nm]]), nm = nm))) if (this[["fixed"]]) { ## Peter Dalgaard claims the next line should be this[["fixed"]][[3]][[1]] != "1" ## but the current version seems to work ok. if (fixed[[nm]][[3]] != "1") { as1F <- asOneSidedFormula(fixed[[nm]][[3]]) this[["fixed"]] <- model.matrix(as1F, model.frame(as1F, dataMix)) auxContr <- attr(this[["fixed"]], "contrasts") contr <- c(contr, auxContr[is.na(match(names(auxContr), names(contr)))]) } } if (any(unlist(this[["random"]]))) { for(i in 1:Q) { wch <- which(!is.na(match(rnames[[i]], nm))) if (length(wch) == 1) { # only one formula for nm at level i if ((rF.i <- ranForm[[i]][[nm]][[3]]) != "1") { this[["random"]][[i]] <- model.matrix(asOneSidedFormula(rF.i), model.frame(asOneSidedFormula(rF.i), dataMix)) auxContr <- attr(this[["random"]][[i]], "contrasts") contr <- c(contr, auxContr[is.na(match(names(auxContr), names(contr)))]) } } else if (length(wch) > 0) { # multiple formulas this[["random"]][[i]] <- th.ran.i <- as.list(lapply(ranForm[[i]][wch], function(el, data) { if (el[[3]] == "1") TRUE else model.matrix(asOneSidedFormula(el[[3]]), model.frame(asOneSidedFormula(el[[3]]), data)) }, data = dataMix)) for(j in seq_along(th.ran.i)) { if (is.matrix(th.ran.i[[j]])) { auxContr <- attr(th.ran.i[[j]], "contrasts") contr <- c(contr, auxContr[is.na(match(names(auxContr), names(contr)))]) } } } } } plist[[nm]] <- this } ## Ensure that all elements of are matrices contrMat <- function(nm, contr, data) { ## nm is a term in a formula, and can be a call x <- eval(parse(text=nm), data) levs <- levels(x) val <- do.call(contr[[nm]], list(n = length(levs))) rownames(val) <- levs val } nms <- names(contr)[sapply(contr, is.character)] contr[nms] <- lapply(nms, contrMat, contr = contr, data = dataMix) if (is.null(sfix <- start$fixed)) stop ("'start' must have a component called 'fixed'") ## ## Fixed effects names ## fn <- character(0) currPos <- 0 fixAssign <- list() for(nm in fnames) { if (is.logical(f <- plist[[nm]]$fixed)) { currPos <- currPos + 1 currVal <- list(currPos) if (all(unlist(lapply(plist[[nm]]$random, is.logical)))) { fn <- c(fn, nm) names(currVal) <- nm } else { aux <- paste(nm, "(Intercept)", sep=".") fn <- c(fn, aux) names(currVal) <- aux } fixAssign <- c(fixAssign, currVal) } else { currVal <- attr(f, "assign") fTerms <- terms(asOneSidedFormula(fixed[[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 = ".") fixAssign <- c(fixAssign, lapply(currVal, function(el) el + currPos )) currPos <- currPos + length(unlist(currVal)) fn <- c(fn, paste(nm, colnames(f), sep = ".")) } } fLen <- length(fn) if (length(sfix) != fLen) stop ("starting values for the 'fixed' component are not the correct length") names(sfix) <- fn ## ## Random effects names ## rn <- wchRnames <- vector("list", Q) names(rn) <- names(wchRnames) <- namGrp for(i in 1:Q) { rn[[i]] <- character(0) uRnames <- unique(rnames[[i]]) wchRnames[[i]] <- integer(length(uRnames)) names(wchRnames[[i]]) <- uRnames for(j in seq_along(rnames[[i]])) { nm <- rnames[[i]][j] wchRnames[[i]][nm] <- wchRnames[[i]][nm] + 1 r <- plist[[nm]]$random[[i]] if (data.class(r) == "list") r <- r[[wchRnames[[i]][nm]]] if (is.logical(r)) { if (r) { rn[[i]] <- c(rn[[i]], if (is.logical(plist[[nm]]$fixed)) nm else paste(nm,"(Intercept)",sep=".")) } # else: keep rn[[i]] } else { rn[[i]] <- c(rn[[i]], paste(nm, colnames(r), sep = ".")) } } } Names(nlmeSt$reStruct) <- rn rNam <- unlist(rn) # unlisted names of random effects rlength <- lengths(rn) # number of random effects per stratum rLen <- sum(rlength) # total number of random effects pLen <- rLen + fLen # total number of parameters ncols <- c(rlength, fLen, 1) Dims <- MEdims(grpShrunk, ncols) if (max(Dims$ZXlen[[1]]) < Dims$qvec[1]) { warning(gettextf("fewer observations than random effects in all level %s groups", Q), domain = NA) } sran <- vector("list", Q) names(sran) <- namGrp if (!is.null(sran0 <- start$random)) { if (inherits(sran0, "data.frame")) { sran0 <- list(as.matrix(sran0)) } else { if (!is.list(sran0)) { if (!is.matrix(sran0)) { stop("starting values for random effects should be a list, or a matrix") } sran0 <- list(as.matrix(sran0)) } } if (is.null(namSran <- names(sran0))) { if (length(sran) != Q) { stop(gettextf("list with starting values for random effects must have names or be of length %d", Q), domain = NA) } names(sran0) <- rev(namGrp) # assume given in outer-inner order } else { if (any(noMatch <- is.na(match(namSran, namGrp)))) { stop(sprintf(ngettext(sum(noMatch), "group name not matched in starting values for random effects: %s", "group names not matched in starting values for random effects: %s"), paste(namSran[noMatch], collapse=", ")), domain = NA) } } } for(i in 1:Q) { if (is.null(sran[[i]] <- sran0[[namGrp[i]]])) { sran[[i]] <- array(0, c(rlength[i], Dims$ngrps[i]), list(rn[[i]], unique(as.character(grps[, Q-i+1])))) } else { if (!is.matrix(sran[[i]])) stop("starting values for the random components should be a list of matrices") dimsran <- dim(sran[[i]]) if (dimsran[1] != Dims$ngrps[i]) { stop(gettextf("number of rows in starting values for random component at level %s should be %d", namGrp[i], Dims$ngrps[i]), domain = NA) } if (dimsran[2] != rlength[i]) { stop(gettextf("number of columns in starting values for random component at level %s should be %d", namGrp[i], rlength[i]), domain = NA) } dnamesran <- dimnames(sran[[i]]) if (is.null(dnamesran[[1]])) { stop("starting values for random effects must include group levels") } else { levGrps <- unique(as.character(grps[, Q-i+1])) if(!all(sort(dnamesran[[1]]) == sort(levGrps))) { stop(gettextf("groups levels mismatch in 'random' and starting values for 'random' at level %s", namGrp[i]), domain = NA) } sran[[i]] <- sran[[i]][levGrps, , drop = FALSE] } if (!is.null(dnamesran[[2]])) { if(!all(sort(dnamesran[[2]]) == sort(rn[[i]]))) { ## first try to resolve it for(j in seq_len(rlength[i])) { if (is.na(match(dnamesran[[2]][j], rn[[i]]))) { if (!is.na(mDn <- match(paste(dnamesran[[2]][j], "(Intercept)", sep="."), rn[[i]]))) { dnamesran[[2]][j] <- rn[[i]][mDn] } else { if (!is.na(mDn <- match(dnamesran[[2]][j], paste(rn[[i]], "(Intercept)", sep = ".")))) { dnamesran[[2]][j] <- rn[[i]][mDn] } else { stop (gettextf("names mismatch in 'random' and starting values for 'random' at level %s", namGrp[i]), domain = NA) } } } } dimnames(sran[[i]]) <- dnamesran } sran[[i]] <- sran[[i]][, rn[[i]], drop = FALSE] } else { dimnames(sran[[i]])[[2]] <- rn[[i]] } sran[[i]] <- t(sran[[i]]) } } names(sran) <- namGrp nPars <- length(unlist(sran)) + fLen # total number of PNLS parameters ## ## defining values of constants used in calculations ## NReal <- sum(naPat) ## ## Creating the fixed and random effects maps ## fmap <- list() n1 <- 1 for(nm in fnames) { if (is.logical(f <- plist[[nm]]$fixed)) { fmap[[nm]] <- n1 n1 <- n1 + 1 } else { fmap[[nm]] <- n1:(n1+ncol(f) - 1) n1 <- n1 + ncol(f) } } rmap <- rmapRel <- vector("list", Q) names(rmap) <- names(rmapRel) <- namGrp n1 <- 1 startRan <- 0 for(i in 1:Q) { wchRnames[[i]][] <- 0 rmap[[i]] <- rmapRel[[i]] <- list() for(nm in rnames[[i]]) { wchRnames[[i]][nm] <- wchRnames[[i]][nm] + 1 r <- plist[[nm]]$random[[i]] if (data.class(r) == "list") { r <- r[[wchRnames[[i]][nm]]] } if (is.logical(r)) { val <- n1 n1 <- n1 + 1 } else { val <- n1:(n1+ncol(r) - 1) n1 <- n1 + ncol(r) } if (is.null(rmap[[i]][[nm]])) { rmap[[i]][[nm]] <- val rmapRel[[i]][[nm]] <- val - startRan } else { rmap[[i]][[nm]] <- c(rmap[[i]][[nm]], list(val)) rmapRel[[i]][[nm]] <- c(rmapRel[[i]][[nm]], list(val - startRan)) } } startRan <- startRan + ncols[i] } ## ## defining the nlFrame, i.e., nlEnv, an environment in R : ## grpsRev <- rev(lapply(grps, as.character)) bmap <- c(0, cumsum(unlist(lapply(sran, function(el) length(as.vector(el)))))) nlEnv <- list2env( list(model = nlmeModel, data = dataMix, groups = grpsRev, plist = plist, beta = as.vector(sfix), bvec = unlist(sran), b = sran, X = array(0, c(N, fLen), list(NULL, fn)), Z = array(0, c(N, rLen), list(NULL, rNam)), fmap = fmap, rmap = rmap, rmapRel = rmapRel, bmap = bmap, level = Q, N = N, Q = Q, naPat = naPat, .parameters = c("bvec", "beta"), finiteDiffGrad = finiteDiffGrad)) modelExpression <- ~ { pars <- getParsNlme(plist, fmap, rmapRel, bmap, groups, beta, bvec, b, level, N) res <- eval(model, data.frame(data, pars)) if (!length(grad <- attr(res, "gradient"))) { grad <- finiteDiffGrad(model, data, pars) } for (nm in names(plist)) { gradnm <- grad[, nm] if (is.logical(f <- plist[[nm]]$fixed)) { if(f) X[, fmap[[nm]]] <- gradnm # else f == FALSE =^= 0 } else X[, fmap[[nm]]] <- gradnm * f for(i in 1:Q) { if (is.logical(r <- plist[[nm]]$random[[i]])) { if (r) Z[, rmap[[i]][[nm]]] <- gradnm # else r == FALSE =^= 0 } else { rm.i <- rmap[[i]][[nm]] if (data.class(rm.i) != "list") { Z[, rm.i] <- gradnm * r } else { for(j in seq_along(rm.i)) { Z[, rm.i[[j]]] <- if (is.logical(rr <- r[[j]])) gradnm else gradnm * rr } } } } } result <- c(Z[naPat, ], X[naPat, ], res[naPat]) result[is.na(result)] <- 0 result }## {modelExpression} modelResid <- ~ eval(model, data.frame(data, getParsNlme(plist, fmap, rmapRel, bmap, groups, beta, bvec, b, level, N)))[naPat] ww <- eval(modelExpression[[2]], envir = nlEnv) w <- ww[NReal * pLen + (1:NReal)] ZX <- array(ww[1:(NReal*pLen)], c(NReal, pLen), list(row.names(dataMixShrunk), c(rNam, fn))) w <- w + as.vector(ZX[, rLen + (1:fLen), drop = FALSE] %*% sfix) if (!is.null(start$random)) { startRan <- 0 for(i in 1:Q) { w <- w + as.vector((ZX[, startRan + 1:ncols[i], drop = FALSE] * t(sran[[i]])[as.character(grpShrunk[, Q-i+1]),,drop = FALSE]) %*% rep(1, ncols[i])) startRan <- startRan + ncols[i] } } ## creating the condensed linear model attr(nlmeSt, "conLin") <- list(Xy = array(c(ZX, w), dim = c(NReal, sum(ncols)), dimnames = list(row.names(dataMixShrunk), c(colnames(ZX), deparse(model[[2]])))), dims = Dims, logLik = 0, ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions sigma = controlvals$sigma, auxSigma = 0) ## additional attributes of nlmeSt attr(nlmeSt, "resp") <- yShrunk attr(nlmeSt, "model") <- modelResid attr(nlmeSt, "local") <- nlEnv attr(nlmeSt, "NReal") <- NReal ## initialization nlmeSt <- Initialize(nlmeSt, dataMixShrunk, grpShrunk, control = controlvals) parMap <- attr(nlmeSt, "pmap") decomp <- length(coef(nlmeSt)) == length(coef(nlmeSt$reStruct)) && !needUpdate(nlmeSt) if(decomp) { # can do one decomposition ## need to save conLin for calculating updating between steps oldConLin <- attr(nlmeSt, "conLin") } numIter <- 0 # number of iterations pnlsSettings <- c(controlvals$pnlsMaxIter, controlvals$minScale, controlvals$pnlsTol, 0, 0, 0) nlModel <- nonlinModel(modelExpression, nlEnv) repeat { ## alternating algorithm numIter <- numIter + 1 ## LME step if (needUpdate(nlmeSt)) { # updating varying weights nlmeSt <- update(nlmeSt, dataMixShrunk) } if (decomp) { attr(nlmeSt, "conLin") <- MEdecomp(oldConLin) } oldPars <- coef(nlmeSt) if (controlvals$opt == "nlminb") { control <- list(trace = controlvals$msVerbose, iter.max = controlvals$msMaxIter) keep <- c("eval.max", "abs.tol", "rel.tol", "x.tol", "xf.tol", "step.min", "step.max", "sing.tol", "scale.init", "diff.g") control <- c(control, controlvals[names(controlvals) %in% keep]) optRes <- nlminb(c(coef(nlmeSt)), function(nlmePars) -logLik(nlmeSt, nlmePars), control = control) aConv <- coef(nlmeSt) <- optRes$par convIter <- optRes$iterations } else { ## nlm(.) aNlm <- nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars), p = c(coef(nlmeSt)), hessian = TRUE, print.level = controlvals$msVerbose, gradtol = if(numIter == 1) controlvals$msTol else 100*.Machine$double.eps, iterlim = if(numIter < 10) 10 else controlvals$msMaxIter, check.analyticals = FALSE) aConv <- coef(nlmeSt) <- aNlm$estimate convIter <- aNlm$iterations } nlmeFit <- attr(nlmeSt, "lmeFit") <- MEestimate(nlmeSt, grpShrunk) if (verbose) { cat("\n**Iteration", numIter) cat(sprintf("\nLME step: Loglik: %s, %s iterations: %d\n", format(nlmeFit$logLik), controlvals$opt, convIter)) print(nlmeSt) } ## PNLS step if (is.null(correlation)) { cF <- 1.0 cD <- 1 } else { cF <- corFactor(nlmeSt$corStruct) cD <- Dim(nlmeSt$corStruct) } vW <- if(is.null(weights)) 1.0 else varWeights(nlmeSt$varStruct) work <- .C(fit_nlme, thetaPNLS = as.double(c(as.vector(unlist(sran)), sfix)), pdFactor = as.double(pdFactor(nlmeSt$reStruct)), as.integer(unlist(rev(grpShrunk))), as.integer(unlist(Dims)), as.integer(attr(nlmeSt$reStruct, "settings"))[-(1:3)], as.double(cF), as.double(vW), as.integer(unlist(cD)), settings = as.double(pnlsSettings), additional = double(NReal * (pLen + 1)), as.integer(!is.null(correlation)), as.integer(!is.null(weights)), ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions as.double(controlvals$sigma), nlModel, NAOK = TRUE) if (work$settings[4] == 1) { ## convResult <- 2 msg <- "step halving factor reduced below minimum in PNLS step" if (controlvals$returnObject) warning(msg) else stop(msg) } ## dim(work$pdFactor) <- dim(pdMatrix(nlmeSt$reStruct[[1]])) ## matrix(nlmeSt$reStruct[[1]]) <- crossprod(work$pdFactor) ## fix from Setzer.Woodrow@epamail.epa.gov for nested grouping factors i.pdF <- 1 for (i in seq_along(nlmeSt$reStruct)) { d.i <- dim(pdMatrix(nlmeSt$reStruct[[i]])) ni.pdF <- i.pdF + prod(d.i) pdF <- array(work$pdFactor[i.pdF:(ni.pdF -1)], dim = d.i) matrix(nlmeSt$reStruct[[i]]) <- crossprod(pdF) i.pdF <- ni.pdF } oldPars <- c(sfix, oldPars) for(i in 1:Q) sran[[i]][] <- work$thetaPNLS[(bmap[i]+1):bmap[i+1]] sfix[] <- work$thetaPNLS[nPars + 1 - (fLen:1)] if (verbose) { cat("PNLS step: RSS = ", format(work$settings[6]), "\n fixed effects: ") for (i in 1:fLen) cat(format(sfix[i])," ") cat("\n iterations:", work$settings[5],"\n") } aConv <- c(sfix, coef(nlmeSt)) # 2nd part added by SDR 04/19/2002 w[] <- work$additional[(NReal * pLen) + 1:NReal] ZX[] <- work$additional[1:(NReal * pLen)] w <- w + as.vector(ZX[, rLen + (1:fLen), drop = FALSE] %*% sfix) startRan <- 0 for(i in 1:Q) { w <- w + as.vector((ZX[, startRan + 1:ncols[i], drop = FALSE] * t(sran[[i]])[as.character(grpShrunk[, Q-i+1]),,drop = FALSE]) %*% rep(1, ncols[i])) startRan <- startRan + ncols[i] } if (decomp) { oldConLin$Xy[] <- c(ZX, w) oldConLin$logLik <- 0 } else { attr(nlmeSt, "conLin")$Xy[] <- c(ZX, w) attr(nlmeSt, "conLin")$logLik <- 0 } conv <- abs((oldPars - aConv)/ ifelse(abs(aConv) < controlvals$tolerance, 1, aConv)) aConv <- c(fixed = max(conv[1:fLen])) conv <- conv[-(1:fLen)] for(i in names(nlmeSt)) { if (any(parMap[,i])) { aConv <- c(aConv, max(conv[parMap[,i]])) names(aConv)[length(aConv)] <- i } } if (verbose) { cat(sprintf("Convergence crit. (must all become <= tolerance = %g):\n", controlvals$tolerance)) print(aConv) } if ((max(aConv) <= controlvals$tolerance) || (aConv["fixed"] <= controlvals$tolerance && convIter == 1)) { ## convResult <- 0 break } if (numIter >= controlvals$maxIter) { ## convResult <- 1 msg <- gettextf( "maximum number of iterations (maxIter = %d) reached without convergence", controlvals$maxIter) if (controlvals$returnObject) { warning(msg, domain=NA) break } else { stop(msg, domain=NA) } } } ## wrapping up nlmeFit <- if (decomp) MEestimate(nlmeSt, grpShrunk, oldConLin) else MEestimate(nlmeSt, grpShrunk) ## degrees of freedom for fixed effects tests fixDF <- getFixDF(ZX[, rLen + (1:fLen), drop = FALSE], grpShrunk, attr(nlmeSt, "conLin")$dims$ngrps, fixAssign) ## 17-11-2015; Fixed sigma patch; SH Heisterkamp; Quantitative Solutions attr(fixDF, "varFixFact") <- varFix <- nlmeFit$sigma * nlmeFit$varFix varFix <- crossprod(varFix) dimnames(varFix) <- list(fn, fn) ## ## fitted.values and residuals (in original order) ## Resid <- if (decomp) resid(nlmeSt, level = 0:Q, oldConLin)[revOrderShrunk, ] else resid(nlmeSt, level = 0:Q)[revOrderShrunk, ] Fitted <- yShrunk[revOrderShrunk] - Resid rownames(Resid) <- rownames(Fitted) <- origOrderShrunk grpShrunk <- grpShrunk[revOrderShrunk, , drop = FALSE] attr(Resid, "std") <- nlmeFit$sigma/(varWeights(nlmeSt)[revOrderShrunk]) ## inverting back reStruct nlmeSt$reStruct <- solve(nlmeSt$reStruct) attr(nlmeSt, "fixedSigma") <- (controlvals$sigma > 0) ## saving part of dims dims <- attr(nlmeSt, "conLin")$dims[c("N", "Q", "qvec", "ngrps", "ncol")] ## getting the approximate var-cov of the parameters apVar <- if (controlvals$apVar) lmeApVar(nlmeSt, nlmeFit$sigma, .relStep = controlvals[[".relStep"]], minAbsPar = controlvals[["minAbsParApVar"]], natural = controlvals[["natural"]]) else "Approximate variance-covariance matrix not available" ## putting sran in the right format sran <- lapply(sran, t) ## getting rid of condensed linear model, fit, and other attributes ###- oClass <- class(nlmeSt) attributes(nlmeSt) <- attributes(nlmeSt)[ c("names", "class", "pmap", "fixedSigma")] ##- class(nlmeSt) <- oClass ## ## creating the nlme object ## isGrpd <- inherits(data, "groupedData") structure(class = c("nlme", "lme"), list(modelStruct = nlmeSt, dims = dims, contrasts = contr, coefficients = list(fixed = sfix, random = rev(sran)), varFix = varFix, sigma = nlmeFit$sigma, apVar = apVar, logLik = nlmeFit$logLik, numIter = numIter, groups = grpShrunk, call = Call, method = method, fitted = Fitted, residuals = Resid, plist = plist, map = list(fmap=fmap,rmap=rmap,rmapRel=rmapRel,bmap=bmap), fixDF = fixDF), ## saving labels and units for plots units = if(isGrpd) attr(data, "units"), labels = if(isGrpd) attr(data, "labels")) } ## {nlme.formula} ### ##' Calculate the parameters from the fixed and random effects getParsNlme <- function(plist, fmap, rmapRel, bmap, groups, beta, bvec, b, level, N) { pars <- array(0, c(N, length(plist)), list(NULL, names(plist))) ## for random effects below iQ <- if (level > 0) { Q <- length(groups) (Q - level + 1L):Q } else integer() # empty for (nm in names(plist)) { ## 1) Fixed effects if (is.logical(f <- plist[[nm]]$fixed)) { if (f) pars[, nm] <- beta[fmap[[nm]]] ## else pars[, nm] <- 0 (as f == FALSE) } else pars[, nm] <- f %*% beta[fmap[[nm]]] ## 2) Random effects for(i in iQ) if(!is.null(rm.i. <- rmapRel[[i]][[nm]])) { b.i <- b[[i]] b.i[] <- bvec[(bmap[i] + 1):bmap[i+1]] ## NB: some groups[[i]] may be *new* levels, i.e. non-matching: gr.i <- match(groups[[i]], colnames(b.i)) # column numbers + NA if (is.logical(r <- plist[[nm]]$random[[i]])) { if (r) pars[, nm] <- pars[, nm] + b.i[rm.i., gr.i] ## else r == FALSE =^= 0 } else if (data.class(r) != "list") { pars[, nm] <- pars[, nm] + (r * t(b.i)[gr.i, rm.i., drop=FALSE]) %*% rep(1, ncol(r)) } else { b.i.gi <- b.i[, gr.i, drop=FALSE] for(j in seq_along(rm.i.)) { if (is.logical(rr <- r[[j]])) { if(rr) pars[, nm] <- pars[, nm] + b.i.gi[rm.i.[[j]], ] ## else rr == FALSE =^= 0 } else pars[, nm] <- pars[, nm] + (rr * t(b.i.gi[rm.i.[[j]], , drop=FALSE])) %*% rep(1, ncol(rr)) } } } # for( i ) if(!is.null(rm.i. ..)) } pars } ### ### Methods for standard generics ### formula.nlme <- function(x, ...) eval(x$call[["model"]]) predict.nlme <- function(object, newdata, level = Q, asList = FALSE, na.action = na.fail, naPattern = NULL, ...) { ## ## method for predict() designed for objects inheriting from class nlme ## Q <- object$dims$Q if (missing(newdata)) { # will return fitted values val <- fitted(object, level, asList) if (length(level) == 1) return(val) return(data.frame(object[["groups"]][,level[level != 0], drop = FALSE], predict = val)) } maxQ <- max(level) # maximum level for predictions nlev <- length(level) newdata <- as.data.frame(newdata) if (maxQ > 0) { # predictions with random effects whichQ <- Q - (maxQ-1):0 reSt <- object$modelStruct$reStruct[whichQ] ## nlmeSt <- nlmeStruct(reStruct = reSt) groups <- getGroupsFormula(reSt) if (any(is.na(match(all.vars(groups), names(newdata))))) { ## groups cannot be evaluated in newdata stop("cannot evaluate groups for desired levels on 'newdata'") } } else { reSt <- NULL } mfArgs <- list(formula = asOneFormula( formula(object), object$call$fixed, formula(reSt), naPattern, omit = c(names(object$plist), "pi", deparse(getResponseFormula(object)[[2]]))), data = newdata, na.action = na.action, drop.unused.levels = TRUE) dataMix <- do.call(model.frame, mfArgs) origOrder <- row.names(dataMix) # preserve the original order whichRows <- match(origOrder, row.names(newdata)) if (maxQ > 0) { ## sort the model.frame by groups and get the matrices and parameters ## used in the estimation procedures grps <- getGroups(newdata, eval(substitute(~ 1 | GRPS, list(GRPS = groups[[2]])))) ## ordering data by groups if (inherits(grps, "factor")) { # single level grps <- grps[whichRows, drop = TRUE] oGrps <- data.frame(grps) ## checking if there are missing groups if (any(naGrps <- is.na(grps))) { grps[naGrps] <- levels(grps)[1L] # input with existing level } ord <- order(grps) #"order" treats a single named argument peculiarly grps <- data.frame(grps) row.names(grps) <- origOrder names(grps) <- names(oGrps) <- as.character(deparse((groups[[2L]]))) } else { grps <- oGrps <- do.call(data.frame, ## FIXME? better lapply(*, drop) ?? lapply(grps[whichRows, ], function(x) x[drop = TRUE])) ## checking for missing groups if (any(naGrps <- is.na(grps))) { ## need to input missing groups for(i in names(grps)) { grps[naGrps[, i], i] <- levels(grps[,i])[1L] } naGrps <- t(apply(naGrps, 1, cumsum)) # propagating NAs } ord <- do.call(order, grps) ## making group levels unique grps[, 1] <- grps[, 1][drop = TRUE] for(i in 2:ncol(grps)) { grps[, i] <- as.factor(paste(as.character(grps[, i-1]), as.character(grps[, i ]), sep = "/")) } } ## if (match(0, level, nomatch = 0)) { ## naGrps <- cbind(FALSE, naGrps) ## } ## naGrps <- as.matrix(naGrps)[ord, , drop = FALSE] naGrps <- cbind(FALSE, naGrps)[ord, , drop = FALSE] grps <- grps[ord, , drop = FALSE] dataMix <- dataMix[ord, ,drop = FALSE] revOrder <- match(origOrder, row.names(dataMix)) # putting in orig. order } ## making sure factor levels are the same as in contrasts contr <- object$contrasts for(i in names(dataMix)) { if (inherits(dataMix[,i], "factor") && !is.null(contr[[i]])) { levs <- levels(dataMix[,i]) levsC <- dimnames(contr[[i]])[[1L]] 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) } ## if (length(levs) < length(levsC)) { ## if (inherits(dataMix[,i], "ordered")) { ## dataMix[,i] <- ordered(as.character(dataMix[,i]), levels = levsC) ## } else { ## dataMix[,i] <- factor(as.character(dataMix[,i]), levels = levsC) ## } ## } attr(dataMix[,i], "contrasts") <- contr[[i]][levs, , drop = FALSE] } } N <- nrow(dataMix) ## ## evaluating the naPattern expression, if any ## naPat <- if(is.null(naPattern)) rep(TRUE, N) else as.logical(eval(asOneSidedFormula(naPattern)[[2]], dataMix)) ## ## Getting the plist for the new data frame ## plist <- object$plist fixed <- eval(object$call$fixed) if (!is.list(fixed)) fixed <- list(fixed) fixed <- do.call(c, lapply(fixed, function(fix.i) { if (is.name(fix.i[[2]])) list(fix.i) else ## multiple parameters on left hand side eval(parse(text = paste0("list(", paste(paste(all.vars(fix.i[[2]]), deparse (fix.i[[3]]), sep = "~"), collapse = ","), ")"))) })) fnames <- unlist(lapply(fixed, function(el) deparse(el[[2]]))) names(fixed) <- fnames fix <- fixef(object) ## fn <- names(fix) for(nm in fnames) { if (!is.logical(plist[[nm]]$fixed)) { oSform <- asOneSidedFormula(fixed[[nm]][[3]]) plist[[nm]]$fixed <- model.matrix(oSform, model.frame(oSform, dataMix)) } } if (maxQ > 0) { grpsRev <- lapply(rev(grps), as.character) ranForm <- formula(reSt)[whichQ] namGrp <- names(ranForm) rnames <- lapply(ranForm, function(el) unlist(lapply(el, function(el1) deparse(el1[[2]])))) for(i in seq_along(ranForm)) { names(ranForm[[i]]) <- rnames[[i]] } ran <- ranef(object) ran <- if(is.data.frame(ran)) list(ran) else rev(ran) ## rn <- lapply(ran[whichQ], names) ran <- lapply(ran, t) ranVec <- unlist(ran) for(nm in names(plist)) { for(i in namGrp) { if (!is.logical(plist[[nm]]$random[[i]])) { wch <- which(!is.na(match(rnames[[i]], nm))) plist[[nm]]$random[[i]] <- if (length(wch) == 1) { # only one formula for nm oSform <- asOneSidedFormula(ranForm[[i]][[nm]][[3]]) model.matrix(oSform, model.frame(oSform, dataMix)) } else { # multiple formulae lapply(ranForm[[i]][wch], function(el) { if (el[[3]] == "1") { TRUE } else { oSform <- asOneSidedFormula(el[[3]]) model.matrix(oSform, model.frame(oSform, dataMix)) } }) } } } } } else { namGrp <- "" grpsRev <- ranVec <- ran <- NULL } val <- vector("list", nlev) modForm <- getCovariateFormula(object)[[2]] omap <- object$map for(i in 1:nlev) { val[[i]] <- eval(modForm, data.frame(dataMix, getParsNlme(plist, omap$fmap, omap$rmapRel, omap$bmap, grpsRev, fix, ranVec, ran, level[i], N)))[naPat] } names(val) <- c("fixed", rev(namGrp))[level + 1] val <- as.data.frame(val) if (maxQ > 0) { val <- val[revOrder, , drop = FALSE] if (any(naGrps)) { val[naGrps] <- NA } } ## putting back in original order if (maxQ > 1) { # making groups unique for(i in 2:maxQ) oGrps[, i] <- as.factor(paste(as.character(oGrps[,i-1]), as.character(oGrps[,i ]), sep = "/")) } if (length(level) == 1) { val <- val[,1] ## ?? not in predict.lme() if (level > 0) { # otherwise 'oGrps' are typically undefined grps <- as.character(oGrps[, level]) if (asList) { val <- split(val, ordered(grps, levels = unique(grps))) } else { names(val) <- grps } } lab <- "Predicted values" if (!is.null(aux <- attr(object, "units")$y)) { lab <- paste(lab, aux) } attr(val, "label") <- lab val } else { data.frame(oGrps, predict = val) } } ## based on R's update.default update.nlme <- function (object, model., ..., evaluate = TRUE) { if (is.null(call <- object$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 <- c(as.list(call), extras[!existing]) call <- as.call(call) } } if(evaluate) eval(call, parent.frame()) else call } ## update.nlme <- ## function(object, model, data, fixed, random, groups, start, correlation, ## weights, subset, method, na.action, naPattern, control, ## verbose, ...) ## { ## thisCall <- as.list(match.call())[-(1:2)] ## if (!is.null(thisCall$start) && is.numeric(start)) { ## thisCall$start <- list(fixed = start) ## } ## if (!is.null(nextCall <- object$origCall) && ## (is.null(thisCall$fixed) && !is.null(thisCall$random))) { ## nextCall <- as.list(nextCall)[-1] ## } else { ## nextCall <- as.list(object$call)[-1] ## if (is.null(thisCall$fixed)) { # no changes in fixef model ## if (is.null(thisCall$start)) { ## thisCall$start <- list(fixed = fixef(object)) ## } else { ## if (is.null(thisCall$start$fixed)) { ## thisCall$start$fixed <- fixef(object) ## } ## } ## } ## if (!is.null(thisCall$start$random)) { # making start random NULL ## thisCall$start$random <- NULL ## } ## if (is.null(thisCall$random) && is.null(thisCall$subset)) { ## ## no changes in ranef model and no subsetting ## thisCall$random <- object$modelStruct$reStruct ## } ## } ## if (!is.null(thisCall$model)) { ## thisCall$model <- update(formula(object), model) ## } ## 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(nlme, nextCall) ## } ###*### nlmeStruct - a model structure for nlme fits nlmeStruct <- ## constructor for nlmeStruct objects function(reStruct, corStruct = NULL, varStruct = NULL)#, resp = NULL, ## model = NULL, local = NULL, N = NULL, naPat = NULL) { val <- list(reStruct = reStruct, corStruct = corStruct, varStruct = varStruct) structure(val[!vapply(val, is.null, NA)], # removing NULL components 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 = c("nlmeStruct", "lmeStruct", "modelStruct")) } ##*## nlmeStruct methods for standard generics fitted.nlmeStruct <- function(object, level = Q, conLin = attr(object, "conLin"), ...) { Q <- attr(object, "conLin")$dims[["Q"]] attr(object, "resp") - resid(object, level, conLin) } residuals.nlmeStruct <- function(object, level = Q, conLin = attr(object, "conLin"), ...) { Q <- conLin$dims[["Q"]] stopifnot(length(level) >= 1) loc <- attr(object, "local") oLev <- get("level", envir = loc) on.exit(assign("level", oLev, envir = loc)) dn <- c("fixed", rev(names(object$reStruct)))[level + 1] val <- array(0, c(attr(object, "NReal"), length(level)), list(dimnames(conLin$Xy)[[1]], dn)) for(i in seq_along(level)) { assign("level", level[i], envir = loc, immediate = TRUE) val[, i] <- c(eval(attr(object, "model")[[2]], envir = loc)) } val } nlmeControl <- ## Set control values for iterations within nlme function(maxIter = 50, pnlsMaxIter = 7, msMaxIter = 50, minScale = 0.001, tolerance = 1e-5, niterEM = 25, pnlsTol = 0.001, msTol = 0.000001, returnObject = FALSE, msVerbose = FALSE, gradHess = TRUE, apVar = TRUE, .relStep = .Machine$double.eps^(1/3), minAbsParApVar = 0.05, opt = c("nlminb", "nlm"), natural = TRUE, sigma = NULL, ...) { 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, pnlsMaxIter = pnlsMaxIter, msMaxIter = msMaxIter, minScale = minScale, tolerance = tolerance, niterEM = niterEM, pnlsTol = pnlsTol, msTol = msTol, returnObject = returnObject, msVerbose = msVerbose, gradHess = gradHess, apVar = apVar, .relStep = .relStep, minAbsParApVar = minAbsParApVar, opt = match.arg(opt), natural = natural, sigma=sigma, ...) } nonlinModel <- function(modelExpression, env, paramNames = get(".parameters", envir = env)) { modelExpression <- modelExpression[[2]] thisEnv <- environment() offset <- 0 ind <- vector("list", length(paramNames)) names(ind) <- paramNames for( i in paramNames ) { v.i <- get(i, envir = env) ind[[ i ]] <- offset + seq_along(v.i) offset <- offset + length(v.i) } modelValue <- eval(modelExpression, env) rm(i, offset, paramNames) function (newPars) { if(!missing(newPars)) { for(i in names(ind)) assign( i, unname(newPars[ ind[[i]] ]), envir = env) assign("modelValue", eval(modelExpression, env), envir = thisEnv) } modelValue } } ## Local Variables: ## ess-indent-offset: 2 ## End: