Parcial 2 Modelos Estadísticos I

library(car)
library(nortest)
library(MASS)
library(latex2exp)
library(ggplot2)
library(tibble)
library(GGally)
library(MVN)
library(pracma)
library(dplyr)
library(ggpubr) # Para Realizar QQplots con función ggqqplot(model$resi)
library(QuantPsyc) # Para verificar multinormalidad
library(tidyverse)
library(devtools) # Instalar paquetes que no están en el CRAN
library(alr3)

1. Considera la base de datos LED.csv. Se cree que la esperanza de vida esta en esta en función del año en cuestión, el porcentaje de ingesta de alcohol y el estado del país(desarrollado o en desarrollo). Ademas se sospecha que existe interacción entre el año y el estado del país. Este tipo de interacción, digamos que entre dos regresores \(x_1\) y \(x_2\), se ingresa en el modelo de regresión considerando el producto de tales variables \(x_1x_2\).

# Carga del conjunto de datos
LED<-read.csv("/Users/danimathud/Downloads/Led.csv",na.strings = "NA")
LED[,1]<-as.factor(LED$Country)
LED<- LED %>% filter(LED$Alcohol!="NA")
LED<- LED %>% filter(LED$Life.expectancy!="NA")
for(i in 1:nrow(LED)){
  if (LED[i,3]=="Developing"){
    LED[i,3]<-1
  }
   else{
     LED[i,3]<-0
   } 
}
LED$Status<-as.numeric(LED$Status)

a)

Ajusta el modelo con la interacción entre el año y el estado del país.

Para poder visualizar, es conveniente convertir la variable estatus a numérica, donde “1” hace referencia a estado Developing y “0” a estado No Developing.

muestreo1<-sample(nrow(LED),size=300)
LedModel<-lm(LED$Life.expectancy~LED$Alcohol+LED$Year+LED$Status+LED$Status*LED$Year)
plot(LED$Year,LED$Life.expectancy,pch=19,main="Grafico Dispersión Año vs Esperanza Vida",xlab="Año",ylab="Esperanza de Vida") 

plot(LED$Status, LED$Life.expectancy,pch=19, main="Dispersión Estatus vs Esperanza de vida",xlab="Estatus",ylab="Esperanza de Vida")

plot(LED$Alcohol,LED$Life.expectancy,pch=19,main="Dispersión Alcohol vs Esperanza de Vida",xlab="Alcohol",ylab="Esperanza de Vida")

summary(LedModel)

Call:
lm(formula = LED$Life.expectancy ~ LED$Alcohol + LED$Year + LED$Status + 
    LED$Status * LED$Year)

Residuals:
    Min      1Q  Median      3Q     Max 
-32.936  -4.910   1.353   6.206  17.976 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)         -607.02598  169.26080  -3.586 0.000341 ***
LED$Alcohol            0.45921    0.04768   9.631  < 2e-16 ***
LED$Year               0.33962    0.08432   4.028 5.79e-05 ***
LED$Status          -110.14343  186.68031  -0.590 0.555232    
LED$Year:LED$Status    0.05031    0.09301   0.541 0.588582    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.085 on 2730 degrees of freedom
Multiple R-squared:  0.2844,    Adjusted R-squared:  0.2834 
F-statistic: 271.3 on 4 and 2730 DF,  p-value: < 2.2e-16

Viendo los diagramas dispersión en cualquiera de los tres podemos observar poca relación lineal entre las variables predictora y la respuesta. Sin embargo analicemos una manera adecuada para ajustar un modelo lineal.

plot(LedModel,which = 1,pch=19,sub.caption = "")

El gráfico de residuales nos da indicio de dos posibles subconjuntos, puesto que visualmente se viola el supuesto de homocedasticidad. Estos subconjuntos corresponden a países no desarrollados y desarrollados. Puesto que la esperanza de vida en países desarrollados son mayores, con lo cual podríamos pensar en que quizás resulte mas conveniente separar el conjunto de datos original en dos subconjuntos independientes. Por otro lado del gráfico QQ, vemos que por lo menos considerando ambos estatus la normalidad no se esta cumpliendo tampoco como se ve en la siguiente gráfica

ggqqplot(LedModel$residuals,title = "QQ Plot")

Claramente se puede ver la no normalidad en los errores, ya que gran parte de los residuos se sale de las bandas de confianza.

Puesto que separamos por estatus, ya no tiene mucho sentido considerarla como regresora, de modo que consideremos alcohol y años como regresores. Ajustando un modelo para cada estatus obtenemos lo siguiente:

EstatusDev<-LED %>% filter(LED$Status==1)
EstatusNoDev <- LED %>% filter(LED$Status==0)
LedModelDev <- lm(EstatusDev$Life.expectancy~EstatusDev$Year+EstatusDev$Alcohol)
LedModelNoDev<-lm(EstatusNoDev$Life.expectancy~EstatusNoDev$Year+EstatusNoDev$Alcohol)

Para estatus desarrollado, observamos un gráfico de residuales poco homogéneo, pues lo residuales se agrupan en una región y por otro lado del QQ vemos que tampoco hay normalidad, de hecho los residuales se quedan fuera de la banda de confianza, y este hecho se confirma aun mas con la prueba de Shapiro la cual rechaza la normalidad

plot(LedModelDev,which=1,pch=19,sub.caption = "")

ggqqplot(LedModelDev$residuals)

shapiro.test(LedModelDev$residuals)

    Shapiro-Wilk normality test

data:  LedModelDev$residuals
W = 0.94791, p-value < 2.2e-16

Ahora considerando estatus no desarrollado, obtenemos también la no normalidad y no homocedasticidad en los errores, puesto que al igual que ocurrió en estatus desarrollado, los residuos se están concentrando de manera no homogénea sobre una franja horizontal. Podemos concluir que el modelo lineal no es el adecuado para este conjunto de datos, pues en ambos casos no se satisfacen ninguno de los supuestos.

plot(LedModelNoDev,which=1,pch=19,sub.caption = "")

ggqqplot(LedModelNoDev$residuals)

shapiro.test(LedModelNoDev$residuals)

    Shapiro-Wilk normality test

data:  LedModelNoDev$residuals
W = 0.96997, p-value = 2.114e-08

b)

Ajusta un nuevo modelo sin considerar los regresores que resultaron no significativos en el inciso anterior.

En el resumen del inciso anterior se puede ver que las variables estatus e interacción tenian p-valores grande, con lo cual es plausible que los sus respectivos \(\beta\)’s sean cero, de modo que podemos considerar un modelo sin estas predictoras, sin embargo también pudimos observar que el conjunto de datos en general si tendía a separarse en dos grupos distintos, lo cual seguiremos observando a pesar de no tener en cuenta al estatus. De modo que es de esperar obtener resultados similares a los obtenidos en el inciso anterior, es decir que es posible que de nuevo el modelo lineal no sea adecuado para este conjunto en general retirando las predictoras ya mencionadas.

LedModelGen<-lm(LED$Life.expectancy~LED$Year+LED$Alcohol)
plot(LedModelGen,which=1,pch=19,sub.caption = "")

ggqqplot(LedModelGen$residuals)

shapiro.test(LedModelGen$residuals)

    Shapiro-Wilk normality test

data:  LedModelGen$residuals
W = 0.95025, p-value < 2.2e-16

Nuevamente al igual que se observo en el inciso anterior, el comportamiento de los residuales es no homogéneo, se esta concentrando y no se esta dispersando de manera uniforme, con lo cual homocedastidad no se esta cumpliendo, por otro lado del QQ, podeos ver que gran parte de los residuales se están quedando fuera de las bandas de confianza. En conclusión considerando o no los regresores mencionados, hay violación de supuestos y por lo tanto el modelo lineal tampoco es adecuado en este caso.

# pureErrorAnova(modelo1)
# ggcorr(maria[,1:3],label=T)
# confidenceEllipse(fit)

c)

#Residuales Estudentizados
plot(LedModelGen$fitted.values,rstudent(LedModelGen),pch=19,main="Estudentizados vs Ajustados")
abline(h=0,lty=2,col="red")

plot(LED$Year,rstudent(LedModelGen),pch=19,col="blue",main="Regresor(Año) vs Estudentizados",ylab="Residuales estudentizados",xlab="Año")

plot(LED$Alcohol,rstudent(LedModelGen),pch=19,main="Regresor(Alcohol) vs Estudentizados",xlab = "Alcohol",ylab="Residuales estudentizados")

# Regresión Parcial
M1<-lm(LED$Life.expectancy~LED$Alcohol)
M2<-lm(LED$Alcohol~LED$Year)
plot(M1$residuals,M2$residuals,pch=19,main="Regresión Parcial- (Alc~Año) vs (EspVida~Alc)",xlab="Residuales Alcohol en función del año",ylab="Residuales Esp. de vida en función de Alcohol")

De la primera gráfica de estudentizados versus ajustados, podemos observar de entrada un comportamiento poco homogéneo de los residuales, pues como podemos observar, estos se están concentrando hacia la izquierda, como en forma de megáfono, lo que indica que igualdad de varianza no se esta cumpliendo. Por lo tanto empezando ya hay indicios de que el modelo lineal que intentamos ajustar no es adecuado para el conjunto de datos, o por lo menos no para describir esperanza de vida en función del año y consumo de alcohol. Veamos los demás gráficos de residuales, para afirmar con mas seguridad la no adecuación del modelo.

De la segunda gráfica, podemos observar el comportamiento de los residuales versus el regresor año. Evidenciamos nuevamente un comportamiento no homogéneo, pues de hecho es casi inmediato por por las tendencia vertical de los datos para cada año, de igual forma aunque el modelo poco se ajusta, recordemos que si el modelo se ajustara bien en principio el gráfico que esperaríamos tener debería no diferir mucho del original, lo cual claramente aquí no se tiene. En conclusión para esta gráfica se sigue viendo que el modelo quizás no sea el adecuado.

En la tercera gráfica comparamos ahora el regresor alcohol versus los residuales estudentizados, aunque aquí también se evidencia la no homocedasticidad, sin embargo si podemos notar que el gráfico tiene mucha similitud al del modelo original observado en el gráfico 1, lo que al menos si nos permite suponer que este regresor si tiene influencia en el modelo general a pesar de que no sea lineal.

Por ultimo comparando el gráfico de residuales para regresión parcial, observamos que el modelo no es adecuado, puesto que de serlo, deberíamos observar los residuales con un comportamiento lineal, es decir una dispersión a la cual podría ajustar un modelo lineal y en lo que vemos esto no pareciera.

ResiPress<-LedModelGen$residuals/(1-influence(LedModelGen)$hat)
syy <- sum((LED$Life.expectancy-mean(LED$Life.expectancy))^2)
EstPress<-sum(ResiPress^2)
p_2<-1-EstPress/syy

En conclusión el modelo lineal no es adecuado. Adicional, su poder predictivo es bastante pobre, pues observamos un \(P^2\) de 0.1966, lo cual esta muy bajo.

d)

Ajusta el modelo con el año como único predictor y realiza una prueba de bondad de ajuste para dicho modelo

M3<-lm(LED$Life.expectancy~LED$Year)
pure.error.anova(M3)

Warning: 'pure.error.anova' is deprecated.
Use 'pureErrorAnova' instead.
See help("Deprecated") and help("alr3-deprecated").

Analysis of Variance Table

Response: LED$Life.expectancy
               Df Sum Sq Mean Sq F value Pr(>F)    
LED$Year        1   6799  6799.3 76.3004 <2e-16 ***
Residuals    2733 242606    88.8                   
 Lack of fit   14    309    22.1  0.2478 0.9979    
 Pure Error  2719 242297    89.1                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Podemos ver del resumen la salida el p-valor (Lack of fit) de 0.9979, lo que por lo menos con prueba de bondad de ajuste, no se rechaza que el modelo sea lineal.

e)

Determina si tiene sentido pensar que el tiempo de vida esperado en 2022 estará entre 80 y 85 años

ResiPress2<-M3$residuals/(1-influence(M3)$hat)
syy_2 <- sum((LED$Life.expectancy-mean(LED$Life.expectancy))^2)
EstPress_2<-sum(ResiPress2^2)
p_cuadrado<-1-EstPress_2/syy_2

Tenemos que el poder predictivo es demasiado bajo (0.0254), con lo cual de entrada esperaríamos que no tuviera tanto sentido pensar que la esperanza de vida este entre 80 y 85 años, por otro lado, tengamos en cuenta que nuestro conjunto de datos considera años como máximo de los años considerados en nuestro ajuste es 2015, con lo cual considerar valores por encima de esto, generaría mas error de predicción. Sin embargo veamos si la intuición es correcta viendo un intervalo de confianza

valort<-qt(0.025,df=nrow(LED)-2)
MSS_res<-sum((M3$residuals)^2)/(nrow(LED)-2)
x0<-2022
y0<--661.23+0.36392*x0
sxx<-sum((LED$Year-mean(LED$Year))^2)
z<-(x0-mean(LED$Year))^2/sxx


ColInf<-y0+valort*sqrt(MSS_res*(1+1/nrow(LED)+z))
ColInf

[1] 56.09807

ColSup<-y0-valort*sqrt(MSS_res*(1+1/nrow(LED)+z))
ColSup

[1] 93.13441

Realizando el calculo de los intervalos de predicción, vemos que la intuición no fue correcta y de hecho las esperanza de vida de entre 80 y 85 son valores en el intervalo de predicción, y por lo tanto si es plausible que para 2022 la esperanza de vida este entre 80 y 85 a pesar de que la intuición pareciera indicar lo contrario.

2. La base de datos ``Salarios.csv’’ contiene información anual del porcentaje de mujeres empleadas en cierto país. Con esta base realiza lo siguiente:

Salarios<-read.csv("/Users/danimathud/Downloads/Salarios.csv")

a)

Ajusta un modelo de regresión simple utilizando el año como variable predictora (incluye la validación de supuestos necesaria).

SalModel<-lm(Salarios$Porcentaje.de.trabajadoras.con.pago.por.horas~Salarios$Anio)
summary(SalModel)

Call:
lm(formula = Salarios$Porcentaje.de.trabajadoras.con.pago.por.horas ~ 
    Salarios$Anio)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.39847 -0.12824 -0.06107  0.12481  0.44198 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -266.88626   38.01103  -7.021 6.18e-05 ***
Salarios$Anio    0.16565    0.01915   8.648 1.18e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2473 on 9 degrees of freedom
Multiple R-squared:  0.8926,    Adjusted R-squared:  0.8807 
F-statistic: 74.79 on 1 and 9 DF,  p-value: 1.182e-05

Viendo el gráfico de dispersión, podemos observar un poco de comportamiento lineal, verifiquemos los supuestos para confirmar si el modelo lineal se ajusta bien al conjunto de datos.

plot(SalModel$fitted.values,rstudent(SalModel),pch=19,cex=1.2,panel.smooth(SalModel$fitted.values,rstudent(SalModel),col.smooth = "blue"),xlab="Ajustados",ylab="Residuales Estudentizados",main = "Ajustados vs Gráfico Residuales Estud.  ")
abline(h=0,lty=2,col="red")

# Separación de residuales
Aleatorios<-sample(nrow(Salarios),size =floor(nrow(Salarios)*0.6))
Grupo1<-SalModel$residuals[Aleatorios]
Grupo2<-SalModel$residuals[-Aleatorios]
Etiquetas1<-rep("Grupo1",floor(nrow(Salarios)*0.6))
Etiquetas2<-rep("Grupo2",nrow(Salarios)-floor(nrow(Salarios)*0.6))
Etiquetas<-c(Etiquetas1,Etiquetas2)
ResidualesAgrup<-c(Grupo1,Grupo2)
ResidualesAgrup<-data.frame(ResidualesAgrup,Etiquetas)
as.factor(ResidualesAgrup$Etiquetas)

 [1] Grupo1 Grupo1 Grupo1 Grupo1 Grupo1 Grupo1 Grupo2 Grupo2 Grupo2 Grupo2
[11] Grupo2
Levels: Grupo1 Grupo2

# Verificación de Homocedasticidad
leveneTest(ResidualesAgrup$ResidualesAgrup,group=ResidualesAgrup$Etiquetas)

Warning: ResidualesAgrup$Etiquetas coerced to factor.

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  1.8869 0.2028
       9               

Viendo el resumen del modelo, rechazamos que \(\beta_0,\beta_1\) son cero por tener p-valores inferiores a \(0.05\).

Por otro lado, haciendo la validación de supuestos podemos notar primeramente que del gráfico de residuos, parece indicar homocedasticidad en los errores, si observamos la curva de ajuste(color azul) podemos notar una curva un “poco horizontal”, y además muy próximos a la linea punteada, con lo cual al menos visualmente observamos media cero en los residuales. Complementando mas la validación de homocedasticidad, usando la prueba de Levene, podemos observar un p-valor de 0.2837, con lo cual considerando un nivel de significancia \(\alpha=0.05\), tenemos que no hay evidencia estadística suficiente para concluir que los residuales no tienen igual varianza, lo cual es consecuente con lo observado previamente en el gráfico de residuales.

# Verificación de Normalidad en los errores
shapiro.test(SalModel$residuals)

    Shapiro-Wilk normality test

data:  SalModel$residuals
W = 0.96828, p-value = 0.8686

ggqqplot(SalModel$residuals,color = "blue",ylab="Muestrale",xlab="Teoricos",title = "Grafico QQ")

Por ultimo, verifiquemos la normalidad en los errores. Del gráfico QQ, observamos que los residuales se mantienen en las bandas de confianza, y están muy cercanos a la recta, con lo cual de inicio; visualmente podemos verificar que los residuos siguen la misma distribución y además con la prueba de Shapiro-Wilk obtenemos un p-valor de 0.8686 con lo cual no rechazamos la normalidad en los residuales. En conclusión vemos que se satisfacen todos los supuestos de regresión simple, y por tanto resulta adecuado ajustar una recta a los datos. El modelo ajustado es:

plot(Salarios$Anio,Salarios$Porcentaje.de.trabajadoras.con.pago.por.horas,pch=19,cex=1.2,xlab="Años",ylab="Porcentaje de Trabajadoras")
abline(SalModel,col="red")

b)

Realiza un análisis de residuales para determinar si existen outliers.

# Residuales estudentizados
plot(SalModel$fitted.values,rstudent(SalModel),pch=19,cex=1.2,ylab="Estudentizados",xlab="Ajustados",main="Ajustados vs Residuales estudentizados",ylim=c(-2.5,2.5))
points(c(SalModel$fitted.values[2],SalModel$fitted.values[9]),c(rstudent(SalModel)[2],rstudent(SalModel)[9]),col="red",pch=19,cex=1.2)

abline(h=0,col="red",lty=2)
legend(62.45,2.5,legend = c("Posibles Atipicos"),fill="red")

Como podemos observar del gráfico de residuales estudentizados, tenemos dos residuales que están más alejados con respecto al resto(los de color rojo), así que retiremoslos del modelo con con el fin de observar, si el modelo varía considerablemente y de esta manera podemos concluir que posiblemente se trate de valores atípicos.

SalariosSinOutliers<-Salarios[-c(which.min(rstudent(SalModel)),which.max(rstudent(SalModel))),]
SalModel2<-lm(SalariosSinOutliers$Porcentaje.de.trabajadoras.con.pago.por.horas~SalariosSinOutliers$Anio)

plot(SalModel2$fitted.values,rstudent(SalModel2),pch=19,cex=1.2,panel.smooth(SalModel2$fitted.values,rstudent(SalModel2),col.smooth = "blue"),xlab="Ajustados",ylab="Residuales Estudentizados",main = " Ajust. vs Res_stud. Quitando 'Outliers' ")
abline(h=0,lty=2,col="red")

#legend(62,0.5,legend = c("azul","rojo"), fill=c("red","blue"))

Como podemos ver del gráfico de residuales estudentizados retirando las observaciones mencionadas, podemos ver que los residuales se dispersan de una manera un poco mas homogénea que como lo hacían los residuales estudentizados originales en los cuales si considerábamos esas observaciones, de modo que es posible que estas observaciones se traten de datos atípicos.

c)

Determina si los posibles valores atípicos hallados en el inciso anterior son valores influyentes.

Al igual que el inciso a) validemos los supuestos. En el inciso b) verificamos que tanto varío el modelo mediante el gráfico de residuales estudentizados, concluimos que había mejoría en la dispersión de los residuales pues se veía una dispersión mas homogénea el modelo retirando estos Outliers, sin embargo veamos que tanta variabilidad tienen estas observaciones en el modelo general, verificando el resumen y el QQ plot.

ggqqplot(SalModel2$residuals,color = "blue",ylab="Muestrales",xlab="Teoricos",title = "Grafico QQ, quitando posibles influyentes")

summary(SalModel2)

Call:
lm(formula = SalariosSinOutliers$Porcentaje.de.trabajadoras.con.pago.por.horas ~ 
    SalariosSinOutliers$Anio)

Residuals:
   Min     1Q Median     3Q    Max 
-0.180 -0.070 -0.035  0.155  0.175 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)              -225.9100    23.4230  -9.645 2.71e-05 ***
SalariosSinOutliers$Anio    0.1450     0.0118  12.286 5.43e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1396 on 7 degrees of freedom
Multiple R-squared:  0.9557,    Adjusted R-squared:  0.9494 
F-statistic: 150.9 on 1 and 7 DF,  p-value: 5.427e-06

summary(SalModel)

Call:
lm(formula = Salarios$Porcentaje.de.trabajadoras.con.pago.por.horas ~ 
    Salarios$Anio)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.39847 -0.12824 -0.06107  0.12481  0.44198 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -266.88626   38.01103  -7.021 6.18e-05 ***
Salarios$Anio    0.16565    0.01915   8.648 1.18e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2473 on 9 degrees of freedom
Multiple R-squared:  0.8926,    Adjusted R-squared:  0.8807 
F-statistic: 74.79 on 1 and 9 DF,  p-value: 1.182e-05

Podemos observar que al igual que en el inciso anterior, se satisfacen los supuestos del modelo lineal, puesto que por un lado del QQplot, los residuales se mantienen dentro de las bandas de confianza, lo que nos indica igual distribución como era de esperar, por otro lado con el gráfico de residuos visto en el inciso b) observamos que los residuales se mantienen en un comportamiento “decente” lo que indica homocedasticidad, y la normalidad se mantiene por la prueba de Shapiro Test donde se obtuvo un p-valor de 0.139. Asi que el modelo se comporta bien extrayendo esas dos observaciones; por otro lado podemos observar del resumen un aumento en el valor del R ajustado para el modelo extrayendo las observaciones, y observando los \(\beta\)’s en ambos modelos vemos una poca variación, con lo cual en conclusión los datos no son suficientemente influyentes en modelo original. Lo cual podemos observar en la grafica donde comparamos modelos muy variación con la obtenida en un inicio.

shapiro.test(SalModel2$residuals)

    Shapiro-Wilk normality test

data:  SalModel2$residuals
W = 0.87498, p-value = 0.139

plot(Salarios$Anio,Salarios$Porcentaje.de.trabajadoras.con.pago.por.horas,pch=19,cex=1.2,xlab="Años",ylab="Porcentaje de Trabajadoras",main="Comparando modelos")
abline(SalModel,col="red")

abline(SalModel2,col="blue")
legend(1989,61.2,legend=c("Con Outliers","Sin Outliers"),fill=c("red","blue"))