SCM

SCM Repository

[matrix] Annotation of /branches/trunk-lme4/R/lmer.R
ViewVC logotype

Annotation of /branches/trunk-lme4/R/lmer.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 423 - (view) (download)

1 : bates 413 setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"),
2 : bates 316 function(x, value)
3 :     {
4 :     if (value$msMaxIter < 1) return(x)
5 :     st <- ccoef(x) # starting values
6 :     nc <- x@nc
7 :     nc <- nc[1:(length(nc) - 2)]
8 :     constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))
9 :     fn <- function(pars) {
10 :     ccoef(x) <- pars
11 :     deviance(x, REML = value$REML)
12 :     }
13 : bates 380 gr <- if (value$analyticGradient)
14 :     function(pars) {
15 :     ccoef(x) <- pars
16 :     grad <- lme4:::gradient(x, REML = value$REML, unconst = TRUE)
17 :     grad[constr] <- -grad[constr]/pars[constr]
18 :     grad
19 :     } else NULL
20 : bates 316 optimRes <- optim(st, fn, gr,
21 :     method = "L-BFGS-B",
22 :     lower = ifelse(constr, 1e-10, -Inf),
23 : bates 362 control = list(maxit = value$msMaxIter,
24 : bates 411 trace = as.integer(value$msVerbose)))
25 : bates 316 if (optimRes$convergence != 0) {
26 :     warning(paste("optim returned message",optimRes$message,"\n"))
27 :     }
28 : bates 411 ccoef(x) <- optimRes$par
29 : bates 316 return(x)
30 :     })
31 :    
32 : bates 413 setMethod("deviance", signature(object = "lmer"),
33 : bates 316 function(object, REML = FALSE, ...) {
34 : bates 413 .Call("lmer_factor", object, PACKAGE = "Matrix")
35 : bates 316 object@deviance[ifelse(REML, 2, 1)]
36 :     })
37 :    
38 : bates 413 setMethod("ranef", signature(object = "lmer"),
39 : bates 316 function(object, ...) {
40 : bates 413 .Call("lmer_ranef", object, PACKAGE = "Matrix")
41 : bates 316 })
42 :    
43 : bates 413 setMethod("fixef", signature(object = "lmer"),
44 : bates 316 function(object, ...) {
45 : bates 413 val = .Call("lmer_fixef", object, PACKAGE = "Matrix")
46 : bates 316 names(val) = object@cnames[[".fixed"]]
47 :     val[-length(val)]
48 :     })
49 :    
50 : bates 413 setMethod("vcov", signature(object = "lmer"),
51 : bates 316 function(object, REML = TRUE, useScale = TRUE,...) {
52 : bates 413 ## force an "lmer_invert"
53 :     sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
54 : bates 316 rr <- object@RXX
55 :     nms <- object@cnames[[".fixed"]]
56 :     dimnames(rr) <- list(nms, nms)
57 :     nr <- nrow(rr)
58 :     rr <- rr[-nr, -nr, drop = FALSE]
59 :     rr <- rr %*% t(rr)
60 :     if (useScale) {
61 :     rr = sc^2 * rr
62 :     }
63 :     rr
64 :     })
65 :    
66 : bates 413 setMethod("VarCorr", signature(x = "lmer"),
67 : bates 316 function(x, REML = TRUE, useScale = TRUE, ...) {
68 : bates 413 val = .Call("lmer_variances", x, PACKAGE = "Matrix")
69 : bates 316 for (i in seq(along = val)) {
70 :     dimnames(val[[i]]) = list(x@cnames[[i]], x@cnames[[i]])
71 :     val[[i]] = as(as(val[[i]], "pdmatrix"), "corrmatrix")
72 :     }
73 :     new("VarCorr",
74 : bates 413 scale = .Call("lmer_sigma", x, REML),
75 : bates 316 reSumry = val,
76 :     useScale = useScale)
77 :     })
78 :    
79 : bates 413 setMethod("gradient", signature(x = "lmer"),
80 : bates 316 function(x, REML, unconst, ...)
81 : bates 413 .Call("lmer_gradient", x, REML, unconst))
82 : bates 316
83 : bates 413 setMethod("summary", "lmer",
84 : bates 316 function(object, REML = TRUE, useScale = TRUE, ...) {
85 :     fcoef <- fixef(object)
86 :     corF <- as(as(vcov(object, REML, useScale), "pdmatrix"),
87 :     "corrmatrix")
88 :     DF <- getFixDF(object)
89 :     coefs <- cbind(fcoef, corF@stdDev, DF)
90 :     nc <- object@nc
91 :     dimnames(coefs) <-
92 :     list(names(fcoef), c("Estimate", "Std. Error", "DF"))
93 :     new("summary.ssclme",
94 :     coefficients = as.matrix(coefs),
95 : bates 413 scale = .Call("lmer_sigma", object, REML),
96 : bates 316 denomDF = as.integer(DF),
97 :     REML = REML,
98 : bates 423 ngrps = unlist(lapply(object@flist,
99 :     function(x) length(levels(x)))),
100 : bates 316 nobs = nc[length(nc)],
101 :     corFixed = corF,
102 :     VarCorr = VarCorr(object, REML, useScale),
103 :     useScale = useScale,
104 :     showCorrelation = FALSE)
105 :     })
106 :    
107 : bates 413 setMethod("show", "lmer",
108 : bates 316 function(object) {
109 :     fcoef <- fixef(object)
110 :     corF <- as(as(vcov(object, REML = TRUE, useScale = TRUE), "pdmatrix"),
111 :     "corrmatrix")
112 :     DF <- getFixDF(object)
113 :     coefs <- cbind(fcoef, corF@stdDev, DF)
114 :     nc <- object@nc
115 :     dimnames(coefs) <-
116 :     list(names(fcoef), c("Estimate", "Std. Error", "DF"))
117 :     new("summary.ssclme",
118 :     coefficients = as.matrix(coefs),
119 : bates 413 scale = .Call("lmer_sigma", object, REML = TRUE),
120 : bates 316 denomDF = as.integer(DF),
121 :     REML = TRUE,
122 : bates 423 ngrps = unlist(lapply(object@flist,
123 :     function(x) length(levels(x)))),
124 : bates 316 nobs = nc[length(nc)],
125 :     corFixed = corF,
126 :     VarCorr = VarCorr(object, REML = TRUE, useScale = TRUE),
127 :     useScale = TRUE,
128 :     showCorrelation = FALSE)
129 :     })
130 :    
131 :     ## calculates degrees of freedom for fixed effects Wald tests
132 :     ## This is a placeholder. The answers are generally wrong. It will
133 :     ## be very tricky to decide what a 'right' answer should be with
134 :     ## crossed random effects.
135 :    
136 : bates 413 setMethod("getFixDF", signature(object="lmer"),
137 : bates 316 function(object, ...)
138 :     {
139 :     nc <- object@nc[-seq(along = object@Omega)]
140 :     p <- nc[1] - 1
141 :     n <- nc[2]
142 :     rep(n-p, p)
143 :     })
144 :    
145 : bates 413 setMethod("anova", signature(object="lmer"),
146 : bates 316 function(object, ...)
147 :     {
148 :     foo <- object
149 :     foo@status["factored"] <- FALSE
150 : bates 413 .Call("lmer_factor", foo, PACKAGE="Matrix")
151 : bates 316 dfr <- getFixDF(foo)
152 :     ss <- foo@RXX[ , ".response"]^2
153 :     ssr <- ss[[".response"]]
154 :     ss <- ss[seq(along = dfr)]
155 :     names(ss) <- object@cnames[[".fixed"]][seq(along = dfr)]
156 :     # FIXME: This only gives single degree of freedom tests
157 :     ms <- ss
158 :     df <- rep(1, length(ms))
159 :     f <- ms/(ssr/dfr)
160 :     P <- pf(f, df, dfr, lower.tail = FALSE)
161 :     table <- data.frame(df, ss, ms, dfr, f, P)
162 :     dimnames(table) <-
163 :     list(names(ss),
164 :     c("Df", "Sum Sq", "Mean Sq", "Denom", "F value", "Pr(>F)"))
165 :     if (any(match(names(ss), "(Intercept)", nomatch = 0)))
166 :     table <- table[-1,]
167 :     attr(table, "heading") <- "Analysis of Variance Table"
168 :     class(table) <- c("anova", "data.frame")
169 :     table
170 :     })

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