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 939, Tue Sep 20 17:23:24 2005 UTC revision 940, Tue Sep 27 00:50:35 2005 UTC
# Line 1066  Line 1066 
1066      attr(ref, "observed") <- obs      attr(ref, "observed") <- obs
1067      ref      ref
1068  }  }
1069    
1070    mer2 <-
1071        function(formula, data, family,
1072                 method = c("REML", "ML", "PQL", "Laplace", "AGQ"),
1073                 control = list(), start,
1074                 subset, weights, na.action, offset,
1075                 model = TRUE, x = FALSE, y = FALSE , ...)
1076              ## x, y : not dealt with at all -- FIXME ? .NotYetImplemented(
1077    {
1078        ## match and check parameters
1079        if (length(formula) < 3) stop("formula must be a two-sided formula")
1080        ## cv <- do.call("Matrix:::lmerControl", control)
1081        cv <- control
1082        ## evaluate glm.fit, a generalized linear fit of fixed effects only
1083        mf <- match.call()
1084        m <- match(c("family", "data", "subset", "weights",
1085                     "na.action", "offset"), names(mf), 0)
1086        mf <- mf[c(1, m)]
1087        frame.form <- subbars(formula) # substitute `+' for `|'
1088        fixed.form <- nobars(formula)  # remove any terms with `|'
1089        if (inherits(fixed.form, "name")) # RHS is empty - use a constant
1090            fixed.form <- substitute(foo ~ 1, list(foo = fixed.form))
1091        environment(fixed.form) <- environment(frame.form) <- environment(formula)
1092        mf$formula <- fixed.form
1093        mf$x <- mf$model <- mf$y <- TRUE
1094        mf[[1]] <- as.name("glm")
1095        glm.fit <- eval(mf, parent.frame())
1096        x <- glm.fit$x
1097        y <- as.double(glm.fit$y)
1098        family <- glm.fit$family
1099        ## check for a linear mixed model
1100        lmm <- family$family == "gaussian" && family$link == "identity"
1101        if (lmm) { # linear mixed model
1102            method <- match.arg(method)
1103            if (method %in% c("PQL", "Laplace", "AGQ")) {
1104                warning(paste('Argument method = "', method,
1105                              '" is not meaningful for a linear mixed model.\n',
1106                              'Using method = "REML".\n', sep = ''))
1107                method <- "REML"
1108            }
1109        } else { # generalized linear mixed model
1110            if (missing(method)) method <- "PQL"
1111            else {
1112                method <- match.arg(method)
1113                if (method == "ML") method <- "PQL"
1114                if (method == "REML")
1115                    warning('Argument method = "REML" is not meaningful ',
1116                            'for a generalized linear mixed model.',
1117                            '\nUsing method = "PQL".\n')
1118            }
1119        }
1120    
1121        ## evaluate a model frame for fixed and random effects
1122        mf$formula <- frame.form
1123        mf$x <- mf$model <- mf$y <- mf$family <- NULL
1124        mf$drop.unused.levels <- TRUE
1125        mf[[1]] <- as.name("model.frame")
1126        frm <- eval(mf, parent.frame())
1127    
1128        ## grouping factors and model matrices for random effects
1129        bars <- findbars(formula[[3]])
1130        random <-
1131            lapply(bars, function(x)
1132                   list(model.matrix(eval(substitute(~ T, list(T = x[[2]]))),
1133                                     frm),
1134                        eval(substitute(as.factor(fac)[,drop = TRUE],
1135                                        list(fac = x[[3]])), frm)))
1136        names(random) <- unlist(lapply(bars, function(x) deparse(x[[3]])))
1137    
1138        ## order factor list by decreasing number of levels
1139        nlev <- sapply(random, function(x) length(levels(x[[2]])))
1140        if (any(diff(nlev) > 0)) {
1141            random <- random[rev(order(nlev))]
1142        }
1143    
1144        ## Create the model matrices and a mixed-effects representation (mer)
1145        mmats <- c(lapply(random, "[[", 1),
1146                   .fixed = list(cbind(glm.fit$x, .response = glm.fit$y)))
1147        mer <- .Call("mer2_create", lapply(random, "[[", 2),
1148                           mmats, method, PACKAGE = "Matrix")
1149        mer
1150    }

Legend:
Removed from v.939  
changed lines
  Added in v.940

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