Estimation and prediction in sparse and unbalanced tables

TL;DR

This paper introduces a hierarchical ANOVA approach for estimating and predicting in multi-way tables with sparse and unbalanced data, effectively borrowing strength from lower-order models to improve estimates.

Contribution

The paper proposes a novel hierarchical ANOVA method that adaptively shrinks estimates towards lower-order models in sparse multi-way tables, enhancing estimation accuracy.

Findings

Hierarchical ANOVA outperforms unshrinked methods in simulations.

The approach effectively handles sparsity and unbalanced data.

Algorithm is simple and computationally efficient.

Abstract

We consider the problem where we have a multi-way table of means, indexed by several factors, where each factor can have a large number of levels. The entry in each cell is the mean of some response, averaged over the observations falling into that cell. Some cells may be very sparsely populated, and in extreme cases, not populated at all. We might still like to estimate an expected response in such cells. We propose here a novel hierarchical ANOVA (HANOVA) representation for such data. Sparse cells will lean more on the lower-order interaction model for the data. These in turn could have components that are poorly represented in the data, in which case they rely on yet lower-order models. Our approach leads to a simple hierarchical algorithm, requiring repeated calculations of sub-table means of modified counts. The algorithm has shown superiority over the unshrinked methods in both simulations and real data sets.

Figures

Equations

$$L(\mu )=\sum _{i,j,k}{n}_{ijk}(y_{ijk}-{\mu }_{ijk}{)}^2+\lambda \sum _{i,j,k}({\mu }_{ijk}-{\gamma }_{ijk}{)}^2,$$

$$E({\mu }_{ijk}|\mathbf{y})=\frac{n_{ijk}{y}_{ijk}}{n_{ijk}+\lambda }+\frac{\lambda {\gamma }_{ijk}}{n_{ijk}+\lambda },$$

$$L(\mu )=\sum _{i,j}(y_{ij}-{\mu }_{ij}{)}^2+\lambda \sum _{i,j}({\mu }_{ij}-{\mu }_{i\cdot }-{\mu }_{\cdot j}+{\mu }_{\cdot \cdot }{)}^2,$$

$$L(\mathbf{\mu })=(\mathbf{y}-\mathbf{\mu }{)}^T(\mathbf{y}-\mathbf{\mu })+\lambda {\mathbf{\mu }}^T(\mathbf{I}-{\mathbf{P}}_A)\mathbf{\mu },$$

$$\hat{\mathbf{\mu }}=(\mathbf{I}+\lambda (\mathbf{I}-{\mathbf{P}}_A){)}^{-1}\mathbf{y}.$$

$$\hat{\mathbf{\mu }}=\frac{\mathbf{y}}{1+\lambda }+\frac{\lambda \tilde{\mathbf{y}}}{1+\lambda },$$

$$\underset{\alpha }{\min }\sum _{i,j,k}{n}_{ijk}(y_{ijk}-({\alpha }_i^1+{\alpha }_j^2+{\alpha }_k^3){)}^2+{\lambda }_1\cdot \sum _{n_{ijk}>0}({\alpha }_i^1+{\alpha }_j^2+{\alpha }_k^3-{\mu }_{ijk}^{(0)}{)}^2.$$

$$\underset{\psi }{\min }$$

$$+{\lambda }_2\cdot \sum _{n_{ijk}>0}({\psi }_{ij}^{12}+{\psi }_{jk}^{23}+{\psi }_{ik}^{13}-{\mu }_{ijk}^{(1)}{)}^2.$$

$$\underset{\mu }{\min }\sum _{i,j,k}{n}_{ijk}(y_{ijk}-{\mu }_{ijk}{)}^2+{\lambda }_3\sum _{n_{ijk}>0}({\mu }_{ijk}-{\mu }_{ijk}^{(2)}{)}^2.$$

$${\mu }_I^{(k)}=f_k(I;{\mathbf{\beta }}^{(k)})=\sum _{J\subset [m],|J|=k}{\beta }_{I_J}^J$$

$$y_{ij}$$

$${\theta }_{ij}$$

$$\alpha ,\beta \sim \mathrm{N}(0,V)$$

$$y_i=\sum _{m=0}^M{\beta }_{j_i^m}^{(m)},$$

$${\beta }_j^{(m)}\sim \mathrm{N}(0,{\sigma }_m^2),\forall j=1,\dots ,J_m,m=1,\dots ,M$$

$$y_{ijkl}=\mu +a_i+b_j+c_k+(ab{)}_{ij}+(ac{)}_{ik}+(bc{)}_{jk}+(abc{)}_{ijk}+{ϵ}_{ijkl}$$

$$-\frac{1}{2}[\Vert \mathbf{y}-{\mathbf{\mu }}^{(k)}{\Vert }^2+{\lambda }_k\Vert {\mathbf{\mu }}^{(k)}-{\mathbf{\mu }}^{(k-1)}{\Vert }^2]$$

$$\mathbf{y}|\mathbf{\mu }$$

$$\mathbf{\mu }|{\mathbf{\beta }}^{(m)}$$

$${\mathbf{\beta }}^{(m)}|{\mathbf{\beta }}^{(m-1)}$$

$$⋮$$

$${\mathbf{\beta }}^{(k)}|{\mathbf{\beta }}^{(k-1)}$$

$$⋮$$

$${\mathbf{\beta }}^{(1)}|{\mathbf{\beta }}^{(0)}$$

$${\lambda }_k=\frac{{\sigma }^2+{\sigma }_m^2+\dots +{\sigma }_k^2}{{\sigma }_{k-1}^2}$$

$$\mathbf{y}|\mathbf{\mu }$$

$$\mathbf{\mu }|{\mathbf{\beta }}^{(m)}$$

$${\mathbf{\beta }}^{(m)}$$

$$\mathbf{\Sigma }$$

$${\tau }_i^2$$

$$\mathrm{E}[{\mathbf{y}}^T{\mathbf{U}}_k{\mathbf{U}}_k^T\mathbf{y}]$$

$$=\mathrm{tr}[{\mathbf{U}}_k{\mathbf{U}}_k^T({\sigma }^2{\mathbf{I}}_n+{\sigma }_m^2{\mathbf{I}}_n+{\mathbf{U}}_m\mathbf{\Sigma }{\mathbf{U}}_m^T)]$$

$$=d_k({\sigma }^2+{\sigma }_m^2)+\sum _{j=0}^{k-1}(d_{j+1}-d_j){\tau }_j^2,$$

$$\forall k=1,\cdots ,m$$

$$\mathrm{E}[{\mathbf{y}}^T\mathbf{U}{\mathbf{U}}^T\mathbf{y}]$$

$$=\mathrm{tr}[({\sigma }^2+{\sigma }_m^2){\mathbf{I}}_n+{\mathbf{U}}_m\mathbf{\Sigma }{\mathbf{U}}_m^T]$$

$$=n({\sigma }^2+{\sigma }_m^2)+\sum _{j=0}^{m-1}{d}_{j+1}{\sigma }_j^2$$

$$\mathrm{E}[{\mathbf{y}}^T{\mathbf{U}}_{k+1}{\mathbf{U}}_{k+1}^T\mathbf{y}-{\mathbf{y}}^T{\mathbf{U}}_k{\mathbf{U}}_k^T\mathbf{y}]=(d_{k+1}-d_k){\tau }_k^2,$$

Full text

Qingyuan Zhao was a Ph.D. candidate at Stanford University when this work completed. Trevor Hastie was partially supported by grant DMS-1407548 from the National Science Foundation, grant 5R01 EB 001988-21 from the National Institutes of Health, and a Google research grant.

1. Introduction

Prediction with factorial features has been studied for a long time in statistics. Recently there are many arising applications of this kind, but with much larger data size than before. For instance, consider data on restaurant ratings, as provided by Zagat, Yelp or other such services. In addition, the restaurants can usually be classified according to various factors, such as zip-code or geographical region (at least 10000 levels), type of cuisine (Italian, French, etc, potentially dozens of levels), and price category (often 1-5 stars). One important task is to estimate the average rating for a particular kind of restaurants defined by these factors, which can be used to answer questions like “What is the most popular cuisine in San Francisco?” or “Do expensive restaurants usually get better ratings?”.

Such data sets and related questions have been studied thoroughly by statisticians from the first formal proposal of analysis of variance (ANOVA) in Fisher (1918). In most of the analysis of variance literature, the focus is how to estimate the relative importance of different variance components and how to test for their statistical significance. The problem studied in this paper differs from classical ANOVA in the following aspects:

• •

The focus is estimating and predicting cell means, whereas the main objective of ANOVA is to find important factors or interactions and test for their significance.

• •

The observed data table is sparse and highly unbalanced. In many applications, we may only have a small proportion of the cells observed with quite different weights. For example, one zip code may have no Ethiopian restaurant, only one unpopular Japanese restaurant with a few dozen ratings, but many American restaurants with thousands of ratings in total.

• •

The size of the data is usually very large and we also have a large number of main effects and interactions to estimate, so it is very important to find an computationally efficient algorithm while maintaining good statistical properties.

In the rest of the paper, we first illustrate our method with the restaurant ratings example in Section 2, then describe the statistical problem we are trying to solve and some previous literature in Section 3. In Section 4.1 we prove some theoretical properties of our algorithm when the observed data is balanced. This gives us an empirical choice of penalty parameter, as described in Section 4.2. The rest of the paper is devoted to algorithm implementation (Section 5) and simulations and some real data results (Section 6).

2. Estimation in Three-way Tables

Our proposed method is general can be applied to any number of factors, but this would require a somewhat technical representation. In this section we will restrict ourselves to three factors ($F_{1}$, $F_{2}$ and $F_{3}$) to demonstrate the key elements in the algorithm.

Suppose these factors have $I$, $J$, and $K$ levels respectively. The restaurant ratings are a collection of $n_{ijk}$ measurements $y_{ijk\ell},\;\ell=1,\ldots,n_{ijk}$ at each cell $(i,j,k)$, and we are interested in predicting the cell mean $\mu_{i,j,k}=E(Y|F_{1}=i,F_{2}=j,F_{3}=k)$. We summarize our data by the cell means $\bar{y}_{ijk}=y_{ijk\cdot}/n_{ijk}$ and a weight $n_{ijk}$. For simplicity, we drop the bar, and refer to our aggregate data as $y_{ijk}$.

For some cells $n_{ijk}$ could be very small, even zero, but we would still like to have a reasonable estimate. In this case we would like to shrink our estimate towards a more stable number that borrows strength from similar cells. This suggests a Bayesian mixed-effects framework Diggle et al. (1994), where we assume the $\mu_{ijk}$ are random, say Gaussian, with distribution $\mu_{ijk}\sim N(\gamma_{ijk},\sigma_{\mu}^{2})$. For the moment assume $\gamma_{ijk}$ is known. If we assume the original measurements $y_{ijk\ell}|\mu_{ijk}\sim N(\mu_{ijk},\sigma^{2})$ and are all independent, then the negative log-posterior likelihood for $\mu_{ijk}$ given data is proportional to

[TABLE]

where $\lambda=\sigma^{2}/\sigma^{2}_{\mu}$. The posterior mode is simple to characterize:

[TABLE]

a simple weighted average of the observed mean $y_{ijk}$ and $\gamma_{ijk}$, with more emphasis on the former when $n_{ijk}$ is large. Now $\gamma_{ijk}$ is not known, and so we can represent it by a parametric model $\gamma_{ijk}(\theta)$. Empirical Bayes amounts to estimating $\theta$ by maximizing the marginal likelihood (integrating out $\mu$), and then using $\gamma_{ijk}(\hat{\theta})$ as the target of shrinkage in (2.2).

Barry (1990), in a balanced two-way layout, proposed a simple main-effects parametrization $\gamma_{ij}=\alpha^{0}+\alpha^{1}_{i}+\alpha^{2}_{j},$ with $\sum_{i}\alpha^{1}_{i}=0$, and $\sum_{j}\alpha^{2}_{j}=0$. One can show in this case that the empirical Bayes procedure amounts to maximizing the following likelihood

[TABLE]

where $\mu_{i\cdot}=\sum_{i=1}^{I}\mu_{ij}/I$ etc. Now one can show two things:

• •

We can write (2.3) in vector notation as

[TABLE]

where $\mathbf{P}_{A}$ is the main-effects ANOVA projection operator. This has solution

[TABLE]

• •

This solution can be shown to have the much simpler form

[TABLE]

where $\tilde{\mathbf{y}}=\mathbf{P}_{A}\mathbf{y}$, which is the main-effects ANOVA fit. This is of the same form as (2.2), with $n_{ij}=1$ and $\gamma_{ij}=\tilde{y}_{ij}$

This nice simplification disappears if we include weights in (2.3) or (2.4), although we do get a closed form expression for the solution along the lines of (2.6) (but involving weights). What we like about (2.6), apart from its simplicity, is that it is easy to compute. The ANOVA fit $\tilde{y}_{ij}$ requires simple marginal means of $y_{ij}$ along the two factors, and the overall mean. There are several problems though. This simplicity, both in representation and especially in computation, goes away when each observation has weights. In addition, it may be that the main-effects ANOVA model is not well estimated for some values of $i$ or $j$, because of sparsity in these margins.

We propose a method that has the simplicity of (2.2) for a weighted model, that has a multi-level hierarchical structure, and is easy to compute.

**Hierarchical penalized ANOVA model

for a three-way table**

Fit $\mu^{(0)}_{ijk}=\bar{y}_{...}$, the overall weighted mean of all the $y$s. 2. 1.

Fit the additive model $\mu^{(1)}_{ijk}=\alpha_{i}^{1}+\alpha_{j}^{2}+\alpha^{3}_{k}$ by solving the following weighted penalized least squares (WPLS) problem:

[TABLE]

This can be solved using a simple backfitting algorithm Hastie and Tibshirani (1990), where each step is a version of (2.2) that only requires table summing. Here we shrink to the overall mean, so if a particular one-way marginal count is low, the shrinkage will still kick in. 3. 2.

Fit the second-order interaction model $\mu^{(2)}_{ijk}=\psi_{ij}^{12}+\psi_{jk}^{23}+\psi^{13}_{ik}$ by solving the WPLS problem:

[TABLE]

This is the case we were concerned about in Section 1, when some two-way tables have sparse counts. In this case, that term is shrunk more towards its parent. 4. 3.

Fit the full third-order model $\mu^{(3)}_{ijk}$ by solving the WPLS problem

[TABLE]

This is where we started, i.e. model (2.2), with $\gamma_{ijk}=\tilde{y}^{(2)}_{ijk}$. The term we shrink to has been regularized to accommodate lower-level sparsity.

We have omitted the details of backfitting in steps 1 and 2. Consider using coordinate-descent to solve (1.), where for example we hold $\alpha^{2}_{j}$ and $\alpha^{3}_{k}$ fixed at their current estimates and solve for $\alpha^{1}_{i}$. By collecting terms, this is again of the form (2.2), and so the update involves simple table averaging of a modified response and shrinkage target.

There are many interesting details that haven’t been discussed yet. For example, we have the same $\lambda$ at each level; these could be different. We can take the empirical Bayes analogy further and use variance components to suggest values for $\lambda$, or at least relative values. We will discuss this in Section 4. Another option is to terminate the hierarchy early. In the restaurant ratings example, three way table is manageable even with many levels for each factor. However, when the number of factors grows, the number of cells grows at a exponential rate. In this case we might truncate the hierarchy early, e.g. at step 2, and use $\tilde{y}^{(2)}_{ijk}=\tilde{\psi}_{ij}^{12}+\tilde{\psi}_{jk}^{23}+\tilde{\psi}^{13}_{ik}$ as the estimate. Note in this case we would only need to store and compute at the level of two-way tables.

3. Extension to General Tables

Let’s extend our algorithm to a general table with $m\geq 1$ factors. All the cell means are represented by a vector $\mathbf{y}$. Each cell mean is indexed by its factor levels $I=(i_{1},i_{2},\cdots,i_{m})$. For example, if $m=3$, then the cell that is in the 2nd, 4th and 1st level of the three factors corresponds to $y_{I}=y_{241}$. Every cell is also associated with a weight, denoted by $n_{I}$.

Let $\Omega$ be the set of observed indices. The general statistical question we are trying to address is: suppose the ratings follow an underlying statistical model $\mathbf{y}=\bm{\mu}+\bm{\epsilon}$ and $\mu_{I}=f(I;\bm{\beta})$ for some deterministic or random function $f$, if we have only observed $\mathbf{y}_{\Omega}$, how do we estimate the underlying cell means $\bm{\mu}$?

3.1. General Algorithm

In this paper, our algorithm assumes a hierarchy of models $f_{0},f_{1},\cdots f_{m}$ increasing in complexity. Each time the algorithm uses fits from the previous model $f_{k}$ as the prior for a new and more complex model $f_{k+1}$. This is a very general hierarchical model and we may choose different sequence of model functions in different applications.

In our problem, the model function $f_{k},~{}0\leq k\leq m$ is a $k$-th order interaction model, i.e.

[TABLE]

Here $J$ is a subset of $[m]=\{1,2,\ldots,m\}$, indicating the factors we are considering; $\beta^{J}_{I_{J}}$ is the ”mean” effect of cells that are in levels $I_{J}$ of the factors specified in $J$, and $\bm{\beta}^{(k)}=\{\bm{\beta}^{J}\mathrel{\mathop{\mathchar 58\relax}}~{}J\subset[m],|J|=k\}$ is the model parameter for $f_{k}$.

Our algorithm is summarized in Algorithm 1. Notice that this idea can be generalized to solve any hierarchical model, not necessarily the full additive model in (3.1).

3.2. Previous Methods

Many statisticians studied this kind of problem in the 1980’s using Bayes models. One typical model (DuMouchel and Harris, 1983) is

[TABLE]

and they also want to extrapolate for some missing cells in the table. The empirical Bayes method is to put a normal prior on $\alpha$ and $\beta$

[TABLE]

and estimate the hyperprior variance parameter $V$ using the data. A more general Bayesian ANOVA model and many more examples are considered in the tutorial paper by Casella (1992). See also Searle et al. (1992, ch. 9) and Searle (2006).

More recently, Beran (2005) considered ANOVA from a pure shrinkage and optimization view. The author considers the whole class of estimators that solve penalized least squares, which may correspond to some Bayes modeling or not. Then among all the possible shrinkage estimators, the paper tries to find the one that minimize estimated risk. This is a very sound approach from a theoretical point of view, but the algorithm proposed in the paper scales poorly and couldn’t be used in real data sets.

In Gelman (2005), the author uses a hierarchical Bayes model to analyze Analysis of Variance in an universal way. The statistical model is

[TABLE]

where $m$ stands for a batch of regression coefficients (or in terms of ANOVA, variance component). Then he assumes normal priors on these coefficients

[TABLE]

and puts hierarchical prior on $\sigma_{m}^{2}$. This bears some similarity to our statistical model described in 3.1. Gelman’s paper focuses mostly on the relationship between this hierarchical Bayes model with classical ANOVA and how to use Gibbs sampler to sample from the posterior distribution of $\sigma_{m}^{2}$.

Following this idea, a more recent paper by Volfovsky and Hoff (2014) considers a three-way factorial model

[TABLE]

where $\epsilon_{ijkl}$ are i.i.d. $\mathrm{N}(0,\sigma^{2})$. However, the effects here are possibly correlated, which makes the problem different from the usual settings. The authors put a normal prior with zero mean and general covariance matrix on the effects and treat the covariance matrices as parameters with invert-Wishart priors. By plugging in empirical Bayes estimate for the hyperpriors, they can run a Gibbs sampler to find the posterior.

However, despite the flexibility of hierarchical Bayes model, it usually needs a Markov Chain Monte Carlo (MCMC) algorithm to sample from posterior distribution. This is unacceptable for more than a few hundred cells. Also, as noticed in Section 1, our focus is on estimating and predicting cell means instead of inference about variance components.

From the perspective of high dimensional estimation, shrinkage is known to be very effective, for example the James-Stein estimator dominates maximum likelihood estimator (Efron, 2010). Recently, in a multi-task averaging problem, Feldman et al. (2012) uses James-Stein estimator with empirically estimated covariance matrix and claims to outperform James-Stein in some occasions. However, the multi-task averaging problem has no associated covariates, and the structural information of these factors are crucial in our algorithm.

In light of the hierarchical Bayes model considered so far, our algorithm can be viewed as an empirical solution to a hierarchical Bayes model.

4. Theory

In this section we develop some theoretical results for HANOVA. These also provide a method to empirically select the regularization parameter $\lambda$ in the algorithm.

4.1. Balanced Table

Let’s first assume all the observed cells have the same weight (i.e. same number of observations).111This is different from what a ”balanced table’ is meant in most of the statistics literature. In our paper, we say a table is balanced if all the observed cells have the same number of observations, but the table could contain many empty cells. Assume we observe $n$ cells and the responses are centered. At each stage, we are maximizing the log ”likelihood” function

[TABLE]

subject to $\bm{\mu}^{(k)}\in\mathcal{S}_{k}$. Here $\mathcal{S}_{k}$ is the linear subspace of $\mathbb{R}^{n}$ generated by all possible $k$-th order effects. This implies the subspaces are nested $\mathcal{S}_{1}\subset\mathcal{S}_{2}\subset\ldots\subset\mathcal{S}_{k}$. Say the dimension of $\mathcal{S}_{k}$ is $d_{k}$ and the projection matrix onto $\mathcal{S}_{k}$ is $\mathbf{P}_{k}$. This means $\mathbf{P}_{k}=\mathbf{U}_{k}\mathbf{U}_{k}^{T}$ where the columns of $\mathbf{U}_{k}\in\mathbb{R}^{n\times d_{k}}$ are orthonormal basis for $\mathcal{S}_{k}$, i.e. $\mathbf{U}_{k}^{T}\mathbf{U}_{k}=\mathbf{I}_{d_{k}}$. Because the subspaces are nested, we can assume the first $d_{k-1}$ columns of $\mathbf{U}_{k}$ are $\mathbf{U}_{k-1}$, i.e. \mathbf{U}_{k}=\left(\begin{array}[]{cc}\mathbf{U}_{k-1}&\mathbf{V}_{k}\end{array}\right). Moreover suppose \mathbf{U}=\left(\begin{array}[]{cc}\mathbf{U}_{K}&\mathbf{V}\end{array}\right)\in\mathbb{R}^{n\times n} is an orthogonal matrix.

Consider a hierarchical Bayes model that is a special case of (3.1).

[TABLE]

As the next theorem indicates, this hierarchical model is closely related to maximizing the log-likelihood (4.1). In the usual ANOVA model, each factor may have totally different main effects (or interactions) and we can test for the significance of them. However, as pointed out in Section 1, we are no longer interested in any individual factor. By using orthogonal matrices in the hierarchical model (4.2), we implicitly treat the effects and interactions of the observable factors as ”randomly chosen” directions in the column space of $\mathbf{U}_{m}$.

Theorem 1**.**

For balanced table, each step of the hierarchical procedure (Algorithm 1) for (4.2) is equivalent to maximizing (4.1) subject to $\bm{\mu}^{(k)}\in\mathcal{S}_{k}$, with

[TABLE]

and $\bm{\beta}^{(k-1)}$ being the fit from last step.

Proof.

See Appendix A.1. ∎

Theorem 1 also indicates we can estimate all the variance parameters $\sigma_{k}^{2}$ from the data and then compute a vector of $\bm{\lambda}$. Let’s assume the ratings are centered (by the global mean) so we put $\bm{\beta}^{(0)}=\mathbf{0}$. Then the multi-level hierarchical model (4.2) is actually a linear random effects model

[TABLE]

Now we can state our main theorem that guarantees the effectiveness of our algorithm for balanced table.

Theorem 2**.**

For balanced table, if we use (4.3) to compute $\bm{\lambda}$, the solution to our algorithm is the posterior mean of $\bm{\beta}$ in the linear random effects model (4.4).

Proof.

See Appendix A.2. ∎

In reality we don’t know what $\sigma_{0}^{2}$, $\sigma_{1}^{2}$,…, $\sigma_{m}^{2}$ are, but we can use various methods developed in ANOVA to estimate them. By doing this, our algorithm is equivalent to an empirical Bayes solution of the multi-level hierarchical model.

Following the general principle of variance component analysis, we can compute the expectation of certain quadratic forms of $\mathbf{y}$, then use these quadratic forms to find an unbiased estimator of $\sigma_{k}^{2}$. For example,

[TABLE]

Also

[TABLE]

Thus we have $m+1$ equations and $m+1$ parameters (recall we assume $\sigma^{2}$ is known), the unbiased estimator of the prior variances can be obtained by solving this linear system.

In fact, if we call $\mathbf{U}_{m+1}=\mathbf{U}$ and $d_{m+1}=n$, we have

[TABLE]

So

[TABLE]

are unbiased for estimating $\sigma_{k}^{2}$. Note that it is also possible to treat $\sigma^{2}$ as a tuning parameter, but usually $\sigma^{2}$ can be estimated quite accurately from the data.

The above formulae is only one specific (perhaps the easiest) choice of estimating equations. More general methods in analysis of variance for unbalanced tables, such as Restricted Maximum Likelihood (REML) can also be used (see, for example, Searle et al. (1992)).

4.2. Unbalanced Table

If the observed table is unbalanced, then we cannot prove any exact conclusions like 1 or 2. In this case, the log likelihood function we maximize is

[TABLE]

where $\mathbf{N}$ is a diagonal matrix with weights.

Similar to the previous section, we may sequentially maximize (4.9) for $k=1,2,\ldots,m$ to estimate $\bm{\mu}$ and $\bm{\beta}$ in the following linear random effects model

[TABLE]

Here $\bm{\Sigma}$ is the same as the one (4.4). This model is simply the unbalanced version of (4.4).

The theory we derived in Section 4.1 can be used to give us a vector of reasonable penalty parameters $\bm{\lambda}$. One way to do this is to pretend all the observed cells have the same weight and use equations in (4.8) to estimate all the $\sigma_{k}^{2}$. Then one can compute ratios between the variance components (i.e. (4.3)) to get $\bm{\lambda}$. The estimate $\hat{\sigma}_{k}^{2}$ is still unbiased, because the estimation equations (4.5) and (4.6) still hold true after replacing $\sigma^{2}$ with $\sigma^{2}\mathrm{tr}(\mathbf{N}^{-1}/n)$.

5. Implementation

In this section we discuss some implementation details of HANOVA.

5.1. Preprocessing

In real applications, many data sets are not summarized in the format of cell means and weights. Usually we may have multiple units in a cell and every unit receive multiple reviews. A mixed-effect model describing this is

[TABLE]

Here $\mu_{c}$ is the cell mean, $\alpha_{i}$ is the unit effect, $\sigma_{r}^{2}$ is the variance of user’s rating. $\sigma_{r}$ may actually depend on the unit, and the procedure below can be slightly modified to this heterogeneous case. $n_{i}$ is the number of reviewers of that restaurant, $\sigma_{u}^{2}$ is the variance of the unit effect.

The preprocessing procedure for this model is described in Algorithm 2. Note that since we use the estimated variance directly as weight in step 3, the $\sigma^{2}$ defined in (4.2) is $1$.

5.2. Normal Equations

In our algorithm, the $k$-th order model is

[TABLE]

Here $\mu^{(k)}_{I}=\sum_{J\subset[m],|J|=k}\beta^{J}_{I_{J}}$. Differentiate the above loss with respect to each $\mu^{J}_{L}$ and equate to zero

[TABLE]

So the coefficient of $\beta^{K}_{M}$ in the $(J,L)$-th equation is

[TABLE]

Now the problem reduces to solve a large system of linear equations

[TABLE]

5.3. Backfitting

We implemented the block coordinate descent algorithm (or backfitting algorithm ) to solve the above equations. Each block contains $\beta^{J}_{\cdot}$, i.e. all the effects for certain margins. For example, when $J$ contains only one factor, this block contains all the main effects associated with that factor. The psuedocode of our algorithm is in Algorithm 3.

Note that the basic building block of every iteration is computing $\sum_{I_{J}=L}(n_{I}+\lambda_{k})\mu_{I}(\mathrm{old})$. For a fixed $J$, this amounts to a weighted table sum over $J^{C}$. This operation can be easily parallelized.

6. Results

6.1. Simulations

We first use simulations to verify the optimality of HANOVA. We simulate $200$ instances from the linear mixed effect model (4.4) with four factors (each has 10 levels) and the true model contains 2-way interactions. The variance parameters are chosen to be $(\sigma_{0},\sigma_{1},\sigma_{2},\sigma)=(2,1,0,0.5),\mathrm{or~{}}(2,1,0,1),\mathrm{or~{}}(1,2,0,1)$ (fairly strong signal). The empirically estimated $\lambda_{1}$ sometimes can be infinity. In this case, we truncate $\lambda_{1}$ to $5$.

The results are shown in the violin plots in Figure 1(a) to 1(c). In all three cases HANOVA with oracle $\bm{\lambda}$ achieves Bayes risk and HANOVA with empirical $\bm{\lambda}$ performs almost as good as the oracle. The unshrinked linear model suffers from increased noise, as we can see from the first and second plots.

The next simulation suggests HANOVA can fit a high order model without overfitting. In this simulation, the true model is of order $3$ but the three way interactions are very weak ($(\sigma_{0},\sigma_{1},\sigma_{2},\sigma_{3},\sigma)=(2,1,\mathbf{0.5},0,1)$). The cells can have up to $10$ times different weights.

We generate $50$ instances and plot the estimation RMSE in Figure 1(d). The unshrinked linear model with all the three-way interactions performs poorly, due to overfitting. The second order HANOVA model is slightly better than second order linear model, and the third order model is able to squeeze a little more.

6.2. Real Data

6.2.1. IMDb Data

We ran our algorithm on an IMDb movie ratings data set. This data set contains all the directors who have more than 20 movies on IMDb record. The data set is a $5\times 20\times 33$ table with 276 observations (about $8\%$ of all the cells). The three factors in the data set are decade, genre and director.

Since the observed cell is too sparse, we will only fit our hierarchical penalized model to first order and compare it with the unshrinked main effects model. The procedure described in Section 4.2 gives a suggested $\lambda_{1}=0.064$. We further ran cross-validation on the training set and choose to use $\lambda_{1}=0.1$. We use the squared root of the number of movies as cell weights.

The first order HANOVA with $\lambda_{1}$ has test set Root Mean Squared Error (RMSE)$=0.798$, while the unshrinked main effects model has $\mathrm{RMSE}=0.804$ and the cross-validated ridge regression main effects model has $\mathrm{RMSE}=0.801$.

6.2.2. Zagat Data

This data set contains survey ratings of 1776 restaurants in New York City in $33$ zip-codes, $25$ different type of cuisines and $7$ price levels. Each rating is from [math] to $30$. $20\%$ of the restaurants are left out as test data.

The standard deviation of reviewer random effect is about $7.51$. This large random effect makes the ratings of many unpopular restaurants not trustable, as pointed out in Section 1. The number of reviews of each restaurants follow a heavy-tailed distribution. One restaurant has $8922$ reviews but $106$ restaurants have less than or equal to $5$ reviews. The average and median number of reviews are $117.1$ and $42.0$.

The standard deviation of restaurant random effect (within cell) is about $2.47$. This is the lower limit of RMSE of any prediction methods based on the three factors. We use the method described in 5.1 to preprocess the data.

The $\sigma_{m}$ estimated by HANOVA (defined in (4.2)) is $(0.423,0.149,0)$ (recall the $\sigma^{2}$ is always $1$ after the preprocessing procedure 2) and the corresponding $\lambda=(5.71,45.0,\inf)$. This is the case that the data has fairly strong main effects and weak interactions. Since the third-order interactions are too weak, we only fit the first three HANOVA models. The RMSEs are $3.010$ (HANOVA grand-mean model), $2.641$ (HANOVA main effects model) and $2.638$ (HANOVA second-order interactions model). The best $\lambda_{2}$ selected by cross-validation is $72.5$ and the corresponding second-order interactions model has RMSE$=2.625$. As comparison, the RMSE of random main effects model is $2.640$ and the RMSE of price-only random effects model is $2.767$.

The RMSE reduction of HANOVA interactions model is more significant in cells that have fewer restaurants. This is one of HANOVA’s the main purposes. Among the $355$ restaurants in the test data set, $32$ of them receive rating adjustment greater than 1 from main effects model to interactions model. There mean absolute prediction error for these $32$ restaurants are reduced from $2.031$ to $1.707$ by fitting a interactions model by HANOVA.

Appendix A Proof of Theorems

A.1. Theorem 1

First look at

[TABLE]

The log of the density function is

[TABLE]

The distribution of $\mathbf{y}$ given $\bm{\beta}_{k}$ is

[TABLE]

Here $\mathbf{D}$ is a diagonal matrix

[TABLE]

Thus

[TABLE]

Notice that \mathbf{U}^{T}(\mathbf{y}-\mathbf{U}_{k}\bm{\beta}^{(k)})=\left(\begin{array}[]{c}\mathbf{U}_{k}^{T}\mathbf{y}-\bm{\beta}^{(k)}\\ \mathbf{U}_{k}^{\perp^{T}}\mathbf{y}\\ \end{array}\right), so

[TABLE]

The claim is immediately followed by the fact that $\|\bm{\beta}^{(k)}-\mathbf{U}_{k}^{T}\mathbf{y}\|^{2}=\|\mathbf{U}_{k}(\bm{\beta}^{(k)}-\mathbf{U}_{k}^{T}\mathbf{y})\|^{2}=\|\bm{\mu}^{(k)}-\mathbf{P}_{k}\mathbf{y}\|^{2}=\|\bm{\mu}^{(k)}-\mathbf{y}\|^{2}-\|\mathbf{y}-\mathbf{P}_{k}\mathbf{y}\|^{2}$. The regularization parameter is

[TABLE]

.

A.2. Theorem 2

Recall our algorithm is just iteratively maximizing likelihood (4.1) to get $\hat{\bm{\beta}}_{m}$ and then find the posterior mean of $\bm{\beta}$ given $\mathbf{y}$ and $\bm{\beta}^{(k)}=\hat{\bm{\beta}}_{m}$.

By empirical Bayes estimator of $\bm{\beta}$, I mean just compute the posterior mean of $\bm{\beta}$ given $\mathbf{y}$ in (4.4) and plug in whatever $\sigma^{2}$ and $\sigma_{k}^{2}$ we used to compute $\lambda_{k}$.

The posterior mean of (4.4) can be computed by Tweedie’s formulaRobbins (1964). The marginal distribution of $\mathbf{y}$ is $\mathrm{N}(\mathbf{0},(\sigma^{2}+\sigma_{m}^{2})\mathbf{I}_{n}+\mathbf{U}_{m}\bm{\Sigma}\mathbf{U}_{m}^{T})$, i.e.

[TABLE]

so the posterior mean is given by

[TABLE]

On the other hand, in the case of a balanced table, each step of our algorithm produces

[TABLE]

and finally

[TABLE]

It suffices to show (A.9) and (A.11) are actually the same. Notice that $\mathbf{U}_{k}\mathbf{U}_{k}^{T}\hat{\bm{\mu}}^{(j)}=\hat{\bm{\mu}}^{(j)}$ for $j<k$ because $\hat{\bm{\mu}}^{(j)}\in\mathcal{S}_{j}\subset\mathcal{S}_{k}$ and $\mathbf{U}_{k}\mathbf{U}_{k}^{T}$ is just the projection matrix onto $\mathcal{S}_{k}$. With some calculation,

[TABLE]

Since $\mathbf{U}_{k}\mathbf{U}_{k}^{T}\mathbf{y}=\sum_{j=1}^{k}\mathbf{V}_{j}\mathbf{V}_{j}^{T}\mathbf{y}$, it is easy to see that

[TABLE]

The theorem is immediately proved if we plug this in (A.11) and compare it to (A.9).