Ejercicios Regresión

Ejempo 1

Tipo: Regresión Lineal Simple
Nombre: Empresa y salarios
Descripción: Vamos a ver si podemor predecir el salario actual de la persona en función del salario con el que entró a trabajar
Biliografia: https://www.youtube.com/watch?v=-p02G7NXlSk&list=PLKc8oG2WuyL9nO8z0ImlZ0SgRWm2NuKBb

Importar y explorar datos 1

url1<-"https://raw.githubusercontent.com/armandovl/datasets_uno/main/empresa_salarios.csv" 
datos1<-read.csv(url(url1)) #llamamos el dataset mediante una url
head(datos1,7) #los primeros 7 registros

##   id sexo    fechnac educ catlab salario_actual salario_inicial antiguedad
## 1  1    1 03/02/1952   15      3          57000           27000         98
## 2  2    1 23/05/1958   16      1          40200           18750         98
## 3  3    0 26/07/1929   12      1          21450           12000         98
## 4  4    0 15/04/1947    8      1          21900           13200         98
## 5  5    1 09/02/1955   15      1          45000           21000         98
## 6  6    1 22/08/1958   15      1          32100           13500         98
## 7  7    1 26/04/1956   15      1          36000           18750         98
##   experiencia minoria directivo
## 1         144       0         1
## 2          36       0         0
## 3         381       0         0
## 4         190       0         0
## 5         138       0         0
## 6          67       0         0
## 7         114       0         0

#estructura de los datos
str(datos1)

## 'data.frame':    474 obs. of  11 variables:
##  $ id             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ sexo           : int  1 1 0 0 1 1 1 0 0 0 ...
##  $ fechnac        : chr  "03/02/1952" "23/05/1958" "26/07/1929" "15/04/1947" ...
##  $ educ           : int  15 16 12 8 15 15 15 12 15 12 ...
##  $ catlab         : int  3 1 1 1 1 1 1 1 1 1 ...
##  $ salario_actual : int  57000 40200 21450 21900 45000 32100 36000 21900 27900 24000 ...
##  $ salario_inicial: int  27000 18750 12000 13200 21000 13500 18750 9750 12750 13500 ...
##  $ antiguedad     : int  98 98 98 98 98 98 98 98 98 98 ...
##  $ experiencia    : int  144 36 381 190 138 67 114 0 115 244 ...
##  $ minoria        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ directivo      : int  1 0 0 0 0 0 0 0 0 0 ...

#estadísticos desrcriptivos
summary(datos1)

##        id             sexo          fechnac               educ      
##  Min.   :  1.0   Min.   :0.0000   Length:474         Min.   : 8.00  
##  1st Qu.:119.2   1st Qu.:0.0000   Class :character   1st Qu.:12.00  
##  Median :237.5   Median :1.0000   Mode  :character   Median :12.00  
##  Mean   :237.5   Mean   :0.5443                      Mean   :13.49  
##  3rd Qu.:355.8   3rd Qu.:1.0000                      3rd Qu.:15.00  
##  Max.   :474.0   Max.   :1.0000                      Max.   :21.00  
##      catlab      salario_actual   salario_inicial   antiguedad   
##  Min.   :1.000   Min.   : 15750   Min.   : 9000   Min.   :63.00  
##  1st Qu.:1.000   1st Qu.: 24000   1st Qu.:12488   1st Qu.:72.00  
##  Median :1.000   Median : 28875   Median :15000   Median :81.00  
##  Mean   :1.411   Mean   : 34420   Mean   :17016   Mean   :81.11  
##  3rd Qu.:1.000   3rd Qu.: 36938   3rd Qu.:17490   3rd Qu.:90.00  
##  Max.   :3.000   Max.   :135000   Max.   :79980   Max.   :98.00  
##   experiencia        minoria         directivo     
##  Min.   :  0.00   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.: 19.25   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median : 55.00   Median :0.0000   Median :0.0000  
##  Mean   : 95.86   Mean   :0.2194   Mean   :0.1772  
##  3rd Qu.:138.75   3rd Qu.:0.0000   3rd Qu.:0.0000  
##  Max.   :476.00   Max.   :1.0000   Max.   :1.0000

vamos a ver la correlación entre esas dos variables

cor.test(datos1$salario_inicial,datos1$salario_actual)

## 
##  Pearson's product-moment correlation
## 
## data:  datos1$salario_inicial and datos1$salario_actual
## t = 40.276, df = 472, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8580696 0.8989267
## sample estimates:
##       cor 
## 0.8801175

#tambien se puede hacer como correlación de kendall o spearman
#cor.test(datos1$salario_inicial,datos1$salario_actual, method = "kendall")

Visualizar la variables 1

Vamos a ver graficamente como se comporta el salario en actual en funcion del salario inicial

#Para ver los parametros, se escribe en la consola ?plot

plot(datos1$salario_inicial,datos1$salario_actual,
     col="green3",
     pch=18,
     xlab="Salario inicial",
     ylab="Salario actual",
     main="Relación entre el salario actual y el salarion inicial",
     sub="relacion lineal")

Parece que sí se ve una forma de relación lineal

Hacer el modelo 1

# alt 126 = ~
modelo1<-lm(salario_actual~salario_inicial,datos1)

summary(modelo1)

## 
## Call:
## lm(formula = salario_actual ~ salario_inicial, data = datos1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -35424  -4031  -1154   2584  49293 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.928e+03  8.887e+02    2.17   0.0305 *  
## salario_inicial 1.909e+00  4.741e-02   40.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8115 on 472 degrees of freedom
## Multiple R-squared:  0.7746, Adjusted R-squared:  0.7741 
## F-statistic:  1622 on 1 and 472 DF,  p-value: < 2.2e-16

Tabla de annova: p-value: < 2.2e-16 ,sí existe una pendiente, e modelo si es significativo
Coeficientes de determinación R2: 0.7741 (se usa el r2 ajustado que es más riguroso a la hora de interpretar la bondad de ajuste del modelo, cuanto del porcentaje de la variabilidad de la variable predicha es explicado por la variable predictora
Coeficientes y constuir la ecuación:

modelo1$coefficients

##     (Intercept) salario_inicial 
##      1928.20576         1.90945

Comprobar los supuestos 1

Linealidad 1

#la línea roja tiene que estar sobre la línea punteada, en este caso no se cumple
plot(modelo1,1)

Media de error 1

residuos=residuals(modelo1)
mean(residuos)

## [1] -1.862293e-13

a pesar de que es negativo, Se acerca mucho al cero

Normalidad resiuos 1

Se ve que tiene una alta curtosis

#hacemos la gráfica de histograma
hist(residuos, breaks=50 , col="lightblue", freq = FALSE)

#obtenemos la línea de densidad
lines(density(residuos), lwd = 2, col = 'blue')

#tendencia de curva normal
x <- seq(min(residuos), max(residuos), length = 40)
f <- dnorm(x, mean = mean(residuos), sd = sd(residuos))
lines(x, f, col = "red", lwd = 2)

#todos se deben acercar a la línea punteada, pero hay varios que no
plot(modelo1,2)

shapiro.test(residuos)# este test solo es para menos de 50 observaciones

## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.85688, p-value < 2.2e-16

nos dice que no hay normalidad, pero vamos a hacer otro tipo de prueba

# Alternativamente conviene usar la modificación de Lilliefors a este test.
# Esta corrección considera que los parámetros son estimados, a diferencia
# del «ks» «a secas»:

library(nortest)
lillie.test(residuos)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuos
## D = 0.12727, p-value < 2.2e-16

# realizar la prueba Anderson-Darling para comprobar la normalidad >50 observaciones
ad.test (residuos)

## 
##  Anderson-Darling normality test
## 
## data:  residuos
## A = 16.125, p-value < 2.2e-16

En el caso de las dos pruebas anteriores, nos salen valores menores a 0.05 entonces no hay normalidad, H0= Normalidad Ha= No normalidad

par(mfrow=c(2,1)) #para dividir el gráfico en dos
hist(residuos, col="yellow",breaks=60) #color y cortes
boxplot(residuos, bty="l", range=1.5, col="yellow", horizontal=T,xlab="residuos") #boxplot

par(mfrow=c(1,1)) #devolver a los parametros iniciales

En el gráfico anterior s eve que hay muchos valores atípicos, sobrepasan 1 y medio el rango intercuartil

#
plot(modelo1,5)

Homogeneidad de varianzas 1

Graficamente se ve con el grafico 1 del modelo, los puntos no tienen que tener forma de embudo

plot(modelo1,1)

predichos=fitted(modelo1)
plot(residuos~predichos)

plot(modelo1,3)

EN el grafico de escalas la linea roja debe ser horizontal y no debe mostrar ningun patron

parece que si hay una forma de embudo, pero vamos a hacer la prueba Breusch-Pagan
h0=no hay heterocedasticidad, ha= si hay heterocedasticidad

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(modelo1)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo1
## BP = 59.224, df = 1, p-value = 1.407e-14

el pvalor es menor a 0.05 es decir que se rechaza la h0, si hay heterocedasticidad

Autocorelación 1

Multicolineidad 1

Hacer gráfico del modelo 1

plot(datos1$salario_inicial,datos1$salario_actual,
     col="green3",
     pch=18,
     xlab="Salario inicial",
     ylab="Salario actual",
     main="Relación entre el salario actual y el salario inicial",
     sub="relacion lineal")
abline(modelo1, col="red")

plot(datos1$salario_inicial,residuos,)

cor(datos1$salario_inicial,residuos)

## [1] 1.729793e-17

Ejempo 2

Tipo: Regresión multiple
Nombre: Empresa y salarios
Descripción: Vamos a ver si podemor predecir el salario actual de la persona en función del salario con el que entró a trabajar y también la variable de años de educacion
Biliografia: https://www.youtube.com/watch?v=-p02G7NXlSk&list=PLKc8oG2WuyL9nO8z0ImlZ0SgRWm2NuKBb

Importar y explorar datos 1

library(rio)
datos2<-import("https://raw.githubusercontent.com/armandovl/datasets_uno/main/empresa_salarios.csv")
head(datos2,7) #los primeros 7 registros

##   id sexo    fechnac educ catlab salario_actual salario_inicial antiguedad
## 1  1    1 03/02/1952   15      3          57000           27000         98
## 2  2    1 23/05/1958   16      1          40200           18750         98
## 3  3    0 26/07/1929   12      1          21450           12000         98
## 4  4    0 15/04/1947    8      1          21900           13200         98
## 5  5    1 09/02/1955   15      1          45000           21000         98
## 6  6    1 22/08/1958   15      1          32100           13500         98
## 7  7    1 26/04/1956   15      1          36000           18750         98
##   experiencia minoria directivo
## 1         144       0         1
## 2          36       0         0
## 3         381       0         0
## 4         190       0         0
## 5         138       0         0
## 6          67       0         0
## 7         114       0         0

#estructura de los datos
str(datos2)

## 'data.frame':    474 obs. of  11 variables:
##  $ id             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ sexo           : int  1 1 0 0 1 1 1 0 0 0 ...
##  $ fechnac        : chr  "03/02/1952" "23/05/1958" "26/07/1929" "15/04/1947" ...
##  $ educ           : int  15 16 12 8 15 15 15 12 15 12 ...
##  $ catlab         : int  3 1 1 1 1 1 1 1 1 1 ...
##  $ salario_actual : int  57000 40200 21450 21900 45000 32100 36000 21900 27900 24000 ...
##  $ salario_inicial: int  27000 18750 12000 13200 21000 13500 18750 9750 12750 13500 ...
##  $ antiguedad     : int  98 98 98 98 98 98 98 98 98 98 ...
##  $ experiencia    : int  144 36 381 190 138 67 114 0 115 244 ...
##  $ minoria        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ directivo      : int  1 0 0 0 0 0 0 0 0 0 ...

#estadísticos desrcriptivos
summary(datos2)

##        id             sexo          fechnac               educ      
##  Min.   :  1.0   Min.   :0.0000   Length:474         Min.   : 8.00  
##  1st Qu.:119.2   1st Qu.:0.0000   Class :character   1st Qu.:12.00  
##  Median :237.5   Median :1.0000   Mode  :character   Median :12.00  
##  Mean   :237.5   Mean   :0.5443                      Mean   :13.49  
##  3rd Qu.:355.8   3rd Qu.:1.0000                      3rd Qu.:15.00  
##  Max.   :474.0   Max.   :1.0000                      Max.   :21.00  
##      catlab      salario_actual   salario_inicial   antiguedad   
##  Min.   :1.000   Min.   : 15750   Min.   : 9000   Min.   :63.00  
##  1st Qu.:1.000   1st Qu.: 24000   1st Qu.:12488   1st Qu.:72.00  
##  Median :1.000   Median : 28875   Median :15000   Median :81.00  
##  Mean   :1.411   Mean   : 34420   Mean   :17016   Mean   :81.11  
##  3rd Qu.:1.000   3rd Qu.: 36938   3rd Qu.:17490   3rd Qu.:90.00  
##  Max.   :3.000   Max.   :135000   Max.   :79980   Max.   :98.00  
##   experiencia        minoria         directivo     
##  Min.   :  0.00   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.: 19.25   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median : 55.00   Median :0.0000   Median :0.0000  
##  Mean   : 95.86   Mean   :0.2194   Mean   :0.1772  
##  3rd Qu.:138.75   3rd Qu.:0.0000   3rd Qu.:0.0000  
##  Max.   :476.00   Max.   :1.0000   Max.   :1.0000

Ver la relación entre variables cuantitativa

library(PerformanceAnalytics)
chart.Correlation(datos2[,4:11], histogram = F )

vamos a ver la correlación entre todas las variables

cor(datos2[,c(4,6:7)])

##                      educ salario_actual salario_inicial
## educ            1.0000000      0.6605589       0.6331956
## salario_actual  0.6605589      1.0000000       0.8801175
## salario_inicial 0.6331956      0.8801175       1.0000000

#tambien se puede hacer como correlación de kendall o spearman
#cor.test(datos1$salario_inicial,datos1$salario_actual, method = "kendall")

parece que hay multi colinealidad entre las x (salario inicial, educación)

Visualizar la variables 1

Vamos a ver graficamente como se comporta el salario en actual en funcion del salario inicial

#Para ver los parametros, se escribe en la consola ?plot

plot(datos2$salario_inicial,datos2$salario_actual,
     col="green3",
     pch=18,
     xlab="Salario inicial",
     ylab="Salario actual",
     main="Relación entre el salario actual y el salarion inicial",
     sub="relacion lineal")

plot(datos2$salario_inicial,datos2$educ,
     col="green3",
     pch=18,
     xlab="Salario inicial",
     ylab="Años educacion",
     main="Relación entre el salario actual y el años de educación",
     sub="relacion lineal")

no parece que haya una relación lineal

Hacer el modelo 1

# alt 126 = ~
modelo2<-lm(salario_actual~salario_inicial+educ,datos2)

summary(modelo2)

## 
## Call:
## lm(formula = salario_actual ~ salario_inicial + educ, data = datos2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -29581  -4093   -741   2917  49218 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -7.809e+03  1.754e+03  -4.452 1.06e-05 ***
## salario_inicial  1.673e+00  5.885e-02  28.423  < 2e-16 ***
## educ             1.020e+03  1.606e+02   6.356 4.91e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7797 on 471 degrees of freedom
## Multiple R-squared:  0.7924, Adjusted R-squared:  0.7915 
## F-statistic: 898.9 on 2 and 471 DF,  p-value: < 2.2e-16

Tabla de annova: p-value: < 2.2e-16, y todas la variables son significativas ,sí existe una pendiente, e modelo si es significativo
Coeficientes de determinación R2: 0.7915 (se usa el r2 ajustado que es más riguroso a la hora de interpretar la bondad de ajuste del modelo, cuanto del porcentaje de la variabilidad de la variable predicha es explicado por la variable predictora
Coeficientes y constuir la ecuación:

modelo2$coefficients

##     (Intercept) salario_inicial            educ 
##    -7808.714157        1.672631     1020.390142

Comprobar los supuestos 1

Linealidad 1

#la línea roja tiene que estar sobre la línea punteada, en este caso no se cumple
plot(modelo2,1)

la línea roja se aleja de la punteada, pero menos de lo que era con una sola variable

Media de error 1

residuos=residuals(modelo2)
mean(residuos)

## [1] -6.146089e-13

a pesar de que es negativo, Se acerca mucho al cero

Normalidad resiuos 1

Se ve que tiene una alta curtosis

#hacemos la gráfica de histograma
hist(residuos, breaks=50 , col="lightblue", freq = FALSE)

#obtenemos la línea de densidad
lines(density(residuos), lwd = 2, col = 'blue')

#tendencia de curva normal
x <- seq(min(residuos), max(residuos), length = 40)
f <- dnorm(x, mean = mean(residuos), sd = sd(residuos))
lines(x, f, col = "red", lwd = 2)

#todos se deben acercar a la línea punteada, pero hay varios que no
plot(modelo2,2)

shapiro.test(residuos)# este test solo es para menos de 50 observaciones

## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.87555, p-value < 2.2e-16

nos dice que no hay normalidad, pero vamos a hacer otro tipo de prueba

# Alternativamente conviene usar la modificación de Lilliefors a este test.
# Esta corrección considera que los parámetros son estimados, a diferencia
# del «ks» «a secas»:

library(nortest)
lillie.test(residuos)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuos
## D = 0.11123, p-value = 2.23e-15

# realizar la prueba Anderson-Darling para comprobar la normalidad >50 observaciones
ad.test (residuos)

## 
##  Anderson-Darling normality test
## 
## data:  residuos
## A = 12.136, p-value < 2.2e-16

En el caso de las dos pruebas anteriores, nos salen valores menores a 0.05 entonces no hay normalidad, H0= Normalidad Ha= No normalidad

par(mfrow=c(2,1)) #para dividir el gráfico en dos
hist(residuos, col="yellow",breaks=60) #color y cortes
boxplot(residuos, bty="l", range=1.5, col="yellow", horizontal=T,xlab="residuos") #boxplot

par(mfrow=c(1,1)) #devolver a los parametros iniciales

En el gráfico anterior s eve que hay muchos valores atípicos, sobrepasan 1 y medio el rango intercuartil

vamos a ver si hay valores influyentes (valores con mayor distancia de cook), los que están fuera de la línea punteada, estos se pueden quitar y ver como sale el modelo

#
plot(modelo2,5)

te nemos que ver que estos valores no sean iguales o mayores que 1

plot(modelo2,4)

#calcularlos manualmente
#datosCook<-cooks.distance(modelo2) #obtener los datos de cook
#which(ddatosCook>1) #solo me de los valores mayores a 1

Homogeneidad de varianzas 1

Graficamente se ve con el grafico 1 del modelo, los puntos no tienen que tener forma de embudo

plot(modelo2,1)

predichos=fitted(modelo2)
plot(residuos~predichos)

plot(modelo2,3)

EN el grafico de escalas la linea roja debe ser horizontal y no debe mostrar ningun patron

Multicolineidad 1

parece que si hay una forma de embudo, pero vamos a hacer la prueba Breusch-Pagan
h0=no hay heterocedasticidad, ha= si hay heterocedasticidad

library(lmtest)
bptest(modelo2)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo2
## BP = 49.22, df = 2, p-value = 2.052e-11

el pvalor es menor a 0.05 es decir que se rechaza la h0, si hay heterocedasticidad

Hacermos la prueba de factores de inflacion de varianza
donde el valor no tiene que estar entre 5 y diez

library(DescTools)
VIF(modelo2)

## salario_inicial            educ 
##        1.669273        1.669273

Autocorelación 1

Hacer gráfico del modelo 1

library(ggplot2)

str(datos2)

## 'data.frame':    474 obs. of  11 variables:
##  $ id             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ sexo           : int  1 1 0 0 1 1 1 0 0 0 ...
##  $ fechnac        : chr  "03/02/1952" "23/05/1958" "26/07/1929" "15/04/1947" ...
##  $ educ           : int  15 16 12 8 15 15 15 12 15 12 ...
##  $ catlab         : int  3 1 1 1 1 1 1 1 1 1 ...
##  $ salario_actual : int  57000 40200 21450 21900 45000 32100 36000 21900 27900 24000 ...
##  $ salario_inicial: int  27000 18750 12000 13200 21000 13500 18750 9750 12750 13500 ...
##  $ antiguedad     : int  98 98 98 98 98 98 98 98 98 98 ...
##  $ experiencia    : int  144 36 381 190 138 67 114 0 115 244 ...
##  $ minoria        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ directivo      : int  1 0 0 0 0 0 0 0 0 0 ...

plot<-ggplot(datos2,aes(x =salario_inicial, y =salario_actual)) + 
          
 geom_point(aes(size=educ, alpha=0.3))+ #, shape=island, color = sex 
  

#           xlab('Largo del pico') +
#           ylab('Profundidad del pico') +
#           ggtitle('Tres variables cuanti, dos cuali') +
#           theme_minimal()
# 
# plot$labels$colour = "Sexo"
# plot$labels$size = "Peso"
# plot$labels$shape = "Isla"

#geom_smooth(method="lm",color="red") +
#geom_smooth(data=datos2,method=lm, formula=y~(x^2) ,size=1)

stat_smooth(data=datos2,method="lm", formula="y~x^2" ,size=1,color="purple")

#poner otro modelo cuadratico
#+ stat_smooth(data=datos2,method="lm", formula="y~x+I(x^2)" ,size=1,color="red", )

plot

library(ggExtra)

plot<-ggplot(datos2,aes(x =salario_inicial, y =salario_actual)) + 
          
 geom_point(aes(size=educ, alpha=0.3))+ #, shape=island, color = sex  size=
  

#           xlab('Largo del pico') +
#           ylab('Profundidad del pico') +
#           ggtitle('Tres variables cuanti, dos cuali') +
#           theme_minimal()
# 
# plot$labels$colour = "Sexo"
# plot$labels$size = "Peso"
# plot$labels$shape = "Isla"

#geom_smooth(method="lm",color="red") +
#geom_smooth(data=datos2,method=lm, formula=y~(x^2) ,size=1)

stat_smooth(data=datos2,method="lm", formula="y~x^2" ,size=1,color="purple")

#poner otro modelo cuadratico
#+ stat_smooth(data=datos2,method="lm", formula="y~x+I(x^2)" ,size=1,color="red", )

#adicional boxplot
ggMarginal(plot,type="boxplot")

Ejempo 3

Tipo: Regresión multiple
Nombre: Empresa y salarios
Descripción: Vamos a ver si podemor predecir el salario actual de la persona en función varias varables, tanto continuas, como dummies
Biliografia: https://www.youtube.com/watch?v=-p02G7NXlSk&list=PLKc8oG2WuyL9nO8z0ImlZ0SgRWm2NuKBb

Importar y explorar datos 1

library(rio)
datos3<-import("https://raw.githubusercontent.com/armandovl/datasets_uno/main/empresa_salarios.csv")
head(datos3,7) #los primeros 7 registros

##   id sexo    fechnac educ catlab salario_actual salario_inicial antiguedad
## 1  1    1 03/02/1952   15      3          57000           27000         98
## 2  2    1 23/05/1958   16      1          40200           18750         98
## 3  3    0 26/07/1929   12      1          21450           12000         98
## 4  4    0 15/04/1947    8      1          21900           13200         98
## 5  5    1 09/02/1955   15      1          45000           21000         98
## 6  6    1 22/08/1958   15      1          32100           13500         98
## 7  7    1 26/04/1956   15      1          36000           18750         98
##   experiencia minoria directivo
## 1         144       0         1
## 2          36       0         0
## 3         381       0         0
## 4         190       0         0
## 5         138       0         0
## 6          67       0         0
## 7         114       0         0

#estructura de los datos
str(datos3)

## 'data.frame':    474 obs. of  11 variables:
##  $ id             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ sexo           : int  1 1 0 0 1 1 1 0 0 0 ...
##  $ fechnac        : chr  "03/02/1952" "23/05/1958" "26/07/1929" "15/04/1947" ...
##  $ educ           : int  15 16 12 8 15 15 15 12 15 12 ...
##  $ catlab         : int  3 1 1 1 1 1 1 1 1 1 ...
##  $ salario_actual : int  57000 40200 21450 21900 45000 32100 36000 21900 27900 24000 ...
##  $ salario_inicial: int  27000 18750 12000 13200 21000 13500 18750 9750 12750 13500 ...
##  $ antiguedad     : int  98 98 98 98 98 98 98 98 98 98 ...
##  $ experiencia    : int  144 36 381 190 138 67 114 0 115 244 ...
##  $ minoria        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ directivo      : int  1 0 0 0 0 0 0 0 0 0 ...

#estadísticos desrcriptivos
summary(datos3)

##        id             sexo          fechnac               educ      
##  Min.   :  1.0   Min.   :0.0000   Length:474         Min.   : 8.00  
##  1st Qu.:119.2   1st Qu.:0.0000   Class :character   1st Qu.:12.00  
##  Median :237.5   Median :1.0000   Mode  :character   Median :12.00  
##  Mean   :237.5   Mean   :0.5443                      Mean   :13.49  
##  3rd Qu.:355.8   3rd Qu.:1.0000                      3rd Qu.:15.00  
##  Max.   :474.0   Max.   :1.0000                      Max.   :21.00  
##      catlab      salario_actual   salario_inicial   antiguedad   
##  Min.   :1.000   Min.   : 15750   Min.   : 9000   Min.   :63.00  
##  1st Qu.:1.000   1st Qu.: 24000   1st Qu.:12488   1st Qu.:72.00  
##  Median :1.000   Median : 28875   Median :15000   Median :81.00  
##  Mean   :1.411   Mean   : 34420   Mean   :17016   Mean   :81.11  
##  3rd Qu.:1.000   3rd Qu.: 36938   3rd Qu.:17490   3rd Qu.:90.00  
##  Max.   :3.000   Max.   :135000   Max.   :79980   Max.   :98.00  
##   experiencia        minoria         directivo     
##  Min.   :  0.00   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.: 19.25   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median : 55.00   Median :0.0000   Median :0.0000  
##  Mean   : 95.86   Mean   :0.2194   Mean   :0.1772  
##  3rd Qu.:138.75   3rd Qu.:0.0000   3rd Qu.:0.0000  
##  Max.   :476.00   Max.   :1.0000   Max.   :1.0000

#convertir a factor varias palabras

# names<-c(1:3,5)
# mydata[,names]<-lapply(mydata[,names], factor)

names<-c("sexo","minoria","directivo","catlab")
datos3[,names]<-lapply(datos3[,names], factor)
str(datos3)

## 'data.frame':    474 obs. of  11 variables:
##  $ id             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ sexo           : Factor w/ 2 levels "0","1": 2 2 1 1 2 2 2 1 1 1 ...
##  $ fechnac        : chr  "03/02/1952" "23/05/1958" "26/07/1929" "15/04/1947" ...
##  $ educ           : int  15 16 12 8 15 15 15 12 15 12 ...
##  $ catlab         : Factor w/ 3 levels "1","2","3": 3 1 1 1 1 1 1 1 1 1 ...
##  $ salario_actual : int  57000 40200 21450 21900 45000 32100 36000 21900 27900 24000 ...
##  $ salario_inicial: int  27000 18750 12000 13200 21000 13500 18750 9750 12750 13500 ...
##  $ antiguedad     : int  98 98 98 98 98 98 98 98 98 98 ...
##  $ experiencia    : int  144 36 381 190 138 67 114 0 115 244 ...
##  $ minoria        : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ directivo      : Factor w/ 2 levels "0","1": 2 1 1 1 1 1 1 1 1 1 ...

table(datos3$catlab)

## 
##   1   2   3 
## 363  27  84

Ver la relación entre variables cuantitativa

library(PerformanceAnalytics)
chart.Correlation(datos3[,c(4,6:9)], histogram = F )

vamos a ver la correlación entre todas las variables

cor(datos3[,c(4,6:9)])

##                        educ salario_actual salario_inicial   antiguedad
## educ             1.00000000     0.66055891      0.63319565  0.047378777
## salario_actual   0.66055891     1.00000000      0.88011747  0.084092267
## salario_inicial  0.63319565     0.88011747      1.00000000 -0.019753475
## antiguedad       0.04737878     0.08409227     -0.01975347  1.000000000
## experiencia     -0.25235252    -0.09746693      0.04513563  0.002978134
##                  experiencia
## educ            -0.252352521
## salario_actual  -0.097466926
## salario_inicial  0.045135627
## antiguedad       0.002978134
## experiencia      1.000000000

#tambien se puede hacer como correlación de kendall o spearman
#cor.test(datos1$salario_inicial,datos1$salario_actual, method = "kendall")

library(corrplot)

## corrplot 0.92 loaded

corrplot(cor(datos3[,c(4,6:9)]),method="circle")

parece que hay multi colinealidad entre variables

Visualizar la variables 1

Hacer el modelo 1

# alt 126 = ~
modelo3<-lm(salario_actual~salario_inicial+educ+minoria+directivo+sexo+catlab,datos3)

summary(modelo3)

## 
## Call:
## lm(formula = salario_actual ~ salario_inicial + educ + minoria + 
##     directivo + sexo + catlab, data = datos3)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -23447  -3513   -859   2548  47093 
## 
## Coefficients: (1 not defined because of singularities)
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       151.2521  2045.6150   0.074   0.9411    
## salario_inicial     1.2100     0.0753  16.070  < 2e-16 ***
## educ              785.6276   164.0812   4.788 2.26e-06 ***
## minoria1        -1276.7892   831.0418  -1.536   0.1251    
## directivo1      12037.1406  1443.0242   8.342 8.32e-16 ***
## sexo1            1913.8399   807.4902   2.370   0.0182 *  
## catlab2          3242.6141  1608.0999   2.016   0.0443 *  
## catlab3                 NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7216 on 467 degrees of freedom
## Multiple R-squared:  0.8237, Adjusted R-squared:  0.8214 
## F-statistic: 363.6 on 6 and 467 DF,  p-value: < 2.2e-16

Tabla de annova: p-value: < 2.2e-16, si hay pendiente pero no todas la variables son significativas
Coeficientes de determinación R2: 0.8214 (se usa el r2 ajustado que es más riguroso a la hora de interpretar la bondad de ajuste del modelo, cuanto del porcentaje de la variabilidad de la variable predicha es explicado por la variable predictora
Coeficientes y constuir la ecuación:

#va otra vez el modelo sin minoria
# alt 126 = ~
modelo3<-lm(salario_actual~salario_inicial+educ+directivo+sexo+catlab,datos3)

summary(modelo3)

## 
## Call:
## lm(formula = salario_actual ~ salario_inicial + educ + directivo + 
##     sexo + catlab, data = datos3)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -23521  -3434   -815   2761  47474 
## 
## Coefficients: (1 not defined because of singularities)
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -1.974e+02  2.036e+03  -0.097   0.9228    
## salario_inicial  1.217e+00  7.527e-02  16.168  < 2e-16 ***
## educ             7.866e+02  1.643e+02   4.787 2.27e-06 ***
## directivo1       1.224e+04  1.439e+03   8.501 2.53e-16 ***
## sexo1            1.757e+03  8.022e+02   2.190   0.0290 *  
## catlab2          3.019e+03  1.604e+03   1.882   0.0604 .  
## catlab3                 NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7227 on 468 degrees of freedom
## Multiple R-squared:  0.8228, Adjusted R-squared:  0.8209 
## F-statistic: 434.5 on 5 and 468 DF,  p-value: < 2.2e-16

#quitar catlab que tambien genera problemas
modelo3<-lm(salario_actual~salario_inicial+educ+directivo+sexo,datos3)
summary(modelo3)

## 
## Call:
## lm(formula = salario_actual ~ salario_inicial + educ + directivo + 
##     sexo, data = datos3)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -23841  -3552   -811   2808  47243 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.067e+03  1.927e+03   0.554  0.57994    
## salario_inicial 1.223e+00  7.541e-02  16.211  < 2e-16 ***
## educ            6.785e+02  1.544e+02   4.395 1.37e-05 ***
## directivo1      1.223e+04  1.443e+03   8.472 3.15e-16 ***
## sexo1           2.258e+03  7.588e+02   2.975  0.00308 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7246 on 469 degrees of freedom
## Multiple R-squared:  0.8214, Adjusted R-squared:  0.8199 
## F-statistic: 539.4 on 4 and 469 DF,  p-value: < 2.2e-16

ya mas o menos quedó el modelo, vamos a sacar los coeficientes

modelo3$coefficients

##     (Intercept) salario_inicial            educ      directivo1           sexo1 
##     1067.324880        1.222557      678.468029    12226.010592     2257.708452

Comprobar los supuestos 1

Linealidad 1

#la línea roja tiene que estar sobre la línea punteada, en este caso no se cumple
plot(modelo3,1)

la línea roja se aleja de la punteada, pero menos de lo que era con una sola variable

Media de error 1

residuos=residuals(modelo3)
mean(residuos)

## [1] 3.295545e-13

Se acerca mucho al cero

Normalidad resiuos 1

Se ve que tiene una alta curtosis

#hacemos la gráfica de histograma
hist(residuos, breaks=50 , col="lightblue", freq = FALSE)

#obtenemos la línea de densidad
lines(density(residuos), lwd = 2, col = 'blue')

#tendencia de curva normal
x <- seq(min(residuos), max(residuos), length = 40)
f <- dnorm(x, mean = mean(residuos), sd = sd(residuos))
lines(x, f, col = "red", lwd = 2)

#todos se deben acercar a la línea punteada, pero hay varios que no
plot(modelo3,2)

shapiro.test(residuos)# este test solo es para menos de 50 observaciones

## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.88391, p-value < 2.2e-16

nos dice que no hay normalidad, pero vamos a hacer otro tipo de prueba

# Alternativamente conviene usar la modificación de Lilliefors a este test.
# Esta corrección considera que los parámetros son estimados, a diferencia
# del «ks» «a secas»:

library(nortest)
lillie.test(residuos)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuos
## D = 0.10022, p-value = 2.429e-12

# realizar la prueba Anderson-Darling para comprobar la normalidad >50 observaciones
ad.test (residuos)

## 
##  Anderson-Darling normality test
## 
## data:  residuos
## A = 11.281, p-value < 2.2e-16

En el caso de las dos pruebas anteriores, nos salen valores menores a 0.05 entonces no hay normalidad, H0= Normalidad Ha= No normalidad

par(mfrow=c(2,1)) #para dividir el gráfico en dos
hist(residuos, col="yellow",breaks=60) #color y cortes
boxplot(residuos, bty="l", range=1.5, col="yellow", horizontal=T,xlab="residuos") #boxplot

par(mfrow=c(1,1)) #devolver a los parametros iniciales

En el gráfico anterior s eve que hay muchos valores atípicos, sobrepasan 1 y medio el rango intercuartil

vamos a ver si hay valores influyentes (valores con mayor distancia de cook), los que están fuera de la línea punteada, estos se pueden quitar y ver como sale el modelo

#
plot(modelo3,5)

te nemos que ver que estos valores no sean iguales o mayores que 1

plot(modelo3,4)

#calcularlos manualmente
#datosCook<-cooks.distance(modelo2) #obtener los datos de cook
#which(ddatosCook>1) #solo me de los valores mayores a 1

Homogeneidad de varianzas 1

Graficamente se ve con el grafico 1 del modelo, los puntos no tienen que tener forma de embudo

plot(modelo3,1)

predichos=fitted(modelo3)
plot(residuos~predichos)

plot(modelo3,3)

EN el grafico de escalas la linea roja debe ser horizontal y no debe mostrar ningun patron

Multicolineidad 1

parece que si hay una forma de embudo, pero vamos a hacer la prueba Breusch-Pagan
h0=no hay heterocedasticidad, ha= si hay heterocedasticidad

library(lmtest)
bptest(modelo3)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo3
## BP = 55.904, df = 4, p-value = 2.1e-11

el pvalor es menor a 0.05 es decir que se rechaza la h0, si hay heterocedasticidad

Hacermos la prueba de factores de inflacion de varianza
donde el valor no tiene que estar entre 5 y diez

library(DescTools)
VIF(modelo3)

## salario_inicial            educ       directivo            sexo 
##        3.173415        1.786305        2.741302        1.289283

Autocorelación 1

Es la prueba de durwin watson , solo se hace cuando sean modelos de series de tiempo

H 0 (hipótesis nula): No existe correlación entre los residuos.

H A (hipótesis alternativa): Los residuos están autocorrelacionados.

dwtest(modelo3) #si hay autocorrelación

## 
##  Durbin-Watson test
## 
## data:  modelo3
## DW = 1.7954, p-value = 0.01193
## alternative hypothesis: true autocorrelation is greater than 0

Hacer gráfico del modelo 1 hasta 3 cuai y 3 cuanti

library(ggplot2)

str(datos3)

## 'data.frame':    474 obs. of  11 variables:
##  $ id             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ sexo           : Factor w/ 2 levels "0","1": 2 2 1 1 2 2 2 1 1 1 ...
##  $ fechnac        : chr  "03/02/1952" "23/05/1958" "26/07/1929" "15/04/1947" ...
##  $ educ           : int  15 16 12 8 15 15 15 12 15 12 ...
##  $ catlab         : Factor w/ 3 levels "1","2","3": 3 1 1 1 1 1 1 1 1 1 ...
##  $ salario_actual : int  57000 40200 21450 21900 45000 32100 36000 21900 27900 24000 ...
##  $ salario_inicial: int  27000 18750 12000 13200 21000 13500 18750 9750 12750 13500 ...
##  $ antiguedad     : int  98 98 98 98 98 98 98 98 98 98 ...
##  $ experiencia    : int  144 36 381 190 138 67 114 0 115 244 ...
##  $ minoria        : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ directivo      : Factor w/ 2 levels "0","1": 2 1 1 1 1 1 1 1 1 1 ...

plot<-ggplot(datos3,aes(x =salario_inicial, y =salario_actual))+ #empezamos por dos simples
geom_point()
plot

plot<-ggplot(datos3,aes(x =salario_inicial, y =salario_actual))+ #empezamos por dos simples
geom_point(aes(size=educ, alpha=0.3)) #agregamos otra cuanti en tamaño
plot

plot<-ggplot(datos3,aes(x =salario_inicial, y =salario_actual))+ 
geom_point(aes(size=educ, alpha=0.3, color=sexo)) #agregamos una cuali en color
plot

#las tres cualitativas van a estar en sexo

plot<-ggplot(datos3,aes(x =salario_inicial, y =salario_actual))+ 
geom_point(aes(size=educ, alpha=0.3, color=sexo,shape=sexo))+ 
  facet_grid(sexo~.)
plot

#linea de regresion

plot<-ggplot(datos3,aes(x =salario_inicial, y =salario_actual))+ 
geom_point(aes(size=educ, alpha=0.3, color=sexo,shape=sexo))+ 
  facet_grid(sexo~.)+

geom_smooth(method=lm, se=FALSE, fullrange=TRUE)

plot

## `geom_smooth()` using formula 'y ~ x'

#hacer un zoom al grafico
#linea de regresion

plot<-ggplot(datos3,aes(x =salario_inicial, y =salario_actual))+ 
geom_point(aes(size=educ, alpha=0.3, color=sexo,shape=sexo))+ 
  facet_grid(sexo~.)+

geom_smooth(method=lm, se=FALSE, fullrange=TRUE)+
xlim(10000, 20000)

plot

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 95 rows containing non-finite values (stat_smooth).

## Warning: Removed 95 rows containing missing values (geom_point).

Ejercicios Regresión

Armando Valdés-López

24/12/2021

Ejempo 1

Importar y explorar datos 1

Visualizar la variables 1

Hacer el modelo 1

Comprobar los supuestos 1

Linealidad 1

Media de error 1

Normalidad resiuos 1

Homogeneidad de varianzas 1

Autocorelación 1

Multicolineidad 1

Hacer gráfico del modelo 1

Ejempo 2

Importar y explorar datos 1

Visualizar la variables 1

Hacer el modelo 1

Comprobar los supuestos 1

Linealidad 1

Media de error 1

Normalidad resiuos 1

Homogeneidad de varianzas 1

Multicolineidad 1

Autocorelación 1

Hacer gráfico del modelo 1

Ejempo 3

Importar y explorar datos 1

Visualizar la variables 1

Hacer el modelo 1

Comprobar los supuestos 1

Linealidad 1

Media de error 1

Normalidad resiuos 1

Homogeneidad de varianzas 1

Multicolineidad 1

Autocorelación 1

Hacer gráfico del modelo 1 hasta 3 cuai y 3 cuanti