massage tests and add template approach to formula parsing

Author: Steve Walker

DESCRIPTION

@@ -10,13 +10,11 @@ Description: The penalized least squares (PLS) and penalized iteratively
     in pure R. The purpose is to clarify how PLS and PIRLS work without having
     to read through C++ code, and as a sandbox for trying out modified versions
     of PLS and PIRLS.
+Depends:
+    Matrix
 Suggests:
     minqa,
     mlmRev
 Imports:
-    lme4,
-    Matrix
+    lme4
 License: GPL-2
-Collate:
-    'pirls.R'
-    'pls.R'

NAMESPACE

@@ -1,17 +1,26 @@
+export(Zsection)
 export(blockLambdat)
+export(mkLambdat)
+export(mkLind)
 export(mkRanefRepresentation)
+export(mkRanefStructures)
+export(mkTemplate)
+export(mkTemplates)
+export(mkTheta)
+export(mkZt)
+export(mkZtSection)
 export(pirls)
 export(pls)
 export(plsform)
-export(Zsection)
-importFrom(lme4,findbars)
-importFrom(lme4,nobars)
-importFrom(lme4,subbars)
-importFrom(Matrix,bdiag)
+export(rLmer)
 importFrom(Matrix,Cholesky)
 importFrom(Matrix,Diagonal)
+importFrom(Matrix,bdiag)
 importFrom(Matrix,rBind)
 importFrom(Matrix,sparse.model.matrix)
+importFrom(lme4,findbars)
+importFrom(lme4,nobars)
+importFrom(lme4,subbars)
 importMethodsFrom(Matrix,"%*%")
 importMethodsFrom(Matrix,crossprod)
 importMethodsFrom(Matrix,determinant)

R/pls.R

@@ -12,7 +12,7 @@ NULL
 ##' to read through C++ code, and as a sandbox for
 ##' trying out modifications to PLS.
 ##'
-##' @param X fixed effects design (model) matrix
+##' @param X fixed effects model matrix
 ##' @param y response
 ##' @param Zt transpose of the sparse model matrix for the random effects
 ##' @param Lambdat upper triangular sparse Cholesky factor of the
@@ -109,14 +109,14 @@ pls <- function(X,y,Zt,Lambdat,thfun,weights,
 ##'
 ##' @return a \code{list} with:
 ##' \itemize{
-##' \item \code{X} Fixed effects design (model) matrix
+##' \item \code{X} Fixed effects model matrix
 ##' \item \code{y} Observed response vector
 ##' \item \code{fr} Model frame
 ##' \item \code{call} Matched call
 ##' \item \code{REML} Logical indicating REML or not
 ##' \item \code{weights} Prior weights or \code{NULL}
 ##' \item \code{offset} Prior offset term or \code{NULL}
-##' \item \code{Zt} Transposed random effects design matrix
+##' \item \code{Zt} Transposed random effects model matrix
 ##' \item \code{Lambdat} Transposed relative covariance factor
 ##' \item \code{theta} Vector of covariance parameters
 ##' \item \code{lower} Vector of lower bounds for \code{theta}
@@ -142,8 +142,9 @@ plsform <- function(formula, data, REML=TRUE, weights, offset, sparseX = FALSE,
            y = rho$y,
            fr = fr, call = mc,
            REML = as.logical(REML)[1]),
-      if (is.null(wts <- model.weights(fr))) wts else list(weights=wts),
-      if (is.null(off <- model.offset(fr))) off else list(offset=off),
+      #if (is.null(wts <- model.weights(fr))) wts else list(weights=wts),
+      #if (is.null(off <- model.offset(fr))) off else list(offset=off),
+      list(weights = rho$weights), list(offset = rho$offset),
       mkRanefRepresentation(lapply(rr, function(t) as.factor(eval(t[[3]], fr))),
                 lapply(rr, function(t)
                        model.matrix(eval(substitute( ~ foo, list(foo = t[[2]]))), fr))))
@@ -155,8 +156,8 @@ initializeResp <- function(fr, REML, family){
     if(!inherits(family,"family")) family <- family()
     y <- model.response(fr)
 ### Why consider there here?  They are handled in plsform.
-    #offset <- model.offset(fr)
-    #weights <- model.weights(fr)
+    # offset <- model.offset(fr)
+    # weights <- model.weights(fr)
     n <- nrow(fr)
     etastart_update <- model.extract(fr, "etastart")
     if(length(dim(y)) == 1) {
@@ -172,15 +173,15 @@ initializeResp <- function(fr, REML, family){
     if (!is.null(REML)) rho$REML <- REML
     rho$etastart <- fr$etastart
     rho$mustart <- fr$mustart
-    #if (!is.null(offset)) {
-    #    if (length(offset) == 1L) offset <- rep.int(offset, n)
-    #    stopifnot(length(offset) == n)
-    #    rho$offset <- unname(offset)
-    #} else rho$offset <- rep.int(0, n)
-    #if (!is.null(weights)) {
-    #    stopifnot(length(weights) == n, all(weights >= 0))
-    #    rho$weights <- unname(weights)
-    #} else rho$weights <- rep.int(1, n)
+    if (!is.null(offset <- model.offset(fr))) {
+       if (length(offset) == 1L) offset <- rep.int(offset, n)
+       stopifnot(length(offset) == n)
+       rho$offset <- unname(offset)
+    } else rho$offset <- rep.int(0, n)
+    if (!is.null(weights <- model.weights(fr))) {
+       stopifnot(length(weights) == n, all(weights >= 0))
+       rho$weights <- unname(weights)
+    } else rho$weights <- rep.int(1, n)
     stopifnot(inherits(family, "family"))
     rho$nobs <- n
     eval(family$initialize, rho)
@@ -202,11 +203,11 @@ mkMFCall <- function(mc, form, nobars=FALSE) {
     mc
 }
 
-##' Create a section of a transposed random effects design matrix
+##' Create a section of a transposed random effects model matrix
 ##'
 ##' @param grp Grouping factor for a particular random effects term.
 ##' @param mm Dense model matrix for a particular random effects term.
-##' @return Section of a random effects design matrix corresponding to a
+##' @return Section of a random effects model matrix corresponding to a
 ##' particular term.
 ##' @export
 ##' @examples
@@ -342,6 +343,8 @@ blockLambdat <- function(nl, nc) {
 }
 
 
+
+
 ##' Make random effects representation
 ##'
 ##' Create all of the elements required to specify the random-effects
@@ -363,7 +366,7 @@ blockLambdat <- function(nl, nc) {
 ##' @return A \code{list} with:
 ##' \itemize{
 ##' \item \code{Lambdat} Transformed relative covariance factor
-##' \item \code{Zt} Transformed random effects design matrix
+##' \item \code{Zt} Transformed random effects model matrix
 ##' \item \code{theta} Vector of covariance parameters
 ##' \item \code{lower} Vector of lower bounds on the covariance parameters
 ##' \item \code{upper} Vector of upper bounds on the covariance parameters

R/templateApproach.R

@@ -1,9 +1,39 @@
+##' Random lmer-type model simulation
+##' @export
+rLmer <- function(grp, mmRE, mmFE) message("not yet written")
+
+
+##' Make random effects structures
+##' @export
+mkRanefStructures <- function(grp, mm){
+                                        # checking things are OK
+    if(!is.list(grp)) grp <- list(grp)
+    if(!is.list(mm)) mm <- list(mm)
+    grp <- lapply(grp, as.factor)
+    
+
+    nl <- sapply(grp, nlevels)   # number of grouping factor levels per term 
+    nc <- sapply(mm, ncol)       # number of model matrix columns per term
+    templates <- mkTemplates(nc) # templates for relative covariance factor
+    theta <- mkTheta(templates)
+
+    list(     Zt = mkZt(grp, mm),
+         Lambdat = mkLambdat(templates, nl),
+            Lind = mkLind(nl, nc),
+           theta = theta,
+           lower = ifelse(theta, 0, -Inf),  # lower and
+           upper = rep(Inf, length(theta))  # upper bounds on theta parameters
+         )
+}
+
+
 ##' Create a section of a transposed random effects model matrix
 ##'
 ##' @param grp Grouping factor for a particular random effects term.
 ##' @param mm Dense model matrix for a particular random effects term.
 ##' @return Section of a random effects design matrix corresponding to a
 ##' particular term.
+##' @export
 ##' @examples
 ##' ## consider a term (x | g) with:
 ##' ## number of observations, n = 6
@@ -23,16 +53,17 @@ mkZtSection <- function(grp,mm) {
 }
 
 ##' Make transposed random-effects model matrix
+##' @export
 mkZt <- function(grp,mm){
-    ZtSections <- mapply(mkZtSection, grp, mm)
-    rBind(ZtSections)
+    ZtSections <- mapply(mkZtSection, grp, mm, SIMPLIFY=FALSE)
+    do.call(rBind, ZtSections)
 }
 
-
 ##' Make a single template for a relative covariance factor
 ##'
 ##' @param nc Number of columns in a dense model matrix for a particular
 ##' random effects term
+##' @export
 ##' @examples
 ##' mkTemplate(5)
 mkTemplate <- function(nc){
@@ -48,9 +79,11 @@ mkTemplate <- function(nc){
 }
 
 ##' Make list of templates for relative covariance factor
+##' @export
 mkTemplates <- function(nc) lapply(nc, mkTemplate)
 
 ##' Make vector of indices giving the mapping from theta to Lambdat
+##' @export
 mkLind <- function(nl, nc){
                                         # number of thetas per term (i.e. number
                                         # of elements in the upper triangle of
@@ -69,17 +102,18 @@ mkLind <- function(nl, nc){
 }
 
 ##' Make initial relative covariance factor from list of templates
+##' @export
 mkLambdat <- function(templates, nl){
                                         # repeat templates once for each level
-    repTemplates <- rep(templates, nl)
+    templateList <- rep(templates, nl)
                                         # return Lambdat by putting blocks
                                         # together along the diagonal
-    .bdiag(repTemplates)
+    .bdiag(templateList)
 }    
 
 ##' Make initial theta from list of templates
+##' @export
 mkTheta <- function(templates){
     thetas <- lapply(templates, slot, "x")
     unlist(thetas)
 }
-

man/Zsection.Rd

@@ -1,8 +1,8 @@
 \name{Zsection}
 \alias{Zsection}
-\title{Create a section of a transposed random effects design matrix}
+\title{Create a section of a transposed random effects model matrix}
 \usage{
-  Zsection(grp, mm)
+Zsection(grp, mm)
 }
 \arguments{
   \item{grp}{Grouping factor for a particular random
@@ -12,11 +12,12 @@
   effects term.}
 }
 \value{
-  Dense section of a random effects design matrix
+Section of a random effects model matrix corresponding to a
+particular term.
 }
 \description{
-  Create a section of a transposed random effects design
-  matrix
+Create a section of a transposed random effects model
+matrix
 }
 \examples{
 ## consider a term (x | g) with:

man/blockLambdat.Rd

@@ -2,7 +2,7 @@
 \alias{blockLambdat}
 \title{Create digonal block on transposed relative covariance factor}
 \usage{
-  blockLambdat(nl, nc)
+blockLambdat(nl, nc)
 }
 \arguments{
   \item{nl}{Number of levels in a grouping factor for a
@@ -14,17 +14,16 @@
   the \code{mm} argument in \code{\link{Zsection}}).}
 }
 \value{
-  A \code{list} with: \itemize{ \item the block \item
-  ititial values of theta for the block \item lower bounds
-  on these initial theta values \item a function that
-  updates the block given the section of theta for this
-  block }
+A \code{list} with: \itemize{ \item the block \item ititial
+values of theta for the block \item lower bounds on these
+initial theta values \item a function that updates the
+block given the section of theta for this block }
 }
 \description{
-  Each random-effects term is represented by diagonal block
-  on the transposed relative covariance factor.
-  \code{blockLambdat} creates such a block, and returns
-  related information along with it.
+Each random-effects term is represented by diagonal block
+on the transposed relative covariance factor.
+\code{blockLambdat} creates such a block, and returns
+related information along with it.
 }
 \examples{
 (l <- blockLambdat(2, 3))

man/mkRanefRepresentation.Rd

@@ -2,7 +2,7 @@
 \alias{mkRanefRepresentation}
 \title{Make random effects representation}
 \usage{
-  mkRanefRepresentation(grps, mms)
+mkRanefRepresentation(grps, mms)
 }
 \arguments{
   \item{grps}{List of factor vectors of length n indicating
@@ -13,29 +13,28 @@
   corresponds to a random effects term.}
 }
 \value{
-  A \code{list} with: \itemize{ \item \code{Lambdat}
-  Transformed relative covariance factor \item \code{Zt}
-  Transformed random effects design matrix \item
-  \code{theta} Vector of covariance parameters \item
-  \code{lower} Vector of lower bounds on the covariance
-  parameters \item \code{upper} Vector of upper bounds on
-  the covariance parameters \item \code{thfun} A function
-  that maps \code{theta} onto the non-zero elements of
-  \code{Lambdat} }
+A \code{list} with: \itemize{ \item \code{Lambdat}
+Transformed relative covariance factor \item \code{Zt}
+Transformed random effects model matrix \item \code{theta}
+Vector of covariance parameters \item \code{lower} Vector
+of lower bounds on the covariance parameters \item
+\code{upper} Vector of upper bounds on the covariance
+parameters \item \code{thfun} A function that maps
+\code{theta} onto the non-zero elements of \code{Lambdat} }
 }
 \description{
-  Create all of the elements required to specify the
-  random-effects structure of a mixed effects model.
+Create all of the elements required to specify the
+random-effects structure of a mixed effects model.
 }
 \details{
-  The basic idea of this function is to call
-  \code{\link{Zsection}} and \code{\link{blockLambdat}}
-  once for each random effects term (ie. each list element
-  in \code{grps} and \code{mms}). The results of
-  \code{\link{Zsection}} for each term are \code{rBind}ed
-  together. The results of \code{\link{blockLambdat}} are
-  \code{bdiag}ed together, unless all terms have only a
-  single column ('predictor') in which case a diagonal
-  matrix is created directly.
+The basic idea of this function is to call
+\code{\link{Zsection}} and \code{\link{blockLambdat}} once
+for each random effects term (ie. each list element in
+\code{grps} and \code{mms}). The results of
+\code{\link{Zsection}} for each term are \code{rBind}ed
+together. The results of \code{\link{blockLambdat}} are
+\code{bdiag}ed together, unless all terms have only a
+single column ('predictor') in which case a diagonal matrix
+is created directly.
 }