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 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:
Vienna University of Economics and Business Powered By FusionForge