\documentclass[11pt]{article}

%\VignettePackage{doBy}
%\VignetteIndexEntry{population means}
%\VignetteIndexEntry{LSMEANS}
%\VignetteIndexEntry{contrasts}
%\VignetteIndexEntry{estimable functions}


\usepackage{a4wide,hyperref}
\usepackage[T1]{fontenc}
\usepackage{url,a4}
\usepackage{boxedminipage,color,xcolor}

\RequirePackage{color,fancyvrb,amsmath,amsfonts}
\DeclareMathOperator{\EE}{\mathbb{E}}

\usepackage{framed}
\usepackage{comment}
\definecolor{shadecolor}{gray}{0.91}
\def\pkg#1{\texttt{#1}}
\def\R{\texttt{R}}
\def\code#1{\texttt{#1}}
\def\popmeans{\code{popMeans()}}
\def\popmatrix{\code{popMatrix()}}
\def\linmeans{\code{linMeans()}}
\def\linmatrix{\code{linMatrix()}}
\def\esticon{\code{esticon()}}

<<echo=FALSE,print=FALSE>>=
require( doBy )
prettyVersion <- packageDescription("doBy")$Version
prettyDate <- format(Sys.Date())
@


\title{Population means (also called marginal means or LSMEANS),
  contrasts and estimable functions
  in the \texttt{doBy} package}
\author{S{\o}ren H{\o}jsgaard and Ulrich Halekoh}
\date{\pkg{doBy} version \Sexpr{prettyVersion} as of \Sexpr{prettyDate}}

\begin{document}
\renewenvironment{Schunk}{\begin{shaded}\small}{\end{shaded}}

\maketitle

\hrule
\tableofcontents

\parindent0pt
\parskip5pt

%\tableofcontents
\setkeys{Gin}{height=3in}
\SweaveOpts{keep.source=T}

\section{Introduction}
\label{sec:xxx}

This is a working document; please feel free to suggest improvements.

\section{A simulated dataset}
\label{sec:simulated-dataset}

Consider these data:
@
<<>>=
library(doBy)
dd <- expand.grid(A=factor(1:3),B=factor(1:3),C=factor(1:2))
dd$y <- rnorm(nrow(dd))
dd$x <- rnorm(nrow(dd))^2
dd$z <- rnorm(nrow(dd))
head(dd,10)
@ %def

Consider the additive model
\begin{equation}
  \label{eq:1}
  y_i = \beta_0 + \beta^1_{A(i)}+\beta^2_{B(i)} + \beta^3_{C(i)} + e_i
\end{equation}
where $e_i \sim N(0,\sigma^2)$. We fit this model:

@
<<>>=
mm <- lm(y~A+B+C, data=dd)
coef(mm)
@ %def

Notice that the parameters corresponding to the factor levels
\code{A1}, \code{B1} and \code{C2} are set to zero to ensure
identifiability of the remaining parameters.


\section{Linear functions of parameters, contrasts}
\label{sec:line-funct-param}

For a regression model with parameters $\beta=(\beta^1, \beta^2,\dots,
\beta^P)$ we shall refer to a weighted sum of the form
\begin{displaymath}
  \sum_j w_j \beta^j
\end{displaymath}
as a contrast. Notice that it is common in the litterature to require
that $sum_j w_j=0$ for the sum $  \sum_j w_j \beta^j$ to be called a
contrast but we do not follow this tradition here.

The effect of changing the factor $A$ from \code{A2} to \code{A3} can
be found as
@
<<>>=
w <- c(0,-1,1,0,0,0)
sum(coef(mm)*w)
@ %def

The \esticon\ function provides this estimate, the standard error
etc.\ as follows:
@
<<>>=
esticon(mm, w)
@ %def


\section{Population means}
\label{sec:xxx}

Population means (sometimes also called marginal means)
are in some sciences much used for reporting marginal effects (to be
described below). Population means are known as lsmeans in SAS
jargon. Population means is a special kind of contrasts as defined in
Section~\ref{sec:line-funct-param}.

The model (\ref{eq:1})
is a model for the conditional mean $\EE(y|A,B,C)$.  Sometimes one is
interested in quantities like $\EE(y|A)$. This quantity can not
formally be found unless $B$ and $C$ are random variables such that we
may find $\EE(y|A)$ by integration.

However, suppose that $A$ is a treatment of main interest, $B$ is a
blocking factor and $C$ represents days on which the experiment was
carried out. Then it is tempting to average $\EE(y|A,B,C)$ over $B$
and $C$ (average over block and day) and think of this average as
$\EE(y|A)$.

\subsection{A brute--force calculation}
\label{sec:xxx}

The population mean for $A=1$ is
\begin{equation}
  \label{eq:2}
  \beta^0 + \beta^1_{A1} + \frac{1}{3} (\beta^2_{B1}+\beta^2_{B2}+\beta^2_{B3})
  + \frac{1}{2}(\beta^3_{C1}+\beta^3_{C2})
\end{equation}

Recall that the
parameters corresponding to the factor levels
\code{A1}, \code{B1} and \code{C2} are set to zero to ensure
identifiability of the remaining parameters. Therefore we may also
write the population mean for $A=1$ as
\begin{equation}
  \label{eq:3}
  \beta^0 + \frac{1}{3} (\beta^2_{B2}+\beta^2_{B3})
  + \frac{1}{2}(\beta^3_{C2})
\end{equation}


This quantity can be estimated as:

@
<<>>=
w <- c(1, 0, 0, 1/3, 1/3, 1/2)
coef(mm)*w
sum(coef(mm)*w)
@ %def


We may find the population mean for all three levels of $A$ as
@
<<>>=
W <- matrix(c(1, 0, 0, 1/3, 1/3, 1/2,
              1, 1, 0, 1/3, 1/3, 1/2,
              1, 0, 1, 1/3, 1/3, 1/2),nr=3, byrow=TRUE)
W
W %*% coef(mm)
@ %def

Notice that the matrix W is based on that the first level of $A$ is
set as the reference level. If the reference level is changed then so
must $W$ be.

\subsection{Using \esticon}
\label{sec:xxx}

Given that one has specified $W$, the \esticon\ function in the
\code{doBy} package be used for the calculations above and the
function also provides standard errors, confidence
limits etc:
@
<<>>=
esticon(mm, W)
@ %def

\section{Using \popmatrix\  and \popmeans}
\label{sec:xxx}

Writing the matrix $W$ is somewhat tedious and hence error prone. In
addition, there is a potential risk of getting the wrong answer if the
the reference level of a factor has been changed.
The \popmatrix\ function provides an automated way of generating
such matrices.
The above \verb+W+ matrix is  constructed by
@
<<>>=
pma <- popMatrix(mm,effect='A')
summary(pma)
@ %def



The \popmeans\ function is simply a wrapper around first a call
to \popmatrix\ followed by a call to (by default) \esticon:
@
<<>>=
pme <- popMeans(mm, effect='A')
pme
@ %def

More details about how the matrix was constructed is provided by the
\code{summary()} function:
@
<<>>=
summary(pme)
@ %def


The \verb+effect+ argument requires  to calculate the population means
for each level of
$A$ aggregating across the levels of the other variables in the data.

Likewise we may do:
@
<<>>=
popMatrix(mm,effect=c('A','C'))
@ %def
This gives the matrix for calculating the estimate for each
combination of \code{A} and \code{C} when averaging over \code{B}.
Consequently
@
<<>>=
popMeans(mm)
@ %def
gives the ``total average''.


\subsection{Using the \code{at} argument}

We may be interested in finding the population means
at all levels of  $A$
but only at $C=1$. This is obtained by using the \code{at} argument:
@
<<>>=
popMatrix(mm,effect='A', at=list(C='1'))
@ %def
Notice here that average is only taken over $B$. Another way of
creating the population means
at  all levels of $(A,C)$ is therefore
<<>>=
popMatrix(mm,effect='A', at=list(C=c('1','2')))
@ %def


We may have several variables in the \code{at} argument:
@
<<>>=
popMatrix(mm,effect='A', at=list(C=c('1','2'), B='1'))
@ %def

\subsection{Ambiguous specification when using the \texttt{effect} and
  \texttt{at} arguments}

There is room for an ambiguous specification if a variable appears in
both the \code{effect} and the \code{at} argument, such as
@
<<>>=
popMatrix(mm,effect=c('A','C'), at=list(C='1'))
@ %def

This ambiguity is due to the fact that the \verb+effect+ argument asks
for the populations means at all levels of the variables but the
\verb+at+ chooses only specific levels.

This ambiguity is resolved as follows: Any variable in the \code{at}
argument is removed from the \code{effect} argument such as the
statement above is equivalent to
@
<<eval=F>>=
popMatrix(mm,effect='A', at=list(C='1'))
@ %def

\subsection{Using covariates}

Next consider the model where a covariate is included:
@
<<>>=
mm2 <- lm(y~A+B+C+C:x, data=dd)
coef(mm2)
@ %def

In this case we get
<<>>=
popMatrix(mm2,effect='A', at=list(C='1'))
@ %def

Above, $x$ has been replaced by its average and that is the general
rule for models including covariates. However we may use the \code{at}
argument to ask for calculation of the population mean at some
user-specified value of $x$, say 12:
<<>>=
popMatrix(mm2,effect='A', at=list(C='1',x=12))
@ %def


\subsection{Using transformed covariates}

Next consider the model where a  transformation of a covariate is included:
@
<<>>=
mm3 <- lm(y~A+B+C+C:log(x), data=dd)
coef(mm3)
@ %def

In this case we can not use \popmatrix\ (and hence
\popmeans\ directly.  Instead we have first to
generate a new variable, say \verb+log.x+, with
\verb+log.x+$=\log(x)$, in the data and then proceed as

<<>>=
dd <- transform(dd, log.x = log(x))
mm3 <- lm(y~A+B+C+C:log.x, data=dd)
popMatrix(mm3,effect='A', at=list(C='1'))
@ %def

\section{The \code{engine} argument of \popmeans}

The \popmatrix is a function to generate a linear tranformation matrix of the model
parameters with emphasis on  constructing such matrices for population
means.
\popmeans\ invokes by default the \esticon\ function on this
linear transformation matrix for calculating parameter estimates and
confidecne intervals.
A similar function to \esticon\ is the \verb+glht+ function of the \verb+multcomp+
 package.

 The \code{glht()} function
 can be chosen via the \verb+engine+ argument of \popmeans:
<<>>=
 library(multcomp)
g<-popMeans(mm,effect='A', at=list(C='1'),engine="glht")
g
@ %def

This allows to apply the methods available on the \verb+glht+ object like
<<>>=
summary(g,test=univariate())
confint(g,calpha=univariate_calpha())
@
which yield the same results as the \esticon\ function.

By default the functions will adjust the tests  and confidence intervals for multiplicity
<<>>=
summary(g)
confint(g)
@


\end{document}

















