# SCM Repository

[matrix] Diff of /branches/trunk-lme4/R/lmer.R
 [matrix] / branches / trunk-lme4 / R / lmer.R

# Diff of /branches/trunk-lme4/R/lmer.R

revision 446, Wed Jan 19 19:08:01 2005 UTC revision 449, Fri Jan 21 22:01:21 2005 UTC
# Line 1  Line 1
1  contr.SAS <- function(n, contrasts = TRUE)  contr.SAS <- function(n, contrasts = TRUE)
2  {  {
3      if (is.numeric(n) && length(n) == 1) contr.treatment(n, n, contrasts)      contr.treatment(n, if (is.numeric(n) && length(n) == 1) n else length(n), contrasts)
else contr.treatment(n, length(n), contrasts)
4  }  }
5
6  lmerControl <-                            # Control parameters for lmer  lmerControl <-                            # Control parameters for lmer
# Line 38  Line 37
37                     model = TRUE, x = FALSE, y = FALSE, ...)                     model = TRUE, x = FALSE, y = FALSE, ...)
38        {        {
39                                          # match and check parameters                                          # match and check parameters
40            method <- match.arg(method)            REML <- match.arg(method) == "REML"
41            controlvals <- do.call("lmerControl", control)            controlvals <- do.call("lmerControl", control)
42            controlvals\$REML <- method == "REML"            controlvals\$REML <- REML
43            if (length(formula) < 3) stop("formula must be a two-sided formula")            if (length(formula) < 3) stop("formula must be a two-sided formula")
44                                          # create the model frame as frm
45            mf <- match.call()            mf <- match.call()           # create the model frame as frm
46            m <- match(c("data", "subset", "weights", "na.action", "offset"),            m <- match(c("data", "subset", "weights", "na.action", "offset"),
47                       names(mf), 0)                       names(mf), 0)
48            mf <- mf[c(1, m)]            mf <- mf[c(1, m)]
# Line 66  Line 65
65            names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))            names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
66
67            ## order factor list by decreasing number of levels            ## order factor list by decreasing number of levels
68            ford <- rev(order(sapply(random, function(x) length(levels(x[[2]])))))            nlev <- sapply(random, function(x) length(levels(x[[2]])))
69            if (any(ford != seq(a = random))) { # re-order both facs and random            if (any(diff(nlev) < 0)) {
70                random <- random[ford]                random <- random[rev(order(nlev))]
71            }            }
72            mmats <- c(lapply(random, "[[", 1),            mmats <- c(lapply(random, "[[", 1),
73                       .fixed = list(cbind(model.matrix(nobars(formula), frm),                       .fixed = list(cbind(model.matrix(nobars(formula), frm),
74                       .response = model.response(frm))))                       .response = model.response(frm))))
75            obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")            obj <- .Call("lmer_create", lapply(random, "[[", 2), mmats, PACKAGE = "Matrix")
76            obj@call <- match.call()            obj@call <- match.call()
77            obj@REML <- method == "REML"            obj@REML <- REML
78            .Call("lmer_initial", obj, PACKAGE="Matrix")            .Call("lmer_initial", obj, PACKAGE="Matrix")
79            .Call("lmer_ECMEsteps", obj,            .Call("lmer_ECMEsteps", obj,
80                  controlvals\$niterEM,                  controlvals\$niterEM,
# Line 83  Line 82
82                  controlvals\$EMverbose,                  controlvals\$EMverbose,
83                  PACKAGE = "Matrix")                  PACKAGE = "Matrix")
84            LMEoptimize(obj) <- controlvals            LMEoptimize(obj) <- controlvals
85            #fitted = .Call("ssclme_fitted", obj, facs, mmats, TRUE, PACKAGE = "Matrix")            #fitted <- .Call("ssclme_fitted", obj, facs, mmats, TRUE, PACKAGE = "Matrix")
86            #residuals = mmats\$.Xy[,".response"] - fitted            #residuals <- mmats\$.Xy[,".response"] - fitted
87            #if (as.logical(x)[1]) x = mmats else x = list()            #if (as.logical(x)[1]) x <- mmats else x = list()
88            #rm(mmats)            #rm(mmats)
89            obj            obj
90        })        })
# Line 121  Line 120
120                   return(x)                   return(x)
121               })               })
122
setMethod("deviance", signature(object = "lmer"),
function(object, REML = FALSE, ...) {
.Call("lmer_factor", object, PACKAGE = "Matrix")
object@deviance[ifelse(REML, 2, 1)]
})

123  setMethod("ranef", signature(object = "lmer"),  setMethod("ranef", signature(object = "lmer"),
124            function(object, ...) {            function(object, ...) {
125                .Call("lmer_ranef", object, PACKAGE = "Matrix")                .Call("lmer_ranef", object, PACKAGE = "Matrix")
# Line 139  Line 132
132                val[-length(val)]                val[-length(val)]
133            })            })
134
setMethod("vcov", signature(object = "lmer"),
function(object, REML = TRUE, useScale = TRUE,...) {
## force an "lmer_invert"
sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
rr <- object@RXX
nms <- object@cnames[[".fixed"]]
dimnames(rr) <- list(nms, nms)
nr <- nrow(rr)
rr <- rr[-nr, -nr, drop = FALSE]
rr <- rr %*% t(rr)
if (useScale) {
rr = sc^2 * rr
}
rr
})

135  setMethod("VarCorr", signature(x = "lmer"),  setMethod("VarCorr", signature(x = "lmer"),
136            function(x, REML = TRUE, useScale = TRUE, ...) {            function(x, REML = TRUE, useScale = TRUE, ...) {
137                val = .Call("lmer_variances", x, PACKAGE = "Matrix")                val = .Call("lmer_variances", x, PACKAGE = "Matrix")
# Line 163  Line 140
140                    val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")                    val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")
141                }                }
142                new("VarCorr",                new("VarCorr",
143                    scale = .Call("lmer_sigma", x, REML),                    scale = .Call("lmer_sigma", x, REML, PACKAGE = "Matrix"),
144                    reSumry = val,                    reSumry = val,
145                    useScale = useScale)                    useScale = useScale)
146            })            })
147
148  setMethod("gradient", signature(x = "lmer"),  setMethod("gradient", signature(x = "lmer"),
149            function(x, REML, unconst, ...)            function(x, REML, unconst, ...)
150            .Call("lmer_gradient", x, REML, unconst))            .Call("lmer_gradient", x, REML, unconst, PACKAGE = "Matrix"))
151
152  setMethod("summary", "lmer",  setMethod("summary", signature(object = "lmer"),
153            function(object, REML = TRUE, useScale = TRUE, ...) {            function(object, ...)
154                fcoef <- fixef(object)            new("summary.lmer", object, useScale = TRUE, showCorrelation = TRUE))
155                corF <- as(as(vcov(object, REML, useScale), "pdmatrix"),
156                           "corrmatrix")  setMethod("show", signature(object = "lmer"),
157                DF <- getFixDF(object)            function(object)
158                coefs <- cbind(fcoef, corF@stdDev, DF)            show(new("summary.lmer", object, useScale = TRUE, showCorrelation = FALSE))
159                nc <- object@nc            )
dimnames(coefs) <-
list(names(fcoef), c("Estimate", "Std. Error", "DF"))
new("summary.lmer",
coefficients = as.matrix(coefs),
scale = .Call("lmer_sigma", object, REML),
denomDF = as.integer(DF),
REML = REML,
ngrps = unlist(lapply(object@flist,
function(x) length(levels(x)))),
nobs = nc[length(nc)],
corFixed = corF,
VarCorr = VarCorr(object, REML, useScale),
useScale = useScale,
showCorrelation = FALSE)
})
160
161  setMethod("show", "lmer",  setMethod("show", "summary.lmer",
162            function(object) {            function(object) {
163                fcoef <- fixef(object)                fcoef <- fixef(object)
164                corF <- as(as(vcov(object, REML = TRUE, useScale = TRUE), "pdmatrix"),                useScale <- object@useScale
165                  corF <- as(as(vcov(object, useScale = useScale), "pdmatrix"),
166                           "corrmatrix")                           "corrmatrix")
167                DF <- getFixDF(object)                DF <- getFixDF(object)
168                coefs <- cbind(fcoef, corF@stdDev, DF)                coefs <- cbind(fcoef, corF@stdDev, DF)
169                nc <- object@nc                nc <- object@nc
170                dimnames(coefs) <-                dimnames(coefs) <-
171                    list(names(fcoef), c("Estimate", "Std. Error", "DF"))                    list(names(fcoef), c("Estimate", "Std. Error", "DF"))
172                new("summary.ssclme",                              digits <- max(3, getOption("digits") - 2)
173                    coefficients = as.matrix(coefs),                REML <- length(object@REML) > 0 && object@REML[1]
174                    scale = .Call("lmer_sigma", object, REML = TRUE),                llik <- logLik(object)
175                    denomDF = as.integer(DF),                dev <- object@deviance
176                    REML = TRUE,
177                    ngrps = unlist(lapply(object@flist,                rdig <- 5
178                                          function(x) length(levels(x)))),                cat("Linear mixed-effects model fit by ")
179                    nobs = nc[length(nc)],                cat(ifelse(object@REML, "REML\n", "maximum likelihood\n") )
180                    corFixed = corF,                if (!is.null(object@call\$formula)) {
181                    VarCorr = VarCorr(object, REML = TRUE, useScale = TRUE),                    cat("Formula:", deparse(object@call\$formula),"\n")
182                    useScale = TRUE,                }
183                    showCorrelation = FALSE)                if (!is.null(object@call\$data)) {
184                      cat("   Data:", deparse(object@call\$data), "\n")
185                  }
186                  if (!is.null(object@call\$subset)) {
187                      cat(" Subset:",
188                          deparse(asOneSidedFormula(object@call\$subset)[[2]]),"\n")
189                  }
190                  print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
191                                   logLik = c(llik),
192                                   MLdeviance = dev["ML"],
193                                   REMLdeviance = dev["REML"],
194                                   row.names = ""))
195                  cat("Random effects:\n")
196                  show(VarCorr(object))
197                  ngrps <- lapply(object@flist, function(x) length(levels(x)))
198                  cat(sprintf("# of obs: %d, groups: ", object@nc[length(object@nc)]))
199                  cat(paste(paste(names(ngrps), ngrps, sep = ", "), collapse = "; "))
200                  cat("\n")
201                  if (!useScale)
202                      cat("\nEstimated scale (compare to 1) ",
203                          .Call("lmer_sigma", object, object@REML, PACKAGE = "Matrix"),
204                          "\n")
205                  if (nrow(coefs) > 0) {
206                      if (useScale) {
207                          stat <- coefs[,1]/coefs[,2]
208                          pval <- 2*pt(abs(stat), coefs[,3], lower = FALSE)
209                          nms <- colnames(coefs)
210                          coefs <- cbind(coefs, stat, pval)
211                          colnames(coefs) <- c(nms, "t value", "Pr(>|t|)")
212                      } else {
213                          coefs <- coefs[, 1:2, drop = FALSE]
214                          stat <- coefs[,1]/coefs[,2]
215                          pval <- 2*pnorm(abs(stat), lower = FALSE)
216                          nms <- colnames(coefs)
217                          coefs <- cbind(coefs, stat, pval)
218                          colnames(coefs) <- c(nms, "z value", "Pr(>|z|)")
219                      }
220                      cat("\nFixed effects:\n")
221                      printCoefmat(coefs, tst.ind = 4, zap.ind = 3)
222                      if (length(object@showCorrelation) > 0 && object@showCorrelation[1]) {
223                          rn <- rownames(coefs)
224                          dimnames(corF) <- list(
225                                                   abbreviate(rn, minlen=11),
226                                                   abbreviate(rn, minlen=6))
227                          if (!is.null(corF)) {
228                              p <- NCOL(corF)
229                              if (p > 1) {
230                                  cat("\nCorrelation of Fixed Effects:\n")
231                                  corF <- format(round(corF, 3), nsmall = 3)
232                                  corF[!lower.tri(corF)] <- ""
233                                  print(corF[-1, -p, drop=FALSE], quote = FALSE)
234                              }
235                          }
236                      }
237                  }
238                  invisible(object)
239            })            })
240
241  ## calculates degrees of freedom for fixed effects Wald tests  ## calculates degrees of freedom for fixed effects Wald tests
# Line 234  Line 252
252            rep(n-p, p)            rep(n-p, p)
253        })        })
254
setMethod("fitted", signature(object="lmer"),
function(object, ...)
{
object@fitted
})

setMethod("residuals", signature(object="lmer"),
function(object, ...) object@residuals )

255  setMethod("logLik", signature(object="lmer"),  setMethod("logLik", signature(object="lmer"),
256            function(object, REML = object@REML, ...) {            function(object, REML = object@REML, ...) {
257                val <- -deviance(object, REML = REML)/2                val <- -deviance(object, REML = REML)/2
# Line 254  Line 263
263                val                val
264            })            })
265
# setMethod("summary", signature(object="lme"),
#           function(object, ...) {
#               llik <- logLik(object)
#               resd <- residuals(object, type="pearson")
#               if (length(resd) > 5) {
#                   resd <- quantile(resd)
#                   names(resd) <- c("Min","Q1","Med","Q3","Max")
#               }
#               new("summary.lme",
#                   call = object@call,
#                   logLik = llik,
#                   re = summary(object@rep, REML = object@REML,
#                                useScale = TRUE),
#                   residuals = resd)
#           })

# setMethod("show", signature(object = "summary.lme"),
#           function(object)
#       {
#           rdig <- 5
#           cat("Linear mixed-effects model fit by ")
#           cat(ifelse(object@re@REML, "REML\n", "maximum likelihood\n") )
#           if (!is.null(object@call\$formula)) {
#               cat("Fixed:", deparse(object@call\$formula),"\n")
#           }
#           if (!is.null(object@call\$data)) {
#               cat(" Data:", deparse(object@call\$data), "\n")
#           }
#           if (!is.null(object@call\$subset)) {
#               cat(" Subset:",
#                   deparse(asOneSidedFormula(object@call\$subset)[[2]]),"\n")
#           }
#           llik <- object@logLik
#           print(data.frame(AIC = AIC(llik), BIC = BIC(llik),
#                            logLik = c(object@logLik), row.names = ""))
#           cat("\n")
#           object@re@useScale <- TRUE
#           object@re@showCorrelation <- TRUE
#           show(object@re)
#           invisible(object)
#       })

266  setMethod("anova", signature(object = "lmer"),  setMethod("anova", signature(object = "lmer"),
267            function(object, ...)            function(object, ...)
268        {        {
# Line 344  Line 311
311            foo <- object            foo <- object
312            foo@status["factored"] <- FALSE            foo@status["factored"] <- FALSE
313            .Call("lmer_factor", foo, PACKAGE="Matrix")            .Call("lmer_factor", foo, PACKAGE="Matrix")
314            dfr <- lme4:::getFixDF(foo)            dfr <- getFixDF(foo)
315            rcol <- ncol(foo@RXX)            rcol <- ncol(foo@RXX)
316            ss <- foo@RXX[ , rcol]^2            ss <- foo@RXX[ , rcol]^2
317            ssr <- ss[[rcol]]            ssr <- ss[[rcol]]
# Line 368  Line 335
335            }            }
336        })        })
337
setMethod("formula", "lmer", function(x, ...) x@call\$formula)

setMethod("plot", signature(x = "lmer"),
function(x, y, ...)
cat("plot method for lmer not yet implemented\n"))

338  setMethod("update", signature(object = "lmer"),  setMethod("update", signature(object = "lmer"),
339            function(object, formula., ..., evaluate = TRUE)            function(object, formula., ..., evaluate = TRUE)
340        {        {
# Line 412  Line 373
373            ci <- array(NA, dim = c(length(parm), 2),            ci <- array(NA, dim = c(length(parm), 2),
374                        dimnames = list(pnames[parm], pct))                        dimnames = list(pnames[parm], pct))
375            ses <- sqrt(diag(vcov(object)))[parm]            ses <- sqrt(diag(vcov(object)))[parm]
376            ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object@rep)[parm], qt))            ci[] <- cf[parm] + ses * t(outer(a, getFixDF(object)[parm], qt))
377            ci            ci
378        })        })
379
380    setReplaceMethod("coef", signature(object = "lmer", value = "numeric"),
381                     function(object, unconst = FALSE, ..., value)
382                     .Call("lmer_coefGets", object, as.double(value),
383                           unconst, PACKAGE = "Matrix"))
384
385    setMethod("coef", signature(object = "lmer"),
386              function(object, unconst = FALSE, ...) {
387                  .Call("lmer_coef", object, unconst, PACKAGE = "Matrix")
388              })
389
390    setMethod("deviance", "lmer",
391              function(object, REML = NULL, ...) {
392                  .Call("lmer_factor", object, PACKAGE = "Matrix")
393                  if (is.null(REML))
394                      REML <- if (length(oR <- object@REML)) oR else FALSE
395                  object@deviance[[ifelse(REML, "REML", "ML")]]
396              })
397
398    setMethod("chol", signature(x = "lmer"),
399              function(x, pivot = FALSE, LINPACK = pivot) {
400                  x@status["factored"] <- FALSE # force a decomposition
401                  .Call("lmer_factor", x, PACKAGE = "Matrix")
402              })
403
404    setMethod("solve", signature(a = "lmer", b = "missing"),
405              function(a, b, ...)
406              .Call("lmer_invert", a)
407              )
408
409    setMethod("formula", "lmer", function(x, ...) x@call\$formula)
410
411    setMethod("vcov", signature(object = "lmer"),
412              function(object, REML = object@REML, useScale = TRUE,...) {
413                  ## force an "lmer_invert"
414                  sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
415                  rr <- object@RXX
416                  nms <- object@cnames[[".fixed"]]
417                  dimnames(rr) <- list(nms, nms)
418                  nr <- nrow(rr)
419                  rr <- rr[-nr, -nr, drop = FALSE]
420                  rr <- rr %*% t(rr)
421                  if (useScale) {
422                      rr = sc^2 * rr
423                  }
424                  rr
425              })
426

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

 root@r-forge.r-project.org ViewVC Help Powered by ViewVC 1.0.0
Thanks to: