[R-sig-ME] modified lmer2 code problem
Douglas Bates
bates at stat.wisc.edu
Tue Feb 13 19:54:14 CET 2007
On 2/13/07, Colin Beale <c.beale at macaulay.ac.uk> wrote:
> I'm trying to replicate some code I have that works in GenStat
> (http://www.vsni.co.uk/products/genstat/) using a modification of the
> lmer2 function. This modification allows a matrix to be passed directly
> to the random effects, rather than being built up from the factor
> levels. I've written a little wrapper that passes things to the modified
> version of lmer2 and called it lmerWrap (code for both the modified
> lmer2 and lmerWrap functions are at the end). I've tested it by passing
> model matrices to the function and comparing this to the lmer2 output
> and all works well. For example:
> fm1 <- lmer(Reaction ~ 1 + (1|Days) + (1|Subject), data = sleepstudy)
> fm2 <- lmerWrap(Reaction ~ 1 + (1|Days) , data = sleepstudy, mat =
> model.matrix ( ~ -1 + sleepstudy$Subject ))
> give apparently identical parameter estimates, as they should do. I've
> then added a small random number to model.matrix(~ -1 +
> sleepstudy$Subject ):
lmer or lmer2? You mentioned lmer2 above but seem to be calling lmer here
> mat2 <- model.matrix ( ~ -1 + sleepstudy$Subject ) + runif (length
> (prod (dim (model.matrix ( ~ -1 + sleepstudy$Subject)))))
>
> and fitted a third model:
>
> fm3 <- lmerWrap (Reaction ~ 1 + (1|Days), data = sleepstudy,
> mat = mat2)
>
> I export this "smudged" dataset to GenStat and run the code I have
> there and I get results identical to 4 significant figures, so all seems
> to be working well.
>
> However, I then try to run the code using the real dataset for which I
> want to replicate the analysis, and I get completely different results
> to that when run in GenStat, as the variance is estimated to be 0 in R:
> here's my code and output (if anyone wants the data too let me know).
>
> try2 <- lmerWrap(BirthWeight ~ m_sum + Motherrs + MotherAge
> + MotherAge2 + Sex + birthday + (1|Mother)
> + (1|BirthYear), data = deer, mat = rand_mat, control =
> list(msVerbose = T,
> niterEM = 0, gradient = FALSE))
> 0 3564.56: 0.916339 0.235294 0.436794
> 1 3563.77: 0.823463 0.263738 0.435422
> 2 3563.44: 0.838527 0.253850 0.434686
> 3 3563.39: 0.844450 0.232247 0.406417
> 4 3563.31: 0.852666 0.247116 0.405113
> 5 3563.31: 0.853871 0.244576 0.403450
> 6 3563.31: 0.852928 0.245936 0.400634
> 7 3563.30: 0.854648 0.244316 0.394544
> 8 3563.24: 0.871616 0.260905 0.198572
> 9 3562.99: 0.859845 0.253159 0.00167238
> 10 3562.96: 0.852482 0.245837 0.00164461
> 11 3562.96: 0.853410 0.247334 0.00000
> 12 3562.96: 0.853197 0.246962 0.00000
> 13 3562.96: 0.853194 0.246967 0.00000
> As I know the GenStat code recovers the correct parameters for
> simulated data, and I get similar results in the sleepstudy example
> above, I am at a loss as to what is happening here. Does anyone have any
> suggestions as to where I go from here? Is lmer struggling to find
> sensible starting values or does anyone have any other suggestions?
I'm not sure I understand the model. As I see it there are two
variance components in the model, (1|Mother) and (1|BirthYear) so I
don't understand the 3 parameters in the optimization.
> Thanks very much,
>
> Colin
>
>
> lmerWrap <- function (formula, data = NULL, mat, ...)
> {
> colsums <- colSums(mat)
> matAdd <- as.factor(rep(1:(dim(mat)[2]), length = dim(mat)[1]))
> formula <- as.formula(formula)
> formula <- update(formula, .~. + (1|matAdd))
> if (!is.null(data)) dat <- as.data.frame(cbind(data, matAdd))
> smoothMat <- Matrix(mat, sparse = T)
> run <- lmer2a(formula = formula, data = dat, matDimI = dim(mat),
> rand_matI = smoothMat, ...)
> }
>
>
> lmer2a <- 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,
> rand_matI, ...)
> {
> print(t(rand_matI))
> 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)
> 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)
> 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))
> nFacs <- list(length(Ztl)) #My insert
> for (i in 1:length(Ztl)) { #My insert
> nFacs[[i]] <- dim(Ztl[[i]]) #My insert
> } #My insert
> id <- unlist(lapply (nFacs, all.equal, dim(t(rand_matI)))) #My
> insert
> id <- which(id == "TRUE") #My
> insert
> Ztl[[id]] <- 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(mer2_create, fl, Zt, t(X), as.double(Y),
> method == "REML", nc, cnames, fr$offset, fr$weights)
> if (!is.null(start))
> mer <- setST(mer, start)
> const <- unlist(lapply(mer at nc, function(n) rep(1:0, c(n,
> (n * (n - 1))/2))))
> optimRes <- nlminb(.Call(mer2_getPars, mer), function(x)
> .Call(mer2_deviance,
> .Call(mer2_setPars, mer, x), as.integer(0)), lower =
> ifelse(const,
> 0, -Inf), control = list(trace = cv$msVerbose, iter.max =
> cv$msMaxIter))
> if (optimRes$convergence)
> warning(paste("nlminb failed to converge:",
> optimRes$message))
> .Call(mer2_setPars, mer, optimRes$par)
> .Call(mer2_update_effects, mer)
> return(new("lmer2", mer, frame = if (model) fr$mf else
> data.frame(),
> terms = mt, call = mc))
> }
> mer
> }
> environment(lmer2a) <- environment(lmer2)
>
>
> --
> 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.
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list