Hi, Dimitris: The change you suggested sounds constructive. Unfortunately, it did NOT solve the problem, at least for the modification of the example from the 'lme' help page I tested.
However, one other similar change (and adding 'nlme:::' to the calls for functions hidden in an nlme namespace, identified partly with 'traceback()') produced a version of 'lme.formula' that actually returned an answer for my test case: > fm1 <- lme(distance ~ age, data = Orthodont, + control=lmeControl(msMaxIter=1, + returnObject=TRUE)) Warning message: nlminb problem, convergence error code = 1 ; message = iteration limit reached without convergence (9) in: lme.formula(distance ~ age, data = Orthodont, control = lmeControl(msMaxIter = 1, Below please find (1) the 2 changes and (2) a complete version of 'lme.formula' that worked for this example. Thanks for your contribution. Spencer Graves ### change 1 ################## if (!needUpdate(lmeSt)) { if(optRes$convergence){ msg <- paste(controlvals$opt, "problem, convergence error code =", optRes$convergence, "; message =", optRes$message) { if(!controlvals$returnObject) stop(msg) else warning(msg) } break } } ### change 2 ################# if (numIter > controlvals$maxIter) { msg <- paste("Maximum number of iterations", "(lmeControl(maxIter)) reached without convergence.") if (controlvals$returnObject) { warning(msg) break } else { stop(msg) } } ### lme.formula revised ########### lme.formula <- function(fixed, data = sys.frame(sys.parent()), random = pdSymm( eval( as.call( fixed[ -2 ] ) ) ), correlation = NULL, weights = NULL, subset, method = c("REML", "ML"), na.action = na.fail, control = list(), contrasts = NULL, keep.data = TRUE) { Call <- match.call() miss.data <- missing(data) || !is.data.frame(data) ## control parameters controlvals <- lmeControl() if (!missing(control)) { if(!is.null(control$nlmStepMax) && control$nlmStepMax < 0) { warning("Negative control$nlmStepMax - using default value") control$nlmStepMax <- NULL } controlvals[names(control)] <- control } ## ## checking arguments ## if (!inherits(fixed, "formula") || length(fixed) != 3) { stop("\nFixed-effects model must be a formula of the form \"resp ~ pred\"") } method <- match.arg(method) REML <- method == "REML" reSt <- reStruct(random, REML = REML, data = NULL) 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 randL <- as.list(randL) reSt <- reStruct(randL, REML = REML, data = NULL) } else { names(reSt) <- namGrp } } else { ## will assume single group groups <- ~ 1 names(reSt) <- "1" } } ## 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(paste("Incompatible formulas for groups in \"random\"", "and \"correlation\"")) } if (corQ < lmeQ) { warning(paste("Cannot use smaller level of grouping for", "\"correlation\" than for \"random\". Replacing", "the former with the latter.")) attr(correlation, "formula") <- eval(parse(text = paste("~", deparse(getCovariateFormula(formula(correlation))[[2]]), "|", deparse(groups[[2]])))) } } else { if (any(lmeGrpsForm != corGrpsForm[1:lmeQ])) { stop(paste("Incompatible formulas for groups in \"random\"", "and \"correlation\"")) } } } else { ## using the same grouping as in random attr(correlation, "formula") <- eval(parse(text = paste("~", deparse(getCovariateFormula(formula(correlation))[[2]]), "|", deparse(groups[[2]])))) corQ <- lmeQ <- 1 } } else { corQ <- lmeQ <- 1 } ## create an lme structure containing the random effects model and plug-ins lmeSt <- lmeStruct(reStruct = reSt, corStruct = correlation, varStruct = varFunc(weights)) ## extract a data frame with enough information to evaluate ## fixed, groups, reStruct, corStruct, and varStruct mfArgs <- list(formula = asOneFormula(formula(lmeSt), fixed, groups), 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 for(i in names(contrasts)) # handle contrasts statement contrasts(dataMix[[i]]) = contrasts[[i]] ## sort the model.frame by groups and get the matrices and parameters ## used in the estimation procedures grps <- getGroups(dataMix, groups) ## 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 = "/")) NULL } } 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 ## obtaining basic model matrices N <- nrow(grps) Z <- model.matrix(reSt, dataMix) ncols <- attr(Z, "ncols") Names(lmeSt$reStruct) <- attr(Z, "nams") ## keeping the contrasts for later use in predict contr <- attr(Z, "contr") X <- model.frame(fixed, dataMix) Terms <- attr(X, "terms") auxContr <- lapply(X, function(el) if (inherits(el, "factor") && length(levels(el)) > 1) contrasts(el)) contr <- c(contr, auxContr[is.na(match(names(auxContr), names(contr)))]) contr <- contr[!unlist(lapply(contr, is.null))] X <- model.matrix(fixed, data=X) y <- eval(fixed[[2]], dataMix) ncols <- c(ncols, dim(X)[2], 1) Q <- ncol(grps) ## creating the condensed linear model attr(lmeSt, "conLin") <- list(Xy = array(c(Z, X, y), c(N, sum(ncols)), list(row.names(dataMix), c(colnames(Z), colnames(X), deparse(fixed[[2]])))), dims = nlme:::MEdims(grps, ncols), logLik = 0) ## checking if enough observations per group to estimate ranef tmpDims <- attr(lmeSt, "conLin")$dims if (max(tmpDims$ZXlen[[1]]) < tmpDims$qvec[1]) { warning(paste("Fewer observations than random effects in all level", Q,"groups")) } ## degrees of freedom for testing fixed effects fixDF <- nlme:::getFixDF(X, grps, attr(lmeSt, "conLin")$dims$ngrps, terms = Terms) ## initialization lmeSt <- Initialize(lmeSt, dataMix, grps, control = controlvals) parMap <- attr(lmeSt, "pmap") ## Checking possibility of single decomposition if (length(lmeSt) == 1) { # reStruct only, can do one decomposition ## need to save conLin for calculating fitted values and residuals oldConLin <- attr(lmeSt, "conLin") decomp <- TRUE attr(lmeSt, "conLin") <- nlme:::MEdecomp(attr(lmeSt, "conLin")) } else decomp <- FALSE ## ## getting the linear mixed effects fit object, ## possibly iterating for variance functions ## numIter <- 0 repeat { oldPars <- coef(lmeSt) optRes <- if (controlvals$opt == "nlminb") { nlminb(c(coef(lmeSt)), function(lmePars) -logLik(lmeSt, lmePars), control = list(iter.max = controlvals$msMaxIter, trace = controlvals$msVerbose)) } else { optim(c(coef(lmeSt)), function(lmePars) -logLik(lmeSt, lmePars), control = list(trace = controlvals$msVerbose, maxit = controlvals$msMaxIter, reltol = if(numIter == 0) controlvals$msTol else 100*.Machine$double.eps), method = controlvals$optimMethod) } numIter0 <- NULL coef(lmeSt) <- optRes$par attr(lmeSt, "lmeFit") <- nlme:::MEestimate(lmeSt, grps) ## checking if any updating is needed if (!needUpdate(lmeSt)) { if(optRes$convergence){ msg <- paste(controlvals$opt, "problem, convergence error code =", optRes$convergence, "; message =", optRes$message) { if(!controlvals$returnObject) stop(msg) else warning(msg) } break } } ## updating the fit information numIter <- numIter + 1 lmeSt <- update(lmeSt, dataMix) ## calculating the convergence criterion aConv <- coef(lmeSt) conv <- abs((oldPars - aConv)/ifelse(aConv == 0, 1, aConv)) aConv <- NULL for(i in names(lmeSt)) { if (any(parMap[,i])) { aConv <- c(aConv, max(conv[parMap[,i]])) names(aConv)[length(aConv)] <- i } } if (max(aConv) <= controlvals$tolerance) { break } if (numIter > controlvals$maxIter) { msg <- paste("Maximum number of iterations", "(lmeControl(maxIter)) reached without convergence.") if (controlvals$returnObject) { warning(msg) break } else { stop(msg) } } } ## wrapping up lmeFit <- attr(lmeSt, "lmeFit") names(lmeFit$beta) <- namBeta <- colnames(X) attr(fixDF, "varFixFact") <- varFix <- lmeFit$sigma * lmeFit$varFix varFix <- crossprod(varFix) dimnames(varFix) <- list(namBeta, namBeta) ## ## fitted.values and residuals (in original order) ## Fitted <- fitted(lmeSt, level = 0:Q, conLin = if (decomp) { oldConLin } else { attr(lmeSt, "conLin") })[revOrder, , drop = FALSE] Resid <- y[revOrder] - Fitted attr(Resid, "std") <- lmeFit$sigma/(varWeights(lmeSt)[revOrder]) ## putting groups back in original order grps <- grps[revOrder, , drop = FALSE] ## making random effects estimates consistently ordered # for(i in names(lmeSt$reStruct)) { # lmeFit$b[[i]] <- lmeFit$b[[i]][unique(as.character(grps[, i])),, drop = F] # NULL # } ## inverting back reStruct lmeSt$reStruct <- solve(lmeSt$reStruct) ## saving part of dims dims <- attr(lmeSt, "conLin")$dims[c("N", "Q", "qvec", "ngrps", "ncol")] ## getting the approximate var-cov of the parameters if (controlvals$apVar) { apVar <- nlme:::lmeApVar(lmeSt, lmeFit$sigma, .relStep = controlvals[[".relStep"]], minAbsPar = controlvals[["minAbsParApVar"]], natural = controlvals[["natural"]]) } else { apVar <- "Approximate variance-covariance matrix not available" } ## getting rid of condensed linear model and fit attr(lmeSt, "conLin") <- NULL attr(lmeSt, "lmeFit") <- NULL ## ## creating the lme object ## estOut <- list(modelStruct = lmeSt, dims = dims, contrasts = contr, coefficients = list( fixed = lmeFit$beta, random = lmeFit$b), varFix = varFix, sigma = lmeFit$sigma, apVar = apVar, logLik = lmeFit$logLik, numIter = if (needUpdate(lmeSt)) numIter else numIter0, groups = grps, call = Call, terms = Terms, method = method, fitted = Fitted, residuals = Resid, fixDF = fixDF) if (keep.data && !miss.data) estOut$data <- data if (inherits(data, "groupedData")) { ## saving labels and units for plots attr(estOut, "units") <- attr(data, "units") attr(estOut, "labels") <- attr(data, "labels") } class(estOut) <- "lme" estOut } ################################### Dimitris Rizopoulos wrote: > I think the following part of lme.formula > > if (numIter > controlvals$maxIter) { > stop("Maximum number of iterations reached without convergence.") > } > > > should be something like > > > if (numIter > controlvals$maxIter) { > if (controlvals$returnObject) { > warning("Maximum number of iterations reached without > convergence.") > break > } else { > stop("Maximum number of iterations reached without > convergence.") > } > } > > > Best, > Dimitris > > ---- > Dimitris Rizopoulos > Ph.D. Student > Biostatistical Centre > School of Public Health > Catholic University of Leuven > > Address: Kapucijnenvoer 35, Leuven, Belgium > Tel: +32/(0)16/336899 > Fax: +32/(0)16/337015 > Web: http://med.kuleuven.be/biostat/ > http://www.student.kuleuven.be/~m0390867/dimitris.htm > > > ----- Original Message ----- > From: "Spencer Graves" <[EMAIL PROTECTED]> > To: "Ravi Varadhan" <[EMAIL PROTECTED]> > Cc: "'Pryseley Assam'" <[EMAIL PROTECTED]>; "'R-Users'" > <R-help@stat.math.ethz.ch>; "Douglas Bates" <[EMAIL PROTECTED]> > Sent: Friday, June 30, 2006 4:08 AM > Subject: Re: [R] lme convergence > > >> Does anyone know how to obtain the 'returnObject' from an 'lme' >> run >> that fails to converge? An argument of this name is described on >> the >> 'lmeControl' help page as, "a logical value indicating whether the >> fitted object should be returned when the maximum number of >> iterations >> is reached without convergence of the algorithm. Default is >> 'FALSE'." >> >> Unfortunately, I've so far been unable to get it to work, as >> witnessed by the following modification of an example from the >> '?lme' >> help page: >> >>> library(nlme) >>> fm1 <- lme(distance ~ age, data = Orthodont, >> + control=lmeControl(msMaxIter=1)) >> Error in lme.formula(distance ~ age, data = Orthodont, control = >> lmeControl(msMaxIter = 1)) : >> iteration limit reached without convergence (9) >>> fm1 >> Error: object "fm1" not found >>> fm1 <- lme(distance ~ age, data = Orthodont, >> + control=lmeControl(msMaxIter=1, >> + returnObject=TRUE)) >> Error in lme.formula(distance ~ age, data = Orthodont, control = >> lmeControl(msMaxIter = 1, : >> iteration limit reached without convergence (9) >>> fm1 >> Error: object "fm1" not found >> >> I might be able to fix the problem myself, working through the >> 'lme' >> code line by line, e.g., using 'debug'. However, I'm not ready to >> do >> that just now. >> >> Best Wishes, >> Spencer Graves >> >> Ravi Varadhan wrote: >>> Use "try" to capture error messages without breaking the loop. >>> ?try >>> >>> -------------------------------------------------------------------------- >>> Ravi Varadhan, Ph.D. >>> Assistant Professor, The Center on Aging and Health >>> Division of Geriatric Medicine and Gerontology >>> Johns Hopkins University >>> Ph: (410) 502-2619 >>> Fax: (410) 614-9625 >>> Email: [EMAIL PROTECTED] >>> Webpage: >>> http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html >>> -------------------------------------------------------------------------- >>> >>>> -----Original Message----- >>>> From: [EMAIL PROTECTED] [mailto:r-help- >>>> [EMAIL PROTECTED] On Behalf Of Pryseley Assam >>>> Sent: Wednesday, June 28, 2006 12:18 PM >>>> To: R-Users >>>> Subject: [R] lme convergence >>>> >>>> Dear R-Users, >>>> >>>> Is it possible to get the covariance matrix from an lme model >>>> that did >>>> not converge ? >>>> >>>> I am doing a simulation which entails fitting linear mixed >>>> models, using >>>> a "for loop". >>>> Within each loop, i generate a new data set and analyze it using >>>> a mixed >>>> model. The loop stops When the "lme function" does not converge >>>> for a >>>> simulated dataset. I want to inquire if there is a method to >>>> suppress the >>>> error message from the lme function, or better still, a way of >>>> going about >>>> this issue of the loop ending once the lme function does not >>>> converge. >>>> >>>> Thanks in advance, >>>> Pryseley >>>> >>>> >>>> --------------------------------- >>>> >>>> >>>> [[alternative HTML version deleted]] >>>> >>>> ______________________________________________ >>>> R-help@stat.math.ethz.ch mailing list >>>> https://stat.ethz.ch/mailman/listinfo/r-help >>>> PLEASE do read the posting guide! >>>> http://www.R-project.org/posting- >>>> guide.html >>> ______________________________________________ >>> R-help@stat.math.ethz.ch mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-help >>> PLEASE do read the posting guide! >>> http://www.R-project.org/posting-guide.html >> ______________________________________________ >> R-help@stat.math.ethz.ch mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide! >> http://www.R-project.org/posting-guide.html >> > > > Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm > > ______________________________________________ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html