---
title: "Group Means, Group Mean Centered Observations"
subtitle: "Multilevel recoding chores #1"
author:
 - name: Paul Johnson
   affiliation: Center for Research Methods and Data Analysis, University of Kansas
   email: pauljohn@ku.edu
 - name: "CRMDA Guide #43. Please visit [http://crmda.ku.edu/guides](http://crmda.ku.edu/guides) for updates."
 - name: "Tags: guides, R, multilevel models"

abstract: The most common recoding chore that comes with
    multilevel models is creating the group mean and individual
    deviations within group variables. This describes how that can be
    done with tools provided in the base of R, as well as in some
    addon packages.
checked_by: "First Last"
Note to Authors: please_dont_change_the_next 4 lines!
date: "`r format(Sys.time(), '%Y %B %d')`"
output:
  html_document:
    highlight: haddock
---

```{r setup, include=FALSE}
##This Invisible Chunk is required in all CRMDA documents
outdir <- paste0("tmpout")
if (!file.exists(outdir)) dir.create(outdir, recursive = TRUE)
knitr::opts_chunk$set(echo=TRUE, comment=NA, fig.path=paste0(outdir, "/p-"))
options(width = 70)
```

# Multilevel modelers need group means

Multilevel modelers might need group-level summary information for
various reasons. They often want to change a multilevel regression
from this

\begin{equation}
y_{ji} = \beta_0 + \beta_1 x_{ji} + \epsilon_{ji}
\label{eq:eq1}
\end{equation}

to

$$
y_{ji} = \beta_0 + \beta_1 \bar{x}_{j.} + \beta_2 \widetilde{x}_{ji} +
\epsilon_{ji}
\label{eq:eq2}
$$

where $\bar{x}_{j.}$ is the mean of $x_{ji}$ within group $j$ and
$\widetilde{x}_{ji} = x_{ji} -\bar{x}_{j.}$. The variable $\widetilde{x}_{ji}$ is
the "group mean-centered" value of $x_{ji}$.

Changing the specification from the first equation \ref{eq:eq1} to the
second equation \ref{eq:eq2} is a way of checking for the model's
specification.  If the effect of changing the within group mean is
markedly different from changing the value of the individual's score,
then the ecological fallacy might be in evidence.  This specification
is widely known as "Mundlak's suggestion" for exploring the effect of
a numeric predictor at the contextual and individual levels.

There is an extended discussion of this coding in Snijders and
Bosker(). It turns out that the preceeding model \ref{eq:eq2} is
statistically equivalent to this,

$$
y_{ji} = \beta_0 + \beta_3 \bar{x}_{j.} + \beta_4 x_{ji} +
\epsilon_{ji}
\label{eq:eq3}
$$

This model appears to be more simple because it includes the
individual variable, $x_{ji}$, as a predictor, rather than
$\widetilde{x}_{ji}$. However, this fit is statistically
equivalent. The fitted regression lines are identical and the
difference between the fits is a bit of an illusion. All of the
information contained in one is actually present in the other.

\begin{align}
y_{ji} & = \beta_0 + \beta_1 \bar{x}_{j.} + \beta_2 \widetilde{x}_{ji} + \epsilon_{ji}\\
       & = \beta_0 + \beta_1 \bar{x}_{j.} + \beta_2( x_{ji} -\bar{x}_{j.}) + \epsilon_{ji}\\
	   & = \beta_0 + (\beta_1-\beta_2)\bar{x}_{j.} + \beta_2(x_{ji}) + \epsilon_{ji}
\end{align}

Clearly, the "individual level effect" is the same in both models,
$\beta_2 = \beta_4$. And, equally clearly, the superficial difference we see
in the parameter estimates is easily understood as $\beta_3 = (\beta_1-\beta_2)$.

While these are equivalent models, estimating a model with
$\widetilde{x}_{ji}$, the individual deviations within groups, is
considered to be suggestive by some authors. It is a way of explicitly
estimationg a "big fish in a small pond" effect, as if to say that
individual variations about the group mean are somehow more
interesting that individual variations themselves. 

In Raudenbush and Bryk (), the variable $\widetilde{x}_{ji}$ is given
a random intercept. That is an entirely different specification, of
course.

In any case, the group means are necessary for the analysis and, for
some purposes, the "group mean centered" individual data is
thought-provoking.

## Example data

The High School and Beyond data set summarizes mathematical
and English language achievement for students. I have instructions
about downloading and importing that data on this web page:

http://pj.freefaculty.org/guides/stat/DataSets/HSB/00-README.txt

The following will download a copy of the data if one is not
already available in the working directory.

```{r}
if (file.exists("hsb.rds")){
    hsb <- readRDS("hsb.rds")
} else {
    library(foreign)
    hsb.url <- "http://www.stata-press.com/data/mlmus3/hsb.dta"
    hsb <- read.dta(hsb.url)
    saveRDS(hsb, "hsb.rds")
}
```

```{r}
hsb$schoolidf <- factor(hsb$schoolid)
```


## Stata offers egen

If a variable is named ```ses``` and we want to calculate
the mean for each school, we run 

```
egen ses_mn = mean(ses), by(schoolid)
```

The deviation within group is easy to calculate after that.

```
gen ses_dev = ses - ses_mn 
```

In Stata, one is accustomed to the idea that the data frame will be
altered by the addition of a new variable.

In the Rabe-Hesketh and Skrondal book, a method is describe to use a
for loop to iterate through many variables and create group mean
columns for them.

Here we can practice by getting the group means for the variables ses
and female.

```
foreach var of varlist ses female {
    egen `var'_mn = mean(`var'), by(schoolid)
}

foreach var of varlist ses  female {
    gen `var'_dev = `var' - `var'_mn
}
```


## How to get the same work done with R

R offers many functions in the base that can produce the group means.
The more difficult challenge is to "duplicate" the means back 
onto the data frame so that we can put them to use.

### Base functions

```aggregate``` and ```by``` are convenience functions in the R base that are intended
to simplify the use of ```tapply```.


```{r}
ses_mn1 <- aggregate(hsb[, "ses", drop=FALSE], by = list(schoolid = hsb$schoolidf), mean, na.rm = TRUE)
head(ses_mn1)
str(ses_mn1)
class(ses_mn1)
```

```{r}
ses_mn2 <- by(hsb$ses, hsb$schoolidf, mean, na.rm = TRUE)
head(ses_mn2)
str(ses_mn2)
class(ses_mn2)
```

```{r}
ses_mn3 <- tapply(hsb$ses, hsb$schoolidf, mean, na.rm = TRUE)
head(ses_mn3)
str(ses_mn3)
class(ses_mn3)
```

Using any of these methods, we are able to retrieve the means of the
groups.

The challenge is to put those means back onto the data frame. It is
easy to get this wrong because the results from tapply do not come
back in the expected order. Here is a toy example to demonstrate that.

```{r}
toy <- data.frame(x1 = rnorm(13),
                  x2 = c(rep("Z", 3), rep("P", 2), rep("A", 3), rep("C", 5)),
                  stringsAsFactors = FALSE)
toy
x1_mn <- tapply(toy$x1, toy$x2, mean)
x1_mn
```

Note that the ordering from top to bottom in the data frame has groups
Z, P, A and C, the ordering of tapply is different.  tapply creates a
factor variable, which, by default, uses the available values in 
alphabetical order as the levels.  Hence, a decision to manufacture
a vector by the obvious sort of trick will result in a ghastly error.

<!-- ```{r} -->
<!-- x1_mn <- aggregate(toy, by = list(toy$x2), function(d){ -->
<!--     mean(d[ , "x1"]) -->
<!--     ##rep(mean(d[ , "x1"], NROW(d))) -->
<!-- }) -->
<!-- ``` -->

To put the data "back onto" the data frame requires the careful use of
index subscripts or the merge function. 

```{r}
x1_mn <- aggregate(toy[ , "x1", drop = FALSE], by = list("x2" = toy$x2), mean, na.rm = TRUE)
x1_mn
toy2 <- merge(toy, x1_mn, by = "x2", suffix = c("", "_mn"), sort = FALSE)
toy2
```

#### You usually want more than one variable {.bs-callout .bs-callout-red}

Variable "mmses" is the mean of ses, within schools, which we verify here
Create a vector of variable names for which we need school level means

```{r}
vars <- c("ses", "female")
ses_mn <- aggregate(hsb[ , vars],
                          by = list("schoolidf" = hsb$schoolidf),  mean, na.rm = TRUE)

hsb2 <- merge(hsb, ses_mn, by = "schoolidf", suffix = c("", "_mn"))
for(i in vars){
    hsb2[ , paste0(i, "_dev")] <- hsb2[ , i] - hsb2[ , paste0(i, "_mn")]
}
```


#### An R function like Stata's egen {.bs-callout .bs-callout-green}

```{r}

##' Generate group summaries and individual deviations within groups
##'
##' Similar to Stata egen, except more versatile and fun! Will
##' create a new data frame with 2 columns. Rows match
##' the rows of the original data frame. 
##'
##' Now I return only the new columns
##' @param dframe a data frame
##' @param x Variable names or a vector of variable names
##' @param by A grouping variable name or a vector of grouping names
##' @param FUN Defaults to the mean, have not tested alternatives
##' @param suffix The suffixes to be added to column 1 and column 2
##' @return new data frame with "x_mn" and "x_dev" added as variables
##' @author Paul Johnson
##' @examples
##' Suppose you get the MLMUS hsb data frame, somehow
##' xx1 <- gmd(hsb, "ses", "schoolidf")
##' xx2 <- gmd(hsb, c("ses", "female"), "schoolidf")
##' xx3 <- gmd(hsb, c("ses", "female"), c("schoolidf", "sector"))
##' xx4 <- gmd(hsb, c("ses", "female"),
##'                    c("schoolidf"), FUN = median)
gmd <- function(dframe, x, by, FUN = mean, suffix = c("_mn", "_dev")){
    xmean <- aggregate(dframe[ , x, drop = FALSE],
                       dframe[ , by, drop = FALSE], FUN,
                       na.rm = TRUE)
    df2 <- merge(dframe, xmean, by = by, suffix = c("", suffix[1]))
    for(i in x){
        df2[ , paste0(i, suffix[2])] <- df2[ , i] - df2[ , paste0(i, suffix[1])]
    }
    df2[ , colnames(df2)[!colnames(df2) %in% colnames(dframe)]]
}


xx1 <- gmd(hsb, "ses", "schoolidf")
head(xx1)
## Could combine with hsb
## hsb <- cbind(hsb, xx1)
## head(hsb[ , c("schoolidf", "ses", "ses_mn", "ses_dev")])
xx2 <- gmd(hsb, c("ses", "female"), "schoolidf")
head(xx2)
hsb <- cbind(hsb, xx2)
head(hsb[ , c("schoolidf", "ses", "ses_mn", "ses_dev", "female", "female_mn", "female_dev")])
```

####  Why do I bother with base functions?{.bs-callout .bs-callout-red}

This function, which uses only R base tools, is not
as fast as possible. The \emph{hsb} data set is is only 7185 rows with 
26 columns, not big data by any sense of the imagination.

If we could simply insert a column into a data.frame without causing a
full memory re-write of the data.frame within R, it would be faster.
The fear of excess data copying is what fuels innovations during the
past 5-10 years, such as the packages ```data.table``` ```plyr``` and
```dplyr```.

All of these packages sit on top of the famous R package "Rcpp", which
is rightly praised as a revolutionary advance in R for high
performance problems.  Unfortunately, it is also one of the sources of
periodic project failures because certain data sets expose bugs (which
are fixed eventually, usually quickly).

Simply put, relying on more contributed packages can cause project
failures.  Recently, we have had 2 projects come to a halt because
of little bugs in packages that were relied upon by other packages
that we were using.  These problems get fixed, eventually, but why
complicate the workflow unnecessarily.

#### If you are willing to use data.table {.bs-callout .bs-callout-blue}

When data sets are truly immense, the ```data.table``` package offers
the only tools I've found to quickly sort, merge, and
calculate-by-subgroup. ```data.table``` introduces an entirely new
style of R programming, touted by its advocates as more "elegant" and
"understandable". The CRMDA has a ```data.table``` guide available in
our collection.

Here is an example of data.table in action.


```{r}
library(data.table)
hsbdt <- as.data.table(hsb)
setkey(hsbdt, schoolidf)
hsbdt[ , ses_mn2 := mean(ses, na.rm = TRUE), by = schoolidf]
hsbdt[ , ses_dev2 := ses - ses_mn2]
hsbnew <- as.data.frame(hsbdt)
head(hsbnew[ , c("schoolidf", "ses", "ses_mn2", "ses_dev2")])
```

The beauty of the data.table is that memory is used efficiently. The
new columns, ses_mn2 and ses_dev2, are added without re-copying the
rest of the data.  In fact, we could snug these 2 column creations
into one command, thus reducing the memory allocation to a single step.

 <!--TODO: write out example  -->

```{r, include=F}
detach(package:data.table)
```

#### The peril of transform()

The base R offers a function called ```transform``` which seems, at
least in theory, to offer an avenue to efficiently append the column
results onto the original data frame.  This one is actually not more
efficient, it is simply different syntax.  It seems to some to offer
the best of both worlds.

```{r}
xxx <- transform(hsb, ses_mn = xx1$ses_mn, ses_dev = xx1$ses_dev)
```

What's the problem? Well, it is not faster, and it is dangerous. There
is a stern warning about the transform function.  It is among the most
clear cautions in all of the R documentation.  Run "?transform", look
for "Warning". The help page says

> "This is a convenience function
> intended for use interactively.  For programming it is better to use
> the standard subsetting arithmetic functions, and in particular the
> non-standard evaluation of argument ‘transform’ can have unanticipated
> consequences."

Speaking sarcastically, except for that warning, the transform function
looks like a good way to go at this.

#### Enter mutate

Hadley Wickham has offered a number of R packages and new idioms. Like
```data.table```, his ````plyr``` package, and its successor
```dplyr```, are intended to make R code both more "elegant" (in the
eyes of the author, at least), faster, and less error prone.  The ```plyr```
package implements an approach that is described as the "split, apply,
combine" paradigm.  It is summarized in this article:
 
Wickham, H. (2011). The Split-Apply-Combine Strategy for Data
Analysis. Journal of Statistical Software, 40(1),
1 - 29. doi:http://dx.doi.org/10.18637/jss.v040.i01

(https://www.jstatsoft.org/article/view/v040i01/v40i01.pdf)

The plyr package supplies a function called ```mutate```, a
replacement for transform.

```{r}
library(plyr)
hsb2 <- ddply(hsb, "schoolidf", mutate, ses_mn = mean(ses, na.rm = TRUE), ses_dev = ses - ses_mn)
head(hsb2)
```

```{r, include=F}
detach(package:plyr)
```

If one enjoys the "pipe" notation that is supported by the newer dplyr
(_which I do not_), the same is written as

```{r}
library(dplyr)
hsb3 <- hsb %>% group_by(schoolidf) %>% mutate(ses_mn = mean(ses, na.rm=TRUE), ses_dev = ses - ses_mn)
hsb3[ , c('schoolidf', 'ses', 'ses_mn', 'ses_dev')]
```
<!-- ##hsb %>% group_by(schoolidf) %>% summarize(ses_mn = mean(ses, na.rm = TRUE)) %>% right_join(hsb, by = "schoolidf") -->

```{r}
table(hsb3$ses_mn == hsb$ses_mn)
```

[//]: (All guides must have this as the final stanza)

```{r sessionInfo, echo = FALSE}
sessionInfo()
```

Available under
[Created Commons license 3.0 <img src="http://crmda.dept.ku.edu/images/cc-3.0.png" alt="CC BY"
style="width: 75px;height: 20px;"/>](http://creativecommons.org/licenses/by/3.0/)

