¿Cuál es la mejor manera de encontrar el error estándar *entre* los ajustes de regresión lineal?

Question

Así que tengo un escenario donde hay $n = 8$ sujetos, que se observan en 20 puntos de tiempo y que tienen heteroscedasticidad en su respuesta. Por ejemplo, considere lo siguiente:

num_datasets = 8;

x = [1:20]';

%define matrix for the response for 8 different datasets
Y = repmat(x,1,8) * nan;

for i = 1:size(X,2)
    Y(:,i) = 2*x + unifrnd(3,8)*randn(size(x));
end

Así que claramente cada observación/sujeto tiene el mismo modelo lineal que relaciona su respuesta ( $y$ ) al regresor ( $x$ ), pero las cantidades/fuentes de ruido varían según el tema. Ahora, sé que el error estándar para el ajuste de regresión lineal tiene la forma

$$\sigma\sqrt{\frac{1}{n}+ \frac{(x^*-\bar x)^2}{\sum_{i=1}^n (x_i-\bar{x})^2} }$$

donde $\sigma$ representa la desviación estándar de los residuos del ajuste, $n$ representa el número de muestras en la observación (en mi ejemplo anterior sería 20, no 8), $(x^* - \bar x)$ representa la distancia de cada $x_i$ muestra de la media (por lo que el error estándar aumenta hiperbólicamente a medida que se desvía de la media), y luego ${\sum_{i=1}^n (x_i-\bar{x})^2}$ es simplemente la varianza en $x$ .

Sin embargo, si interpreto correctamente esta ecuación, creo que esto da el error estándar a través de la dimensión de $x$ y no me dice directamente el error estándar entre sujetos. En otras palabras, sospecho que no sería una buena idea utilizar esta fórmula para cada sujeto y luego tomar el error estándar medio (por favor, corríjanme si me equivoco). Así que tengo dos preguntas:

1. ¿Cuál sería la mejor manera de calcular el error estándar entre sujetos? ¿Sería simplemente realizar el ajuste para cada sujeto, y luego tomar la desviación estándar de los ajustes?

2. ¿Cómo sería la forma del error estándar del ajuste, y cuál es la intuición detrás de eso? ¿Seguiría siendo hiperbólico? Creo que no, pero en realidad no estoy seguro.

Answer 1

Lo mejor sería pensar en esta situación en términos de meta-análisis : reunir la información de varios estudios para estimar el modelo de la población subyacente. Los estudios se combinan ponderándolos según la información que aportan, normalmente ponderando de forma inversa cada estudio por la varianza de sus estimaciones.

Puede pensar en su caso como si representara 8 "estudios" diferentes (8 sujetos diferentes), con cada uno de ellos con valores de $y$ medido en 20 valores de $x$ . Suponemos que dentro de cada tema la norma supuestos de la regresión lineal en particular, que las observaciones no están correlacionadas y que la varianza de $y$ sobre la regresión es independiente del valor de $x$ . A diferencia de muchos meta-análisis prácticos que dependen de los resúmenes de los resultados de cada uno de los diversos estudios, usted sigue teniendo los datos individuales de cada "estudio".

Por lo tanto, si se desea un modelo para la población subyacente, una forma sencilla de proceder sería hacer cada una de las 8 regresiones individuales y determinar el valor de la varianza residual estimada $\hat\sigma_j^2$ para cada sujeto $j$ A continuación, reponga cada punto de datos individual de forma inversa a la varianza estimada para el sujeto correspondiente, y realice un regresión por mínimos cuadrados ponderados sobre los 160 puntos de datos.

Lo que usted llama la forma "hiperbólica" del error en $\hat y$ para nuevas predicciones en función de $x$ será el mismo. Proviene de la incertidumbre en la estimación de la pendiente en la regresión. El error es el más pequeño ( $\sigma/\sqrt{n}$ ) en el valor medio de $x$ y luego aumenta con esa forma simplemente porque no se sabe con qué rapidez $y$ cambios con $x$ a medida que se aleja de $\bar x$ . La regresión lineal con observaciones no correlacionadas ponderadas por sus varianzas proporciona los mejores estimadores lineales insesgados (BLUE) del coeficiente de regresión,** y con la regresión ponderada que combina todos los casos se tiene ahora un $n$ valor de 160. Así que la anchura de esa zona de incertidumbre tenderá a minimizarse.

Hay que tener cierta precaución, ya que el Página del NIST lo pone:

La mayor desventaja de los mínimos cuadrados ponderados, de la que mucha gente no es consciente, es probablemente el hecho de que la teoría que sustenta este método se basa en la suposición de que los pesos se conocen con exactitud. Por supuesto, esto casi nunca ocurre en las aplicaciones reales, por lo que hay que utilizar pesos estimados.

Y como dijo @cardinal:

Aprender una variante es difícil.

Para una distribución normal con varianza $\sigma^2$ El varianza de una estimación de varianza $\hat\sigma^2$ de $n$ observaciones es $2\sigma^4/(n-1)$ . Por lo tanto, a menos que tenga muchos puntos de datos y una razón para creer que hay diferencias sustanciales en la verdadera $\sigma_j^2$ valores entre los sujetos $j$ puede que este enfoque de ponderación no sea muy beneficioso.

Lo anterior supone que todos los sujetos tienen las mismas pendientes e interceptos para la relación entre $y$ y $x$ . Uno podría interpretar su sugerencia de "tomar la desviación estándar de los ajustes" como que usted espera verdaderas diferencias entre los sujetos en estos valores de los parámetros. En ese caso, podría obtener estimaciones de las varianzas de los interceptos y las pendientes entre los sujetos con un modelo mixto . Todavía se puede hacer una ponderación de los puntos de datos individuales.

---

*Creo que hay una forma de estimar todas las varianzas intra-sujeto y los coeficientes de regresión compartidos en un solo modelo, pero no recuerdo inmediatamente cuál es. Probablemente requeriría un enfoque iterativo o de máxima verosimilitud. Este es un enfoque simple que llega a la esencia de su pregunta.

**Véase la página enlazada sobre la regresión ponderada. Esto supone que las varianzas son conocidas.

Answer 2

Pregunta 1. ¿Cuál sería la mejor manera de calcular el error estándar entre los sujetos? ¿Sería simplemente realizar el ajuste para cada sujeto, y luego tomar la desviación estándar de los ajustes?

Opción 1: Utilizar los mínimos cuadrados ponderados. El teorema de Gauss Markov nos dice que el estimador del error estándar ponderado de la varianza inversa será el mejor estimador lineal insesgado (BLUE). Obsérvese que, aunque el modelo de la media es correcto, y en consecuencia la estimación no ponderada es insesgada, se añade la eficiencia de utilizar el estimador iterativo de mínimos cuadrados generalizados para proporcionar una mejor estimación de los residuos. Ayuda a identificar los grados de libertad apropiados para la varianza intracluster. Como referencia, he incluido la estimación en dos etapas, pero tengo problemas para identificar la corrección correcta de los grados de libertad.

Un resultado interesante en el que estoy trabajando es la idea de que los programas informáticos con opciones para la correlación intraclúster pueden proporcionar estimaciones coherentes de la heteroscedasticidad. Es decir, independientemente de si una muestra está muy intracorrelacionada o es muy variable, el efecto neto es una ponderación a la baja de esa muestra, por lo que se puede obtener el mismo error estándar óptimo en ambos casos.

Utilizando residuos no ponderados para estimar la varianza del cluster, estoy encontrando que es difícil identificar el grado de libertad apropiado para la estimación de la varianza intracluster. Estoy añadiendo mi código a continuación para que otros lo verifiquen. $n-1$ es demasiado conservador, y $n-2$ es demasiado conservador.

Opción 2: Utilizar el estimador de la varianza sándwich (consistente con la heteroscedasticidad) o el bootstrap.

Pregunta 2: ¿Cómo sería la forma del error estándar del ajuste, y cuál es la intuición que hay detrás? ¿Seguiría siendo hiperbólico? Creo que no, pero en realidad no estoy seguro.

La distribución límite de la distribución de errores sigue siendo normal siempre que la muestra "crezca más rápido" en términos de número de puntos temporales que de número de sujetos, o al menos de forma que la heteroscedasticidad a nivel de sujeto esté algo acotada. La intuición es que se trata de un resultado del teorema del límite central de Lyapunov.

require(gee)
`%covers%` <- function(x, y) x[1] < y & y < x[2]
sse.df <-function(x, df=1) {
  sum({x-mean(x)}^2)/{length(x)-df}
}
confint.gee <- function (object, parm, level = 0.95, ...) 
{
  cf <- coef(object)
  pnames <- names(cf)
  if (missing(parm)) 
    parm <- pnames
  else if (is.numeric(parm)) 
    parm <- pnames[parm]
  a <- (1 - level)/2
  a <- c(a, 1 - a)
  # pct <- format.perc(a, 3)
  pct <- paste0(formatC(100*a, format='f', digits=1), '%')
  fac <- qnorm(a)
  ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm, 
                                                             pct))
  # ses <- sqrt(diag(vcov(object)))[parm]
  ses <- sqrt(diag(object$robust.variance))[parm]
  ci[] <- cf[parm] + ses %o% fac
  ci
}

do.one <- function() {
  s1 <- 1
  s2 <- 1
  nc <- 8
  nt <- 20
  i <- rep(1:8, each=nt)

  e <- rnorm(nc, 0, s1)[i] + rnorm(nc*nt, 0, s2)

  x <- rep(seq(-3, 3, length.out = nt), times=nc)
  y <- 2*x + e

  r <- lm.fit(cbind(1,x), y)$residuals

  wls <- lm(y ~ x, weights=rep(1/tapply(r^2, i, sse.df, df=1), each=nt))
  gls <- gls(y ~ x, correlation=corCompSymm(form=~1|i))

  gee <- gee(y ~ x, id = i)

  c( ## coverage of 80% CIs
    confint(wls, parm='x', level = .8) %covers% 2,
    confint(gee, parm='x', level = .8) %covers% 2,
    confint(gls, parm='x', level= 0.8) %covers% 2,
    vcov(wls)[2,2]^.5,
    gee$robust.variance[2,2]^.5,
    vcov(gls)[2,2]^.5
  )
}

set.seed(123)
out <- replicate(500, do.one())

## 80% coverage of CIs
rowMeans(out[1:3, ])

par(mfrow=c(1,3))
hist(out[4, ], xlab='Sigma two-pass', main='')
hist(out[5, ], xlab='Sigma GEE', main='')
hist(out[6, ], xlab='Sigma GLS', main='')

Nos da una cobertura del 70% para el WLS de 2 grados de libertad y del 74% para el GEE. y del 82,54% para el GLS. Los histogramas de las estimaciones del error estándar muestran una distribución muy normal en todos los casos.