Linear regression with multiple y values per x and their arithmetic means I noticed that if there are multiple values of $y$ for each value of $x$, I can replace the values of $y$ with their arithmetic mean at each value of $x$ and still get the same regression model (obtained via least squares). This seems to be true only if the number of $y$ values is constant across all values of $x$.
Here is some R code to illustrate this question:
set.seed(12)

df <- data.frame(x=1:10, y=1:10 + runif(100, max=5))
df_agg <- aggregate(y ~ x, data=df, FUN=mean)

with(df, plot(x, y, pch=19, col=rgb(0, 0, 0, 0.25)))
with(df_agg, points(x, y, pch=17, col="blue"))

(ab <- lm(y ~ x, df_agg))
lm(y ~ x, df)

abline(ab, col="blue")


How can I show mathematically that this relationship is true? I assume it has something to do with the properties of the mean and the sum of residual squares, but I can't seem to be able to figure out why exactly this is the case.
 A: I will use here the classical matrix notation for the linear
regression
$$
  \tag{1}
   \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} 
$$
where $\mathbf{y}$ is a vector of length $n$ and $\boldsymbol{\beta}$
is a vector of length $p$. The $n \times p$ design matrix $\mathbf{X}$
has as its $i$-th row the transpose of the vector $\mathbf{x}_i$ of
covariates. The error vector $\boldsymbol{\varepsilon}$ is assumed to
be normal with mean zero and covariance matrix $\sigma^2
\mathbf{I}_n$.
Assume that there are only $m$ distinct values $\mathbf{x}^\star_k$
for the vector of covariates $\mathbf{x}_i$ with $m<n$, and that the
value $\mathbf{x}^\star_k$ corresponds to $n_k$ observations, so that
$n = n_1 + \dots + n_m$. We can further assume then that the
observations are given by groups corresponding to the distinct values
of the covariates: $\mathbf{x}_1^\star$, $\dots$,
$\mathbf{x}_m^\star$. So, the observations $1$ to $n_1$ correspond to
$\mathbf{x}^\star_1$, then the $n_1 + 1$ to $n_1 + n_2$ correspond to
$\mathbf{x}^\star_2$ and so on.
We can consider then two matrices $\mathbf{A}$ and $\mathbf{G}$. The
"aggregating matrix" $\mathbf{A}$ has size $m \times n$ and transforms
a vector with length $n$ into the vector of the group averages, with
length $m$. The matrix $\mathbf{G}$ has size $n \times m$ and "picks
out" the group $k_i$ corresponding to the observation $i$. These are
$$ \mathbf{A} = \begin{bmatrix} n_1^{-1}
\mathbf{1}^\top_{n_1} & & & \\
& n_2^{-1} \mathbf{1}^\top_{n_2} & & \\
& & \ddots & \\  
&  & &  n_m^{-1} \mathbf{1}^\top_{n_m}
\end{bmatrix} \qquad
\mathbf{G} = \begin{bmatrix} 
\mathbf{1}_{n_1} & & & \\
& \mathbf{1}_{n_2} & & \\
& & \ddots & \\  
&  & &  \mathbf{1}_{n_m}
\end{bmatrix} 
$$
where $\mathbf{1}_r$ stands for a vector of $r$ ones.
Note that $\mathbf{A} \mathbf{G}= \mathbf{I}_m$ and that
$\mathbf{G}^\top\mathbf{G}$ is the diagonal matrix
$\mathbf{D} := \text{diag}_i \{n_i \}$ with size $m$. Note also
that $\mathbf{A} = \mathbf{D}^{-1} \mathbf{G}^\top$.
Now consider the "aggregated regression": Its response vector is
$\mathbf{y}^\star := \mathbf{A} \mathbf{y}$ with length $m$ and its design
matrix is $\mathbf{X}^\star = \mathbf{A} \mathbf{X}$ with dimension
$m \times p$. It is clear that the aggregated regression takes the form
$$
\tag{2}
\mathbf{y}^\star = \mathbf{X}^\star \boldsymbol{\beta} +
\boldsymbol{\varepsilon}^\star
$$
where $\boldsymbol{\varepsilon}^\star := \mathbf{A}
\boldsymbol{\varepsilon}$. In other words, the aggregated regression
is obtained by left multiplying (1) by $\mathbf{A}$, the parameter
$\boldsymbol{\beta}$ being the same in both cases. The question is
whether the OLS estimates say $\widehat{\boldsymbol{\beta}}$ and
$\widehat{\boldsymbol{\beta}}_{\text{ag}}$ corresponding to (1) and
(2) are the same. Note that the error $\boldsymbol{\varepsilon}^\star$
has mean zero and covariance $\sigma^2 \mathbf{D}^{-1}$, so the
Maximum-Likelihood estimate for the aggregated regression is given by
Weighted Least Squares (WLS) in the general case where the $n_i$ are
different.
Since $\mathbf{X} = \mathbf{G} \mathbf{X}^\star$, we have 
$$
   \widehat{\boldsymbol{\beta}} = [\mathbf{X}^\top \mathbf{X}]^{-1}
   \mathbf{X}^\top \mathbf{y} =
   [\mathbf{X}^{\star\top}\mathbf{G}^\top \mathbf{G} \mathbf{X}^\star]^{-1}
   \mathbf{X}^{\star\top} \mathbf{G}^\top \mathbf{y} =
   [\mathbf{X}^{\star\top}\mathbf{D} \mathbf{X}^\star]^{-1}
   \mathbf{X}^{\star\top} \mathbf{G}^\top \mathbf{y}.
$$
If all the $n_i$ are equal we have $\mathbf{D} = n_1
\mathbf{I}_m$ and $\mathbf{G}^\top \mathbf{y} = n_1 \mathbf{y}^\star$,
so
$$
  \widehat{\boldsymbol{\beta}} = [\mathbf{X}^{\star\top}\mathbf{X}^\star]^{-1}
   \mathbf{X}^{\star\top} \mathbf{y}^\star,
$$
which is the OLS estimate for the aggregated regression. In the general case,
 WLS should be used for the aggregate regression,
 with weights proportional to the inverse group sizes. The estimates
 of both regressions will then continue to be the same.
A: are you saying that the regression of the means is the same as regression of the whole data set if N is same for each x?
this seems counter-intuitive to me. it seems that the error on beta should be dependent upon the dispersion of the original data set whereas this is not apparent in using the means. this is important for hypothesis testing when comparing betas from two data sets/curves or the difference of beta from zero.
