Regression model with aggregated targets Similar as in this self-answered question, I want to ask about possible approaches for modelling data with aggregated targets, i.e. things like
$$
\bar y_{j[i]} = \alpha + \beta x_i + \varepsilon_i
$$
where $j[i]$ is the $j$-th group, where $i$-th observation belongs, and for each $j$-th group of size $|j|$, we are predicting the target that is an average value of all the $y_i$ observations within the group, $\bar y_{j[i]} = |j|^{-1} \sum_{i \in j[i]} y_i$. Of course, the means are given, and cannot be disaggregated, this is the data we have.
Additional assumption that can be made in here, is that there is clustering within the $j[i]$ groups, so the group assignment is not completely random, the subjects within each group share some characteristics.
For example, imagine that you have data on average test score per class (something to predict), and features on both student level, e.g. individual IQ scores (that should be highly predictive, but not perfect, for exam scores), class level features, and features on higher level of aggregation (school level). I am interested in finding factors that contributed to each individual test score, and predict them. The data is a random sample of classes, the final predictions will be made for students from classes that were not observed in the training data.
Can we use such data to learn anything (approximately) about the unobserved individual-level targets?
What are the approaches used for modelling such data? Can you give some references? Obviously with aggregated data we loose precision, and the variance of the means $\bar y_{j[i]}$ is smaller then of the individual observations $y_i$, so predicting the average target is not the same as predicting individual values. Is there any way how to translate the predictions of the group averages to possible variability between subjects? 
 A: To verify the solution suggested in the great answer by @user20160 I prepared a toy example that demonstrates it. As suggested by @user20160, I am posting the code as a supplement to the answer. For explanations of this approach, check the other answer.
First, let's generate the independent variable and append the column of ones to it, to use matrix formulation of the model.
set.seed(42)
n <- 5000; k <- 50; m <- n/k

x <- rnorm(n, mean = (1:n)*0.01, sd = 10)
X <- cbind(Intercept=1, x)

Next, let's generate the individual predictions $y = X\beta + \varepsilon$.
beta <- rbind(3, 0.75)
sigma <- 10
y <- rnorm(n, X %*% beta, sigma)

To aggregate the results, we use the matrix $W$ of zeros and ones to indicate group membership of size $k \times n$. To estimate group means, we take $\bar y = \tfrac{1}{m}W y$ (same results as tapply(y, grp, mean)).
grp <- factor(rep(1:k, each=m))
W <- t(model.matrix(~grp-1))
ybar <- as.vector((W/m) %*% y)

What leads to the following results, where as expected, the conditional variability of $\bar y$ is much smaller then $y$.

lm_loss <- function(pars) mean((mu_rep - as.vector(X %*% pars))^2)
aggr_loss <- function(pars) mean((mu - as.vector((W/m) %*% (X %*% pars)))^2)

Results from the regular regression model are pretty poor.
init <- rbind(0, 0)
(est1 <- optim(init, lm_loss))$par
##          [,1]
## [1,] 9.058655
## [2,] 0.502987

The "aggregated" model gives results that are really close to the true values of $\beta$.
(est2 <- optim(init, aggr_loss))$par
##           [,1]
## [1,] 3.1029468
## [2,] 0.7424815

You can also see on the plot below, that besides that the input data was aggregated, if we use the "aggregated" model, we are able to recover the true regression line almost perfectly.

Also if we compare mean squared error of predictions for the individual values given the estimated parameters, the "aggregated" model has smaller squared error.
mean((y - as.vector(X %*% est1$par))^2)
## [1] 119.4491
mean((y - as.vector(X %*% est2$par))^2)
## [1] 101.4573

Same thing happens if we minimize the negative log-likelihood. Additionally, this lets us to estimate $\sigma$, and also gives much better result (43.95 for linear regression vs 8.02 for the "aggregated" model).
lm_llik <- function(pars) -1 * sum(dnorm(mu_rep, as.vector(X %*% pars[1:2]), pars[3]/sqrt(k), log=TRUE))
aggr_llik <- function(pars) -1 * sum(dnorm(mu, as.vector((W/m) %*% (X %*% pars[1:2])), pars[3]/sqrt(k), log=TRUE))

A: Here's an approach for solving this type of problem using latent variable models. It's not a specific model, but a general way to formulate a model by breaking the description of the system into two parts: the relationship between individual inputs and (unobserved) individual outputs, and the relationship between individual outputs and (observed) aggregate group outputs. This gives a natural way to think about the problem that (hopefully somewhat) mirrors the data generating process, and makes assumptions explicit. Linear or nonlinear relationships can be accommodated, as well as various types of noise model. There's well-developed, general-purpose machinery for performing inference in latent variable models (mentioned below). Finally, explicitly including individual outputs in the model gives a principled way to make predictions about them. But, of course there's no free lunch--aggregating data destroys information.
General approach
The central idea is to treat the individual outputs as latent variables, since they're not directly observed.
Suppose the individual inputs are $\{x_1, \dots, x_n\}$, where each $x_i \in \mathbb{R}^d$ contains both individual and group-level features for the $i$th individual (group-level features would be duplicated across individuals). Inputs are stored on the rows of matrix $X \in \mathbb{R}^{n \times d}$. The corresponding individual outputs are represented by $y = [y_1, \dots, y_n]^T$ where $y_i \in \mathbb{R}$.
The first step is to postulate a relationship between the individual inputs and outputs, even though the individual outputs are not directly observed in the training data. This takes the form of a joint conditional distribution $p(y \mid X, \theta)$ where $\theta$ is a parameter vector. Of course, it factorizes as $\prod_{i=1}^n p(y_i \mid x_i, \theta)$ if the outputs are conditionally independent, given the inputs (e.g. if error terms are independent).
Next, we relate the unobserved individual outputs to the observed aggregate group outputs $\bar{y} = [\bar{y}_1, \dots, \bar{y}_k]^T$ (for $k$ groups). In general, this takes the form of another conditional distribution $p(\bar{y} \mid y, \phi)$, since the observed group outputs may be a noisy function of the individual outputs (with parameters $\phi$). Note that $\bar{y}$ is conditionally independent of $X$, given $y$. If group outputs are a deterministic function of the individual outputs, then $p(\bar{y} \mid y)$ takes the form of a delta function.
The joint likelihood of the individual and group outputs can then be written as:
$$p(y, \bar{y} \mid X, \theta, \phi) = p(\bar{y} \mid y, \phi) p(y \mid X, \theta)$$
Since the individual outputs are latent variables, they must be integrated out of the joint likelihood to obtain the marginal likelihood for the observed group outputs:
$$p(\bar{y} \mid X, \theta, \phi) =
\int p(\bar{y} \mid y, \phi) p(y \mid X, \theta) dy$$
If group outputs are a known, deterministic function of the individual outputs, the marginal likelihood can be written directly without having to think about this integral (and $\phi$ can be ignored).
Maximum likelihood estimation
Maximum likelihood estimation of the parameters proceeds by maximizing the marginal likelihood:
$$\theta_{ML}, \phi_{ML} \ = \
\arg \max_{\theta,\phi} \ p(\bar{y} \mid X, \theta, \phi)$$
If the above integral can be solved analytically, it's possible to directly optimize the resulting marginal likelihood (either analytically or numerically). However, the integral may be intractable, in which case the expectation maximization algorithm can be used.
The maximum likelihood parameters $\theta_{ML}$ could be studied to learn about the data generating process, or used to predict individual outputs for out-of-sample data. For example, given a new individual input $x_*$, we have the predictive distribution $p(y_* \mid x_*, \theta_{ML})$ (whose form we already chose in the first step above). Note that this distribution doesn't account for uncertainty in estimating the parameters, unlike the Bayesian version below. But, one could construct frequentist prediction intervals (e.g. by bootstrapping).
Care may be needed when making inferences about individuals based on aggregated data (e.g. see various forms of ecological fallacy). It's possible that these issues may be mitigated to some extent here, since individual inputs are known, and only the outputs are aggregated (and parameters are assumed to be common to all individuals). But, I don't want to make any strong statements about this without thinking about it more carefully.
Bayesian inference
Alternatively, we may be interested in the posterior distribution over parameters:
$$p(\theta, \phi \mid \bar{y}, X) =
\frac{1}{Z} p(\bar{y} \mid X, \theta, \phi) p(\theta, \phi)$$
where $Z$ is a normalizing constant. Note that this is based on the marginal likelihood, as above. It also requires that we specify a prior distribution over parameters $p(\theta, \phi)$. In some cases, it may be possible to find a closed form expression for the posterior. This requires an analytical solution to the integral in the marginal likelihood, as well as the integral in the normalizing constant. Otherwise, the posterior can be approximated, e.g. by sampling (as in MCMC) or variational methods.
Given a new individual input $x_*$, we can make predictions about the output $y_*$ using the posterior predictive distribution. This is obtained by averaging the predictive distributions for each possible choice of parameters, weighted by the posterior probability of these parameters given the training data:
$$p(y_* \mid x_*, X, \bar{y}) =
\iint p(y_* \mid x_*, \theta) p(\theta, \phi \mid \bar{y}, X) d\theta d\phi$$
As above, approximations may be necessary.
Example
Here's an example showing how to apply the above approach with a simple, linear model, similar to that described in the question. One could naturally apply the same techniques using nonlinear functions, more complicated noise models, etc.
Generating individual outputs
Let's suppose the unobserved individual outputs are generated as a linear function of the inputs, plus i.i.d. Gaussian noise. Assume the inputs include a constant feature (i.e. $X$ contains a column of ones), so we don't need to worry about an extra intercept term.
$$y_i = \beta \cdot x_i + \epsilon_i \quad \quad
\epsilon_i \sim \mathcal{N}(0, \sigma^2)$$
Therefore, $y = [y_1, \dots, y_n]^T$ has a Gaussian conditional distribution:
$$p(y \mid X, \beta, \sigma^2) = \mathcal{N}(y \mid X \beta, \sigma^2 I)$$
Generating aggregate group outputs
Suppose there are $k$ non-overlapping groups, and the $i$th group contains $n_i$ known points. For simplicity, assume we observe the mean output for each group:
$$\bar{y} = W y$$
where $W$ is a $k \times n$ weight matrix that performs averaging over individuals in each group. $W_{ij} = \frac{1}{n_i}$ if group $i$ contains point $j$, otherwise $0$. Alternatively, we might have assumed that observed group outputs are contaminated with additional noise (which would lead to a different expression for the marginal likelihood below).
Marginal likelihood
Note that $\bar{y}$ is a deterministic, linear transformation $y$, and $y$ has a Gaussian conditional distribution. Therefore, the conditional distribution of $\bar{y}$ (i.e. the marginal likelihood) is also Gaussian, with mean $W X \beta$ and covariance matrix $\sigma^2 W W^T$. Note that $W W^T = \text{diag}(\frac{1}{n_1}, \dots, \frac{1}{n_k})$, which follows from the structure of $W$ above. Let $\bar{X} = W X$ be a matrix whose $i$th row contains the mean of the inputs in the $i$th group. Then, the marginal likelihood can be written as:
$$p(\bar{y} \mid X, \beta, \sigma^2) = \mathcal{N} \left(
  \bar{y} \ \Big| \
  \bar{X} \beta, \
  \sigma^2 \text{diag} \big( \frac{1}{n_1}, \dots, \frac{1}{n_k} \big)
\right)$$
The covariance matrix is diagonal, so the observed outputs are conditionally independent. But, they're not identically distributed; the variances are scaled by the reciprocal of the number of points in each group. This reflects the fact that larger groups average out the noise to a greater extent.
Maximum likelihood estimation
Maximizing the likelihood is equivalent to minimizing the following loss function, which was obtained by writing out the negative log marginal likelihood and then discarding constant terms:
$$\mathcal{L}(\beta, \sigma^2) =
k \log(\sigma^2)
+ \frac{1}{\sigma^2}
(\bar{y} - \bar{X} \beta)^T N (\bar{y} - \bar{X} \beta)$$
where $N = \text{diag}(n_1, \dots, n_k)$. From the loss function, it can be seen that the maximum likelihood weights $\beta_{ML}$ are equivalent to those obtained by a form of weighted least squares. Specifically, by regressing the group-average outputs $\bar{y}$ against the group-average inputs $\bar{X}$, with each group weighted by the number of points it contains.
$$\beta_{ML} = (\bar{X}^T N \bar{X})^{-1} \bar{X}^T N \bar{y}$$
The estimated variance is given by a weighted sum of the squared residuals:
$$\sigma^2_{ML} = \frac{1}{k}
(\bar{y} - \bar{X} \beta_{ML})^T N (\bar{y} - \bar{X} \beta_{ML})$$
Prediction
Given a new input $x_*$, the conditional distribution for the corresponding individual output $y_*$ is:
$$p(y_* \mid x_*, \beta_{ML}, \sigma^2_{ML}) = \mathcal{N}(y_* \mid \beta_{ML} \cdot x_*, \sigma^2_{ML})$$
The conditional mean $\beta_{ML} \cdot x_*$ could be used as a point prediction.
References
Machine learning: A probabilistic perspective (Murphy 2012). I don't recall that it speaks specifically about aggregated data, but, it covers concepts related to latent variable models quite well.
A: Different approaches could be appropriate depending on your goal. I'll describe one approach in case your goal is group-level prediction.
You could use the individual-level features to build a bunch of aggregated features for each group (mean, std, median, max, min, ...). You now have richer features for each group which are likely to perform well on the group level. I've seen this work thousands of times in Kaggle competitions. Also, don't stick to linear regression, gradient boosting works in many cases with tabular data, and can even help you weed out some features (make lots of them, you never know what will work).
As a bonus, this also gives you a way of predicting individual scores by feeding the model a group of one (this feels a little shady though).
