Regresión lineal

Regresión lineal en R

Limpiamos la memoria con rm() y gc(), cargamos los directorios y cargamos las librerías y la base que vamos a utilizar.

rm(list=ls())
gc()

          used (Mb) gc trigger   (Mb)   max used   (Mb)
Ncells 2132921  114    3205452  171.2    3205452  171.2
Vcells 3526159   27  263950304 2013.8 1006890509 7682.0

dir <- paste0(dirname(rstudioapi::getActiveDocumentContext()$path),"/")
bases.dir      <-  paste0(dirname(dir),"/Fuentes/")
resultados.dir <- paste0(dirname(dir),"/Resultados/")
library(tidyverse, warn.conflicts = FALSE)
individual.416 <- read.table(paste0(bases.dir, "usu_individual_t416.txt"), sep=";", dec=",", header = TRUE, fill = TRUE)

El objetivo será hacer inferencia respecto del ingreso de la ocupación principal P21, para ello selecionamos las columnas con las que trabajaremos: variable dependiente y algunas posibles covariables:

Variables a utilizar:

• P21: MONTO DE INGRESO DE LA OCUPACIÓN PRINCIPAL
• NIVEL_ED
• ESTADO
• CAT_OCUP
• CH04: Sexo
• CH06: Edad
• PP04C99: Tamaño del establecimiento, reducido
• PP07H: ¿Por ese trabajo tiene descuento jubilatorio?

creamos un vector para seleccionar estas varibales, y un vectr que especifique qué variables son categóricas, para convertirlas en factores. Es necesario convertir las variables categóricas a la clase factor, para que el algorítmo de regresión lineal las reconozca como lo que son

La función lapply() aplica la función factor() a toda la base[factores]

var <- c('P21','NIVEL_ED','ESTADO','CAT_OCUP', 'CH04','CH06','PP04C99','PP07H')
base <- individual.416 %>% 
  select(var)
factores <- c('NIVEL_ED','ESTADO','CAT_OCUP', 'CH04','PP04C99','PP07H')
base[factores] <- lapply(base[factores], factor)
# Nos quedamos sólo con los valores positivos de la variable explicada
base <- base %>% filter(P21 >0)

Correlación y test de correlación

• Para calcular un test de correlacón de Pearson utilizamos el comando cor().
• Si queremos calcular un test de hipótesis sobre la correlación, o un interavalo de confianza, utilizamos cor.test()
• Si no se cumplen los supestos de Pearson (Normalidad, independencia, linealidad en la relación), podemos utilizar el test de correlación no paramétrico de Spearman o kendall con el parámetro method =
• si queremos que el test sea a una cola, lo podemos especificar con el parámetro alternative =
• La significativiadd por default es 0.95, pero se puede modificar con el parámetro conf.level =

cor(base$P21, base$CH06)

[1] 0.1460676

cor.test(base$P21, base$CH06)

    Pearson's product-moment correlation

data:  base$P21 and base$CH06
t = 20.09, df = 18514, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1319413 0.1601345
sample estimates:
      cor 
0.1460676 

cor.test(base$P21, base$CH06,method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  base$P21 and base$CH06
S = 8.9399e+11, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.1550252 

Regresión lineal

con el comando lm() podemos ajustar el modelo. Algunos parametros importantes de esta función son:

• formula: definimos el modelo como _ Y ~ X _.
  • regresion multiple: y ~ \(X_1 + X_2 + ... X_n\).
  • regresión polinómica: y ~ poly(x = X, degree = k)
  • interacción de variables: y ~ \(X_1*X_2\). Si sólo queremos el término de interacción: \(X_1:X_2\)
• data : especificamos el dataset a utilizar
• subset : si dividimos la muestra en training y testing, podemos indicar el subconjunto de entrenamiento con un vector que indique sus números de fila

Regresión lineal simple

Comenzamos regresando a P21 contra la única variable continua que seleccionamos, la edad.

El modelo propuesto resulta

\[ E(P_{21}/ch_{06})= \beta_0 + \beta_1 ch_{06} \]

modelo.ajustado <- lm(formula = P21 ~ CH06, data = base)

La función lm() genera un objeto que contiene varias propiedades de la regresión que ajustamos.

con summary() podemos acceder los principales resultados

summary(modelo.ajustado)

Call:
lm(formula = P21 ~ CH06, data = base)

Residuals:
   Min     1Q Median     3Q    Max 
-15780  -5630  -1599   3351 386290 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6744.583    222.854   30.27   <2e-16 ***
CH06         105.533      5.253   20.09   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9105 on 18514 degrees of freedom
Multiple R-squared:  0.02134,   Adjusted R-squared:  0.02128 
F-statistic: 403.6 on 1 and 18514 DF,  p-value: < 2.2e-16

Para revisar los supuestos del modelo, así como la posible presencia de outliers, podemos revisar unos gráficos que se generan vía el comando plot().

plot(modelo.ajustado)

• Residuals ~ Fitted values: Aquí buscamos que no exista ningúna estructura en los datos. Si la presentan, según la forma que tome, podría ser un indicador de heterocedasticidad o no linealidad en la relación
• QQ plot: Muestra si los datos se ajustan al supuesto de normalidad. Si los puntos se ubican sobre la linea puntada, esto indica que tienen una distribución normal. En nuestro caso, podemos ver que se ve que se alejan en el extremo derecho de la distribución.
• \(\sqrt{residuos \space estandarizados}\) ~ Fitted Values: Es similar al primer gráfico, con otra escala.
• Leverage ~ residuos estandarizados: Nos permite encontrar observaciones influyentes. En particular, aquellas que queden fuera de los márgenes definidos por la distancia de cook deberían ser revisadas.

Regresión lineal múltiple

\[ E(P_{21}/X)= \beta_0 + \beta_1 ch_{06} + \beta_2 ch_{06}^2 \]

modelo.ajustado <- lm(formula = P21 ~ poly(x = CH06,degree = 2) , data = base)
summary(modelo.ajustado)

Call:
lm(formula = P21 ~ poly(x = CH06, degree = 2), data = base)

Residuals:
   Min     1Q Median     3Q    Max 
-12489  -5442  -1673   3145 388558 

Coefficients:
                              Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   11015.19      66.44  165.79   <2e-16 ***
poly(x = CH06, degree = 2)1  182928.34    9040.73   20.23   <2e-16 ***
poly(x = CH06, degree = 2)2 -147552.65    9040.73  -16.32   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9041 on 18513 degrees of freedom
Multiple R-squared:  0.03522,   Adjusted R-squared:  0.03511 
F-statistic: 337.9 on 2 and 18513 DF,  p-value: < 2.2e-16

Si queremos elegir el mejor modelo polinomico entre los de la forma

\[ E(P_{21}/X)= \beta_0 + \beta_1 ch_{06} + \beta_2 ch_{06}^2 + ... + \beta_n ch_{06}^n \]

podemos hacer uso de los loops, para ajustar varios modelos polinómicos de distinto grado, y luego comparar sus r cuadrados ajustados.

r.ajustado <- c()
for (i in 1:10) {
  modelo.ajustado <- lm(formula = P21 ~ poly(x = CH06,degree = i) , data = base)
  r.ajustado[i] <- summary(modelo.ajustado)$adj.r.squared  
}
plot(r.ajustado, type = "b")

Regresión lineal con variables categóricas

\[ E(P_{21}/X)= \beta_0 + \beta_1 ch_{06} + \beta_2 cat.ocup \]

CAT_OCUP: 1 = Patrón 2 = Cuenta propia 3 = Obrero o empleado

table(base$CAT_OCUP)

    0     1     2     3     4     9 
    0   548  3510 14458     0     0 

Para utilizar dummies, no es necesario que las creemos, simplemente agregamos la variable categórica, y la función se encarga de crear las variables.

modelo.ajustado <- lm(formula = P21 ~ CH06 + CAT_OCUP, data = base)
summary(modelo.ajustado)

Call:
lm(formula = P21 ~ CH06 + CAT_OCUP, data = base)

Residuals:
   Min     1Q Median     3Q    Max 
-19731  -5185  -1593   3060 380342 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11791.53     454.04   25.97   <2e-16 ***
CH06          119.20       5.21   22.88   <2e-16 ***
CAT_OCUP2   -9464.41     408.12  -23.19   <2e-16 ***
CAT_OCUP3   -4873.89     388.90  -12.53   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8875 on 18512 degrees of freedom
Multiple R-squared:  0.07022,   Adjusted R-squared:  0.07007 
F-statistic:   466 on 3 and 18512 DF,  p-value: < 2.2e-16

Best subset selection

Si queremos hacer Best subset selection, podemos utilizar el paquete leaps, aunque no es el único. El comando regsubsets() nos permite hacer selección de modelos exahustiva, así como forward o backward stepwise, entre otros métodos.

library(leaps)
names(base)

[1] "P21"      "NIVEL_ED" "ESTADO"   "CAT_OCUP" "CH04"     "CH06"     "PP04C99"  "PP07H"   

regsubsets.out <-
    regsubsets(P21 ~ NIVEL_ED + ESTADO + CAT_OCUP + CH04 + CH06 + PP04C99 + PP07H,
               data = base,
               nbest = 1,       # 1 best model for each number of predictors
               nvmax = NULL,    # NULL for no limit on number of variables
               force.in = NULL, force.out = NULL,
               method = "exhaustive")

7  linear dependencies found

Reordering variables and trying again:

summary(regsubsets.out)

Subset selection object
Call: regsubsets.formula(P21 ~ NIVEL_ED + ESTADO + CAT_OCUP + CH04 + 
    CH06 + PP04C99 + PP07H, data = base, nbest = 1, nvmax = NULL, 
    force.in = NULL, force.out = NULL, method = "exhaustive")
23 Variables  (and intercept)
          Forced in Forced out
NIVEL_ED2     FALSE      FALSE
NIVEL_ED3     FALSE      FALSE
NIVEL_ED4     FALSE      FALSE
NIVEL_ED5     FALSE      FALSE
NIVEL_ED6     FALSE      FALSE
NIVEL_ED7     FALSE      FALSE
CAT_OCUP1     FALSE      FALSE
CAT_OCUP2     FALSE      FALSE
CH042         FALSE      FALSE
CH06          FALSE      FALSE
PP04C991      FALSE      FALSE
PP04C992      FALSE      FALSE
PP04C993      FALSE      FALSE
PP04C999      FALSE      FALSE
PP07H1        FALSE      FALSE
PP07H2        FALSE      FALSE
ESTADO1       FALSE      FALSE
ESTADO2       FALSE      FALSE
ESTADO3       FALSE      FALSE
ESTADO4       FALSE      FALSE
CAT_OCUP3     FALSE      FALSE
CAT_OCUP4     FALSE      FALSE
CAT_OCUP9     FALSE      FALSE
1 subsets of each size up to 16
Selection Algorithm: exhaustive
          NIVEL_ED2 NIVEL_ED3 NIVEL_ED4 NIVEL_ED5 NIVEL_ED6 NIVEL_ED7 ESTADO1 ESTADO2 ESTADO3 ESTADO4 CAT_OCUP1
1  ( 1 )  " "       " "       " "       " "       " "       " "       " "     " "     " "     " "     " "      
2  ( 1 )  " "       " "       " "       " "       " "       " "       " "     " "     " "     " "     " "      
3  ( 1 )  " "       " "       " "       " "       "*"       " "       " "     " "     " "     " "     " "      
4  ( 1 )  " "       " "       " "       " "       "*"       " "       " "     " "     " "     " "     "*"      
5  ( 1 )  " "       " "       " "       " "       "*"       " "       " "     " "     " "     " "     "*"      
6  ( 1 )  "*"       " "       " "       " "       "*"       " "       " "     " "     " "     " "     "*"      
7  ( 1 )  " "       " "       "*"       "*"       "*"       " "       " "     " "     " "     " "     "*"      
8  ( 1 )  " "       "*"       "*"       "*"       "*"       " "       " "     " "     " "     " "     "*"      
9  ( 1 )  "*"       "*"       "*"       "*"       "*"       " "       " "     " "     " "     " "     "*"      
10  ( 1 ) "*"       "*"       "*"       "*"       "*"       " "       " "     " "     " "     " "     "*"      
11  ( 1 ) "*"       "*"       "*"       "*"       "*"       " "       " "     " "     " "     " "     "*"      
12  ( 1 ) "*"       "*"       "*"       "*"       "*"       " "       " "     " "     " "     " "     "*"      
13  ( 1 ) "*"       "*"       "*"       "*"       "*"       " "       " "     " "     " "     " "     "*"      
14  ( 1 ) "*"       "*"       "*"       "*"       "*"       " "       " "     " "     " "     " "     "*"      
15  ( 1 ) "*"       "*"       "*"       "*"       "*"       " "       " "     " "     " "     " "     "*"      
16  ( 1 ) "*"       "*"       "*"       "*"       "*"       "*"       " "     " "     " "     " "     "*"      
          CAT_OCUP2 CAT_OCUP3 CAT_OCUP4 CAT_OCUP9 CH042 CH06 PP04C991 PP04C992 PP04C993 PP04C999 PP07H1 PP07H2
1  ( 1 )  " "       " "       " "       " "       " "   " "  " "      " "      " "      " "      "*"    " "   
2  ( 1 )  "*"       " "       " "       " "       " "   " "  " "      " "      " "      " "      " "    "*"   
3  ( 1 )  " "       " "       " "       " "       "*"   " "  " "      " "      " "      " "      "*"    " "   
4  ( 1 )  " "       " "       " "       " "       "*"   " "  " "      " "      " "      " "      "*"    " "   
5  ( 1 )  " "       " "       " "       " "       "*"   "*"  " "      " "      " "      " "      "*"    " "   
6  ( 1 )  " "       " "       " "       " "       "*"   "*"  " "      " "      " "      " "      "*"    " "   
7  ( 1 )  " "       " "       " "       " "       "*"   "*"  " "      " "      " "      " "      "*"    " "   
8  ( 1 )  " "       " "       " "       " "       "*"   "*"  " "      " "      " "      " "      "*"    " "   
9  ( 1 )  " "       " "       " "       " "       "*"   "*"  " "      " "      " "      " "      "*"    " "   
10  ( 1 ) " "       " "       " "       " "       "*"   "*"  " "      "*"      " "      " "      "*"    " "   
11  ( 1 ) " "       " "       " "       " "       "*"   "*"  "*"      "*"      " "      " "      "*"    " "   
12  ( 1 ) " "       " "       " "       " "       "*"   "*"  "*"      "*"      " "      "*"      "*"    " "   
13  ( 1 ) " "       " "       " "       " "       "*"   "*"  "*"      "*"      "*"      "*"      "*"    " "   
14  ( 1 ) " "       " "       " "       " "       "*"   "*"  "*"      "*"      "*"      "*"      "*"    "*"   
15  ( 1 ) " "       "*"       " "       " "       "*"   "*"  "*"      "*"      "*"      "*"      "*"    "*"   
16  ( 1 ) "*"       " "       " "       " "       "*"   "*"  "*"      "*"      "*"      "*"      "*"    "*"   

## Adjusted R2
plot(regsubsets.out, scale = "adjr2", main = "Adjusted R^2")

con which.max() podemos ubicar el modelo con el mejor ajuste. En nuestro caso elegimos comparar por \(R^2\).

which.max(summary(regsubsets.out)$adjr2)

[1] 13

Sabiendo que el modelo con 13 covariables es el que mejor ajusta, podemos observar con más detalle el mejor de los modelos con 13 covariables

summary(regsubsets.out)$which[13,]

(Intercept)   NIVEL_ED2   NIVEL_ED3   NIVEL_ED4   NIVEL_ED5   NIVEL_ED6   NIVEL_ED7     ESTADO1     ESTADO2 
       TRUE        TRUE        TRUE        TRUE        TRUE        TRUE       FALSE       FALSE       FALSE 
    ESTADO3     ESTADO4   CAT_OCUP1   CAT_OCUP2   CAT_OCUP3   CAT_OCUP4   CAT_OCUP9       CH042        CH06 
      FALSE       FALSE        TRUE       FALSE       FALSE       FALSE       FALSE        TRUE        TRUE 
   PP04C991    PP04C992    PP04C993    PP04C999      PP07H1      PP07H2 
       TRUE        TRUE        TRUE        TRUE        TRUE       FALSE 

Como las dummies de la misma categorica se deben incluir todas, sólo puedo concluir que la variable ESTADO no es significativa para el modelo, pero todas las demás sí. Como se vé, este método no es muy útil si pensamos trabajar con muchas variables categóricas, ya que no podemos indicarle que todas las dummies pertenecientes a una misma covariable deben ser consideradas en conjunto (por lo menos para esta librería).

modelo.ajustado <- lm(P21 ~ NIVEL_ED + CAT_OCUP + CH04 + CH06 + PP04C99 + PP07H,
               data = base)
summary(modelo.ajustado)

Call:
lm(formula = P21 ~ NIVEL_ED + CAT_OCUP + CH04 + CH06 + PP04C99 + 
    PP07H, data = base)

Residuals:
   Min     1Q Median     3Q    Max 
-19890  -4015  -1105   2313 376024 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9740.395    521.493  18.678  < 2e-16 ***
NIVEL_ED2     1157.565    337.677   3.428 0.000609 ***
NIVEL_ED3     2204.099    337.717   6.526 6.91e-11 ***
NIVEL_ED4     3732.148    328.143  11.374  < 2e-16 ***
NIVEL_ED5     4671.176    351.836  13.277  < 2e-16 ***
NIVEL_ED6     7922.073    338.660  23.392  < 2e-16 ***
NIVEL_ED7      -74.977   1016.887  -0.074 0.941224    
CAT_OCUP2    -8336.445    360.867 -23.101  < 2e-16 ***
CAT_OCUP3   -15104.260   7829.490  -1.929 0.053728 .  
CH042        -3946.775    121.013 -32.614  < 2e-16 ***
CH06            95.655      4.853  19.709  < 2e-16 ***
PP04C991     -1671.932    877.720  -1.905 0.056814 .  
PP04C992     -1194.144    377.407  -3.164 0.001558 ** 
PP04C993      -434.592    327.381  -1.327 0.184366    
PP04C999      -537.287    305.261  -1.760 0.078408 .  
PP07H1       12970.539   7822.940   1.658 0.097332 .  
PP07H2        6810.680   7821.678   0.871 0.383905    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7819 on 18499 degrees of freedom
Multiple R-squared:  0.279, Adjusted R-squared:  0.2783 
F-statistic: 447.3 on 16 and 18499 DF,  p-value: < 2.2e-16

Ejemplo trivial

Por útlimo, utilizamos un dataset de datos de manual, para hacer algunos ejercicios más. Los datos son tomados de Heinz, Peterson, Johnson, y Kerk {2003}, y se encuentran en la base de datos bdims del paquete openintro.

library(openintro,warn.conflicts = FALSE)

Primero realizamos un diagrama de dispersión que muestre la relación entre el peso medido en kilogramos (wgt)y la circunferencia de la cadera medida en centímetros (hip.gi), diferenciando por sexo (sex = 1 para los hombres y sex = 0 para las mujeres)

Primero calculamos el modelo sin diferenciar por sexo

lm.ajustado <- lm(wgt~hip.gi,data = bdims)
summary(lm.ajustado)

Call:
lm(formula = wgt ~ hip.gi, data = bdims)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.449  -7.111  -0.367   7.264  22.353 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -78.21054    5.56902  -14.04   <2e-16 ***
hip.gi        1.52417    0.05747   26.52   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.636 on 505 degrees of freedom
Multiple R-squared:  0.5821,    Adjusted R-squared:  0.5813 
F-statistic: 703.5 on 1 and 505 DF,  p-value: < 2.2e-16

Lo graficamos

plot(bdims$hip.gi,bdims$wgt)
abline(lm.ajustado, col = "red")

Una alternativa en ggplot (no utilizamos el modelo construido, sino que lo volvemos a generar)

ggplot(bdims, aes(hip.gi, wgt)) + 
  geom_point()+ 
  geom_smooth(method = "lm")

Prediccion

Para predecir una nueva observación, podemos calular el intervalo de confianza o el intervalo de predicción con el comando predict.lm()
Por ejemplo, si elegimos una persona con contorno de cadera mide 100 cm.

new <- data.frame(hip.gi = 100)
#intervalo de confianza
predict.lm(lm.ajustado,newdata = new,interval="confidence",level = 0.95)

       fit      lwr    upr
1 74.20646 73.36492 75.048

#intervalo de predicción
predict.lm(lm.ajustado,newdata = new,interval="prediction",level = 0.95)

       fit      lwr      upr
1 74.20646 57.21927 91.19365

anova

So queremos los datos de la tabla anova, lo hacemos con el comando anova()

anova(lm.ajustado)

Analysis of Variance Table

Response: wgt
           Df Sum Sq Mean Sq F value    Pr(>F)    
hip.gi      1  52463   52463  703.49 < 2.2e-16 ***
Residuals 505  37661      75                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Si queremos calcular el modelo con la interacción por sexo.

lm.ajustado <- lm(wgt ~ hip.gi*sex ,data = bdims)
summary(lm.ajustado)

Call:
lm(formula = wgt ~ hip.gi * sex, data = bdims)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.9539  -2.6721  -0.2353   2.6641  21.6628 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -58.97920    3.96073 -14.891  < 2e-16 ***
hip.gi        1.25014    0.04130  30.270  < 2e-16 ***
sex          -7.67232    6.09012  -1.260 0.208326    
hip.gi:sex    0.23095    0.06274   3.681 0.000257 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.613 on 503 degrees of freedom
Multiple R-squared:  0.8812,    Adjusted R-squared:  0.8805 
F-statistic:  1244 on 3 and 503 DF,  p-value: < 2.2e-16

Para graficar el resultado, recupero los valores predichos, que estan guardados en el objeto generado por lm(), y se los agrego a la tabla original.

bdims$yhat <- lm.ajustado$fitted.values
ggplot(bdims, aes(group = sex, color = factor(sex))) + 
  geom_point(aes(hip.gi, wgt))+
  geom_line(aes(hip.gi, yhat))+
  theme_minimal()

La interacción resulta significativa, pero no la dummy. Además, parece que el término de interacción no es realmente necesario, sino más bien que es preferible el supuesto de que ambas poblaciones comparten la misma pendiente, con una diferente ordenada al origen.

Podemos probar qué sucede cuando en lugar de utilizar la interacción, simplemente agregamos la dummy de sexo

lm.ajustado <- lm(wgt ~ hip.gi + sex ,data = bdims)
lm.ajustado <- lm(wgt ~ hip.gi + sex ,data = bdims)
summary(lm.ajustado)

Call:
lm(formula = wgt ~ hip.gi + sex, data = bdims)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.8293  -2.9544  -0.1366   2.4290  19.8266 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -68.55225    3.02439  -22.67   <2e-16 ***
hip.gi        1.35022    0.03147   42.90   <2e-16 ***
sex          14.69455    0.42024   34.97   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.67 on 504 degrees of freedom
Multiple R-squared:  0.878, Adjusted R-squared:  0.8775 
F-statistic:  1814 on 2 and 504 DF,  p-value: < 2.2e-16

La diferencia en el R cuadrado ajustado es menor, por lo que por simplicidad eligiríamos este modelo.

Nuevamente, si queremos graficar ambas rectas, podemos recuperar los valores predichos

Logit

Por útlimo, hacemos un ejercicio simple de una regresión logística, que contrasta el deseo de trabajar más horas respecto a la cantidad de hroas trabajadas.

variables <- c('CH04',        # Sexo
               'PP3E_TOT',    # Horas trabajadas
               'PP03G')       # La semana pasada, ¿quería trabajar más horas? 1-Si 2-No
# Me quedo con los asalariados que respondieron sobre su ingreso
datos <- individual.416 %>% 
  select(variables) %>% 
   filter(!is.na(PP3E_TOT),
          !is.na(PP03G),
          !PP3E_TOT %in% c(0,999),
          PP03G    !=9) %>%
  mutate(PP03G= case_when(PP03G==1 ~ 1,
                          PP03G==2 ~ 0),
         Sexo = case_when(CH04==1 ~'Varon',
                          CH04==2 ~'Mujer')) %>% 
  select(-CH04)
datos

model <- glm(PP03G ~.,family=binomial(link='logit'),data=datos)
summary(model)

Call:
glm(formula = PP03G ~ ., family = binomial(link = "logit"), data = datos)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4542  -0.5422  -0.4012  -0.2565   3.6664  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.089322   0.047374   1.885   0.0594 .  
PP3E_TOT    -0.070682   0.001607 -43.993   <2e-16 ***
SexoVaron    0.612064   0.044361  13.797   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 18255  on 22989  degrees of freedom
Residual deviance: 15865  on 22987  degrees of freedom
AIC: 15871

Number of Fisher Scoring iterations: 5

Si queremos graficarlo. Podemos utilizar lo valores estimados del modelo.

Otra opción es calcular dos modelos logit por separado, según sexo. Esto tiene la ventaja de que incluye los intervalos de confianza del modelo, y la desventaja de que no se trata de un mismo modelo con un control por sexo, sino de dos modelos calculados por separado.

 ggplot(datos, aes(x=PP3E_TOT, y=PP03G, color= Sexo, group=Sexo))+
  geom_jitter(alpha=0.03,width = 0,height = 0.01)+ 
  stat_smooth(method="glm", method.args=list(family="binomial"), se=T)+
  geom_vline(xintercept = 35, linetype= 'dashed')+
  labs(x = "Total de Horas trabajadas en la semana en la ocupación principal",
       y ="¿Quería trabajar más horas?",
       title= "Probabilidad de querer trabajar más horas según género")+
  theme_minimal()+
  theme(legend.position = "bottom",
        text = element_text(size=15))+
  scale_y_continuous(breaks = c(0,.5,1),labels = c("No","50%","Si"))+
  scale_x_continuous(limits = c(0,100))