Varianza de predicción GLM con observaciones medias

Question

Supongamos que tengo un conjunto de datos en el que los valores observados son promedios y no necesariamente puntos de datos individuales. Por ejemplo, supongamos que el registro 1 tiene el valor observado de $Y_1 = 2.0$ . Sin embargo, sé que el número de observaciones de este registro es 3. Es decir, $Y_1$ es una media de tres puntos de datos digamos $1.0, 2.0, 3.0$ o incluso $2.0, 2.0, 2.0$ . Se desconoce la varianza asociada al registro 1.

Supongamos además que tengo 25.000 registros y que puedo construir una matriz de diseño con todas estas medias como respuesta. Tenga en cuenta que el peso observado es correcto.

Creo que un GLM con una matriz de diseño como la descrita produciría medias similares a un GLM en el que la matriz de diseño contiene todos los puntos de datos subyacentes. El problema radica en la varianza de la predicción.

Mi pregunta es: ¿puedo confiar en los errores estándar producidos? Sin embargo, no de la manera convencional. Tendría que hacer algún ajuste porque la varianza de la predicción tiene que ser menor.

¿Qué opinas?

Answer 1

Creo que debería buscar el uso de ponderaciones de tal manera que se minimice una suma de cuadrados residual ponderada. Echa un vistazo a la documentación de ${\tt glm}$ en ${\tt R}$

https://stat.ethz.ch/R-manual/R-patched/library/stats/html/glm.html

o eche un vistazo a la documentación de SAS

http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_glm_sect024.htm

Las ponderaciones que debe utilizar son simplemente el número de observaciones subyacentes. Las ponderaciones serán inversamente proporcionales a las dispersiones de las observaciones. Esto es lo que usted quiere, ya que una media de variables iid tiene una varianza que es inversamente proporcional al número de sumandos.

EDIT (respondiendo al comentario):

Consideremos el modelo normal lineal. Suponemos que $Y_1,...,Y_N$ et $Z_1,...,Z_n$ son variables estocásticas tales que el $Y_i$ son independientes y el $Z_i$ son medios de subconjuntos disjuntos del $Y_i$ 's, $Z_i = \frac{1}{m_i} \sum_{j=1}^{m_i} Y_{ij}$ con una reindexación adecuada del $Y_i$ 's. Para simplificar las cosas, supongamos también que no tenemos covariables. En este caso, la MLE del único parámetro medio será la media del $Y_i$ si estamos modelando el $Y_i$ . En cambio, podemos modelar el $Z_i$ utilizando los pesos correspondientes a $m_1,...,m_n$ . Son estos dos enfoques de modelado los que estamos comparando. Al modelar el $Z_i$ el MLE del parámetro medio será la media ponderada de los $Z_i$ , $\bar{Z}$ . Podemos calcular fácilmente que $\bar{Z}= \bar{Y}$ . ¿Pero qué pasa con el parámetro de la varianza? Para el $Y$ -modelo el estimador ML será

$ \hat{\sigma} = \frac{1}{N}\sum_{i=1}^N (Y_i - \bar{Y})^2 $

y para el $Z$ -modelo

$\bar{\sigma} = \frac{1}{n}\sum_{i=1}^n m_i(Z_i - \bar{Z}).$

Sin embargo, estos estimadores no son iguales, de modo que, dados algunos datos, podemos obtener diferentes estimaciones del parámetro de la varianza utilizando los dos enfoques. Ambos estimarán el parámetro verdadero de forma consistente, pero sabemos que en el $Y$ -el MLE es asintóticamente óptimo en una gran clase de estimadores. Y en realidad el MLE en el $Z$ -está en este modelo ya que es sólo una bonita función del $Y_i$ 's. En términos generales, hay una pérdida de información en la transformación de la $Y_i$ 's a la $Z_i$ 's, tal vez no sea sorprendente.

Sin embargo, normalmente utilizamos los estimadores insesgados de la varianza escalando los MLE. Podemos estudiarlos más a fondo utilizando la simulación. Sin embargo, asintóticamente son equivalentes a los MLE, por lo que sus varianzas se relacionarán entre sí de la misma manera que para los MLE.

La simulación dibuja el $Y$ -vector 10000 veces y para cada una de ellas calcula un $Z$ vectorial. A continuación, se ajustan ambos modelos y se calcula la estimación de la varianza (parámetro de dispersión). Posteriormente se inspeccionan visualmente las distribuciones de estas estimaciones.

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))

Al menos en esta simulación, vemos que utilizando esta $Z$ -modelo da más varianza en la estimación del parámetro de dispersión.

En conclusión, si uno tiene la $Y$ no es aconsejable hacer este tipo de reducción de datos antes de analizarlos.