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 754, Mon May 30 18:40:45 2005 UTC revision 755, Mon Jun 6 15:37:45 2005 UTC
# Line 48  Line 48 
48           analyticGradient=analyticGradient)           analyticGradient=analyticGradient)
49  }  }
50    
51  setMethod("lmer", signature(formula = "formula", family = "missing"),  setMethod("lmer", signature(formula = "formula"),
52            function(formula, data, family,            function(formula, data, family,
53                     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),                     method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
54                     control = list(),                     control = list(),
55                     subset, weights, na.action, offset,                     subset, weights, na.action, offset,
56                     model = TRUE, x = FALSE, y = FALSE, ...)                     model = TRUE, x = FALSE, y = FALSE, ...)
57        {        {                                 ## match and check parameters
58                                          # match and check parameters            if (length(formula) < 3) stop("formula must be a two-sided formula")
59              cv <- do.call("lmerControl", control)
60              if (lmm <- missing(family)) { # linear mixed model
61            method <- match.arg(method)            method <- match.arg(method)
62            if (method %in% c("PQL", "Laplace", "AGQ")) {            if (method %in% c("PQL", "Laplace", "AGQ")) {
63                warning(paste('Argument method = "', method,                warning(paste('Argument method = "', method,
# Line 63  Line 65 
65                              'Using method = "REML".\n', sep = ''))                              'Using method = "REML".\n', sep = ''))
66                method <- "REML"                method <- "REML"
67            }            }
68            controlvals <- do.call("lmerControl", control)            } else {                      # generalized linear mixed model
69            controlvals$REML <- REML <- method == "REML"                method <- if (missing(method)) "PQL" else match.arg(method)
70                  if (method == "ML") method <- "PQL"
71            if (length(formula) < 3) stop("formula must be a two-sided formula")                if (method == "REML")
72                      warning(paste('Argument method = "REML" is not meaningful',
73                                    'for a generalized linear mixed model.',
74                                    '\nUsing method = "PQL".\n'))
75                  if (method %in% c("AGQ"))
76                      stop("AGQ method not yet implemented")
77              }
78    
79            mf <- match.call()           # create the model frame as frm            ## evaluate glm.fit, a generalized linear fit of fixed effects only
80            m <- match(c("data", "subset", "weights", "na.action", "offset"),            mf <- match.call()
81                       names(mf), 0)            m <- match(c("family", "data", "subset", "weights",
82                           "na.action", "offset"), names(mf), 0)
83            mf <- mf[c(1, m)]            mf <- mf[c(1, m)]
84            mf[[1]] <- as.name("model.frame")            frame.form <- subbars(formula) # substitute `+' for `|'
85            frame.form <- subbars(formula)            fixed.form <- nobars(formula)  # remove any terms with `|'
86            environment(frame.form) <- environment(formula)            if (inherits(fixed.form, "name")) # RHS is empty - add a constant
87                  fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
88              environment(fixed.form) <- environment(frame.form) <- environment(formula)
89              mf$formula <- fixed.form
90              mf$x <- mf$model <- mf$y <- TRUE
91              mf[[1]] <- as.name("glm")
92              glm.fit <- eval(mf, parent.frame())
93    
94              ## evaluate a model frame for fixed and random effects
95            mf$formula <- frame.form            mf$formula <- frame.form
96              mf$x <- mf$model <- mf$y <- mf$family <- NULL
97            mf$drop.unused.levels <- TRUE            mf$drop.unused.levels <- TRUE
98              mf[[1]] <- as.name("model.frame")
99            frm <- eval(mf, parent.frame())            frm <- eval(mf, parent.frame())
100    
101            ## grouping factors and model matrices for random effects            ## grouping factors and model matrices for random effects
# Line 95  Line 114 
114            if (any(diff(nlev) > 0)) {            if (any(diff(nlev) > 0)) {
115                random <- random[rev(order(nlev))]                random <- random[rev(order(nlev))]
116            }            }
117            fixed.form <- nobars(formula)            ## Create the model matrices and a mixed-effects representation
           if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula  
           Xmat <- model.matrix(fixed.form, frm)  
118            mmats <- c(lapply(random, "[[", 1),            mmats <- c(lapply(random, "[[", 1),
119                       .fixed = list(cbind(Xmat, .response = model.response(frm))))                       .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
120            obj <- .Call("lmer_create", lapply(random, "[[", 2),            mer <- .Call("lmer_create", lapply(random, "[[", 2),
121                         mmats, PACKAGE = "Matrix")                         mmats, method, PACKAGE = "Matrix")
122            slot(obj, "frame") <- frm            if (lmm) {
123            slot(obj, "terms") <- attr(model.frame(fixed.form, data), "terms")                .Call("lmer_initial", mer, PACKAGE="Matrix")
124            slot(obj, "assign") <- attr(Xmat, "assign")                .Call("lmer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")
125            slot(obj, "call") <- match.call()                LMEoptimize(mer) <- cv
126            slot(obj, "REML") <- REML                fits <- .Call("lmer_fitted", mer, mmats, TRUE, PACKAGE = "Matrix")
127            rm(Xmat)                return(new("lmer", mer, frame = glm.fit$model,
128            .Call("lmer_initial", obj, PACKAGE="Matrix")                           terms = glm.fit$terms,
129            .Call("lmer_ECMEsteps", obj,                           assign = attr(glm.fit$x, "assign"), call = match.call(),
130                  controlvals$niterEM,                           residuals = unname(model.response(frm) - fits),
131                  controlvals$REML,                           fitted = fits))
132                  controlvals$EMverbose,            }
133    
134              ## The rest of the function applies to generalized linear mixed models
135              gVerb <- getOption("verbose")
136              family <- glm.fit$family
137              mmo <- mmats
138              mmats[[1]][1,1] <- mmats[[1]][1,1]
139              conv <- FALSE
140              firstIter <- TRUE
141              msMaxIter.orig <- cv$msMaxIter
142              responseIndex <- ncol(mmats$.fixed)
143              etaold <- eta <- glm.fit$linear.predictors
144              weights <- glm.fit$prior.weights
145              offset <- glm.fit$offset
146              if (is.null(offset)) offset <- numeric(length(eta))
147    
148              for (iter in seq(length = cv$PQLmaxIt))
149              {
150                  mu <- family$linkinv(eta)
151                  dmu.deta <- family$mu.eta(eta)
152                  ## weights (note: weights is already square-rooted)
153                  w <- weights * dmu.deta / sqrt(family$variance(mu))
154                  ## adjusted response (should be comparable to X \beta, not including offset
155                  z <- eta - offset + (mmo$.fixed[, responseIndex] - mu) / dmu.deta
156                  .Call("nlme_weight_matrix_list",
157                        mmo, w, z, mmats, PACKAGE="Matrix")
158                  .Call("lmer_update_mm", mer, mmats, PACKAGE="Matrix")
159                  if (firstIter) {
160                      .Call("lmer_initial", mer, PACKAGE="Matrix")
161                      if (gVerb) cat(" PQL iterations convergence criterion\n")
162                  }
163                  .Call("lmer_ECMEsteps", mer, cv$niterEM, cv$EMverbose, PACKAGE = "Matrix")
164                  LMEoptimize(mer) <- cv
165                  eta[] <- offset + .Call("lmer_fitted", mer,
166                                          mmo, TRUE, PACKAGE = "Matrix")
167                  crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta)))
168                  if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit))
169                  ## use this to determine convergence
170                  if (crit < cv$tolerance) {
171                      conv <- TRUE
172                      break
173                  }
174                  etaold[] <- eta
175                  if (firstIter) {          # Change the number of EM and optimization
176                      cv$niterEM <- 2       # iterations for subsequent PQL iterations.
177                      cv$msMaxIter <- 10
178                      firstIter <- FALSE
179                  }
180              }
181              if (!conv) warning("IRLS iterations for glmm did not converge")
182              cv$msMaxIter <- msMaxIter.orig
183    
184              reducedObj <- .Call("lmer_collapse", mer, PACKAGE = "Matrix")
185              reducedMmats.unadjusted <- mmo
186              reducedMmats.unadjusted$.fixed <-
187                  reducedMmats.unadjusted$.fixed[, responseIndex, drop = FALSE]
188              reducedMmats <- mmats
189              reducedMmats$.fixed <-
190                  reducedMmats$.fixed[, responseIndex, drop = FALSE]
191              fixInd <- seq(along = fixef(mer))
192    
193              bhat <- function(pars = NULL) { # 1:(responseIndex-1) - beta, rest - theta
194                  if (is.null(pars)) {
195                      off <- drop(glm.fit$x %*% fixef(mer)) + offset
196                  } else {
197                      .Call("lmer_coefGets", reducedObj,
198                            as.double(pars[responseIndex:length(pars)]),
199                            2, PACKAGE = "Matrix")
200                      off <- drop(glm.fit$x %*% pars[fixInd]) + offset
201                  }
202                  niter <- 20
203                  conv <- FALSE
204    
205                  eta <- offset + .Call("lmer_fitted", mer, mmo, TRUE, PACKAGE = "Matrix")
206                  etaold <- eta
207    
208                  for (iter in seq(length = niter)) {
209                      mu <- family$linkinv(eta)
210                      dmu.deta <- family$mu.eta(eta)
211                      w <- weights * dmu.deta / sqrt(family$variance(mu))
212                      z <- eta - off + (reducedMmats.unadjusted$.fixed[, 1]
213                                        - mu) / dmu.deta
214                      .Call("nlme_weight_matrix_list",
215                            reducedMmats.unadjusted, w, z, reducedMmats,
216                            PACKAGE="Matrix")
217                      .Call("lmer_update_mm", reducedObj, reducedMmats,
218                  PACKAGE = "Matrix")                  PACKAGE = "Matrix")
219            LMEoptimize(obj) <- controlvals                    eta[] <-
220            slot(obj, "residuals") <-                        off + .Call("lmer_fitted", reducedObj,
221                unname(model.response(frm) -                                    reducedMmats.unadjusted, TRUE,
222                       (slot(obj, "fitted") <-                                    PACKAGE = "Matrix")
223                        .Call("lmer_fitted", obj, mmats, TRUE, PACKAGE = "Matrix")))                    if (max(abs(eta - etaold)) <
224            obj                        (0.1 + max(abs(eta))) * cv$tolerance) {
225                          conv <- TRUE
226                          break
227                      }
228                      etaold[] <- eta
229                  }
230                  if (!conv) warning("iterations for bhat did not converge")
231    
232                  ## bhat doesn't really need to return anything, we
233                  ## just want the side-effect of modifying reducedObj
234                  ## In particular, we are interested in
235                  ## ranef(reducedObj) and reducedObj@bVar (?). But
236                  ## the mu-scale response will be useful for log-lik
237                  ## calculations later, so return them anyway
238    
239                  invisible(family$linkinv(eta))
240              }
241    
242              ## function that calculates log likelihood (the thing that
243              ## needs to get evaluated at each Gauss-Hermite location)
244    
245              ## log scale ? worry about details later, get the pieces in place
246    
247              ## this is for the Laplace approximation only. GH is more
248              ## complicated
249    
250              devLaplace <- function(pars = NULL) {
251                  ## FIXME: This actually returns half the deviance.
252                  mu <- bhat(pars = pars)
253                  ## gets correct values of bhat and bvars. As a side
254                  ## effect, mu now has fitted values
255    
256                  ## GLM family log likelihood (upto constant ?)(log scale)
257                  ## FIXME: need to adjust for sigma^2 for appropriate models (not trivial!)
258    
259                  ## Keep everything on (log) likelihood scale
260    
261                  ## log lik from observations given fixed and random effects
262                  ## get deviance, then multiply by -1/2 (since deviance = -2 log lik)
263                  sum(family$dev.resids(y = glm.fit$y,
264                                        mu = mu, wt = weights(glm.fit)^2))/2 -
265                      .Call("lmer_laplace_devComp", reducedObj, PACKAGE = "Matrix")
266              }
267    
268              if (method == "Laplace") {
269                  nc <- mer@nc
270                  const <- c(rep(FALSE, length(fixef(mer))),
271                             unlist(lapply(nc[1:(length(nc) - 2)],
272                                           function(k) 1:((k*(k+1))/2) <= k)))
273                  RV <- lapply(R.Version()[c("major", "minor")], as.numeric)
274                  if (RV$major == 2 && RV$minor >= 2.0) {
275                      optimRes <-
276                          nlminb(c(fixef(mer),
277                                   .Call("lmer_coef", mer, 2, PACKAGE = "Matrix")),
278                                 devLaplace,
279                                 lower = ifelse(const, 5e-10, -Inf),
280                                 control = list(trace = getOption("verbose"),
281                                 iter.max = cv$msMaxIter))
282                      optpars <- optimRes$par
283                      if (optimRes$convergence != 0)
284                          warning("nlminb failed to converge")
285                  } else {
286                      optimRes <-
287                          optim(c(fixef(mer),
288                                  .Call("lmer_coef", mer, 2, PACKAGE = "Matrix")),
289                                devLaplace, method = "L-BFGS-B",
290                                lower = ifelse(const, 5e-10, -Inf),
291                                control = list(trace = getOption("verbose"),
292                                     reltol = cv$msTol, maxit = cv$msMaxIter))
293                      optpars <- optimRes$par
294                      if (optimRes$convergence != 0)
295                          warning("optim failed to converge")
296                  }
297    
298                  if (getOption("verbose")) {
299                      cat(paste("optim convergence code",
300                                optimRes$convergence, "\n"))
301                      cat("Fixed effects:\n")
302                      print(fixef(mer))
303                      print(optimRes$par[seq(length = responseIndex - 1)])
304                      cat("(Box constrained) variance coefficients:\n")
305                      print(.Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))
306                      print(optimRes$par[responseIndex:length(optimRes$par)] )
307                  }
308              } else {
309                  optpars <- c(fixef(mer),
310                               .Call("lmer_coef", mer, 2, PACKAGE = "Matrix"))
311              }
312    
313              loglik <- -devLaplace(optpars)
314              ff <- optpars[fixInd]
315              names(ff) <- names(fixef(mer))
316              .Call("lmer_coefGets", mer, optpars[-fixInd], 2, PACKAGE = "Matrix")
317              new("glmer", new("lmer", mer, frame = frm, terms = glm.fit$terms,
318                  assign = attr(glm.fit$x, "assign"), call = match.call()),
319                  family = family, glmmll = loglik, fixed = ff)
320        })        })
321    
322  setReplaceMethod("LMEoptimize", signature(x="lmer", value="list"),  setReplaceMethod("LMEoptimize", signature(x="mer", value="list"),
323                   function(x, value)                   function(x, value)
324               {               {
325                   if (value$msMaxIter < 1) return(x)                   if (value$msMaxIter < 1) return(x)
326                   nc <- x@nc                   nc <- x@nc
327                   nc <- nc[1:(length(nc) - 2)]                   constr <- unlist(lapply(nc[1:(length(nc) - 2)],
328                   constr <- unlist(lapply(nc, function(k) 1:((k*(k+1))/2) <= k))                                           function(k) 1:((k*(k+1))/2) <= k))
329                   fn <- function(pars)                   fn <- function(pars)
330                       deviance(.Call("lmer_coefGets", x, pars,                       deviance(.Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix"))
331                                      2, PACKAGE = "Matrix"),                   gr <- NULL
332                                REML = value$REML)                   if (value$analyticGradient)
333                   gr <- if (value$analyticGradient)                       gr <-
334                       function(pars) {                       function(pars) {
335                           if (!isTRUE(all.equal(pars,                           if (!isTRUE(all.equal(pars,
336                                                 .Call("lmer_coef", x, 2))))                                                     .Call("lmer_coef", x,
337                               .Call("lmer_coefGets", x, pars, 2)                                                           2, PACKAGE = "Matrix"))))
338                           .Call("lmer_gradient", x, value$REML, 2)                                   .Call("lmer_coefGets", x, pars, 2, PACKAGE = "Matrix")
339                       } else NULL                               .Call("lmer_gradient", x, 2, PACKAGE = "Matrix")
340                             }
341                   RV <- lapply(R.Version()[c("major", "minor")], as.numeric)                   RV <- lapply(R.Version()[c("major", "minor")], as.numeric)
342                   if ((RV$major > 1 && RV$minor >= 2.0) ||                   if (RV$major == 2 && RV$minor >= 2.0) {
343                       require("port", quietly = TRUE)) {                       optimRes <- nlminb(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
                      optimRes <- nlminb(.Call("lmer_coef", x, 2),  
344                                          fn, gr,                                          fn, gr,
345                                          lower = ifelse(constr, 1e-10, -Inf),                                          lower = ifelse(constr, 5e-10, -Inf),
346                                          control = list(iter.max = value$msMaxIter,                                          control = list(iter.max = value$msMaxIter,
347                                          trace = as.integer(value$msVerbose)))                                          trace = as.integer(value$msVerbose)))
348                   } else {                   } else {
349                       optimRes <- optim(.Call("lmer_coef", x, 2), fn, gr,                       optimRes <- optim(.Call("lmer_coef", x, 2, PACKAGE = "Matrix"),
350                                         method = "L-BFGS-B",                                         fn, gr, method = "L-BFGS-B",
351                                         lower = ifelse(constr, 1e-10, -Inf),                                         lower = ifelse(constr, 5e-10, -Inf),
352                                         control = list(maxit = value$msMaxIter,                                         control = list(maxit = value$msMaxIter,
353                                         trace = as.integer(value$msVerbose)))                                         trace = as.integer(value$msVerbose)))
354                   }                   }
355                   .Call("lmer_coefGets", x, optimRes$par, 2)                   .Call("lmer_coefGets", x, optimRes$par, 2, PACKAGE = "Matrix")
356                   if (optimRes$convergence != 0) {                   if (optimRes$convergence != 0) {
357                       warning(paste("optim returned message",optimRes$message,"\n"))                       warning(paste("optim returned message",optimRes$message,"\n"))
358                   }                   }
# Line 169  Line 366 
366                                  data.frame, check.names = FALSE),                                  data.frame, check.names = FALSE),
367                           varFac = object@bVar,                           varFac = object@bVar,
368                           stdErr = .Call("lmer_sigma", object,                           stdErr = .Call("lmer_sigma", object,
369                           object@REML, PACKAGE = "Matrix"))                           object@method == "REML", PACKAGE = "Matrix"))
370                if (!accumulate || length(val@varFac) == 1) return(val)                if (!accumulate || length(val@varFac) == 1) return(val)
371                ## check for nested factors                ## check for nested factors
372                L <- object@L                L <- object@L
# Line 178  Line 375 
375                val                val
376            })            })
377    
378  setMethod("fixef", signature(object = "lmer"),  setMethod("fixef", signature(object = "mer"),
379            function(object, ...) {            function(object, ...) {
380                val <- .Call("lmer_fixef", object, PACKAGE = "Matrix")                val <- .Call("lmer_fixef", object, PACKAGE = "Matrix")
381                val[-length(val)]                val[-length(val)]
# Line 203  Line 400 
400            })            })
401    
402  setMethod("gradient", signature(x = "lmer"),  setMethod("gradient", signature(x = "lmer"),
403            function(x, REML, unconst, ...)            function(x, unconst, ...)
404            .Call("lmer_gradient", x, REML, unconst, PACKAGE = "Matrix"))            .Call("lmer_gradient", x, unconst, PACKAGE = "Matrix"))
405    
406  setMethod("summary", signature(object = "lmer"),  setMethod("summary", signature(object = "lmer"),
407            function(object, ...)            function(object, ...)
408            new("summary.lmer", object, useScale = TRUE,            new("summary.lmer", object, useScale = TRUE,
409                showCorrelation = TRUE,                showCorrelation = TRUE,
410                method = if (object@REML) "REML" else "ML",                method = object@method,
411                family = gaussian(),                family = gaussian(),
412                logLik = logLik(object),                logLik = logLik(object),
413                fixed = fixef(object)))                fixed = fixef(object)))
# Line 229  Line 426 
426            function(object)            function(object)
427            show(new("summary.lmer", object, useScale = TRUE,            show(new("summary.lmer", object, useScale = TRUE,
428                     showCorrelation = FALSE,                     showCorrelation = FALSE,
429                     method = if (object@REML) "REML" else "ML",                     method = object@method,
430                     family = gaussian(),                     family = gaussian(),
431                     logLik = logLik(object),                     logLik = logLik(object),
432                     fixed = fixef(object)))                     fixed = fixef(object)))
# Line 257  Line 454 
454                dimnames(coefs) <-                dimnames(coefs) <-
455                    list(names(fcoef), c("Estimate", "Std. Error", "DF"))                    list(names(fcoef), c("Estimate", "Std. Error", "DF"))
456                              digits <- max(3, getOption("digits") - 2)                              digits <- max(3, getOption("digits") - 2)
457                REML <- length(object@REML) > 0 && object@REML[1]                REML <- object@method == "REML"
458                llik <- object@logLik                llik <- object@logLik
459                dev <- object@deviance                dev <- object@deviance
460    
# Line 267  Line 464 
464                              object@method, "\n"))                              object@method, "\n"))
465                } else {                } else {
466                    cat("Linear mixed-effects model fit by ")                    cat("Linear mixed-effects model fit by ")
467                    cat(if(object@REML) "REML\n" else "maximum likelihood\n")                    cat(if(REML) "REML\n" else "maximum likelihood\n")
468                }                }
469                if (!is.null(object@call$formula)) {                if (!is.null(object@call$formula)) {
470                    cat("Formula:", deparse(object@call$formula),"\n")                    cat("Formula:", deparse(object@call$formula),"\n")
# Line 301  Line 498 
498                cat("\n")                cat("\n")
499                if (!useScale)                if (!useScale)
500                    cat("\nEstimated scale (compare to 1) ",                    cat("\nEstimated scale (compare to 1) ",
501                        .Call("lmer_sigma", object, object@REML, PACKAGE = "Matrix"),                        .Call("lmer_sigma", object, FALSE, PACKAGE = "Matrix"),
502                        "\n")                        "\n")
503                if (nrow(coefs) > 0) {                if (nrow(coefs) > 0) {
504                    if (useScale) {                    if (useScale) {
# Line 340  Line 537 
537            })            })
538    
539    
 setMethod("lmer", signature(formula = "formula"),  
           function(formula, family, data,  
                    method = c("REML", "ML", "PQL", "Laplace", "AGQ"),  
                    control = list(),  
                    subset, weights, na.action, offset,  
                    model = TRUE, x = FALSE, y = FALSE, ...)  
       {  
           if (missing(method)) {  
               method <- "PQL"  
           } else {  
               method <- match.arg(method)  
               if (method == "ML") method <- "PQL"  
               if (method == "REML")  
                   warning(paste('Argument method = "REML" is not meaningful',  
                                 'for a generalized linear mixed model.',  
                                 '\nUsing method = "PQL".\n'))  
           }  
 ##           if (method %in% c("Laplace", "AGQ"))  
 ##               stop("Laplace and AGQ methods not yet implemented")  
           if (method %in% c("AGQ"))  
               stop("AGQ method not yet implemented")  
           gVerb <- getOption("verbose")  
                                         # match and check parameters  
           controlvals <- do.call("lmerControl", control)  
           controlvals$REML <- FALSE  
           if (length(formula) < 3) stop("formula must be a two-sided formula")  
           ## initial glm fit  
           mf <- match.call()  
           m <- match(c("family", "data", "subset", "weights",  
                        "na.action", "offset"),  
                      names(mf), 0)  
           mf <- mf[c(1, m)]  
           mf[[1]] <- as.name("glm")  
           fixed.form <- nobars(formula)  
           if (!inherits(fixed.form, "formula")) fixed.form <- ~ 1 # default formula  
           environment(fixed.form) <- environment(formula)  
           mf$formula <- fixed.form  
           mf$x <- mf$model <- mf$y <- TRUE  
           glm.fit <- eval(mf, parent.frame())  
           family <- glm.fit$family  
           ## Note: offset is on the linear scale  
           offset <- glm.fit$offset  
           if (is.null(offset)) offset <- 0  
           weights <- sqrt(abs(glm.fit$prior.weights))  
           ## initial 'fitted' values on linear scale  
           etaold <- eta <- glm.fit$linear.predictors  
   
           ## evaluation of model frame  
           mf$x <- mf$model <- mf$y <- mf$family <- NULL  
           mf$drop.unused.levels <- TRUE  
           this.form <- subbars(formula)  
           environment(this.form) <- environment(formula)  
           mf$formula <- this.form  
           mf[[1]] <- as.name("model.frame")  
           frm <- eval(mf, parent.frame())  
   
           ## grouping factors and model matrices for random effects  
           bars <- findbars(formula[[3]])  
           random <-  
               lapply(bars,  
                      function(x) list(model.matrix(eval(substitute(~term,  
                                                                    list(term=x[[2]]))),  
                                                    frm),  
                                       eval(substitute(as.factor(fac)[,drop = TRUE],  
                                                       list(fac = x[[3]])), frm)))  
           names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))  
   
           ## order factor list by decreasing number of levels  
           nlev <- sapply(random, function(x) length(levels(x[[2]])))  
           if (any(diff(nlev) > 0)) {  
               random <- random[rev(order(nlev))]  
           }  
           mmats <- c(lapply(random, "[[", 1),  
                      .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))  
           ## FIXME: Use Xfrm and Xmat to get the terms and assign  
           ## slots, pass these to lmer_create, then destroy Xfrm, Xmat, etc.  
           obj <- .Call("lmer_create", lapply(random, "[[", 2),  
                        mmats, PACKAGE = "Matrix")  
           slot(obj, "frame") <- frm  
           slot(obj, "terms") <- attr(glm.fit$model, "terms")  
           slot(obj, "assign") <- attr(glm.fit$x, "assign")  
           slot(obj, "call") <- match.call()  
           slot(obj, "REML") <- FALSE  
           rm(glm.fit)  
           .Call("lmer_initial", obj, PACKAGE="Matrix")  
           mmats.unadjusted <- mmats  
           mmats[[1]][1,1] <- mmats[[1]][1,1]  
           conv <- FALSE  
           firstIter <- TRUE  
           msMaxIter.orig <- controlvals$msMaxIter  
           responseIndex <- ncol(mmats$.fixed)  
   
           for (iter in seq(length = controlvals$PQLmaxIt))  
           {  
               mu <- family$linkinv(eta)  
               dmu.deta <- family$mu.eta(eta)  
               ## weights (note: weights is already square-rooted)  
               w <- weights * dmu.deta / sqrt(family$variance(mu))  
               ## adjusted response (should be comparable to X \beta, not including offset  
               z <- eta - offset + (mmats.unadjusted$.fixed[, responseIndex] - mu) / dmu.deta  
               .Call("nlme_weight_matrix_list",  
                     mmats.unadjusted, w, z, mmats, PACKAGE="Matrix")  
               .Call("lmer_update_mm", obj, mmats, PACKAGE="Matrix")  
               if (firstIter) {  
                   .Call("lmer_initial", obj, PACKAGE="Matrix")  
                   if (gVerb) cat(" PQL iterations convergence criterion\n")  
               }  
               .Call("lmer_ECMEsteps", obj,  
                     controlvals$niterEM,  
                     FALSE,  
                     controlvals$EMverbose,  
                     PACKAGE = "Matrix")  
               LMEoptimize(obj) <- controlvals  
               eta[] <- offset + ## FIXME: should the offset be here ?  
                   .Call("lmer_fitted", obj,  
                         mmats.unadjusted, TRUE, PACKAGE = "Matrix")  
               crit <- max(abs(eta - etaold)) / (0.1 + max(abs(eta)))  
               if (gVerb) cat(sprintf("%03d: %#11g\n", as.integer(iter), crit))  
               ## use this to determine convergence  
               if (crit < controlvals$tolerance) {  
                   conv <- TRUE  
                   break  
               }  
               etaold[] <- eta  
   
               ## Changing number of iterations on second and  
               ## subsequent iterations.  
               if (firstIter)  
               {  
                   controlvals$niterEM <- 2  
                   controlvals$msMaxIter <- 10  
                   firstIter <- FALSE  
               }  
           }  
           if (!conv) warning("IRLS iterations for glmm did not converge")  
           controlvals$msMaxIter <- msMaxIter.orig  
   
   
 ###          if (TRUE) ## Laplace  
 ###          {  
           ## Need to optimize L(theta, beta) using Laplace approximation  
   
           ## Things needed for that:  
           ##  
           ## 1. reduced ssclme object, offset, weighted model matrices  
           ## 2. facs, reduced model matrices  
   
           ## Of these, those in 2 will be fixed given theta and beta,  
           ## and can be thought of arguments to the L(theta, beta)  
           ## function. However, the ones in 1 will have the same  
           ## structure. So the plan is to pre-allocate them and pass  
           ## them in too so they can be used without creating/copying  
           ## them more than once  
   
   
           ## reduced ssclme  
   
           reducedObj <- .Call("lmer_collapse", obj, PACKAGE = "Matrix")  
           reducedMmats.unadjusted <- mmats.unadjusted  
           reducedMmats.unadjusted$.fixed <-  
               reducedMmats.unadjusted$.fixed[, responseIndex, drop = FALSE]  
           reducedMmats <- mmats  
           reducedMmats$.fixed <-  
               reducedMmats$.fixed[, responseIndex, drop = FALSE]  
   
           ## define function that calculates bhats given theta and beta  
   
           bhat <-  
               function(pars = NULL) # 1:(responseIndex-1) - beta, rest - theta  
               {  
                   if (is.null(pars))  
                   {  
                       off <- drop(mmats.unadjusted$.fixed %*%  
                                   c(fixef(obj), 0)) + offset  
                   }  
                   else  
                   {  
                       .Call("lmer_coefGets",  
                             reducedObj,  
                             as.double(pars[responseIndex:length(pars)]),  
                             TRUE,  
                             PACKAGE = "Matrix")  
                       off <- drop(mmats.unadjusted$.fixed %*%  
                                   c(pars[1:(responseIndex-1)], 0)) + offset  
                   }  
                   niter <- 20  
                   conv <- FALSE  
   
                   eta <- offset +  
                       .Call("lmer_fitted",  
                             obj, mmats.unadjusted, TRUE,  
                             PACKAGE = "Matrix")  
                   etaold <- eta  
   
                   for (iter in seq(length = niter))  
                   {  
                       mu <- family$linkinv(eta)  
                       dmu.deta <- family$mu.eta(eta)  
                       w <- weights * dmu.deta / sqrt(family$variance(mu))  
                       z <- eta - off + (reducedMmats.unadjusted$.fixed[, 1]  
                                         - mu) / dmu.deta  
                       .Call("nlme_weight_matrix_list",  
                             reducedMmats.unadjusted, w, z, reducedMmats,  
                             PACKAGE="Matrix")  
                       .Call("lmer_update_mm",  
                             reducedObj, reducedMmats,  
                             PACKAGE="Matrix")  
                       eta[] <- off +  
                           .Call("lmer_fitted", reducedObj,  
                                 reducedMmats.unadjusted, TRUE,  
                                 PACKAGE = "Matrix")  
                       ##cat(paste("bhat Criterion:", max(abs(eta - etaold)) /  
                       ##          (0.1 + max(abs(eta))), "\n"))  
                       ## use this to determine convergence  
                       if (max(abs(eta - etaold)) <  
                           (0.1 + max(abs(eta))) * controlvals$tolerance)  
                       {  
                           conv <- TRUE  
                           break  
                       }  
                       etaold[] <- eta  
   
                   }  
                   if (!conv) warning("iterations for bhat did not converge")  
   
                   ## bhat doesn't really need to return anything, we  
                   ## just want the side-effect of modifying reducedObj  
                   ## In particular, we are interested in  
                   ## ranef(reducedObj) and reducedObj@bVar (?). But  
                   ## the mu-scale response will be useful for log-lik  
                   ## calculations later, so return them anyway  
   
                   invisible(family$linkinv(eta))  
               }  
   
           ## function that calculates log likelihood (the thing that  
           ## needs to get evaluated at each Gauss-Hermite location)  
   
           ## log scale ? worry about details later, get the pieces in place  
   
           ## this is for the Laplace approximation only. GH is more  
           ## complicated  
   
           devLaplace <- function(pars = NULL)  
           {  
               ## FIXME: This actually returns half the deviance.  
   
               ## gets correct values of bhat and bvars. As a side  
               ## effect, mu now has fitted values  
               mu <- bhat(pars = pars)  
   
               ## GLM family log likelihood (upto constant ?)(log scale)  
               ## FIXME: need to adjust for sigma^2 for appropriate models (not trivial!)  
   
               ## Keep everything on (log) likelihood scale  
   
               ## log lik from observations given fixed and random effects  
               ## get deviance, then multiply by -1/2 (since deviance = -2 log lik)  
               ans <- -sum(family$dev.resids(y = mmats.unadjusted$.fixed[, responseIndex],  
                                             mu = mu,  
                                             wt = weights^2))/2  
   
 ##               if (is.null(getOption("laplaceinR")))  
 ##               {  
               ans <- ans +  
                   .Call("lmer_laplace_devComp", reducedObj,  
                         PACKAGE = "Matrix")  
 ##               }  
 ##               else  
 ##               {  
 ##                   ranefs <- .Call("lmer_ranef", reducedObj, PACKAGE = "Matrix")  
 ##                   ## ans <- ans + reducedObj@devComp[2]/2 # log-determinant of Omega  
   
 ##                   Omega <- reducedObj@Omega  
 ##                   for (i in seq(along = ranefs))  
 ##                   {  
 ##                       ## contribution for random effects (get it working,  
 ##                       ## optimize later)  
 ##                       ## symmetrize RE variance  
 ##                       Omega[[i]] <- Omega[[i]] + t(Omega[[i]])  
 ##                       diag(Omega[[i]]) <- diag(Omega[[i]]) / 2  
   
 ##                       ## want log of `const det(Omega) exp(-1/2 b'  
 ##                       ## Omega b )` i.e., const + log det(Omega) - .5  
 ##                       ## * (b' Omega b)  
   
 ##                       ## FIXME: need to adjust for sigma^2 for appropriate  
 ##                       ## models (easy).  These are all the b'Omega b,  
 ##                       ## summed as they eventually need to be.  Think of  
 ##                       ## this as sum(rowSums((ranefs[[i]] %*% Omega[[i]])  
 ##                       ## * ranefs[[i]]))  
   
 ##                       ranef.loglik.det <- nrow(ranefs[[i]]) *  
 ##                           determinant(Omega[[i]], logarithm = TRUE)$modulus/2  
   
 ## #                      print(ranef.loglik.det)  
   
 ##                       ranef.loglik.re <-  
 ##                           -sum((ranefs[[i]] %*% Omega[[i]]) * ranefs[[i]])/2  
   
 ## #                      print(ranef.loglik.re)  
   
 ##                       ranef.loglik <- ranef.loglik.det + ranef.loglik.re  
   
 ##                       ## Jacobian adjustment  
 ##                       log.jacobian <-  
 ##                           sum(log(abs(apply(reducedObj@bVar[[i]],  
 ##                                             3,  
   
 ##                                             ## next line depends on  
 ##                                             ## whether bVars are variances  
 ##                                             ## or Cholesly factors  
   
 ##                                             ## function(x) sum(diag(x)))  
 ## ## Was this a bug?                          function(x) sum(diag( La.chol( x ) )))  
 ##                                             function(x) prod(diag( La.chol( x ) )))  
 ##                                       )))  
   
 ## #                      print(log.jacobian)  
   
   
 ##                       ## the constant terms from the r.e. and the final  
 ##                       ## Laplacian integral cancel out both being:  
 ##                       ## ranef.loglik.constant <- 0.5 * length(ranefs[[i]]) * log(2 * base::pi)  
   
 ##                       ans <- ans + ranef.loglik + log.jacobian  
 ##                   }  
 ##               }  
               ## ans is (up to some constant) log of the Laplacian  
               ## approximation of the likelihood. Return it's negative  
               ## to be minimized  
   
               ##              cat("Parameters: ")  
               ##              print(pars)  
   
               ##              cat("Value: ")  
               ##              print(as.double(-ans))  
   
               -ans  
           }  
   
           if (method == "Laplace")  
           {  
               RV <- lapply(R.Version()[c("major", "minor")], as.numeric)  
               if ((RV$major > 1 && RV$minor >= 2.0) ||  
                   require("port", quietly = TRUE)) {  
                   optimRes <-  
                       nlminb(c(fixef(obj),  
                                .Call("lmer_coef", obj, TRUE, PACKAGE = "Matrix")),  
                              devLaplace,  
                              control = list(trace = getOption("verbose"),  
                              iter.max = controlvals$msMaxIter))  
                   optpars <- optimRes$par  
                   if (optimRes$convergence != 0)  
                       warning("nlminb failed to converge")  
               } else {  
                   optimRes <-  
                       optim(c(fixef(obj),  
                               .Call("lmer_coef", obj, TRUE, PACKAGE = "Matrix")),  
                             devLaplace,  
                             method = "BFGS", hessian = TRUE,  
                             control = list(trace = getOption("verbose"),  
                             reltol = controlvals$msTol,  
                             maxit = controlvals$msMaxIter))  
                   optpars <- optimRes$par  
                   if (optimRes$convergence != 0)  
                       warning("optim failed to converge")  
                   Hessian <- optimRes$hessian  
               }  
   
               ##fixef(obj) <- optimRes$par[seq(length = responseIndex - 1)]  
               if (getOption("verbose")) {  
                   cat(paste("optim convergence code",  
                             optimRes$convergence, "\n"))  
                   cat("Fixed effects:\n")  
                   print(fixef(obj))  
                   print(optimRes$par[seq(length = responseIndex - 1)])  
                   cat("(Unconstrained) variance coefficients:\n")  
   
                   print(  
                         .Call("lmer_coef",  
                               obj,  
                               TRUE,  
                               PACKAGE = "Matrix"))  
   
                   ##coef(obj, unconst = TRUE) <-  
                   ##    optimRes$par[responseIndex:length(optimRes$par)]  
                   ##print(coef(obj, unconst = TRUE))  
                   print( optimRes$par[responseIndex:length(optimRes$par)] )  
               }  
   
               ## need to calculate likelihood.  also need to store  
               ## new estimates of fixed effects somewhere  
               ## (probably cannot update standard errors)  
 ###Rprof(NULL)  
           }  
           else  
           {  
               optpars <-  
                   c(fixef(obj),  
                     .Call("lmer_coef",  
                           obj,  
                           TRUE,  
                           PACKAGE = "Matrix"))  
               Hessian <- new("matrix")  
           }  
   
   
           ## Before finishing, we need to call devLaplace with the  
           ## optimum pars to get the final log likelihood (still need  
           ## to make sure it's the actual likelihood and not a  
           ## multiple). This would automatically call bhat() and hence  
           ## have the 'correct' random effects in reducedObj.  
   
           loglik <- -devLaplace(optpars)  
           ##print(loglik)  
           ff <- optpars[1:(responseIndex-1)]  
           names(ff) <- names(fixef(obj))  
   
           if (!x) mmats <- list()  
   
 ###      }  
   
           ##obj  
   
           new("glmer", obj,  
               family = family,  
               glmmll = loglik,  
               method = method,  
               fixed = ff)  
       })  
540    
541    
542  ## calculates degrees of freedom for fixed effects Wald tests  ## calculates degrees of freedom for fixed effects Wald tests
# Line 787  Line 553 
553            rep(n-p, p)            rep(n-p, p)
554        })        })
555    
556  setMethod("logLik", signature(object="lmer"),  setMethod("logLik", signature(object="mer"),
557            function(object, REML = object@REML, ...) {            function(object, REML = object@method == "REML", ...) {
558                val <- -deviance(object, REML = REML)/2                val <- -deviance(object, REML = REML)/2
559                nc <- object@nc[-seq(a = object@Omega)]                nc <- object@nc[-seq(a = object@Omega)]
560                attr(val, "nall") <- attr(val, "nobs") <- nc[2]                attr(val, "nall") <- attr(val, "nobs") <- nc[2]
561                attr(val, "df") <-                attr(val, "df") <- nc[1] +
562                    nc[1] + length(.Call("lmer_coef", object, 2))                    length(.Call("lmer_coef", object, 0, PACKAGE = "Matrix"))
563                attr(val, "REML") <- REML                attr(val, "REML") <- REML
564                class(val) <- "logLik"                class(val) <- "logLik"
565                val                val
# Line 804  Line 570 
570                val <- object@glmmll                val <- object@glmmll
571                nc <- object@nc[-seq(a = object@Omega)]                nc <- object@nc[-seq(a = object@Omega)]
572                attr(val, "nall") <- attr(val, "nobs") <- nc[2]                attr(val, "nall") <- attr(val, "nobs") <- nc[2]
573                attr(val, "df") <-                attr(val, "df") <- nc[1] +
574                    nc[1] + length(.Call("lmer_coef", object, 2))                    length(.Call("lmer_coef", object, 0, PACKAGE = "Matrix"))
575                class(val) <- "logLik"                class(val) <- "logLik"
576                val                val
577            })            })
# Line 931  Line 697 
697                .Call("lmer_coef", object, unconst, PACKAGE = "Matrix")                .Call("lmer_coef", object, unconst, PACKAGE = "Matrix")
698            })            })
699    
700  setMethod("deviance", "lmer",  setMethod("deviance", "mer",
701            function(object, REML = NULL, ...) {            function(object, REML = NULL, ...) {
702                .Call("lmer_factor", object, PACKAGE = "Matrix")                .Call("lmer_factor", object, PACKAGE = "Matrix")
703                if (is.null(REML))                if (is.null(REML))
704                    REML <- if (length(oR <- object@REML)) oR else FALSE                    REML <- object@method == "REML"
705                object@deviance[[ifelse(REML, "REML", "ML")]]                object@deviance[[ifelse(REML, "REML", "ML")]]
706            })            })
707    
# Line 959  Line 725 
725  setMethod("formula", "lmer", function(x, ...) x@call$formula)  setMethod("formula", "lmer", function(x, ...) x@call$formula)
726    
727  setMethod("vcov", signature(object = "lmer"),  setMethod("vcov", signature(object = "lmer"),
728            function(object, REML = object@REML, useScale = TRUE,...) {            function(object, REML = object@method == "REML", useScale = TRUE,...) {
729                sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")                sc <- .Call("lmer_sigma", object, REML, PACKAGE = "Matrix")
730                rr <- object@RXX                rr <- object@RXX
731                nms <- object@cnames[[".fixed"]]                nms <- object@cnames[[".fixed"]]

Legend:
Removed from v.754  
changed lines
  Added in v.755

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