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

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