Re: [R] lme convergence

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Sat 01 Jul 2006 - 14:42:55 EST

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" <spencer.graves@pdf.com>
> To: "Ravi Varadhan" <rvaradhan@jhmi.edu>
> Cc: "'Pryseley Assam'" <assampryseley@yahoo.com>; "'R-Users'" 
> <R-help@stat.math.ethz.ch>; "Douglas Bates" <bates@stat.wisc.edu>
> 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:  rvaradhan@jhmi.edu
>>> Webpage: 
>>> http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>>> --------------------------------------------------------------------------
>>>
>>>> -----Original Message-----
>>>> From: r-help-bounces@stat.math.ethz.ch [mailto:r-help-
>>>> bounces@stat.math.ethz.ch] 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 Received on Mon Jul 03 17:55:43 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Mon 03 Jul 2006 - 18:15:04 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.