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 445, Wed Jan 19 19:06:24 2005 UTC revision 446, Wed Jan 19 19:08:01 2005 UTC
# Line 1  Line 1 
1    contr.SAS <- function(n, contrasts = TRUE)
2    {
3        if (is.numeric(n) && length(n) == 1) contr.treatment(n, n, contrasts)
4        else contr.treatment(n, length(n), contrasts)
5    }
6    
7  lmerControl <-                            # Control parameters for lmer  lmerControl <-                            # Control parameters for lmer
8    function(maxIter = 50,    function(maxIter = 50,
9             msMaxIter = 50,             msMaxIter = 50,
# Line 68  Line 74 
74                       .fixed = list(cbind(model.matrix(nobars(formula), frm),                       .fixed = list(cbind(model.matrix(nobars(formula), frm),
75                       .response = model.response(frm))))                       .response = model.response(frm))))
76            obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")            obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")
77              obj@call <- match.call()
78              obj@REML <- method == "REML"
79            .Call("lmer_initial", obj, PACKAGE="Matrix")            .Call("lmer_initial", obj, PACKAGE="Matrix")
80            .Call("lmer_ECMEsteps", obj,            .Call("lmer_ECMEsteps", obj,
81                  controlvals$niterEM,                  controlvals$niterEM,
# Line 75  Line 83 
83                  controlvals$EMverbose,                  controlvals$EMverbose,
84                  PACKAGE = "Matrix")                  PACKAGE = "Matrix")
85            LMEoptimize(obj) <- controlvals            LMEoptimize(obj) <- controlvals
           obj@call <- match.call()  
86            #fitted = .Call("ssclme_fitted", obj, facs, mmats, TRUE, PACKAGE = "Matrix")            #fitted = .Call("ssclme_fitted", obj, facs, mmats, TRUE, PACKAGE = "Matrix")
87            #residuals = mmats$.Xy[,".response"] - fitted            #residuals = mmats$.Xy[,".response"] - fitted
88            #if (as.logical(x)[1]) x = mmats else x = list()            #if (as.logical(x)[1]) x = mmats else x = list()
# Line 175  Line 182 
182                nc <- object@nc                nc <- object@nc
183                dimnames(coefs) <-                dimnames(coefs) <-
184                    list(names(fcoef), c("Estimate", "Std. Error", "DF"))                    list(names(fcoef), c("Estimate", "Std. Error", "DF"))
185                new("summary.ssclme",                new("summary.lmer",
186                    coefficients = as.matrix(coefs),                    coefficients = as.matrix(coefs),
187                    scale = .Call("lmer_sigma", object, REML),                    scale = .Call("lmer_sigma", object, REML),
188                    denomDF = as.integer(DF),                    denomDF = as.integer(DF),
# Line 227  Line 234 
234            rep(n-p, p)            rep(n-p, p)
235        })        })
236    
237    setMethod("fitted", signature(object="lmer"),
238              function(object, ...)
239          {
240              object@fitted
241          })
242    
243    setMethod("residuals", signature(object="lmer"),
244              function(object, ...) object@residuals )
245    
246    setMethod("logLik", signature(object="lmer"),
247              function(object, REML = object@REML, ...) {
248                  val <- -deviance(object, REML = REML)/2
249                  nc <- object@nc[-seq(a = object@Omega)]
250                  attr(val, "nall") <- attr(val, "nobs") <- nc[2]
251                  attr(val, "df") <- nc[1] + length(coef(object))
252                  attr(val, "REML") <- REML
253                  class(val) <- "logLik"
254                  val
255              })
256    
257    # setMethod("summary", signature(object="lme"),
258    #           function(object, ...) {
259    #               llik <- logLik(object)
260    #               resd <- residuals(object, type="pearson")
261    #               if (length(resd) > 5) {
262    #                   resd <- quantile(resd)
263    #                   names(resd) <- c("Min","Q1","Med","Q3","Max")
264    #               }
265    #               new("summary.lme",
266    #                   call = object@call,
267    #                   logLik = llik,
268    #                   re = summary(object@rep, REML = object@REML,
269    #                                useScale = TRUE),
270    #                   residuals = resd)
271    #           })
272    
273    # setMethod("show", signature(object = "summary.lme"),
274    #           function(object)
275    #       {
276    #           rdig <- 5
277    #           cat("Linear mixed-effects model fit by ")
278    #           cat(ifelse(object@re@REML, "REML\n", "maximum likelihood\n") )
279    #           if (!is.null(object@call$formula)) {
280    #               cat("Fixed:", deparse(object@call$formula),"\n")
281    #           }
282    #           if (!is.null(object@call$data)) {
283    #               cat(" Data:", deparse(object@call$data), "\n")
284    #           }
285    #           if (!is.null(object@call$subset)) {
286    #               cat(" Subset:",
287    #                   deparse(asOneSidedFormula(object@call$subset)[[2]]),"\n")
288    #           }
289    #           llik <- object@logLik
290    #           print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
291    #                            logLik = c(object@logLik), row.names = ""))
292    #           cat("\n")
293    #           object@re@useScale <- TRUE
294    #           object@re@showCorrelation <- TRUE
295    #           show(object@re)
296    #           invisible(object)
297    #       })
298    
299  setMethod("anova", signature(object="lmer"),  setMethod("anova", signature(object="lmer"),
300            function(object, ...)            function(object, ...)
301        {        {
302              mCall <- match.call(expand.dots = TRUE)
303              dots <- list(...)
304              modp <- logical(0)
305              if (length(dots))
306                  modp <- sapply(dots, inherits, "lmer") | sapply(dots, inherits, "lm")
307              if (any(modp)) {              # multiple models - form table
308                  opts <- dots[!modp]
309                  mods <- c(list(object), dots[modp])
310                  names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)], as.character)
311                  mods <- mods[order(sapply(lapply(mods, logLik, REML = FALSE), attr, "df"))]
312                  calls <- lapply(mods, slot, "call")
313                  data <- lapply(calls, "[[", "data")
314                  if (any(data != data[[1]])) stop("all models must be fit to the same data object")
315                  header <- paste("Data:", data[[1]])
316                  subset <- lapply(calls, "[[", "subset")
317                  if (any(subset != subset[[1]])) stop("all models must use the same subset")
318                  if (!is.null(subset[[1]]))
319                      header <-
320                          c(header, paste("Subset", deparse(subset[[1]]), sep = ": "))
321                  llks <- lapply(mods, logLik, REML = FALSE)
322                  Df <- sapply(llks, attr, "df")
323                  llk <- unlist(llks)
324                  chisq <- 2 * pmax(0, c(NA, diff(llk)))
325                  dfChisq <- c(NA, diff(Df))
326                  val <- data.frame(Df = Df,
327                                    AIC = sapply(llks, AIC),
328                                    BIC = sapply(llks, BIC),
329                                    logLik = llk,
330                                    "Chisq" = chisq,
331                                    "Chi Df" = dfChisq,
332                                    "Pr(>Chisq)" = pchisq(chisq, dfChisq, lower = FALSE),
333                                    check.names = FALSE)
334                  class(val) <- c("anova", class(val))
335                  attr(val, "heading") <-
336                      c(header, "", "Models:",
337                        paste(names(mods),
338                              unlist(lapply(lapply(calls, "[[", "formula"), deparse)),
339                              sep = ": "),"")
340                  return(val)
341              } else {
342    #               mCall$terms <- object@terms
343    #               mCall$assign <- object@assign
344            foo <- object            foo <- object
345            foo@status["factored"] <- FALSE            foo@status["factored"] <- FALSE
346            .Call("lmer_factor", foo, PACKAGE="Matrix")            .Call("lmer_factor", foo, PACKAGE="Matrix")
347            dfr <- getFixDF(foo)            dfr <- lme4:::getFixDF(foo)
348            ss <- foo@RXX[ , ".response"]^2            rcol <- ncol(foo@RXX)
349            ssr <- ss[[".response"]]            ss <- foo@RXX[ , rcol]^2
350              ssr <- ss[[rcol]]
351            ss <- ss[seq(along = dfr)]            ss <- ss[seq(along = dfr)]
352            names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]            names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]
353            # FIXME: This only gives single degree of freedom tests            # FIXME: This only gives single degree of freedom tests
# Line 252  Line 364 
364            attr(table, "heading") <- "Analysis of Variance Table"            attr(table, "heading") <- "Analysis of Variance Table"
365            class(table) <- c("anova", "data.frame")            class(table) <- c("anova", "data.frame")
366            table            table
367    
368              }
369        })        })
370    
371    setMethod("formula", "lmer", function(x, ...) x@call$formula)
372    
373    setMethod("plot", signature(x = "lmer"),
374              function(x, y, ...)
375              cat("plot method for lmer not yet implemented\n"))
376    
377    setMethod("update", signature(object = "lmer"),
378              function(object, formula., ..., evaluate = TRUE)
379          {
380              call <- object@call
381              if (is.null(call))
382                  stop("need an object with call component")
383              extras <- match.call(expand.dots = FALSE)$...
384              if (!missing(formula.))
385                  call$formula <- update.formula(formula(object), formula.)
386              if (length(extras) > 0) {
387                  existing <- !is.na(match(names(extras), names(call)))
388                  for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
389                  if (any(!existing)) {
390                      call <- c(as.list(call), extras[!existing])
391                      call <- as.call(call)
392                  }
393              }
394              if (evaluate)
395                  eval(call, parent.frame())
396              else call
397          })
398    
399    
400    setMethod("confint", signature(object = "lmer"),
401              function (object, parm, level = 0.95, ...)
402          {
403              cf <- fixef(object)
404              pnames <- names(cf)
405              if (missing(parm))
406                  parm <- seq(along = pnames)
407              else if (is.character(parm))
408                  parm <- match(parm, pnames, nomatch = 0)
409              a <- (1 - level)/2
410              a <- c(a, 1 - a)
411              pct <- paste(round(100 * a, 1), "%")
412              ci <- array(NA, dim = c(length(parm), 2),
413                          dimnames = list(pnames[parm], pct))
414              ses <- sqrt(diag(vcov(object)))[parm]
415              ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object@rep)[parm], qt))
416              ci
417          })
418    
419    
420    
421    

Legend:
Removed from v.445  
changed lines
  Added in v.446

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