\SweaveOpts{engine=R, keep.source=TRUE} \SweaveOpts{eps=FALSE, pdf=TRUE, width=9, height=6, strip.white=TRUE} \setkeys{Gin}{width=\textwidth} <>= options(width=75) library(Matrix) @ \begin{frame}[fragile]\frametitle{Sparse \alert{Model} Matrices} New model matrix classes, %The \pkg{Matrix} package now has model matrix classes, generalizing \Rp's standard \Rfun{model.matrix}:% function: <>= str(dd <- data.frame(a = gl(3,4), b = gl(4,1,12)))# balanced 2-way model.matrix(~ 0+ a + b, dd) @ \end{frame} \begin{frame}[fragile]\frametitle{Sparse \alert{Model} Matrices} The model matrix above \begin{itemize} \item \dots\dots has many zeros, and \item \emph{ratio} ((zeros) : (non-zeros)) increases dramatically with many-level factors \item even more zeros for factor \emph{interactions}: \end{itemize} \smallskip \begin{footnotesize} <>= model.matrix(~ 0+ a * b, dd) @ \end{footnotesize} \end{frame} \begin{frame}[fragile]\frametitle{Sparse \alert{Model} Matrices in 'Matrix'} \begin{itemize} \item These matrices can become very large: Both many rows (large $n$), \emph{and} many columns, large $p$. \item Eg., in Linear Mixed Effects Models, \begin{equation*} \Ew{\bc Y|\bc B=\bm b} = \gX \bm\beta+\rZ\bm b, \end{equation*} \begin{itemize} \item the $\rZ$ matrix is often large and very sparse, and in \pkg{lme4} has always been stored as \code{"sparseMatrix"} (\code{"dgCMatrix"}, specifically). \item Sometimes, $\gX$, (fixed effect matrix) is large, too. \nlQ $\to$ optionally also \code{"sparseMatrix"} in \pkg{lme4}\footnote{the development version of \pkg{lme4}, currently called \pkg{lme4a}.}. \end{itemize} \item We've extended \Rfun{model.\Ul{m}atrix} to \Rfun{model.\Ul{M}atrix} in package \pkg{Matrix} with optional argument \code{sparse = TRUE}. \end{itemize} \end{frame} \begin{frame}[fragile] \frametitle{Sparse Model Matrix \alert{Classes} in 'Matrix'} \begin{footnotesize} \begin{Schunk} \begin{Sinput} setClass("modelMatrix", representation(assign = "integer", contrasts = "list", "VIRTUAL"), contains = "Matrix", validity = function(object) { ........... }) setClass("sparseModelMatrix", representation("VIRTUAL"), contains = c("CsparseMatrix", "modelMatrix")) setClass("denseModelMatrix", representation("VIRTUAL"), contains = c("denseMatrix", "modelMatrix")) ## The actual'' *ModelMatrix classes: setClass("dsparseModelMatrix", contains = c("dgCMatrix", "sparseModelMatrix")) setClass("ddenseModelMatrix", contains = c("dgeMatrix", "ddenseMatrix", "denseModelMatrix")) \end{Sinput} \end{Schunk} \end{footnotesize} \medskip % (adding \code{"ddenseMatrix"} above does \emph{not} influence slots, but yields % consistent superclass ordering.) (\code{"ddenseMatrix"}: \emph{not} for slots, but consistent superclass ordering) \end{frame} \begin{frame}[fragile]\frametitle{model.\alert{M}atrix(*, sparse=TRUE)} Constructing \alert{sparse} models matrices (\pkg{Matrix} package): \begin{footnotesize} <>= model.Matrix(~ 0+ a * b, dd, sparse=TRUE) @ \end{footnotesize} % the above is cut off, i.e. not visible here % identical syntax, just \code{model.\alert{M}atrix(..)}. \end{frame} \section[modelMatrix \Also "glpModel"s]{% modelMatrix \Also General Linear Prediction Models} \begin{frame}[fragile] \frametitle{"modelMatrix" \Also General Linear Prediction Models} Idea: Very \emph{general} setup for \begin{block}{Statistical models based on linear predictors} Class \code{"glpModel"} := General Linear Prediction Models \end{block} \begin{Schunk} \begin{Sinput} setClass("Model", representation(call = "call", fitProps = "list", "VIRTUAL")) setClass("glpModel", representation(resp = "respModule", pred = "predModule"), contains = "Model") \end{Sinput} \end{Schunk} \medskip Two main ingredients: \begin{enumerate} \item Response module \code{"respModule"} \item (Linear) Prediction module \code{"predModule"} \end{enumerate} \end{frame} \begin{frame}[fragile]\frametitle{(1) Response Module} \code{"respModule"}: Response modules for models with a linear predictor, which can include linear models (\code{lm}), generalized linear models (\code{glm}), nonlinear models (\code{nls}) and generalized nonlinear models (\code{nglm}): \begin{footnotesize}\begin{Schunk} \begin{Sinput} setClass("respModule", representation(mu = "numeric", # of length n offset = "numeric", # of length n * k sqrtXwt = "matrix", # of dim(.) == dim(X) == (n, k) sqrtrwt = "numeric", # sqrt(residual weights) weights = "numeric", # prior weights wtres = "numeric", y = "numeric"), validity = function(object) { ....... }) setClass("glmRespMod", representation(family = "family", eta = "numeric", n = "numeric"), # for evaluation of the aic contains = "respModule", validity = function(object) { ..... }) setClass("nlsRespMod", representation(nlenv = "environment", ......), ............) setClass("nglmRespMod", contains = c("glmRespMod", "nlsRespMod")) \end{Sinput} \end{Schunk}\end{footnotesize} \end{frame} \begin{frame}[fragile]\frametitle{(2) Prediction Module} %\item \code{"predModule"}: Linear predictor module consists of \begin{itemize} \item the model matrix \code{X}, \item the coefficient vector \code{coef}, \item a triangular factor of the weighted model matrix \code{fac}, \item (\code{Vtr} $= \bm {V\trans r}$, where $\bm r =$ residuals (typically) \end{itemize} currently in \textcolor{Orange}{dense} and \textcolor{Mulberry}{sparse} flavor: \medskip \begin{footnotesize}\begin{Schunk} \begin{Sinput} setClass("predModule", representation(X = "modelMatrix", coef = "numeric", Vtr = "numeric", fac = "CholeskyFactorization", "VIRTUAL")) ## the sub classes specify more specific classes for the two non-trivial slots: setClass("dPredModule", contains = "predModule", representation(X = "ddenseModelMatrix", fac = "Cholesky")) setClass("sPredModule", contains = "predModule", representation(X = "dsparseModelMatrix", fac = "CHMfactor")) \end{Sinput} \end{Schunk}\end{footnotesize} \end{frame} \begin{frame} \frametitle{Fitting all glpModel''s with One IRLS algorithm} Fitting via IRLS (Iteratively Reweighted Least Squares), where the prediction and response module parts each update themselves''. \bigskip These 3 Steps are iterated till convergence:\\ \begin{enumerate} \item prediction module (\textsc{pm}) only passes \code{X \%*\% coef}$= \bm X \beta$ to the response module (\textsc{rm}) \item from that, the \textsc{rm} \begin{itemize} \item updates its $\bm \mu$, \item then its weighted residuals and X weights'' \end{itemize} \item these two are in turn passed to \textsc{pm} which \begin{itemize} \item reweights itself and \item \code{solve()}s for $\Delta\bm\beta$, the \emph{increment} of $\bm\beta$. \end{itemize} \end{enumerate} Convergence only if Bates-Watts orthgonality criterion is fulfilled. \end{frame} %%% Local Variables: %%% TeX-command-default: "LaTeX PDF" %%% TeX-master: "MaechlerBates.tex" %%% End: