1. Base de datos:

library(moments) #PARA CALCULAR EL SESGO Y OTROS ESTADÍSTICOS
library(dplyr) #MANIPULACIÓN DE DATOS
library(ggplot2) #GRÁFICOS ESTILIZADOS
library(readxl)
empleados <- read_excel("empleados.xls", 
sheet = "Respuestas",  range = "A1:D100")
View(empleados)
colnames(empleados)[c(1,2,3,4)] <- c("Edad","Estatura","Peso","Sexo")
empleados$Sexo <- as.factor(empleados$Sexo)
View(empleados)
str(empleados)
## Classes 'tbl_df', 'tbl' and 'data.frame':    99 obs. of  4 variables:
##  $ Edad    : num  20 18 19 19 21 18 20 18 19 18 ...
##  $ Estatura: num  178 168 194 159 177 180 180 168 190 187 ...
##  $ Peso    : num  82 87 94 62 78 53 62 68 82 79 ...
##  $ Sexo    : Factor w/ 2 levels "Hombre","Mujer": 1 1 1 2 1 1 2 1 1 2 ...
summary(empleados)
##       Edad          Estatura        Peso            Sexo   
##  Min.   :18.00   Min.   :159   Min.   : 51.00   Hombre:55  
##  1st Qu.:18.00   1st Qu.:171   1st Qu.: 67.00   Mujer :44  
##  Median :19.00   Median :178   Median : 72.00              
##  Mean   :20.52   Mean   :177   Mean   : 74.76              
##  3rd Qu.:22.00   3rd Qu.:182   3rd Qu.: 82.00              
##  Max.   :51.00   Max.   :200   Max.   :120.00
## [1] 4.163406
## [1] 8.213975
## [1] 13.11573

Genereración de diagramas de caja para variables continuas, descripción de resultados

TABLA DE FRECUENCIA POR INTERVALOS

##  [1] [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29)
##  [9] [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29)
## [17] [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29)
## [25] [18,29) [18,29) [18,29) [18,29) [18,29) [29,40) [18,29) [18,29)
## [33] [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29)
## [41] [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29)
## [49] [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29)
## [57] [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29)
## [65] [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29)
## [73] [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29)
## [81] [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29) [18,29)
## [89] [18,29) [18,29) [18,29) [18,29) [18,29) [40,51] [18,29) [18,29)
## [97] [18,29) [18,29) [18,29)
## Levels: [18,29) [29,40) [40,51]
## 
## [18,29) [29,40) [40,51] 
##      97       1       1

HISTOGRAMA

Edad media obtenida por sexo:

barplot(prop.table(table(empleados$Edad)),main="Edad")

-Calcule la correlación entre la variable dependiente y cada una de las variables explicativas (numéricas).

La correlación se define como:

En las estadística, la dependencia o asociación es una relación estadística, causal o no, entre dos variables aleatorias o dos conjuntos de datos. La correlación es cualquiera de una amplia clase de relaciones estadísticas que involucran dependencia, aunque en el uso común a menudo se refiere a la medida en que dos variables tienen una relación lineal entre sí.

Habitualmente, al iniciar un estudio de regresión lineal simple se suelen representar los valores de la variable dependiente y de la variable independiente de forma conjunta mediante un diagrama de dispersión para determinar si realmente existe una relación lineal entre ambas. Para realizar un diagrama de dispersión en R utilizaremos la orden plot

Compruebe si hay NA en el marco de datos.

Edad=empleados$Edad
Peso=empleados$Peso
plot(Edad,Peso)

Compruebe si hay NA en el marco de datos.

Estatura=empleados$Estatura
plot(Estatura,Peso)

Instalación y cargar los paquetes requeridos.

library(readr)
library(ggplot2)
library(corrplot) 
library(mlbench) 
library(Amelia) 
library(plotly) 
library(reshape2) 
library(caret)  
library(caTools) 
library(dplyr)

Comprobación si hay NA en el marco de datos.

missmap(empleados,col=c('yellow','black'),y.at=1,y.labels='',legend=TRUE)

colores fríos: correlación lineal directa, R aproximadamente 1 colores cálidos: correlación lineal inversa, R se aproxima a -1 *cuando R tiende a 0 no hay correlación lineal

Diagrama de correlación

corrplot(cor(select(empleados,-Sexo)))

corrplot(cor(select(empleados,-Sexo)), method = "number")

Gráfico de densidad con ggplot # visualizar la distribución de la variable indepeniente

empleados %>% 
    ggplot(aes(Estatura)) +
  stat_density() + 
  theme_bw()

empleados %>% 
    ggplot(aes(Edad)) +
  stat_density() + 
  theme_bw()

empleados %>%
  select(c(Edad, Estatura, Peso)) %>%
  melt(id.vars = "Peso") %>%
  ggplot(aes(x = value, y = Peso, colour = variable)) +
  geom_point(alpha = 0.7) +
  stat_smooth(aes(colour = "black")) +
  facet_wrap(~variable, scales = "free", ncol = 2) +
  labs(x = "Variable Value", y = "Median Peso Empleado (76)") +
  theme_minimal()

Correlación con las variables son: * Peso-Edad

cor(Peso,Edad)
## [1] 0.2553284

Peso-Estatura

cor(Peso,Estatura)
## [1] 0.5439601
cor.test (Peso,Estatura)
## 
##  Pearson's product-moment correlation
## 
## data:  Peso and Estatura
## t = 6.3846, df = 97, p-value = 5.922e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3882390 0.6694803
## sample estimates:
##       cor 
## 0.5439601

El p-valor asociado a este contraste es de 5.922e-09 < 0.05, por lo que rechazamos la hipótesis de que la correlación lineal entre estas dos variables sea 0.

cor.test (Peso,Edad)
## 
##  Pearson's product-moment correlation
## 
## data:  Peso and Edad
## t = 2.6009, df = 97, p-value = 0.01075
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.06099079 0.43101497
## sample estimates:
##       cor 
## 0.2553284

El p-valor asociado a este contraste es de 0.01075 < 0.05, por lo que rechazamos la hipótesis de que la correlación lineal entre estas dos variables sea 0.

2. Análisis ANOVA

Considere una variable categórica y realice un análisis ANOVA (como el revisado en clase), incluya resultados y conclusion al final Esta es la forma de pedir un ANOVA en R: Esta es continua en base a una variable continua y categórica

reg_lin_mul<-lm(empleados$Peso ~ empleados$Sexo+empleados$Edad+empleados$Estatura)
summary(reg_lin_mul)
## 
## Call:
## lm(formula = empleados$Peso ~ empleados$Sexo + empleados$Edad + 
##     empleados$Estatura)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.714  -6.911  -0.846   4.932  40.930 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -81.7705    24.1594  -3.385  0.00104 ** 
## empleados$SexoMujer  -2.3658     2.2160  -1.068  0.28840    
## empleados$Edad        0.6168     0.2643   2.333  0.02173 *  
## empleados$Estatura    0.8188     0.1343   6.099  2.3e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.77 on 95 degrees of freedom
## Multiple R-squared:  0.3458, Adjusted R-squared:  0.3252 
## F-statistic: 16.74 on 3 and 95 DF,  p-value: 8.275e-09

El modelo podría escribirse tal y como sigue: Peso=-81.7705-2.3658Sexo+0.6168Edad+0.8188Altura

Sexo=empleados$Sexo
fm = aov(lm(Peso ~ Sexo+Edad+Estatura))
summary(fm)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Sexo         1    592     592   5.103 0.02618 *  
## Edad         1    920     920   7.923 0.00593 ** 
## Estatura     1   4318    4318  37.192 2.3e-08 ***
## Residuals   95  11028     116                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En este caso, hemos encontrado un efecto significativo de las variables Sexo, Edad y Estatura (el valor de p ha sido menor de .05).

Elementos generados en el ANOVA:

names(fm)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"

¿Cuál es el valor crítico de F bajo la hipótesis nula con un nivel de significación alfa = 0.05? (Este valor nos delimitará la región de aceptación y rechazo)

Bajo la Ho el estadístico de contraste F se distribuye como una F de grados de libertad (I-1), (n-I) donde I es el número de grupos que disponemos y n el tamaño total de la muestral. Así obtenemos el cuantil buscado:

qf(0.05, 2-1, 99-2, lower.tail = F)
## [1] 3.939126

Valores del estadístico > 3.939126 sestarán incluidos en la región de rechazo. En nuestro caso los valors de F obtenidos para las 3 variables (5.103 Sexo,7.923 Edad y 37.192 Estatura) son mayores que el valor crítico obtenido.

¿Qué valor de la tabla ANOVA nos proporciona la varianza muestral común (pooled variance en inglés)? ¿Para qué es útil?

La raíz cuadrada de la media de los cuadrados del error, además de proporcionarnos una estimación de la varianza muestral de todos los datos, se utiliza en la obtención de los intervalos de confianza de las medias en cada uno de los grupos de interés.

Por ejemplo, este sería el intervalo de confianza de la media de los empleados con peso para el sexo Femenino, con un nivel de confianza del 72%:

media <- mean(empleados$Peso[empleados$Sexo =="Mujer"]) 
valor_t <- pt(0.05/2, 99 - 2) # total de elementos y n -1 de grupos
sp <- sqrt(116)  #desviación típica de la varianza muestral común (valor obtenido del anova)
ee  <- valor_t * (sp/ sqrt(99))  #error de estimación nos ayuda a sacar el intervalo de confianza el 99 es el nùmero de elementos
media
## [1] 72.02273

límite superior del intervalo de confianza de la media de empleados cuyo peso para sexo Mujer

media + ee 
## [1] 72.57472

límite inferior del intervalo de confianza de la media de empleados para sexo Mujer

media - ee 
## [1] 71.47073

Si hemos detectado diferencias significativas entre las medias de las poblaciones. ¿Sería posible saber cuáles son los grupos que generan estas diferencias?

a1 <- aov(empleados$Peso ~ empleados$Edad + empleados$Estatura+ empleados$Sexo)
intervals <- TukeyHSD(x=a1, 'empleados$Sexo', conf.level=0.95)
intervals
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = empleados$Peso ~ empleados$Edad + empleados$Estatura + empleados$Sexo)
## 
## $`empleados$Sexo`
##                   diff       lwr      upr    p adj
## Mujer-Hombre -2.288006 -6.614343 2.038332 0.296423
plot(intervals)

Validación del modelo ANOVA

A partir de los residuos del modelo comprobaremos si el modelo ANOVA es adecuado. Los supuestos que se deben cumplir son tres: independencia, homocedasticidad y normalidad.

Independencia

Mientras mas dispersos estan los residuos, menos correlacionados estan.

plot(fm$residuals)

Normalidad

Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos:

summary(fm$residuals)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -23.7142  -6.9114  -0.8464   0.0000   4.9322  40.9302

homostesidad = varianza

boxplot(fm$residuals)

hist(fm$residuals)

qqnorm(fm$residuals) 
qqline(fm$residuals)

El test de Shapiro-Wilk indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (normalidad de los residuos)

Si p>=0.05 podemos afiramar que son normales

shapiro.test(fm$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fm$residuals
## W = 0.9109, p-value = 5.267e-06
Homocedasticidad

Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos:

En este grafico no nos fijamos en la linea sino en los anchos. Varianza se obtiene a travez de los residuos

boxplot(fm$residuals~empleados$Sexo+empleados$Edad, col = c("yellow", "blue"))  

desviaciones <- tapply(fm$residuals, empleados$Sexo, sd)

Comparando la desviación máxima con la mínima obtenemos una orientación sobre la falta de homocedasticidad (>2 aproximadamente)

Si la diferencia entre el maximo y el minimo es >2 la varianza no es igual.

max(desviaciones) / min(desviaciones)    
## [1] 1.253551

El test de Bartlett indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (las varianzas son iguales).

bartlett.test(fm$residuals ~  empleados$Sexo)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  fm$residuals by empleados$Sexo
## Bartlett's K-squared = 2.3595, df = 1, p-value = 0.1245
Kruskal-Wallis y pruebas post-hoc

Es un anova multivariante ¿Qué hipótesis contrasta el test ANOVA?

Ho: las medias son iguales en todas las poblaciones

Ha: hay alguna media distinta

¿Qué hipótesis contrasta la prueba de Kruskal-Wallis?

Ho: la variable respuesta es la misma en todas las poblaciones valoradas

Ha: la variable respuesta es mayor en alguna de las poblaciones

Cuando no se cumplen las hipótesis exigidas por el modelo ANOVA, es posible utilizar la prueba no paramétrica Kruskal-Wallis: ¿hay diferencias significativas entre las poblaciones?

kruskal.test(empleados$Peso, empleados$Sexo)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  empleados$Peso and empleados$Sexo
## Kruskal-Wallis chi-squared = 3.5706, df = 1, p-value = 0.05881

Indica cuál es el estadístico de contraste, los grados de libertad, el p-valor correspondiente y cuál sería el valor crítico que definiría las regiones de aceptación y rechazo con un nivel de significación alfa = 0.05.

El chicuadrado es el obtenido 3.5706

Bajo la Ho el estadístico de contraste H del test de Kruskal-Wallis se distribuye como una Chi-cuadrado de grados de libertad (I-1) (donde I es el número de grupos que disponemos). Así obtenemos el cuantil buscado: El chi teorico es 5.99

qchisq(0.05, 2-1, lower.tail = F)
## [1] 3.841459

Entonces se resta el chi-obtenido -chi teorico entonces se rechza la nula y se acepta la alternativa.

Valores del estadístico > 3.841459 estarán incluidos en la región de rechazo. En nusetro caso 3.5706 es menor que el valor crítico obtenido.

Si transformáramos los datos de la variable respuesta, utilizando logaritmos y después aplicáramos el test de KrusKal-Wallis al logaritmo del número de insectos atrapados, ¿variarían los resultados del test estadístico?

kruskal.test(log(empleados$Peso), empleados$Sexo) 
## 
##  Kruskal-Wallis rank sum test
## 
## data:  log(empleados$Peso) and empleados$Sexo
## Kruskal-Wallis chi-squared = 3.5706, df = 1, p-value = 0.05881

Los resultados son exactamente los mismos. No se producen variaciones porque el test de Kruskal-Wallis trabaja sobre rangos, es decir, sobre ordenaciones de los valores de la variable en cada uno de los grupos. Aunque realicemos una transformación logarítmica, el orden entre los valores de la variable se mantiene y por lo tanto la transformación no afecta a los resultados del test.

Si hemos detectado diferencias significativas en la variable respuesta para las distintas poblaciones. ¿Sería posible saber cuáles son los grupos que generan estas diferencias?

La libreria PCMC resume todo lo realizado anteriormente de anova y llega a la conclusion.

library(PMCMR)
posthoc.kruskal.nemenyi.test(empleados$Peso, empleados$Sexo, method = "Chisq")
## 
##  Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
##                    with Tukey-Dist approximation for independent samples 
## 
## data:  empleados$Peso and empleados$Sexo 
## 
##       Hombre
## Mujer 0.059 
## 
## P value adjustment method: none

3. Construcción del modelo y predicción

3.1 Considerando los cálculos anteriores genere el modelo de regresión lineal múltiple que mejor se ajuste a los datos.

El modelo de regresión lineal general en R:

Train Y Test Data

Permite dividir los datos en entrenamento y prueba de datos utilizando la biblioteca caTools. Lets split the data into train and test data using caTools library.

# establecer una semilla
set.seed(123)

#Seccionar los datos , `split ()` asigna un booleano a una nueva columna basada en el SplitRatio especificado.

split <- sample.split(empleados,SplitRatio =0.75)

trainemp <- subset(empleados,split==TRUE)
testemp <- subset(empleados,split==FALSE)


# train <- select(train,-b)
# test <- select(test,-b)
Entrenando nuestro modelo

Vamos a construir nuestro modelo teniendo en cuenta que Edad, Altura y Sexo son influyentes en la variable objetivo (Peso).

modelemp <- lm(empleados$Peso ~ empleados$Edad + empleados$Estatura+ empleados$Sexo, data =trainemp)
summary(modelemp)
## 
## Call:
## lm(formula = empleados$Peso ~ empleados$Edad + empleados$Estatura + 
##     empleados$Sexo, data = trainemp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.714  -6.911  -0.846   4.932  40.930 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -81.7705    24.1594  -3.385  0.00104 ** 
## empleados$Edad        0.6168     0.2643   2.333  0.02173 *  
## empleados$Estatura    0.8188     0.1343   6.099  2.3e-08 ***
## empleados$SexoMujer  -2.3658     2.2160  -1.068  0.28840    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.77 on 95 degrees of freedom
## Multiple R-squared:  0.3458, Adjusted R-squared:  0.3252 
## F-statistic: 16.74 on 3 and 95 DF,  p-value: 8.275e-09

->Un modelo es aceptable si: a) Las variables son estadisticamente significativas b) R al cuadrado ajustado >= 0.5 (mide la calidad del modelo) c) Residuo N(0,1)ajusten a una distribuciòn normal.(media=0, Varianza=1, tiene que pareceer una campana de Gaus) Residuos no correlacionados (NO SIGUEN UN PATRON)

Se procede a retirar la variable Sexo:

empleados$Sexo(PR(>t))=0.2884 indica la variable que es NO estadisticamente signitifcativa para el modelo entre mas pequeña es este valor es mas signitictiva se excluye.

# para mejorar el modelo una vez que vinmos que sexo no es significativa para el modelo
# esto se lo realiza al finalizar toda la evaluacion 
modelemp <- lm(empleados$Peso~ empleados$Edad + empleados$Estatura -1 , data=trainemp)
summary(modelemp)
## 
## Call:
## lm(formula = empleados$Peso ~ empleados$Edad + empleados$Estatura - 
##     1, data = trainemp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.720  -6.665  -2.443   6.349  41.953 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## empleados$Edad      0.52015    0.27656   1.881    0.063 .  
## empleados$Estatura  0.36310    0.03267  11.115   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.45 on 97 degrees of freedom
## Multiple R-squared:  0.9777, Adjusted R-squared:  0.9772 
## F-statistic:  2125 on 2 and 97 DF,  p-value: < 2.2e-16
La ecuaciòn del modelo encontrado
##     empleados$Edad empleados$Estatura 
##          0.5201547          0.3630955
Visualizando nuestro modelo

Permite visualizar nuestro modelo de regresión lineal trazando los residuos. La diferencia entre el valor observado de la variable dependiente (y) y el valor predicho (y) se denomina residual (e).

res <- residuals(modelemp)

# Convertir residuos en un DataFrame
# 
res <- as.data.frame(res)
ggplot(res,aes(res)) +  geom_histogram(fill='blue',alpha=0.5)

Eje x Son los datos ajustados y eje y la relacion. en la mitad 373m 3700 son valores atìpicos

plot(modelemp)

##### Predicciones Probemos nuestro modelo prediciendo en nuestro conjunto de datos de prueba. predict ajusta el modelo con los datos test

testemp$predicted.Peso<- tail(predict(modelemp,testemp), n = 24) 
pl1 <-testemp %>% 
  ggplot(aes(Peso,predicted.Peso)) +
  geom_point(alpha=0.5) + 
  stat_smooth(aes(colour='black')) +
  xlab('Actual value of peso') +
  ylab('Predicted value of peso')+
  theme_bw()

#ggplotly(pl1)
Evaluemos nuestro modelo

usando Root Mean Square Error, una medida estandarizada de cuán lejos estábamos de los valores reales con nuestros valores predichos. Si rmse es cercano a 1 se puede consdierar bueno.

error <- testemp$Peso-testemp$predicted.Peso
rmse <- sqrt(mean(error)^2)
rmse
## [1] 2.762652

El Root Mean Square Error (RMSE) para nuestro modelo es 2.762652 lo cual no es considerado como bueno ya que es > 1 los resultados pueden mejorarse aún más utilizando la extracción de variables y entrenando el modelo.