¿Es correcta la figura 3.6 de Elementos de aprendizaje estadístico?

Question

Aquí está la figura del libro de texto:

Muestra una relación decreciente entre el tamaño del subconjunto $k$ y el error cuadrático medio (ECM) de los parámetros verdaderos, $\beta$ y las estimaciones $\hat{\beta}(k)$ . Está claro que no debería ser así: añadir más variables a un modelo lineal no implica mejores estimaciones de los parámetros verdaderos. Lo que sí implica añadir más variables es un menor error de entrenamiento, es decir, una menor suma residual de cuadrados.

¿Es el $y$ -¿eje etiquetado incorrectamente? En concreto, ¿es posible que el $y$ muestra, por ejemplo Suma de cuadrados residuales en lugar de $\mathbb{E}|| \hat{\beta}(k) - \beta||^2$ ?

EDITAR :

Las discusiones y los múltiples intentos de reproducción revelaron que es probable que el eje esté etiquetado correctamente. En concreto, no se trata de RSS, ya que estará en una escala completamente diferente.

La pregunta del título sigue en pie: "¿Es correcta la figura 3.6 de ESL?". Mi intuición es que el MSE debería ser más bajo alrededor del óptimo $k$ (La respuesta de @SextusEmpiricus sugiere que es así, pero la correlación es menor). Si echamos un vistazo a la figura 3.6, vemos que el MSE sigue bajando más allá del $k=10$ .

En concreto, espero ver curvas similares a las de la Figura 3.16:

Se muestran los procedimientos adicionales debido a que se encuentra en una diferente $x$ -También utiliza un número diferente de muestras (300 frente a 100). Lo relevante aquí es la forma de, por ejemplo, "Forward stepwise" (común en ambos gráficos: naranja en el primero, negro en el segundo), que muestra un comportamiento bastante diferente en las dos figuras.

Edición final

Aquí puedes encontrar mi intento de replicar Fig3.6; el gráfico muestra diferentes niveles de correlación y número de parámetros distintos de cero. Código fuente aquí .

Answer 1

Muestra una relación decreciente entre el tamaño del subconjunto $k$ y el error cuadrático medio (ECM) de los parámetros verdaderos, $\beta$ y las estimaciones $\hat{\beta}(k)$ .

El gráfico muestra los resultados de las alternativas selección de subconjuntos métodos. El pie de foto explica el diseño experimental: hay 10 elementos de $\beta$ que son distintos de cero. Los 21 elementos restantes son cero. El método ideal de selección de subconjuntos informará correctamente de qué $\beta$ son distintos de cero y que $\beta$ son cero; en otras palabras, no se incluye ni se excluye ningún elemento incorrectamente.

Sesgo de la variable omitida se produce cuando se omiten una o más características en el proceso de generación de datos. Sesgado las estimaciones de los parámetros tienen valores esperados que no son iguales a sus valores reales (ésta es la definición de sesgo ), por lo que la elección de trazar $\mathbb{E}\|\beta -\hat{\beta}(k) \|^2$ tiene sentido. (Nótese que la definición de sesgo no coincide exactamente con este escenario experimental porque $\beta$ también es aleatorio). En otras palabras, el gráfico muestra lo incorrectas que son las estimaciones para varios $k$ para varios métodos de selección de subconjuntos. Cuando $k$ es demasiado pequeño (en este caso, cuando $k<10$ ) las estimaciones de los parámetros están sesgadas, razón por la cual el gráfico muestra grandes valores de $\mathbb{E}\|\beta -\hat{\beta}(k) \|^2$ para pequeños $k$ .

Está claro que no debería ser así: añadir más variables a un modelo lineal no implica mejores estimaciones de los parámetros reales.

Afortunadamente, eso no es lo que muestra la trama. En su lugar, la trama muestra que emplear métodos de selección de subconjuntos puede producir resultados correctos o incorrectos dependiendo de la elección de $k$ .

Sin embargo, este gráfico muestra un caso especial cuando se añaden características adicionales hace mejorar las estimaciones de los parámetros. Si se construye un modelo que presenta un sesgo de variable omitida, el modelo que incluya esas variables logrará un error de estimación de los parámetros menor porque el sesgo de variable omitida no está presente.

Lo que sí implica añadir más variables es un menor error de entrenamiento, es decir, una menor suma residual de cuadrados.

Estás confundiendo la demostración en este pasaje con una alternativa que no emplea selección de subconjuntos . En general, estimar una regresión con una base mayor disminuye el error residual medido con los datos de entrenamiento; eso no es lo que ocurre aquí.

¿Es el $y$ -¿eje etiquetado incorrectamente? En concreto, ¿es posible que el $y$ eje muestra Suma de cuadrados residuales en lugar de $\mathbb{E}\|\beta -\hat{\beta}(k) \|^2$ ?

No lo creo; la línea de razonamiento planteada en el post original no establece por sí misma que la etiqueta sea incorrecta. Los experimentos de Sextus encuentran un patrón similar; no es idéntico, pero la forma de la curva es bastante parecida.

Como apunte, creo que ya que esta trama muestra empírico resultados de un experimento, sería más claro escribir el estimador utilizado para la expectativa, según la sugerencia de Cagdas Ozgenc.

¿Es correcta la figura 3.6 de ESL?

El único definitivo forma de responder a esta pregunta es obtener el código utilizado para generar el gráfico. El código no está a disposición del público ni lo distribuyen los autores.

Sin acceso al código utilizado en el procedimiento, siempre es posible que hubiera algún error en el etiquetado del gráfico, o en la escala/ubicación de los datos o coeficientes; el hecho de que Sextus haya tenido problemas para recrear el gráfico utilizando el procedimiento descrito en el pie de foto proporciona algunas pruebas circunstanciales de que el pie de foto podría no ser del todo exacto. Se podría argumentar que estos problemas de reproducibilidad apoyan la hipótesis de que las propias etiquetas o los puntos graficados pueden ser incorrectos. Por otro lado, es posible que la descripción sea incorrecta, pero que la etiqueta en sí sea correcta.

Otra edición del libro publica una imagen diferente. Pero la existencia de una imagen diferente no implica que ninguna de las dos sea correcta.

Answer 2

añadir más variables a un modelo lineal no implica mejores estimaciones de los verdaderos parámetros

No se trata sólo de estimar variables, sino également selección variable. Cuando sólo subselecciona <10 variables, entonces inevitablemente va a cometer un error.

• Por eso el error disminuye cuando se elige un tamaño mayor para el subconjunto. Porque se están estimando más coeficientes, que probablemente sean coeficientes del modelo verdadero (en lugar de dejarlos igual a cero).

• La disminución del error va un poco más allá de $k=10$ debido a la alta correlación entre las variables.

  La mayor mejora se produce antes de k=10. Pero con $k=10$ aún no estás ahí, y vas a seleccionar ocasionalmente los coeficientes erróneos del modelo verdadero.

  Además, las variables adicionales pueden tener algún efecto regularizador .

• Tenga en cuenta que después de un cierto punto, alrededor de $k=16$ el error es arriba al añadir más variables.

Reproducción del gráfico

En el código R al final estoy tratando de reproducir el gráfico para el caso paso a paso hacia adelante. (esta es también la pregunta aquí: Recreación de la figura de Elementos de aprendizaje estadístico )

Puedo hacer que la figura se parezca

Pero, tuve que hacer algunos ajustes a la generación, utilizando $\beta \sim N(1,0.4)$ en lugar de $\beta \sim N(0,0.4)$ (y aún así no obtengo lo mismo que en la figura, que empieza en 0,95 y baja hasta 0,65, mientras que el MSE calculado con el código aquí es mucho más bajo en cambio). Aun así, la forma es cualitativamente la misma.

El error de este gráfico no se debe tanto a un sesgo: Quería dividir el error cuadrático medio en sesgo y varianza (calculando el error medio del coeficiente y la varianza del error). Sin embargo, ¡el sesgo es muy bajo! Esto se debe a la alta correlación entre los parámetros. Cuando se tiene un subconjunto con sólo 1 parámetro, entonces el parámetro seleccionado en ese subconjunto compensará los parámetros que faltan (puede hacerlo porque está altamente correlacionado). La cantidad en que los otros parámetros son demasiado bajos será más o menos la cantidad en que el parámetro seleccionado será demasiado alto. Así que, por término medio, un parámetro será tanto demasiado alto como demasiado bajo.

• El gráfico anterior está hecho con una correlación 0,15 en lugar de 0,85.
• Además, utilicé un $X$ y $\beta$ (De lo contrario, el sesgo sería igual a cero, como se explicará más adelante).

Distribución del error de estimación del parámetro

A continuación se muestra cómo el error en la estimación del parámetro $\hat\beta_1- \beta_1$ se distribuye en función del tamaño del subconjunto. Esto facilita ver por qué el cambio en el error cuadrático medio se comporta como lo hace.

Tenga en cuenta las siguientes características

• Hay un único pico para subconjuntos de pequeño tamaño. Esto se debe a que el parámetro a menudo no se incluye en el subconjunto y la estimación $\hat\beta$ será cero haciendo que el error $\hat\beta - \beta$ igual a $-\beta$ . Este pico disminuye de tamaño a medida que aumenta el tamaño del subconjunto y aumenta la probabilidad de que se incluya el parámetro.
• Existe un componente distribuido de forma más o menos gaussiana que aumenta de tamaño cuando el pico único disminuye de tamaño. Se trata del error cuando el parámetro se incluye en el subconjunto. Para subconjuntos de pequeño tamaño, el error de este componente no se centra en cero. La razón es que el parámetro necesita compensar la omisión del otro parámetro (con el que está altamente correlacionado). Esto hace que el cálculo del sesgo sea en realidad muy bajo. Lo que es alto es la varianza.

El ejemplo anterior es para $\beta$ y $X$ . Si desea cambiar el $\beta$ para cada simulación, entonces el sesgo sería cada vez diferente. Si se calcula el sesgo como $\mathbb{E}(\hat \beta - \beta)$ entonces te acercas mucho a cero.

library(MASS)

### function to do stepforward regression
### adding variables with best increase in RSS
stepforward <- function(Y,X, intercept) {
  kl <- length(X[1,])  ### number of columns
  inset <- c()
  outset <- 1:kl

  best_RSS <- sum(Y^2)
  ### outer loop increasing subset size
  for (k in 1:kl) {
    beststep_RSS <- best_RSS ### RSS to beat
    beststep_par <- 0
    ### inner looping trying all variables that can be added
    for (par in outset) {
      ### create a subset to test
      step_set <- c(inset,par)
      step_data <- data.frame(Y=Y,X=X[,step_set])
      ### perform model with subset
      if (intercept) {
        step_mod <- lm(Y ~ . + 1, data = step_data)
      }
      else {
        step_mod <- lm(Y ~ . + 0, data = step_data)
      }
      step_RSS <- sum(step_mod$residuals^2)
      ### compare if it is an improvement
      if (step_RSS <= beststep_RSS) {
        beststep_RSS <- step_RSS
        beststep_par <- par
      }
    }
    bestRSS <- beststep_RSS
    inset <- c(inset,beststep_par)
    outset[-which(outset == beststep_par)] 
  }
  return(inset)
}

get_error <- function(X = NULL, beta = NULL, intercept = 0) {
  ### 31 random X variables, standard normal 
  if (is.null(X)) {
    X <- mvrnorm(300,rep(0,31), M)
  }
  ### 10 random beta coefficients 21 zero coefficients
  if (is.null(beta)) {
    beta <- c(rnorm(10,1,0.4^0.5),rep(0,21))
  }
  ### Y with added noise
  Y <- (X %*% beta) + rnorm(300,0,6.25^0.5)

  ### get step order
  step_order <- stepforward(Y,X, intercept)

  ### error computation
  l <- 10
  error <- matrix(rep(0,31*31),31) ### this variable will store error for 31 submodel sizes
  for (l in 1:31) {

    ### subdata
    Z <- X[,step_order[1:l]]
    sub_data <- data.frame(Y=Y,Z=Z)

    ### compute model
    if (intercept) {
      sub_mod <- lm(Y ~ . + 1, data = sub_data)
    }
    else {
      sub_mod <- lm(Y ~ . + 0, data = sub_data)    
    }
    ### compute error in coefficients
    coef <- rep(0,31)
    if (intercept) {
      coef[step_order[1:l]] <- sub_mod$coefficients[-1]
    }
    else {
      coef[step_order[1:l]] <- sub_mod$coefficients[]
    }   
    error[l,] <- (coef - beta)
  }
  return(error)
}

### correlation matrix for X
M <- matrix(rep(0.15,31^2),31)
for (i in 1:31) {
  M[i,i] = 1
}

### perform 50 times the model 
set.seed(1)
X <- mvrnorm(300,rep(0,31), M)           
beta <- c(rnorm(10,1,0.4^0.5),rep(0,21)) 
nrep <- 500
me <- replicate(nrep,get_error(X,beta, intercept = 1)) ### this line uses fixed X and beta
###me <- replicate(nrep,get_error(X,beta, intercept = 1)) ### this line uses random X and fixed beta
###me <- replicate(nrep,get_error(X,beta, intercept = 1)) ### random X and beta each replicate

### storage for error statistics per coefficient and per k
mean_error <- matrix(rep(0,31^2),31)
mean_MSE <- matrix(rep(0,31^2),31)
mean_var <- matrix(rep(0,31^2),31)

### compute error statistics
### MSE, and bias + variance for each coefficient seperately
### k relates to the subset size 
### i refers to the coefficient
### averaging is done over the multiple simulations
for (i in 1:31) {
  mean_error[i,] <- sapply(1:31, FUN = function(k) mean(me[k,i,]))
  mean_MSE[i,] <- sapply(1:31, FUN = function(k) mean(me[k,i,]^2))
  mean_var[i,] <- mean_MSE[i,] - mean_error[i,]^2
}

### plotting curves
### colMeans averages over the multiple coefficients
layout(matrix(1))
plot(1:31,colMeans(mean_MSE[1:31,]), ylim = c(0,0.4), xlim = c(1,31), type = "l", lwd = 2,
     xlab = "Subset size k", ylab = "mean square error of parameters",
     xaxs = "i", yaxs = "i")
points(1:31,colMeans(mean_MSE[1:31,]), pch = 21 , col = 1, bg = 0, cex = 0.7)
lines(1:31,colMeans(mean_var[1:31,]), lty = 2)
lines(1:31,colMeans(mean_error[1:31,]^2), lty = 3)

legend(31,0.4, c("MSE", "variance component", "bias component"),
       lty = c(1,2,3), lwd = c(2,1,1), pch = c(21,NA,NA), col = 1, pt.bg = 0, xjust = 1,
       cex = 0.7)

### plotting histogram
layout(matrix(1:5,5))
par(mar = c(4,4,2,1))
xpar = 1
for (col in c(1,4,7,10,13)) {
  hist(me[col,xpar,], breaks = seq(-7,7,0.05), 
       xlim = c(-1,1), ylim = c(0,500),
       xlab = "", ylab = "",         main=paste0("error in parameter ",xpar," for subset size ",col),
       )
}

Answer 3

Intento dar una respuesta intuitiva sin comprobar e intentar reproducir el código. Ni idea de si el gráfico está mal, pero voy a explicar cómo se corresponde con mi intuición.

La pregunta tiene: "Creo que muestra una relación decreciente entre el tamaño del subconjunto k y el error cuadrático medio (ECM) de los parámetros verdaderos, β y las estimaciones β^(k). Está claro que no debería ser así: añadir más variables a un modelo lineal no implica mejores estimaciones de los parámetros verdaderos. (...) Mi intuición es que el MSE debería ser más bajo alrededor del k óptimo (en algún lugar entre 5-10 debido a las correlaciones)."

Lo que creo que está pasando es lo siguiente. Se trata de la selección de variables. El MSE de las betas estimadas debería ser menor si se seleccionan exactamente las 10 variables correctas . Debería ser sustancialmente mayor si se omite al menos una de estas variables. Tenga en cuenta que la correlación agrava este problema, ya que si se omite una de las variables beta correctas distintas de cero, su contribución se atribuirá a las que ya están en el modelo debido a la correlación. Esto hará que sus estimadores sean peores, además del hecho de que hay un error por la falta de una variable beta. $\beta$ sí mismo. Es no cierto que el efecto de la correlación es que podemos hacerlo bien con menos variables que las 10 correctas, respecto al MSE de los estimadores. Puede ser cierto para la predicción, ya que la información de una variable que falta se compensa con otras variables correlacionadas que ya están en el modelo. Pero el gráfico no trata de esto. El mismo efecto que puede ser útil para la predicción será perjudicial para la estimación, porque el efecto de las variables beta correctas no nulas que faltan se dividirá entre las que ya están en el modelo, afectando a su estimación.

Esto significa que el mínimo debería producirse en 10 sólo si siempre o casi siempre se seleccionan exactamente las 10 variables correctas. Pero esto es muy poco probable, porque la correlación en realidad hace que sea muy difícil encontrar las variables correctas. Lo más probable es que si el procedimiento selecciona 11, 12, incluso 15 variables, todavía no es demasiado improbable que se pase por alto una variable beta verdadera distinta de cero. Por otra parte, las verdaderas variables beta nulas probablemente tendrán coeficientes estimados bastante bajos, por lo que no perjudicarán tanto al MSE del estimador como una variable beta correcta no nula. En mi opinión, esto explica que el MSE del estimador sólo aumente a partir de aproximadamente $k=16$ o incluso $k=27$ más o menos por etapas. Todo esto me parece bien. Lo que muestra es cuántas variables en esta configuración deben seleccionarse para encontrar todos los no ceros verdaderos con una probabilidad suficientemente grande. 16 parece realista para mí, y también está claro que stagewise tiene un tiempo difícil en este problema, ya que necesitará muchos pasos para bajar los parámetros inicialmente sobreestimados.