SCM

SCM Repository

[matrix] View of /www/slides/2010-07-21-Gaithersburg/sparse-model-matrices.Rnw
ViewVC logotype

View of /www/slides/2010-07-21-Gaithersburg/sparse-model-matrices.Rnw

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2575 - (download) (annotate)
Wed Jul 21 09:18:06 2010 UTC (9 years, 3 months ago) by mmaechler
File size: 7581 byte(s)
update; notably IRLS
\SweaveOpts{engine=R, keep.source=TRUE}
\SweaveOpts{eps=FALSE, pdf=TRUE, width=9, height=6, strip.white=TRUE}
\setkeys{Gin}{width=\textwidth}

<<preliminaries,echo=FALSE,results=hide>>=
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:
<<model.matrix>>=
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>>=
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}
<<modelMat-ex>>=
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:


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