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 766, Tue Jun 7 23:03:34 2005 UTC revision 767, Tue Jun 7 23:05:00 2005 UTC
# Line 1  Line 1 
1  ## Methods for lmer and for the objects that it produces  # Methods for lmer and for the objects that it produces
2    
3  ## Some utilities  ## Some utilities
4    
# Line 7  Line 7 
7      ((i - 1) * i)/2 + j      ((i - 1) * i)/2 + j
8  }  }
9    
10  Dhalf <- function(from) {  ##Dhalf <- function(from) {
11      D <- from@D  ##    D <- from@D
12      nf <- length(D)  ##    nf <- length(D)
13      Gp <- from@Gp  ##    Gp <- from@Gp
14      res <- array(0, rep(Gp[nf+1],2))  ##    res <- array(0, rep(Gp[nf+1],2))
15      for (i in 1:nf) {  ##    for (i in 1:nf) {
16          DD <- D[[i]]  ##        DD <- D[[i]]
17          dd <- dim(DD)  ##        dd <- dim(DD)
18          for (k in 1:dd[3]) {  ##        for (k in 1:dd[3]) {
19              mm <- array(DD[ , , k], dd[1:2])  ##            mm <- array(DD[ , , k], dd[1:2])
20              base <- Gp[i] + (k - 1)*dd[1]  ##            base <- Gp[i] + (k - 1)*dd[1]
21              res[cbind(c(base + row(mm)), c(base + col(mm)))] <- c(mm)  ##            res[cbind(c(base + row(mm)), c(base + col(mm)))] <- c(mm)
22          }  ##        }
23      }  ##    }
24      res  ##    res
25  }  ##}
26    
27  lmerControl <-                            # Control parameters for lmer  lmerControl <-                            # Control parameters for lmer
28    function(maxIter = 50,    function(maxIter = 50,
# Line 83  Line 83 
83            mf <- mf[c(1, m)]            mf <- mf[c(1, m)]
84            frame.form <- subbars(formula) # substitute `+' for `|'            frame.form <- subbars(formula) # substitute `+' for `|'
85            fixed.form <- nobars(formula)  # remove any terms with `|'            fixed.form <- nobars(formula)  # remove any terms with `|'
86            if (inherits(fixed.form, "name")) # RHS is empty - add a constant            if (inherits(fixed.form, "name")) # RHS is empty - use a constant
87                fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))                fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
88            environment(fixed.form) <- environment(frame.form) <- environment(formula)            environment(fixed.form) <- environment(frame.form) <- environment(formula)
89            mf$formula <- fixed.form            mf$formula <- fixed.form
90            mf$x <- mf$model <- mf$y <- TRUE            mf$x <- mf$model <- mf$y <- TRUE
91            mf[[1]] <- as.name("glm")            mf[[1]] <- as.name("glm")
92            glm.fit <- eval(mf, parent.frame())            glm.fit <- eval(mf, parent.frame())
93              family <- glm.fit$family
94              x <- glm.fit$x
95              y <- as.double(glm.fit$y)
96    
97            ## evaluate a model frame for fixed and random effects            ## evaluate a model frame for fixed and random effects
98            mf$formula <- frame.form            mf$formula <- frame.form
# Line 114  Line 117 
117            if (any(diff(nlev) > 0)) {            if (any(diff(nlev) > 0)) {
118                random <- random[rev(order(nlev))]                random <- random[rev(order(nlev))]
119            }            }
120            ## Create the model matrices and a mixed-effects representation  
121              ## Create the model matrices and a mixed-effects representation (mer)
122            mmats <- c(lapply(random, "[[", 1),            mmats <- c(lapply(random, "[[", 1),
123                       .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))                       .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
124            mer <- .Call("lmer_create", lapply(random, "[[", 2),            mer <- .Call("lmer_create", lapply(random, "[[", 2),
125                         mmats, method, PACKAGE = "Matrix")                         mmats, method, PACKAGE = "Matrix")
126            if (lmm) {            if (lmm) {                    ## linear mixed model
127                .Call("lmer_initial", mer, PACKAGE="Matrix")                .Call("lmer_initial", mer, PACKAGE="Matrix")
128                .Call("lmer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")                .Call("lmer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")
129                LMEoptimize(mer) <- cv                LMEoptimize(mer) <- cv
130                fits <- .Call("lmer_fitted", mer, mmats, TRUE, PACKAGE = "Matrix")                fits <- .Call("lmer_fitted", mer, mmats, TRUE, PACKAGE = "Matrix")
131                return(new("lmer", mer, frame = glm.fit$model,                return(new("lmer",
132                             mer, frame = if (frame) frm else data.frame(),
133                           terms = glm.fit$terms,                           terms = glm.fit$terms,
134                           assign = attr(glm.fit$x, "assign"), call = match.call(),                           assign = attr(x, "assign"),
135                             call = match.call(),
136                           residuals = unname(model.response(frm) - fits),                           residuals = unname(model.response(frm) - fits),
137                           fitted = fits))                           fitted = fits))
138            }            }
139    
140            ## The rest of the function applies to generalized linear mixed models            ## The rest of the function applies to generalized linear mixed models
141            gVerb <- getOption("verbose")            gVerb <- getOption("verbose")
142            family <- glm.fit$family            etaold <- eta <- glm.fit$linear.predictors
143              wts <- glm.fit$prior.weights
144              offset <- glm.fit$offset
145              if (is.null(offset)) offset <- numeric(length(eta))
146    
147              dev.resids <- quote(family$dev.resids(y, mu, wts))
148              linkinv <- quote(family$linkinv(eta))
149              mu.eta <- quote(family$mu.eta(eta))
150              variance <- quote(family$variance(mu))
151    
152            mmo <- mmats            mmo <- mmats
153            mmats[[1]][1,1] <- mmats[[1]][1,1]            mmats[[1]][1,1] <- mmats[[1]][1,1]
154            conv <- FALSE            conv <- FALSE
155            firstIter <- TRUE            firstIter <- TRUE
156            msMaxIter.orig <- cv$msMaxIter            msMaxIter.orig <- cv$msMaxIter
           responseIndex <- ncol(mmats$.fixed)  
           etaold <- eta <- glm.fit$linear.predictors  
           weights <- glm.fit$prior.weights  
           offset <- glm.fit$offset  
           if (is.null(offset)) offset <- numeric(length(eta))  
157    
158            for (iter in seq(length = cv$PQLmaxIt))            for (iter in seq(length = cv$PQLmaxIt))
159            {            {
160                mu <- family$linkinv(eta)                mu <- eval(linkinv) # family$linkinv(eta)
161                dmu.deta <- family$mu.eta(eta)                dmu.deta <- eval(mu.eta) # family$mu.eta(eta)
162                ## weights (note: weights is already square-rooted)                ## weights (note: wts is already square-rooted)
163                w <- weights * dmu.deta / sqrt(family$variance(mu))                .Call("glmer_weight_matrix_list", mmo,
164                ## adjusted response (should be comparable to X \beta, not including offset                      wts * dmu.deta / sqrt(eval(variance)), ## weights
165                z <- eta - offset + (mmo$.fixed[, responseIndex] - mu) / dmu.deta                      eta - offset + (y - mu) / dmu.deta, ## working residual
166                .Call("nlme_weight_matrix_list",                      mmats, PACKAGE="Matrix")
                     mmo, w, z, mmats, PACKAGE="Matrix")  
167                .Call("lmer_update_mm", mer, mmats, PACKAGE="Matrix")                .Call("lmer_update_mm", mer, mmats, PACKAGE="Matrix")
168                if (firstIter) {                if (firstIter) {
169                    .Call("lmer_initial", mer, PACKAGE="Matrix")                    .Call("lmer_initial", mer, PACKAGE="Matrix")
170                    if (gVerb) cat(" PQL iterations convergence criterion\n")                    if (gVerb) cat(" PQL iterations convergence criterion\n")
171                }                }
172                .Call("lmer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")                .Call("lmer_ECMEsteps", mer, cv$niterEM, cv$EMverbose,
173                        PACKAGE = "Matrix")
174                LMEoptimize(mer) <- cv                LMEoptimize(mer) <- cv
175                eta[] <- offset + .Call("lmer_fitted", mer,                eta <- offset + .Call("lmer_fitted", mer, mmo, TRUE,
176                                        mmo, TRUE, PACKAGE = "Matrix")                                      PACKAGE = "Matrix")
177                crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta)))                crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta)))
178                if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit))                if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit))
179                ## use this to determine convergence                ## use this to determine convergence
# Line 181  Line 191 
191            if (!conv) warning("IRLS iterations for glmm did not converge")            if (!conv) warning("IRLS iterations for glmm did not converge")
192            cv$msMaxIter <- msMaxIter.orig            cv$msMaxIter <- msMaxIter.orig
193    
194            reducedObj <- .Call("lmer_collapse", mer, PACKAGE = "Matrix")            rdobj <- .Call("lmer_collapse", mer, PACKAGE = "Matrix")
195            reducedMmats.unadjusted <- mmo            unwtd <- mmo
196            reducedMmats.unadjusted$.fixed <-            unwtd$.fixed <- as.matrix(y)
197                reducedMmats.unadjusted$.fixed[, responseIndex, drop = FALSE]            wtd <- unwtd
198            reducedMmats <- mmats            wtd[[1]][1,1] <- wtd[[1]][1,1]
           reducedMmats$.fixed <-  
               reducedMmats$.fixed[, responseIndex, drop = FALSE]  
199            fixInd <- seq(along = fixef(mer))            fixInd <- seq(along = fixef(mer))
200              env <- environment()
201            bhat <- function(pars = NULL) { # 1:(responseIndex-1) - beta, rest - theta            bhat <- function(pars) { #pars[fixInd] == beta, pars[-fixInd] == theta
202                if (is.null(pars)) {                if (cv$analyticHessian) {
203                    off <- drop(glm.fit$x %*% fixef(mer)) + offset                    .Call("glmer_bhat_iterate", pars, cv$tolerance,
204                            env, PACKAGE = "Matrix")
205                } else {                } else {
206                    .Call("lmer_coefGets", reducedObj,                    etaold <- off <- eta <- drop(glm.fit$x %*% pars[fixInd]) + offset
207                          as.double(pars[responseIndex:length(pars)]),                    .Call("lmer_coefGets", rdobj, as.double(pars[-fixInd]),
208                          2, PACKAGE = "Matrix")                          2, PACKAGE = "Matrix")
                   off <- drop(glm.fit$x %*% pars[fixInd]) + offset  
               }  
209                niter <- 20                niter <- 20
210                conv <- FALSE                conv <- FALSE
   
               eta <- offset + .Call("lmer_fitted", mer, mmo, TRUE, PACKAGE = "Matrix")  
               etaold <- eta  
   
211                for (iter in seq(length = niter)) {                for (iter in seq(length = niter)) {
212                    mu <- family$linkinv(eta)                        mu <- eval(linkinv)
213                    dmu.deta <- family$mu.eta(eta)                        dmu.deta <- eval(mu.eta)
214                    w <- weights * dmu.deta / sqrt(family$variance(mu))                        .Call("glmer_weight_matrix_list", unwtd,
215                    z <- eta - off + (reducedMmats.unadjusted$.fixed[, 1]                              wts * dmu.deta / sqrt(eval(variance)),
216                                      - mu) / dmu.deta                              eta - off + (y - mu)/dmu.deta, wtd,
                   .Call("nlme_weight_matrix_list",  
                         reducedMmats.unadjusted, w, z, reducedMmats,  
217                          PACKAGE="Matrix")                          PACKAGE="Matrix")
218                    .Call("lmer_update_mm", reducedObj, reducedMmats,                        .Call("lmer_update_mm", rdobj, wtd,
                         PACKAGE="Matrix")  
                   eta[] <-  
                       off + .Call("lmer_fitted", reducedObj,  
                                   reducedMmats.unadjusted, TRUE,  
219                                    PACKAGE = "Matrix")                                    PACKAGE = "Matrix")
220                          eta <- off + .Call("lmer_fitted", rdobj,
221                                             unwtd, TRUE, PACKAGE = "Matrix")
222                    if (max(abs(eta - etaold)) <                    if (max(abs(eta - etaold)) <
223                        (0.1 + max(abs(eta))) * cv$tolerance) {                        (0.1 + max(abs(eta))) * cv$tolerance) {
224                        conv <- TRUE                        conv <- TRUE
225                        break                        break
226                    }                    }
227                    etaold[] <- eta                        etaold <- eta
228                }                }
229                if (!conv) warning("iterations for bhat did not converge")                if (!conv) warning("iterations for bhat did not converge")
230                  }
231                ## bhat doesn't really need to return anything, we                ## bhat doesn't really need to return anything, we
232                ## just want the side-effect of modifying reducedObj                ## just want the side-effect of modifying rdobj
233                ## In particular, we are interested in                ## In particular, we are interested in
234                ## ranef(reducedObj) and reducedObj@bVar (?). But                ## ranef(rdobj) and rdobj@bVar (?). But
235                ## the mu-scale response will be useful for log-lik                ## the mu-scale response will be useful for log-lik
236                ## calculations later, so return them anyway                ## calculations later, so return them anyway
237    
238                invisible(family$linkinv(eta))                invisible(eval(linkinv))
239            }            }
240    
241    
242            ## function that calculates log likelihood (the thing that            ## function that calculates log likelihood (the thing that
243            ## needs to get evaluated at each Gauss-Hermite location)            ## needs to get evaluated at each Gauss-Hermite location)
244    
# Line 247  Line 247 
247            ## this is for the Laplace approximation only. GH is more            ## this is for the Laplace approximation only. GH is more
248            ## complicated            ## complicated
249    
250            devLaplace <- function(pars = NULL) {            devLaplace <- function(pars) {
251                ## FIXME: This actually returns half the deviance.                mu <- bhat(pars)
               mu <- bhat(pars = pars)  
               ## gets correct values of bhat and bvars. As a side  
               ## effect, mu now has fitted values  
252    
253                ## GLM family log likelihood (upto constant ?)(log scale)                ## GLM family log likelihood (upto constant ?)(log scale)
254                ## FIXME: need to adjust for sigma^2 for appropriate models (not trivial!)                ## FIXME: need to adjust for sigma^2 for appropriate models (not trivial!)
# Line 260  Line 257 
257    
258                ## log lik from observations given fixed and random effects                ## log lik from observations given fixed and random effects
259                ## get deviance, then multiply by -1/2 (since deviance = -2 log lik)                ## get deviance, then multiply by -1/2 (since deviance = -2 log lik)
260                sum(family$dev.resids(y = glm.fit$y,                ## FIXME: This actually returns half the deviance.
261                                      mu = mu, wt = weights(glm.fit)^2))/2 -                dev <- sum(family$dev.resids(y = glm.fit$y,
262                    .Call("lmer_laplace_devComp", reducedObj, PACKAGE = "Matrix")                                             mu = mu, wt = wts^2))/2
263                  lap <- .Call("glmer_Laplace_devComp", rdobj, PACKAGE = "Matrix")
264                  if (gVerb) cat(sprintf("  %#11g %#11g\n", dev, lap))
265                  dev - lap
266            }            }
267    
268            if (method == "Laplace") {            if (method == "Laplace") {
# Line 299  Line 299 
299                    cat(paste("optim convergence code",                    cat(paste("optim convergence code",
300                              optimRes$convergence, "\n"))                              optimRes$convergence, "\n"))
301                    cat("Fixed effects:\n")                    cat("Fixed effects:\n")
302                    print(fixef(mer))                    #print(fixef(mer))  ## no longer correct values
303                    print(optimRes$par[seq(length = responseIndex - 1)])                    print(optimRes$par[fixInd])
304                    cat("(Box constrained) variance coefficients:\n")                    cat("(box constrained) variance coefficients:\n")
305                    print(.Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))                    #print(.Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))
306                    print(optimRes$par[responseIndex:length(optimRes$par)] )                    print(optimRes$par[-fixInd])
307                }                }
308            } else {            } else {
309                optpars <- c(fixef(mer),                optpars <- c(fixef(mer),
# Line 868  Line 868 
868    
869  setMethod("terms", signature(x = "lmer"),  setMethod("terms", signature(x = "lmer"),
870            function(x, ...) x@terms)            function(x, ...) x@terms)
871    
872    setMethod("show", signature(object="VarCorr"),
873              function(object)
874          {
875              digits <- max(3, getOption("digits") - 2)
876              useScale <- length(object@useScale) > 0 && object@useScale[1]
877              sc <- ifelse(useScale, object@scale,  1.)
878              reStdDev <- c(lapply(object@reSumry,
879                                   function(x, sc)
880                                   sc*x@stdDev,
881                                   sc = sc), list(Residual = sc))
882              reLens <- unlist(c(lapply(reStdDev, length)))
883              reMat <- array('', c(sum(reLens), 4),
884              list(rep('', sum(reLens)),
885                   c("Groups", "Name", "Variance", "Std.Dev.")))
886              reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens)
887              reMat[,2] <- c(unlist(lapply(reStdDev, names)), "")
888              reMat[,3] <- format(unlist(reStdDev)^2, digits = digits)
889              reMat[,4] <- format(unlist(reStdDev), digits = digits)
890              if (any(reLens > 1)) {
891                  maxlen <- max(reLens)
892                  corr <-
893                      do.call("rbind",
894                              lapply(object@reSumry,
895                                     function(x, maxlen) {
896                                         cc <- format(round(x, 3), nsmall = 3)
897                                         cc[!lower.tri(cc)] <- ""
898                                         nr <- dim(cc)[1]
899                                         if (nr >= maxlen) return(cc)
900                                         cbind(cc, matrix("", nr, maxlen-nr))
901                                     }, maxlen))
902                  colnames(corr) <- c("Corr", rep("", maxlen - 1))
903                  reMat <- cbind(reMat, rbind(corr, rep("", ncol(corr))))
904              }
905              if (!useScale) reMat <- reMat[-nrow(reMat),]
906              print(reMat, quote = FALSE)
907          })

Legend:
Removed from v.766  
changed lines
  Added in v.767

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