SCM

SCM Repository

[matrix] Diff of /pkg/R/lmer.R
ViewVC logotype

Diff of /pkg/R/lmer.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1279, Tue May 23 08:01:11 2006 UTC revision 1280, Tue May 23 08:01:44 2006 UTC
# Line 110  Line 110 
110  }  }
111    
112  rWishart <- function(n, df, invScal)  rWishart <- function(n, df, invScal)
113      .Call("Matrix_rWishart", n, df, invScal, PACKAGE = "Matrix")      .Call(Matrix_rWishart, n, df, invScal)
114    
115  setMethod("coef", signature(object = "mer"),  setMethod("coef", signature(object = "mer"),
116            function(object, ...)            function(object, ...)
# Line 298  Line 298 
298                                         mf)))                                         mf)))
299            if (lmm) {            if (lmm) {
300                ## Create the mixed-effects representation (mer) object                ## Create the mixed-effects representation (mer) object
301                mer <- .Call("mer_create", fl,                mer <- .Call(mer_create, fl,
302                             .Call("Zt_create", fl, Ztl, PACKAGE = "Matrix"),                             .Call(Zt_create, fl, Ztl),
303                             X, Y, method, sapply(Ztl, nrow),                             X, Y, method, sapply(Ztl, nrow),
304                             c(lapply(Ztl, rownames), list(.fixed = colnames(X))),                             c(lapply(Ztl, rownames), list(.fixed = colnames(X))),
305                             !(family$family %in% c("binomial", "poisson")),                             !(family$family %in% c("binomial", "poisson")),
306                             match.call(), family,                             match.call(), family,
307                             PACKAGE = "Matrix")                             PACKAGE = "Matrix")
308                .Call("mer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")                .Call(mer_ECMEsteps, mer, cv$niterEM, cv$EMverbose)
309                LMEoptimize(mer) <- cv                LMEoptimize(mer) <- cv
310                return(new("lmer", mer,                return(new("lmer", mer,
311                           frame = if (model) mf else data.frame(),                           frame = if (model) mf else data.frame(),
# Line 327  Line 327 
327            dev.resids <- quote(family$dev.resids(Y, mu, wtssqr))            dev.resids <- quote(family$dev.resids(Y, mu, wtssqr))
328            LMEopt <- get("LMEoptimize<-")            LMEopt <- get("LMEoptimize<-")
329            doLMEopt <- quote(LMEopt(x = mer, value = cv))            doLMEopt <- quote(LMEopt(x = mer, value = cv))
330            mer <- .Call("mer_create", fl,            mer <- .Call(mer_create, fl,
331                         .Call("Zt_create", fl, Ztl, PACKAGE = "Matrix"),                         .Call(Zt_create, fl, Ztl),
332                         X, Y, method, sapply(Ztl, nrow),                         X, Y, method, sapply(Ztl, nrow),
333                         c(lapply(Ztl, rownames), list(.fixed = colnames(X))),                         c(lapply(Ztl, rownames), list(.fixed = colnames(X))),
334                         !(family$family %in% c("binomial", "poisson")),                         !(family$family %in% c("binomial", "poisson")),
335                         match.call(), family,                         match.call(), family,
336                         PACKAGE = "Matrix")                         PACKAGE = "Matrix")
337    
338            GSpt <- .Call("glmer_init", environment(), PACKAGE = "Matrix")            GSpt <- .Call(glmer_init, environment())
339            if (cv$usePQL) {            if (cv$usePQL) {
340                .Call("glmer_PQL", GSpt, PACKAGE = "Matrix")  # obtain PQL estimates                .Call(glmer_PQL, GSpt)  # obtain PQL estimates
341                PQLpars <- c(fixef(mer),                PQLpars <- c(fixef(mer),
342                             .Call("mer_coef", mer, 2, PACKAGE = "Matrix"))                             .Call(mer_coef, mer, 2))
343            } else {            } else {
344                PQLpars <- c(coef(glmFit),                PQLpars <- c(coef(glmFit),
345                             .Call("mer_coef", mer, 2, PACKAGE = "Matrix"))                             .Call(mer_coef, mer, 2))
346            }            }
347            if (method == "PQL") {            if (method == "PQL") {
348                .Call("glmer_devLaplace", PQLpars, GSpt, PACKAGE = "Matrix")                .Call(glmer_devLaplace, PQLpars, GSpt)
349                .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")                .Call(glmer_finalize, GSpt)
350                return(new("lmer", mer,                return(new("lmer", mer,
351                           frame = if (model) mf else data.frame(),                           frame = if (model) mf else data.frame(),
352                           terms = mt))                           terms = mt))
# Line 360  Line 360 
360                                     function(k) 1:((k*(k+1))/2) <= k)                                     function(k) 1:((k*(k+1))/2) <= k)
361                              ))                              ))
362            devLaplace <- function(pars)            devLaplace <- function(pars)
363                .Call("glmer_devLaplace", pars, GSpt, PACKAGE = "Matrix")                .Call(glmer_devLaplace, pars, GSpt)
364    
365            optimRes <-            optimRes <-
366                nlminb(PQLpars, devLaplace,                nlminb(PQLpars, devLaplace,
367                       lower = ifelse(const, 5e-10, -Inf),                       lower = ifelse(const, 5e-10, -Inf),
368                       control = list(trace = cv$msVerbose,                       control = list(trace = cv$msVerbose,
369                       iter.max = cv$msMaxIter))                       iter.max = cv$msMaxIter))
370            .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")            .Call(glmer_finalize, GSpt)
371            return(new("lmer", mer,            return(new("lmer", mer,
372                       frame = if (model) mf else data.frame(),                       frame = if (model) mf else data.frame(),
373                       terms = mt))                       terms = mt))
# Line 376  Line 376 
376    
377  ## Extract the L matrix  ## Extract the L matrix
378  setAs("mer", "dtCMatrix", function(from)  setAs("mer", "dtCMatrix", function(from)
379        .Call("mer_dtCMatrix", from, PACKAGE = "Matrix"))        .Call(mer_dtCMatrix, from))
380    
381  ## Extract the fixed effects  ## Extract the fixed effects
382  setMethod("fixef", signature(object = "mer"),  setMethod("fixef", signature(object = "mer"),
383            function(object, ...)            function(object, ...)
384            .Call("mer_fixef", object, PACKAGE = "Matrix"))            .Call(mer_fixef, object))
385    
386  ## Extract the random effects  ## Extract the random effects
387  setMethod("ranef", signature(object = "mer"),  setMethod("ranef", signature(object = "mer"),
388            function(object, postVar = FALSE, ...) {            function(object, postVar = FALSE, ...) {
389                ans <- new("ranef.lmer",                ans <- new("ranef.lmer",
390                           lapply(.Call("mer_ranef", object, PACKAGE = "Matrix"),                           lapply(.Call(mer_ranef, object),
391                                  data.frame, check.names = FALSE))                                  data.frame, check.names = FALSE))
392                names(ans) <- names(object@flist)                names(ans) <- names(object@flist)
393                if (postVar) {                if (postVar) {
394                    pV <- .Call("mer_postVar", object, PACKAGE = "Matrix")                    pV <- .Call(mer_postVar, object)
395                    for (i in seq(along = ans))                    for (i in seq(along = ans))
396                        attr(ans[[i]], "postVar") <- pV[[i]]                        attr(ans[[i]], "postVar") <- pV[[i]]
397                }                }
# Line 405  Line 405 
405                   nc <- x@nc                   nc <- x@nc
406                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
407                   fn <- function(pars)                   fn <- function(pars)
408                       deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix"))                       deviance(.Call(mer_coefGets, x, pars, 2))
409                   gr <- if (value$gradient)                   gr <- if (value$gradient)
410                       function(pars) {                       function(pars) {
411                           if (!isTRUE(all.equal(pars,                           if (!isTRUE(all.equal(pars,
412                                                 .Call("mer_coef", x,                                                 .Call(mer_coef, x,
413                                                       2, PACKAGE = "Matrix"))))                                                       2))))
414                               .Call("mer_coefGets", x, pars, 2, PACKAGE = "Matrix")                               .Call(mer_coefGets, x, pars, 2)
415                           .Call("mer_gradient", x, 2, PACKAGE = "Matrix")                           .Call(mer_gradient, x, 2)
416                       }                       }
417                   else NULL                   else NULL
418                   optimRes <- nlminb(.Call("mer_coef", x, 2, PACKAGE = "Matrix"),                   optimRes <- nlminb(.Call(mer_coef, x, 2),
419                                      fn, gr,                                      fn, gr,
420                                      lower = ifelse(constr, 5e-10, -Inf),                                      lower = ifelse(constr, 5e-10, -Inf),
421                                      control = list(iter.max = value$msMaxIter,                                      control = list(iter.max = value$msMaxIter,
422                                      trace = as.integer(value$msVerbose)))                                      trace = as.integer(value$msVerbose)))
423                   .Call("mer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")                   .Call(mer_coefGets, x, optimRes$par, 2)
424                   if (optimRes$convergence != 0) {                   if (optimRes$convergence != 0) {
425                       warning(paste("nlminb returned message",                       warning(paste("nlminb returned message",
426                                     optimRes$message,"\n"))                                     optimRes$message,"\n"))
# Line 476  Line 476 
476    
477  setMethod("deviance", signature(object = "mer"),  setMethod("deviance", signature(object = "mer"),
478            function(object, ...) {            function(object, ...) {
479                .Call("mer_factor", object, PACKAGE = "Matrix")                .Call(mer_factor, object)
480                object@deviance[[ifelse(object@method == "REML", "REML", "ML")]]                object@deviance[[ifelse(object@method == "REML", "REML", "ML")]]
481            })            })
482    
# Line 488  Line 488 
488            lmm <- family$family == "gaussian" && family$link == "identity"            lmm <- family$family == "gaussian" && family$link == "identity"
489            if (!lmm)            if (!lmm)
490                stop("mcmcsamp for GLMMs not yet implemented in supernodal representation")                stop("mcmcsamp for GLMMs not yet implemented in supernodal representation")
491            ans <- t(.Call("mer_MCMCsamp", object, saveb, n,            ans <- t(.Call(mer_MCMCsamp, object, saveb, n,
492                           trans, PACKAGE = "Matrix"))                           trans))
493            attr(ans, "mcpar") <- as.integer(c(1, n, 1))            attr(ans, "mcpar") <- as.integer(c(1, n, 1))
494            class(ans) <- "mcmc"            class(ans) <- "mcmc"
495            glmer <- FALSE            glmer <- FALSE
# Line 540  Line 540 
540                family$link != "identity")                family$link != "identity")
541                stop("simulation of generalized linear mixed models not yet implemented")                stop("simulation of generalized linear mixed models not yet implemented")
542            ## similate the linear predictors            ## similate the linear predictors
543            lpred <- .Call("mer_simulate", object, nsim, PACKAGE = "Matrix")            lpred <- .Call(mer_simulate, object, nsim)
544            sc <- 1            sc <- 1
545            if (object@useScale)            if (object@useScale)
546                sc <- .Call("mer_sigma", object, object@method == "REML",                sc <- .Call(mer_sigma, object, object@method == "REML",
547                            PACKAGE = "Matrix")                            PACKAGE = "Matrix")
548            ## add fixed-effects contribution and per-observation noise term            ## add fixed-effects contribution and per-observation noise term
549            lpred <- as.data.frame(lpred + drop(object@X %*% fixef(object)) +            lpred <- as.data.frame(lpred + drop(object@X %*% fixef(object)) +
# Line 644  Line 644 
644            function(object, REML = object@method == "REML",            function(object, REML = object@method == "REML",
645                     useScale = object@useScale,...) {                     useScale = object@useScale,...) {
646                sc <- if (object@useScale) {                sc <- if (object@useScale) {
647                    .Call("mer_sigma", object, REML, PACKAGE = "Matrix")                    .Call(mer_sigma, object, REML)
648                } else { 1 }                } else { 1 }
649                rr <- as(sc^2 * tcrossprod(solve(object@RXX)), "dpoMatrix")                rr <- as(sc^2 * tcrossprod(solve(object@RXX)), "dpoMatrix")
650                rr@factors$correlation <- as(rr, "correlation")                rr@factors$correlation <- as(rr, "correlation")
# Line 669  Line 669 
669                devc <- as.integer(object@devComp[1:2])                devc <- as.integer(object@devComp[1:2])
670                attr(val, "nall") <- attr(val, "nobs") <- devc[1]                attr(val, "nall") <- attr(val, "nobs") <- devc[1]
671                attr(val, "df") <- abs(devc[2]) +                attr(val, "df") <- abs(devc[2]) +
672                    length(.Call("mer_coef", object, 0, PACKAGE = "Matrix"))                    length(.Call(mer_coef, object, 0))
673                attr(val, "REML") <- REML                attr(val, "REML") <- REML
674                class(val) <- "logLik"                class(val) <- "logLik"
675                val                val
# Line 679  Line 679 
679            function(x, REML = x@method == "REML", useScale = x@useScale, ...)            function(x, REML = x@method == "REML", useScale = x@useScale, ...)
680        {        {
681            sc <- if (useScale)            sc <- if (useScale)
682                .Call("mer_sigma", x, REML, PACKAGE = "Matrix") else 1                .Call(mer_sigma, x, REML) else 1
683            sc2 <- sc * sc            sc2 <- sc * sc
684            cnames <- x@cnames            cnames <- x@cnames
685            ans <- x@Omega            ans <- x@Omega
# Line 742  Line 742 
742            else { ## ------ single model ---------------------            else { ## ------ single model ---------------------
743                foo <- object                foo <- object
744                #foo@status["factored"] <- FALSE                #foo@status["factored"] <- FALSE
745                #.Call("mer_factor", foo, PACKAGE="Matrix")                #.Call(mer_factor, foo)
746                #dfr <- getFixDF(foo)                #dfr <- getFixDF(foo)
747                ss <- foo@rXy^2                ss <- foo@rXy^2
748                ssr <- exp(foo@devComp["logryy2"])                ssr <- exp(foo@devComp["logryy2"])
# Line 779  Line 779 
779    
780  setMethod("fitted", signature(object = "mer"),  setMethod("fitted", signature(object = "mer"),
781            function(object, ...)            function(object, ...)
782            .Call("mer_fitted", object, PACKAGE = "Matrix")            .Call(mer_fitted, object)
783            )            )
784    
785  setMethod("formula", signature(x = "mer"),  setMethod("formula", signature(x = "mer"),
# Line 809  Line 809 
809  setMethod("summary", signature(object = "mer"),  setMethod("summary", signature(object = "mer"),
810            function(object, ...) {            function(object, ...) {
811    
812                fcoef <- .Call("mer_fixef", object, PACKAGE = "Matrix")                fcoef <- .Call(mer_fixef, object)
813                useScale <- object@useScale                useScale <- object@useScale
814                vcov <- vcov(object)                vcov <- vcov(object)
815                corF <- vcov@factors$correlation                corF <- vcov@factors$correlation
# Line 863  Line 863 
863                    methTitle = methTitle,                    methTitle = methTitle,
864                    logLik = llik,                    logLik = llik,
865                    ngrps = sapply(object@flist, function(x) length(levels(x))),                    ngrps = sapply(object@flist, function(x) length(levels(x))),
866                    sigma = .Call("mer_sigma", object, REML, PACKAGE = "Matrix"),                    sigma = .Call(mer_sigma, object, REML),
867                    coefs = coefs,                    coefs = coefs,
868                    vcov = vcov,                    vcov = vcov,
869                    REmat = REmat,                    REmat = REmat,
# Line 909  Line 909 
909      cv <- list(gradient = FALSE, msMaxIter = 200:200,      cv <- list(gradient = FALSE, msMaxIter = 200:200,
910                 msVerbose = 0:0)                 msVerbose = 0:0)
911      sapply(ysim, function(yy) {      sapply(ysim, function(yy) {
912          .Call("mer_update_y", fm0, yy, PACKAGE = "Matrix")          .Call(mer_update_y, fm0, yy)
913          LMEoptimize(fm0) <- cv          LMEoptimize(fm0) <- cv
914          .Call("mer_update_y", fma, yy, PACKAGE = "Matrix")          .Call(mer_update_y, fma, yy)
915          LMEoptimize(fma) <- cv          LMEoptimize(fma) <- cv
916          exp(c(H0 = fm0@devComp[["logryy2"]],          exp(c(H0 = fm0@devComp[["logryy2"]],
917                Ha = fma@devComp[["logryy2"]]))                Ha = fma@devComp[["logryy2"]]))
# Line 938  Line 938 
938                  cat(paste("Using", nAGQ, "quadrature points per column\n"))                  cat(paste("Using", nAGQ, "quadrature points per column\n"))
939          }          }
940          obj <- function(pars)          obj <- function(pars)
941              .Call("glmer_devAGQ", pars, GSpt, nAGQ, PACKAGE = "Matrix")              .Call(glmer_devAGQ, pars, GSpt, nAGQ)
942          optimRes <-          optimRes <-
943              nlminb(PQLpars, obj,              nlminb(PQLpars, obj,
944                     lower = ifelse(const, 5e-10, -Inf),                     lower = ifelse(const, 5e-10, -Inf),
# Line 951  Line 951 
951          if (gVerb)          if (gVerb)
952              cat(paste("convergence message", optimRes$message, "\n"))              cat(paste("convergence message", optimRes$message, "\n"))
953          fxd[] <- optpars[fixInd]  ## preserve the names          fxd[] <- optpars[fixInd]  ## preserve the names
954          .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")          .Call(lmer_coefGets, mer, optpars[-fixInd], 2)
955      }      }
956    
957      .Call("glmer_finalize", GSpt, PACKAGE = "Matrix")      .Call(glmer_finalize, GSpt)
958      loglik[] <- -deviance/2      loglik[] <- -deviance/2
959  }  }

Legend:
Removed from v.1279  
changed lines
  Added in v.1280

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge