SCM

SCM Repository

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

Annotation of /pkg/R/lmer.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 380 - (view) (download)
Original Path: branches/trunk-lme4/R/lmeRep.R

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

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