18 de noviembre de 2016

Modelos lineales 2: Regresiones múltiples

Regresión múltiple

  • Modelos lineales con varias variables predictoras.
  • Pueden ser categóricas, numéricas, o ambas.
  • Complejidad añadida: interacciones

Una variable categórica y una numérica

  • Modelo lineal con:
    • Un intercepto:
      • Nivel de referencia de la variable categórica
      • Valor para variable numérica igual a 0
    • Un coeficiente para la variable numérica (pendiente)
    • Un coeficiente para cada nivel (extra) de la variale categóŕica.

Modelo lineal

\(y= a_2x_2 + a_1x_1 + b\)
  • Variable numérica: \(x_2\)
  • Variable categórica de dos niveles:
    • \(x_1 = 1\): nivel 1
    • \(x_2 = 0\): nivel de referencia

Dummy coding

\(y= a_2x_2 + a_1x_1 + b\)
  • Entonces
    • si \(x_1=0\) y \(x_2=0\), \(y=b\)
      • \(b\) es el tamaño de firma para un hombre de narcicismo 0.
    • si \(x_1=1\) y \(x_2=0\), \(y=a_1+b\)
      • \(a_1+b\) es el tamaño de firma para una mujer de narcicismo 0.
      • \(a_1\) representa la la diferencia entre los tamaños de firmas entre hombres y mujeres de narcicismo 0.

Dummy coding

\(y= a_2x_2 + a_1x_1 + b\)
  • Y…
    • si \(x_2=w\), y \(x_1=0\) \(y=a_2w+b\)
      • \(a_2w+b\) es el tamaño de firma para un hombre de narcicismo \(w\).
      • \(a_2\) es la variación en el tamaño de la firma según el narcicismo.
    • si \(x_2=w\), y \(x_1=1\) \(y=a_2w+a_1+b\)
      • \(a_2w+a_1+b\) es el tamaño de firma para una mujer de narcicismo \(w\).
      • \(a_2\) es la variación en el tamaño de la firma según el narcicismo.

Firmas, narcisismo y sexo

sa.lm=lm(Sign.Area~NPI16TOTAL+Sex,data=datosLimpios)
summary(sa.lm)
## 
## Call:
## lm(formula = Sign.Area ~ NPI16TOTAL + Sex, data = datosLimpios)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.984  -5.743  -1.614   4.576  25.202 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  24.3982     1.0305  23.676   <2e-16 ***
## NPI16TOTAL    0.2086     0.2001   1.042    0.298    
## SexF         -2.0397     0.9472  -2.153    0.032 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.29 on 321 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.02042,    Adjusted R-squared:  0.01432 
## F-statistic: 3.346 on 2 and 321 DF,  p-value: 0.03647

Firmas, narcisismo y sexo

ggplot(datosLimpios,aes(x=NPI16TOTAL,y=Sign.Area))+
  geom_point(aes(col=Sex))+geom_smooth(method="lm",se=F)

Firmas, narcisismo y sexo

ggplot(datosLimpios,aes(x=NPI16TOTAL,y=Sign.Area,col=Sex))+
  geom_point()+geom_smooth(method="lm",se=F)

Interacciones

  • Interacción entre variables
    • El efecto de una variable predictora depende de los valores (niveles) de otra.
  • Las interacciones no están incluidas en un modelo lineal como el que vimos.
    • Es necesario especificar la interacción:
      • Como un producto de las variables que interactúan.
    • Cada interacción debe incluir su propio coeficiente.

Modelo lineal con interacción

\(y= a_{21}x_2x_1 + a_2x_2 + a_1x_1 + b\)
  • Variable categórica de dos niveles:
    • \(x_1 = 1\): nivel 1
    • \(x_1 = 0\): nivel de referencia
  • Variable numérica: \(x_2\)
    • Efecto para el nivel de referencia
  • Interacción: \(x_2x_1\) (doble)
  • Diferencia de efecto para el nivel 1

Dummy coding interacción

  • Efecto del narcicismo es diferente según el sexo:
    • Una pendiente para sexo F.
    • Otra para sexo M.

Modelo con interacción

Nivel sexo \(x_2x_1\) \(x_2\) \(x_1\) intercepto
M 0 narcicismo 0 si
F narcicismo narcicismo 1 si

Dummy coding interacción (por defecto)

\(y= a_{21}x_2x_1 + a_2x_2 + a_1x_1 + b\)
  • Entonces
    • si \(x_2=w\), y \(x_1=0\) \(y=a_2w+b\)
      • \(a_2w+b\) es el tamaño de firma para un hombre de narcicismo \(w\).
      • \(a_2\) es la variación en el tamaño de la firma según el narcicismo, para hombres.
    • si \(x_2=w\), y \(x_1=1\) \(y=a_3w+a_2w+a_1+b\)
      • \(a_{21}w+a_2w+a_1+b\) es el tamaño de firma para una mujer de narcicismo \(w\).
      • \(a_{21}\) es la diferencia del efecto del narcicismo sobre el tamaño de la firma, entre hombres y mujeres.

Firmas, narcisismo y sexo

sa.lmi=lm(Sign.Area~NPI16TOTAL*Sex,data=datosLimpios)
summary(sa.lmi)
## 
## Call:
## lm(formula = Sign.Area ~ NPI16TOTAL * Sex, data = datosLimpios)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.954  -5.892  -1.577   4.846  24.631 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      26.4911     1.2826  20.654  < 2e-16 ***
## NPI16TOTAL       -0.3528     0.2876  -1.227 0.220789    
## SexF             -5.6133     1.6247  -3.455 0.000624 ***
## NPI16TOTAL:SexF   1.0694     0.3969   2.694 0.007428 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.21 on 320 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.04215,    Adjusted R-squared:  0.03317 
## F-statistic: 4.693 on 3 and 320 DF,  p-value: 0.003189

Firmas, narcisismo y sexo

## 
## Call:
## lm(formula = Sign.Area ~ NPI16TOTAL * Sex, data = datosLimpios)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.954  -5.892  -1.577   4.846  24.631 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      26.4911     1.2826  20.654  < 2e-16 ***
## NPI16TOTAL       -0.3528     0.2876  -1.227 0.220789    
## SexF             -5.6133     1.6247  -3.455 0.000624 ***
## NPI16TOTAL:SexF   1.0694     0.3969   2.694 0.007428 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.21 on 320 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.04215,    Adjusted R-squared:  0.03317 
## F-statistic: 4.693 on 3 and 320 DF,  p-value: 0.003189
  • No parece haber efecto de narcicismo en hombres.
  • El efecto del narcicismo difiere entre hombres y mujeres.

Firmas, narcisismo y sexo

  • Dummy coding alternativo
sa.lmi2=lm(Sign.Area~Sex+NPI16TOTAL:Sex,data=datosLimpios)
summary(sa.lmi2)
## 
## Call:
## lm(formula = Sign.Area ~ Sex + NPI16TOTAL:Sex, data = datosLimpios)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.954  -5.892  -1.577   4.846  24.631 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      26.4911     1.2826  20.654  < 2e-16 ***
## SexF             -5.6133     1.6247  -3.455 0.000624 ***
## SexM:NPI16TOTAL  -0.3528     0.2876  -1.227 0.220789    
## SexF:NPI16TOTAL   0.7166     0.2736   2.619 0.009229 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.21 on 320 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.04215,    Adjusted R-squared:  0.03317 
## F-statistic: 4.693 on 3 and 320 DF,  p-value: 0.003189
  • No parece haber efecto de narcicismo en hombres.
  • Sí parece haber efecto de narcicismo en mujeres.

Firmas, narcisismo y sexo

sa.lmi3=lm(Sign.Area~NPI16TOTAL+NPI16TOTAL:Sex,data=datosLimpios)
summary(sa.lmi3)
## 
## Call:
## lm(formula = Sign.Area ~ NPI16TOTAL + NPI16TOTAL:Sex, data = datosLimpios)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.767  -5.757  -1.559   4.116  26.514 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     22.99259    0.80056  28.721   <2e-16 ***
## NPI16TOTAL       0.30290    0.21973   1.379    0.169    
## NPI16TOTAL:SexF -0.05026    0.23308  -0.216    0.829    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.349 on 321 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.006414,   Adjusted R-squared:  0.0002233 
## F-statistic: 1.036 on 2 and 321 DF,  p-value: 0.356
  • No parece haber efecto de narcicismo en hombres.
  • Sí parece haber efecto de narcicismo en mujeres.

Pruebas F

Bondad de ajuste

  • Los residuos son variabilidad que el modelo no explica.
  • Se asume que son normales, de media 0.
    • Se utilizan QQplots para comprobar eso.
  • Se asume que la variabilidad de los residuos es siempre igual (para todos los niveles o valores de la variable independiente).
  • Un buen modelo tiene residuos pequeños, simétricos en 0, de variabilidad constante.
  • Un modelo es robusto si sus parámetros no varían mucho al cambiar levemente las observaciones.
    • Ver Distancia de Cook y Leverage.
  • Un modelo es mejor que otro, si sus residuos son significativamente menores.

Ejemplo: Residuos

plot(sa.lmi2)

Ejemplo: Residuos QQnorm plot

plot(sa.lmi2)

Ejemplo: Residuos estandarizados

plot(sa.lmi2)

Ejemplo: Residuos y apalancamiento (leverage)

plot(sa.lmi2)

Residuos al cuadrado

  • Como vimos, el ajuste lineal minimiza el error cuadrático: esto es, la suma de los residuos al cuadrado.
  • Se usa la notación \(SS\) para denotar sumas de cuadrados.
    • \(RSS\) denota las sumas de residuos al cuadrado (Residual Sum of Squares).
  • Siempre que hay un modelo, hay un \(RSS\).
  • El \(RSS\) dividido los grados de libertad, es una especie de promedio de los residuos al cuadrado.
    • A veces se le denomina con \(MRSS\), por "Mean Residual Sum of Squares" \(MRSS\) o \(MRSE\), etc.
  • La varianza es un caso especial de \(MRSS\).
    • Es la suma de los cuadrados de los "residuos" a la media, dividido los grados de libertad.

Grados de libertad

  • Tengo 10 datos. Calculo la media.
  • ¿Puedo cambiar todos los 10 datos, y volver a tener la misma media?
  • Si, pero:
    • Tengo "libertad total" para poner el elegir cualquier valor en el primer dato.
    • Tengo "libertad total" para poner el elegir cualquier valor en el segundo dato.
    • Tengo "libertad total" para poner el elegir cualquier valor en el tercer dato.
    • Etc…
    • Cuando llego al décimo dato, ya no puedo elegir cualquier valor. Este dato está determinado por los otros 9.

Grados de Libertad

  • Cuando se calcula la varianza, conocemos la media.
  • La varianza depende de las \(n\) observaciones, y además de la media.
  • Por eso, el cálculo de la varianza tiene \(n-1\) grados de libertad.
  • En general, un modelo lineal tiene:
    • \(n\) observaciones.
    • \(p\) parámetros.
    • \(n-p\) grados de libertad.
  • Es decir, dados los valores de los \(p\) parámetros, puedo elegir libremente \(n-p\) valores para las observaciones.

Prueba F

  • Con la Prueba F de Snédecor se puede comparar el ajuste de dos modelos lineales diferentes.
  • Típicamente la conocemos en el contexto de las tablas de ANOVA.
  • El estadístico \(F\) es un cociente, que bajo \(H_0\) sigue la distribución F.

  • El cociente F, en términos genéricos, se puede escribir como:

  • \(F=\frac{ \mathrm{Varianza}\quad \mathrm{explicada}}{\mathrm{Varianza}\quad \mathrm{no}\quad \mathrm{explicada}}\)

  • \(F>1\) : Más varianza explicada que no explicada (modelo ajusta bien)

  • \(F=1\) o \(F<1\) : Varianza explicada es menor que la no explicada (modelo no ajusta bien)

Un ejemplo sencillo

  • Varible categórica, dos niveles, con medias diferentes.
    • Hipótesis nula: una única media general (las medias no son diferentes para cada nivel).
  • \(SS_{total} = SS_{explicada} + SS_{noexplicada}\)
    • \(SS_{explicada}: SS_{total} - SS_{no explicada}\)
  • \(SS_{no explicada}= RSS\)

Partición de varianza

  • \(SS_{total} = (2-5)^2 +(4-5)^2 + (6-5)^2 + (8-5)^2 = 20\)
  • \(SS_{explicada} = (3-5)^2 +(3-5)^2 + (7-5)^2 + (7-5)^2 = 16\)
  • \(SS_{noExplicada} = (2-3)^2 +(4-3)^2 + (6-7)^2 + (8-7)^2 = 4\)

Valor de F

  • \(F=\frac{ \mathrm{Varianza}\quad \mathrm{explicada}}{\mathrm{Varianza}\quad \mathrm{no}\quad \mathrm{explicada}}\)
  • \(F=\frac{SS_{explicada} / \mathrm{gln}}{SS_{noExplicada} / \mathrm{gld}}\)
  • gln: grados de libertad del numerador: \(p-1 = 2-1 = 1\)
  • gld: grados de libertad del denominador: \(n-p = 4-2=1\)
  • \(F=\frac{16 / 1}{4 / 2} = 8\)
  • \(p=0.1056\)

Calculando \(F\)s y ANOVAs

anova(sa.lm)
## Analysis of Variance Table
## 
## Response: Sign.Area
##             Df  Sum Sq Mean Sq F value  Pr(>F)  
## NPI16TOTAL   1   141.2  141.19  2.0546 0.15272  
## Sex          1   318.6  318.61  4.6364 0.03204 *
## Residuals  321 22058.5   68.72                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Calculando \(F\)s y ANOVAs

summary(sa.lm)
## 
## Call:
## lm(formula = Sign.Area ~ NPI16TOTAL + Sex, data = datosLimpios)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.984  -5.743  -1.614   4.576  25.202 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  24.3982     1.0305  23.676   <2e-16 ***
## NPI16TOTAL    0.2086     0.2001   1.042    0.298    
## SexF         -2.0397     0.9472  -2.153    0.032 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.29 on 321 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.02042,    Adjusted R-squared:  0.01432 
## F-statistic: 3.346 on 2 and 321 DF,  p-value: 0.03647

Calculando \(F\)s y ANOVAs

library(car)
Anova(sa.lm,type=3)
## Anova Table (Type III tests)
## 
## Response: Sign.Area
##             Sum Sq  Df  F value  Pr(>F)    
## (Intercept)  38520   1 560.5439 < 2e-16 ***
## NPI16TOTAL      75   1   1.0862 0.29811    
## Sex            319   1   4.6364 0.03204 *  
## Residuals    22059 321                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Calculando \(F\)s y ANOVAs

library(car)
Anova(sa.lm,type=2)
## Anova Table (Type II tests)
## 
## Response: Sign.Area
##             Sum Sq  Df F value  Pr(>F)  
## NPI16TOTAL    74.6   1  1.0862 0.29811  
## Sex          318.6   1  4.6364 0.03204 *
## Residuals  22058.5 321                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F para un modelo lineal

  • Es una \(F\) que compara dos modelos
  • Un modelo con \(p_1\) parámetros (modelo 1)
  • Un modelo con \(p_2\) parámetros (modelo 2)
  • En el ejemplo: \(p_1=1\), una única media general: \(y=b\)
  • En el ejemplo: \(p_2=2\), una media para cada condición: \(y=ax+b\)
  • El modelo 2, con más parámetros siempre tiene menor \(RSS\)
    • La pregunta es si la diferencia es significativa

F para un modelo lineal

  • Recordemos que \(SS_{explicada}: SS_{total} - SS_{no explicada}\)
  • La \(\mathrm{Varianza}\quad \mathrm{explicada}\) en este contexto es la diferencia:
    • \(\frac{RSS_1 - RSS_2}{p_2-p_1}\)
    • esto es, la diferencia de la suma de residuos cuadrados
    • si el ajuste del modelo 2 es mucho mejor, esto será grande
  • La \(\mathrm{Varianza}\quad \mathrm{no} \quad \mathrm{explicada}\) es:
    • \(\frac{RSS_2}{p_2}\)
    • o sea, la varianza que el mejor modelo no puede explicar.

F siempre compara modelos

  • En el ejemplo más simple:
    • Modelo 1: \(y=b\)
      • Una única media general. (1 parámetro)
      • \(RSS_1\) = Varianza de los datos (no explico nada de varianza).
    • Modelo 2: \(y=ax+b\)
      • Una media por nivel. (2 parámetros)
      • \(RSS_2\) = Varianza no explicada por este modelo
    • $  frac{(RSS_1 - RSS_2)/1}{RSS_2/2}$

F siempre compara modelos

  • Cada vez que vean una \(F\)
    • Se contrastan dos modelos
    • Con distinto número de parámetros
  • Si es una F para una variable numérica:
    • Se compara un modelo con vs. sin esa variable/parámetro.
  • Si es una F para una variable categórica:
    • Se compararan modelos con vs. sin las variables/parámetros relevantes.
  • Si es una interacción
    • Se compararan modelos con vs. sin las variables/parámetros relevantes.

La F en summary(modelo.lineal)

summary(sa.lm)
## 
## Call:
## lm(formula = Sign.Area ~ NPI16TOTAL + Sex, data = datosLimpios)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.984  -5.743  -1.614   4.576  25.202 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  24.3982     1.0305  23.676   <2e-16 ***
## NPI16TOTAL    0.2086     0.2001   1.042    0.298    
## SexF         -2.0397     0.9472  -2.153    0.032 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.29 on 321 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.02042,    Adjusted R-squared:  0.01432 
## F-statistic: 3.346 on 2 and 321 DF,  p-value: 0.03647
  • F global. Contraste entre el modelo completo y un modelo \(y=b\) (media general).

ANOVAs tipo II y III

  • El ANOVA genera
    • Una F para cada variable
      • Contraste entre modelos con/sin coeficientes relevantes.
    • Una F para cada interacción
      • Contraste entre modelos con/sin coeficientes relevantes.

El orden importa (tipo II)

## Analysis of Variance Table
## 
## Response: Sign.Area
##                 Df  Sum Sq Mean Sq F value   Pr(>F)   
## NPI16TOTAL       1   141.2  141.19  2.0947 0.148792   
## Sex              1   318.6  318.61  4.7268 0.030428 * 
## NPI16TOTAL:Sex   1   489.3  489.27  7.2587 0.007428 **
## Residuals      320 21569.3   67.40                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Variance Table
## 
## Response: Sign.Area
##                 Df  Sum Sq Mean Sq F value   Pr(>F)   
## Sex              1   385.2  385.16  5.7142 0.017406 * 
## NPI16TOTAL       1    74.6   74.64  1.1073 0.293454   
## Sex:NPI16TOTAL   1   489.3  489.27  7.2587 0.007428 **
## Residuals      320 21569.3   67.40                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA tipo II

  • Tipo secuencial
    • Introduce variables de a una en el modelo
    • Contrasta cada modelo con el anterior
      • Modelo: inter.
      • var 1 + inter. vs inter.
      • var 2 + var 1 + inter. vs var 1 + inter.
      • etc.
    • El orden con que se escriben las variables en el modelo lineal importa
  • Función anova()

ANOVA tipo III

  • Tipo ajustado
    • "Quita" variables de a una en el modelo
    • Contrasta los modelos resultantes contra el modelo completo
      • Modelo completo
      • Modelo - var1 vs Modelo completo
      • Modelo - var2 vs Modelo completo
      • etc.
    • El orden con que se escriben las variables en el modelo lineal no importa
  • Función Anova(,type=3)

En síntesis

  • TODO (o casi) son modelos lineales
  • TODOs (casi) se escriben igual, no importa qué tipos de variables predictoras
  • Cada variable numérica: 1 parámetro extra
  • Cada variable categórica: 1 parámetro extra por cada nivel extra
  • Interacciones: 1 parámetro por cada combinación de niveles (extra)
    • Explosión combinatoria de parámetros
    • Pruebas F al rescate!
  • Las pruebas F siempre comparan 2 modelos
  • ¿Qué tanto mejor ajusta el modelo con más parámetros?

Próximas clases

  • Jueves: Ejercicios…
  • Lunes: "Pantallazo" de "medidas repetidas" (efectos mixtos)