Regresión Lineal

I.

Si uno hace

n <- 1000

# dos clases del mismo tamaño n
x <- c(rep(0, n), rep(1, n))

# mean(y0) = .45, mean(y1) = .55
y0 <- y1 <- rep(0, n)
y0[1:(.45 * n)] <- 1
y1[1:(.55 * n)] <- 1

# mean(y) = .5
y <- c(y0, y1)

summary(lm(y ~ x))

obtiene

Residuals:
   Min     1Q Median     3Q    Max
 -0.55  -0.45   0.00   0.45   0.55

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.45000    0.01574  28.590  < 2e-16 ***
x            0.10000    0.02226   4.492 7.44e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4977 on 1998 degrees of freedom
Multiple R-squared:   0.01,	Adjusted R-squared:  0.009505
F-statistic: 20.18 on 1 and 1998 DF,  p-value: 7.444e-06

donde quiero subrayar que la R² es del 1% o muy pequeña.

Después de publicar Una regresión de Poisson casi trivial con numpyro me riñeron por usar la identidad como función de enlace en la regresión de Poisson. Es decir, por especificarlo como

$$\lambda_t = a + b t$$

en lugar del estándar

$$\lambda_t = \exp(a + b t).$$

Hay varias cosas bastante bien conocidas y una que lo es bastante menos —y que resulta mucho más paradójica— que decir al respecto.

Antes necesito añadir que:

En los libros de texto, imperan las funciones de pérdida simétricas, como el RMSE o el MAE. Pero hay casos —muchos, de hecho, en la práctica— en que las pérdidas son asimétricas: es más oneroso pasarse, p.e., que no llegar. En esta entrada voy a analizar un ejemplo motivado por el siguiente tuit:

El resumen de lo que sigue es el siguiente:

• Voy a bajar datos de producción y consumo eléctrico de REE.
• Voy a dejar en 0 el carbón, el gas y la nuclear.
• Voy a ver por cuánto hay que multiplicar eólica y solar (dejando tal cual el resto de las renovables y cogeneraciones) para alcanzar un óptimo.

Obviamente, en el óptimo:

Esta es la tercera entrada de la serie sobre diagramas causales hiperbásicos, que, como la segunda, no se entenderá sin —y remito a— la primera que define el contexto, objetivo e hipótesis subyacentes de la serie completa. Además, sería conveniente haber leído la segunda.

Esta vez, el diagrama causal es una pequeña modificación del de la anterior:

Ahora, la variable $X$ influye sobre $Y$ por dos vías: directamente y a través de $Z$. Variables como $Z$, conocidas como mediadores son muy habituales. Uno podría pensar que, realmente, ninguna $X$ actúa directamente sobre ninguna $Y$ sino a través de una serie de mecanismos que involucran a variables intermedias $Z_1, \dots, Z_n$ que constituyen una cadena causal. Puede incluso que se desencadenen varias de estas cadenas causales que transmitan a $Y$ la potencia de $X$. Que hablemos de la influencia causal de $X$ sobre $Y$ es casi siempre una hipersimplificación de la realidad.

Esta es la segunda entrada de la serie sobre diagramas causales hiperbásicos. No se entenderá sin —y remito a— la entrada anterior que define el contexto, objetivo e hipótesis subyacentes de la serie completa.

El diagrama causal objeto de esta entrada es apenas una arista más complejo que el de la anterior:

Ahora la variable $Z$ afecta tanto a $Y$ (como en la entrada anterior) como a $X$ (esta es la novedad). Es una situación muy común en el análisis de datos. Algunos ejemplos:

Comienzo hoy una serie de cuatro entradas (¡creo!) sobre diagramas causales supersimples que involucran a tres variables aleatorias: $X$, $Y$ y $Z$. En todos los casos, estaré argumentaré alrededor de en las regresiones lineales Y ~ X e Y ~ X + Z porque nos permiten leer, interpretar y comparar rápida y familiarmente los resultados obtenidos. En particular, me interesará la estimación del efecto (causal, si se quiere) de $X$ sobre $Y$, identificable a través del coeficiente de $X$ en las regresiones. No obstante, quiero dejar claro que:

El otro día publiqué un minihilo en Twitter que terminaba con una encuesta. Proponía el siguiente problema:

1. Quiero, abusando del lenguaje, estimar el efecto de $x$ sobre $y$ usando el modelo lineal clásico $y = a_0 + a_1 x + \epsilon_1$.
2. Pero no puedo medir $x$ con precisión. Solo tengo una medida ruidosa/aproximada de $x$, $z = x + \eta$, donde $\eta$ es normal, independiente de $\epsilon_1$, etc.
3. Uso el modelo $y = b_0 + b_1 z + \epsilon_2$.

La pregunta que planteé consistía en elegir entre las siguientes tres opciones:

Hace unos meses escribí una entrada en defensa (parcial) de una regresión lineal con una R² pequeña. He vuelto a pensar sobre ella y retomo la discusión para esclarecer —sobre todo, para profanos— qué mide la R² y cómo interpretarla según el contexto.

Comienzo por un experimento físico mental. En un laboratorio se realiza un experimento para medir la relación entre dos magnitudes físicas, un efecto $latex y$ y una causa $latex x$. La teoría especifica una relación del tipo $y = a + b x$ y experimentalmente (y en condiciones de laboratorio ideales) se obtienen una serie de datos $latex (x_i, y_i)$. La relación entre ambos es de la consabida forma

Todo esto arranca con el tuit:

En conjunto, como digo, los países con Estados grandes tienden a ser poco progresivos pic.twitter.com/oeI6hkUZwd

— Juan Ramón Rallo (@juanrallo) February 1, 2021

Esa gráfica, extraída de un documento de la OCDE, creo, fue uno de los argumentos esgrimidos por JR Rallo para defender cierta postura que no viene al caso. Lo relevante para estas páginas es que fue contestado y protestado por muchos —de algunos de los cuales, dada su autoproclamada condición de divulgadores científicos, cabría esperar más— en términos exclusivamente de lo pequeño de la R².

Hay cosas tan obvias que ni se plantea la alternativa. Pero luego va R. Gomila y escribe Logistic or Linear? Estimating Causal Effects of Treatments on Binary Outcomes Using Regression Analysis que se resume en lo siguiente: cuando te interese la explicación y no la predicción, aunque tu y sea binaria, usa regresión lineal y pasa de la logística.

Nota: La sección 4.2 de An Introduction to Statistical Learning de se titula precisamente Why Not Linear Regression?

El llamado QualityScore tiene su relevancia en Google Ads. Es un indicador con valores entre 1 y 10 asignado por Google que se basa en tres variables que están descritas por ahí:

• PostClickQualityScore
• SearchPredictedCtr
• CreativeQualityScore

Se trata de variables categóricas con tres niveles: en / por encima de / por debajo de la media.

Haciendo

modelo <- lm(QualityScore ~ PostClickQualityScore +
    SearchPredictedCtr + CreativeQualityScore,
    data = tmp)

summary(modelo)

se obtiene

Call:
lm(formula = QualityScore ~ PostClickQualityScore + SearchPredictedCtr +
    CreativeQualityScore, data = tmp)

Residuals:
        Min       1Q   Median       3Q      Max
-0.25003 -0.07395  0.00775  0.06344  0.86470

Coefficients:
                                    Estimate Std. Error t value Pr(>|t|)
(Intercept)                        1.079603   0.008688   124.3   <2e-16 ***
PostClickQualityScoreAVERAGE       2.114012   0.009037   233.9   <2e-16 ***
PostClickQualityScoreABOVE_AVERAGE 3.856228   0.008448   456.5   <2e-16 ***
SearchPredictedCtrAVERAGE          1.137396   0.003284   346.4   <2e-16 ***
SearchPredictedCtrABOVE_AVERAGE    3.055694   0.004707   649.2   <2e-16 ***
CreativeQualityScoreAVERAGE        0.999580   0.004274   233.9   <2e-16 ***
CreativeQualityScoreABOVE_AVERAGE  2.000725   0.003862   518.1   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1574 on 11426 degrees of freedom
Multiple R-squared:  0.9915,	Adjusted R-squared:  0.9915
F-statistic: 2.212e+05 on 6 and 11426 DF,  p-value: < 2.2e-16

Que no merece mayor explicación. Creo.

En esta entrada voy a crear un conjunto de datos donde dos variables tienen una correlación muy alta, ajustar un modelo de regresión y obtener la siguiente representación de la distribución a posteriori de los coeficientes,

donde se aprecia el efecto de la correlación entre x1 y x2.

El código,

library(mvtnorm)
library(rstan)
library(psych)

n <- 100
corr_coef <- .9

x <- rmvnorm(n, c(0, 0),
  sigma = matrix(c(1, corr_coef, corr_coef, 1), 2, 2))
plot(x)

x1 <- x[,1]
x2 <- x[,2]
x3 <- runif(n) - 0.5

y <- 1 + .4 * x1 - .2 * x2 + .1 * x3 + rnorm(n, 0, .1)

summary(lm(y ~ x1 + x2 + x3))

stan_code <- "
data {
  int N;
  vector[N] y;
  vector[N] x1;
  vector[N] x2;
  vector[N] x3;
}
parameters {
  real a;
  real a1;
  real a2;
  real a3;
  real sigma;
}

model {
  a ~ cauchy(0,10);
  a1 ~ cauchy(0,2.5);
  a2 ~ cauchy(0,2.5);
  a3 ~ cauchy(0,2.5);

  y ~ normal(a + a1 * x1 + a2 * x2 + a3 * x3, sigma);
}"


datos_stan <- list(
    N = n,
    y = y,
    x1 = x1,
    x2 = x2,
    x3 = x3
)

fit2 <- stan(model_code = stan_code,
              data = datos_stan,
              iter = 10000, warmup = 2000,
              chains = 2, thin = 4)

res <- as.data.frame(fit2)
pairs.panels(res[, c("a", "a1", "a2", "a3", "sigma")])

Liberados del estrecho ámbito de nuestra original mentira sugerente gracias a la relación que descubrimos entre residuos y gradientes cuando las pérdidas son cuadráticas podemos adentrarnos en ámbitos más extensos.

Lo que discutimos del gradiente tiene una interpretación fácilmente inteligible en el caso de pérdidas cuadráticas. Pero ni la pérdida de interpretabilidad nos impide extender el razonamiento de la entrada anterior a funciones de pérdida distintas de la cuadrática siempre que podamos calcular un gradiente.