[R] lme convergence
Spencer Graves
spencer.graves at pdf.com
Sat Jul 1 06:42:55 CEST 2006
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 at pdf.com>
> To: "Ravi Varadhan" <rvaradhan at jhmi.edu>
> Cc: "'Pryseley Assam'" <assampryseley at yahoo.com>; "'R-Users'"
> <R-help at stat.math.ethz.ch>; "Douglas Bates" <bates at 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 at jhmi.edu
>>> Webpage:
>>> http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>>> --------------------------------------------------------------------------
>>>
>>>> -----Original Message-----
>>>> From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-
>>>> bounces at 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 at 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 at 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 at 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 at 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
>
More information about the R-help
mailing list