[R] R crash with modified lmer code

From: Colin Beale <c.beale_at_macaulay.ac.uk>
Date: Fri 26 Jan 2007 - 17:01:46 GMT


Hi all,

I've now got a problem with some modified lmer code (function lmer1 pasted at end) - I've made only three changes to the lmer code (marked), and I'm not really looking for comments on this function, but would like to know why execution of the following commands that use it almost invariably (but not quite predictably) leads to the R session terminating.

Here's the command that almost always (if not the first time its run, then certainly the second time) crashes R:

library(lme4)
source("M:/Rcode/lmer2.R") ##(Simply sourcing the modified lmer1 code below)
set.seed(1)
dat <- data.frame(Y = rnorm(400), X1 = rnorm(400), X2 = rnorm(400),

         F1 = as.factor(sample(1:4, 400, replace = T))) forSm <- Matrix(c(runif(1200)), ncol = 3, sparse = T) matDim <- dim(forSm)
smoother <- as.factor(sample(1:4, matDim[1], replace = T)) test <- lmer1 (Y ~ X1 + X2 + (1|F1) + (1|smoother), data = dat, rand_mat = forSm)

the sessionInfo (called after sourcing the code but before the crash) is:

R version 2.4.1 (2006-12-18)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United Kingdom.1252;LC_MONETARY=English_United
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

attached base packages:

[1] "stats"     "graphics"  "grDevices" "datasets"  "tcltk"     "utils"
    "methods"   "base"     

other attached packages:
       lme4      Matrix     lattice    svSocket        svIO      R2HTML
     svMisc       svIDE 
"0.9975-10"  "0.9975-8"   "0.14-16"     "0.9-5"     "0.9-5"      "1.58"
    "0.9-5"     "0.9-5" 



and the modified lmer code (sourced in the above code) is:

lmer1 <- function (formula, data, family = gaussian, method = c("REML",

    "ML", "PQL", "Laplace", "AGQ"), control = list(), start = NULL,     subset, weights, na.action, offset, contrasts = NULL, model = TRUE,  

    matDimI = matDim, rand_matI = rand_mat, ...) {

    method <- match.arg(method)
    formula <- as.formula(formula)
    if (length(formula) < 3)

        stop("formula must be a two-sided formula")

    cv <- do.call("lmerControl", control)
    mc <- match.call()
    fr <- lmerFrames(mc, formula, data, contrasts)
    if (matDimI[1] != length(fr$Y))
        stop("forSmoothing is not a sensible length")                 
#insert

    Y <- fr$Y
    X <- fr$X
    weights <- fr$weights
    offset <- fr$offset
    mf <- fr$mf
    mt <- fr$mt
    if (is.character(family))

        family <- get(family, mode = "function", envir = parent.frame())

    if (is.function(family))

        family <- family()
    if (is.null(family$family)) {

        print(family)
        stop("'family' not recognized")

}

    fltype <- mkFltype(family)
    FL <- lmerFactorList(formula, mf, fltype)
    nFacs <- length(FL)                                          #My
insert

    cnames <- with(FL, c(lapply(Ztl, rownames), list(.fixed = colnames(X))))

    nc <- with(FL, sapply(Ztl, nrow))
    Ztl <- with(FL, .Call(Ztl_sparse, fl, Ztl))

    Ztl[[nFacs]] <- t(rand_matI)                                  #My
insert

    Zt <- if (length(Ztl) == 1)

        Ztl[[1]]
    else do.call("rbind", Ztl)
    fl <- FL$fl
    if (fltype < 0) {

        mer <- .Call(mer_create, fl, Zt, X, as.double(Y), method == 
            "REML", nc, cnames)
        if (!is.null(start)) 
            mer <- setOmega(mer, start)
        .Call(mer_ECMEsteps, mer, cv$niterEM, cv$EMverbose)
        LMEoptimize(mer) <- cv
        return(new("lmer", mer, frame = if (model) fr$mf else
data.frame(), 
            terms = mt, call = mc))

}

    if (method %in% c("ML", "REML"))

        method <- "Laplace"
    if (method == "AGQ")

        stop("method = \"AGQ\" not yet implemented for supernodal representation")

    if (method == "PQL")

        cv$usePQL <- TRUE
    glmFit <- glm.fit(X, Y, weights = weights, offset = offset,

        family = family, intercept = attr(mt, "intercept") > 
            0)
    Y <- as.double(glmFit$y)

    mer <- .Call(mer_create, fl, Zt, X, Y, 0, nc, cnames)     if (!is.null(start))

        mer <- setOmega(mer, start)
    gVerb <- getOption("verbose")
    weights <- glmFit$prior.weights
    eta <- glmFit$linear.predictors
    linkinv <- quote(family$linkinv(eta))     mu.eta <- quote(family$mu.eta(eta))
    mu <- family$linkinv(eta)
    variance <- quote(family$variance(mu))     dev.resids <- quote(family$dev.resids(Y, mu, weights))     LMEopt <- get("LMEoptimize<-")
    doLMEopt <- quote(LMEopt(x = mer, value = cv))     if (family$family %in% c("binomial", "poisson"))

        mer@devComp[8] <- -1
    mer@status["glmm"] <- as.integer(switch(method, PQL = 1,

        Laplace = 2, AGQ = 3))
    GSpt <- .Call(glmer_init, environment(), fltype)     if (cv$usePQL) {

        .Call(glmer_PQL, GSpt)
        PQLpars <- c(fixef(mer), .Call(mer_coef, mer, 2))

}

    else {

        PQLpars <- c(coef(glmFit), .Call(mer_coef, mer, 2))
}

    if (method == "PQL") {

        .Call(glmer_devLaplace, PQLpars, GSpt)
        .Call(glmer_finalize, GSpt)
        return(new("glmer", new("lmer", mer, frame = if (model) mf else
data.frame(), 
            terms = mt, call = match.call()), weights = weights, 
            family = family))

}

    fixInd <- seq(ncol(X))
    const <- c(rep(FALSE, length(fixInd)), unlist(lapply(mer@nc[seq(along = fl)],

        function(k) 1:((k * (k + 1))/2) <= k)))     devLaplace <- function(pars) .Call(glmer_devLaplace, pars,

        GSpt)
    rel.tol <- abs(0.01/devLaplace(PQLpars))     if (cv$msVerbose)

        cat(paste("relative tolerance set to", rel.tol, "\n"))     optimRes <- nlminb(PQLpars, devLaplace, lower = ifelse(const,

        5e-10, -Inf), control = list(trace = cv$msVerbose, iter.max = cv$msMaxIter,

        rel.tol = rel.tol))
    .Call(glmer_finalize, GSpt)
    new("glmer", new("lmer", mer, frame = if (model)

        mf
    else data.frame(), terms = mt, call = match.call()), weights = weights,

        family = family)
}
environment(lmer1) <- environment(lmer)

Thanks,

Colin

Dr. Colin Beale
Spatial Ecologist
The Macaulay Institute
Craigiebuckler
Aberdeen
AB15 8QH
UK

Tel: 01224 498245 ext. 2427
Fax: 01224 311556
Email: c.beale@macaulay.ac.uk

-- 
Please note that the views expressed in this e-mail are those of the
sender and do not necessarily represent the views of the Macaulay
Institute. This email and any attachments are confidential and are
intended solely for the use of the recipient(s) to whom they are
addressed. If you are not the intended recipient, you should not read,
copy, disclose or rely on any information contained in this e-mail, and
we would ask you to contact the sender immediately and delete the email
from your system. Thank you.
Macaulay Institute and Associated Companies, Macaulay Drive,
Craigiebuckler, Aberdeen, AB15 8QH.

______________________________________________
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
and provide commented, minimal, self-contained, reproducible code.
Received on Sat Jan 27 04:36:34 2007

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 Fri 26 Jan 2007 - 18:30:30 GMT.

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