GLM Prediction Variance with Average Observations Suppose I have data set where the observed values are averages and not necessarily individual data points. For example, suppose record 1 has the observed value of $Y_1 = 2.0$. However, I know that the number of observations for this record is 3.  That is, $Y_1$ is an average of three data points say $1.0, 2.0, 3.0$ or even $2.0, 2.0, 2.0$.  It is unknown how much variance is associated with record 1.
Further suppose that I have 25,000 records and I can build a design matrix with all these averages as the response.  Keep in mind that the observed weight is correct.  
I believe a GLM with a design matrix as described would produce similar averages to a GLM where the design matrix contains all the underlying data points.  The problem lies with the variance of the prediction.
My question is -- can I rely on the standard errors produced?  Not in the conventional way, however.  I would need to make some adjustment because the variance of the prediction has got to be smaller.  
Thoughts? 
 A: I think you should look into using weights such that a weighted residual sum of Squares is minimized. Take a look at the documentation for ${\tt glm}$ in ${\tt R}$
https://stat.ethz.ch/R-manual/R-patched/library/stats/html/glm.html
or take a look at the SAS documentation
http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_glm_sect024.htm
The weights you should use is simply the number of underlying observations. The weights will be inversely proportional to the dispersions of the observations. This is what you want as a mean of iid variables has a variance that is inversely proportional to the number of summants.
EDIT (answering the comment):
Let consider the linear normal model. We assume that $Y_1,...,Y_N$ and $Z_1,...,Z_n$ are stochastic variables such that the $Y_i$'s are independent and the $Z_i$'s are means of disjoint subsets of the $Y_i$'s, $Z_i = \frac{1}{m_i} \sum_{j=1}^{m_i} Y_{ij}$ with a suitable re-indexing of the $Y_i$'s. To keep things simple let's also assume that we have no covariates. In this case, the MLE of the lone mean parameter will be the mean of the $Y_i$'s if we are modelling the $Y_i$. Instead we can model the $Z_i$'s using weights corresponding to $m_1,...,m_n$. It's these two modeling approaches that we are comparing. When modeling the $Z_i$'s the MLE of the mean parameter will be the weighted mean of the $Z_i$, $\bar{Z}$. We can easily calculate that $\bar{Z}= \bar{Y}$. But what about the variance parameter? For the $Y$-model the ML estimator will be 
$
\hat{\sigma} = \frac{1}{N}\sum_{i=1}^N (Y_i - \bar{Y})^2
$
and for the $Z$-model 
$\bar{\sigma} = \frac{1}{n}\sum_{i=1}^n m_i(Z_i - \bar{Z}).$
These estimators are not equal however, such that given some data we can get different estimates of the variance parameter using the two approaches. They'll both estimate the true parameter consistently, but we know that in the $Y$-model the MLE is asymptotically optimal in a large class of estimators. And actually the MLE in the $Z$-model is in this model as it is just a nice function of the $Y_i$'s. Loosely speaking, there is a information loss in the transformation from the $Y_i$'s to the $Z_i$'s, maybe not surprisingly. 
However, normally we use the unbiased estimators of the variance by scaling the MLE's. We can study them further using simulation. However, asymptotically they are equivalent with the MLE's, thus their variances will related to each other in the same way as for the MLE's.
The simulation draws the $Y$-vector 10000 times and for each of them calculates a $Z$ vector. Then both models are fitted and the estimate of the variance (dispersion parameter) is calculated. The distributions of these estimates are inspected visually afterwards.
group <- rep(rep(rep(seq(40), length.out = 50), length.out = 75), length.out = 100)
simFunc <- function() {
  m <- aggregate(group, list(group), length)$x
  y <- rnorm(100, mean = 1, sd = 2)
  z <- aggregate(y, list(group), mean)$x
  c(summary(glm(y ~ 1))$dispersion, summary(glm(z ~ 1, weights = m))$dispersion)
}

library(ggplot2)

sim <- replicate(10000, simFunc())
sim <- t(sim)

sim <- c(sim)
sim <- data.frame(type = c(rep("Y", 10000), rep("Z", 10000)), est = sim)

ggplot(data = sim, aes(x = est)) + geom_density(aes(group = type, col = type))


At least in this simulation, we see that using this $Z$-model gives more variance in the estimate of the dispersion parameter.
In conclusion, if one has the $Y$ values, it may not be advisable to do this sort of data reduction before analysing data.
