Examen Final Regresión Avanzada

1 Ejercicio 1

El Departamento de Nutrición de la Ciudad desea estudiar el efecto del tipo de harina y el porcentaje de azúcar sobre el endulzamiento de las galletitas que prepara para los comedores escolares. Los siguientes datos corresponden a la densidad específica de muestras de galletitas preparadas con las distintas combinaciones posibles. Para cada combinación se realizaron tres ensayos independientes.

1. Plantear las hipótesis de interés.
2. Realizar las pruebas correspondientes.
3. Analizar la validez de los resultados de los contrastes.
4. Concluir en el contexto del problema y efectuar recomendaciones.

2 Ejercicio 2

El siguiente conjunto de datos ventascigarrillos.xlsx corresponde a los datos de los 50 estados de EEUU y las variables regitradas son las siguientes:

• Age Edad mediana de las personas que viven en el Estado.
• HS Porcentaje de personas con 25 años o más que tienen su nivel de educación secundaria concluido.
• Income Ingreso personal per cápita.
• Black Porcentaje de individuos de raza negra.
• Female Porcentaje de mujeres.
• Price Precio promedio del paquete de cigarrillos.
• Sales Número de paquetes promedio vendido per cápita.

Se requiere construir una ecuación de regresión que vincule el promedio de paquetes vendidos de cigarrillos per cápita (Sales) con las diversas variables económicas y socio demográficas disponibles en la basse (en caso que sea útiles para predecir a la respuesta de interés). Para ello sugerimos tener en consideración los siguientes lineamientos:

2.1 Dataset

library(readxl)
library(ggplot2)
library(plotly)
library(dplyr)
library(tibble)
library(readxl)
library(dplyr)
library(GGally)
library(Hmisc)
library(corrplot)
library(PerformanceAnalytics)
library(robustbase); library(MASS)
library(olsrr); library(car); library(quantreg)

data <- read.csv("C:/Users/ander/Desktop/1. Ciencias de Datos/4. Regresión avanzada/9. Examen 11042023/ventascigarrillos.csv", sep =";")

head(data)

2.2 Qué nivel de asociación lineal presentan estas variables con la variable respuesta Sales?

Un primer acercamiento para estudiar el nivel de asociación lineal entre las variables, puede ser a través de la matriz de correlación. En ella se observa como Sales tiene una correlación negativa y significativa de -0.30 con Price y positiva con Income de un 0.33.

Para un posible modelo lineal multiple estas serian las variables que mejor explicarian a Sales

data <-data[,-c(1)]
chart.Correlation(data, histogram = F)

2.3 Utilizando un método de selección automática elegir el numero de variables y cuáles son para construir un modelo en el que no se considere interacción.

Se utiliza una validación combinada de modelos, en donde se evaluan dichos modelos en función de la cantidad de regresores y se elige el mejor por el criterio de AIC.

library(olsrr)
lm <- lm(Sales~ ., data=data)
k_best <- ols_step_best_subset(lm)
k_best

El mejor modelo resulto ser el tercero ya que cuenta con el mayor R2 adj y el menor AIC. Por lo que queraría:

modelo1 <- lm(Sales~Age + Income + Price, data=data)
summary(modelo1)

## 
## Call:
## lm(formula = Sales ~ Age + Income + Price, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -50.490 -14.336  -5.164   7.175 128.880 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 64.521307  62.752062   1.028  0.30923   
## Age          4.156318   2.222385   1.870  0.06783 . 
## Income       0.019298   0.006963   2.772  0.00802 **
## Price       -3.407476   1.008475  -3.379  0.00149 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 27.91 on 46 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3017, Adjusted R-squared:  0.2562 
## F-statistic: 6.626 on 3 and 46 DF,  p-value: 0.0008171

Si bien resulto el mejor modelo respecto a las variables disponibles, el R2 solo logra explicar un 30% de la variación de la variable Sales. Las variables mas significativas dentro del modelo resultan ser Income y Price cuyo p-value es menor al nivel de confianza del 0.05.

2.4 Agregue en caso de que tenga sentido alguna interacción al modelo.

Al agregar una interacción entre las variables Female y Price se ve como mejora el R2 adj y esta interaccion resulta significativa para el modelo

modelo2 <- lm(Sales ~ Age + Price*Female + Income, data=data)
summary(modelo2)

## 
## Call:
## lm(formula = Sales ~ Age + Price * Female + Income, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.127 -12.054  -4.804   4.443 126.868 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -6.501e+03  2.043e+03  -3.183  0.00268 **
## Age           5.427e+00  2.638e+00   2.057  0.04563 * 
## Price         1.621e+02  5.196e+01   3.119  0.00320 **
## Female        1.284e+02  3.989e+01   3.219  0.00242 **
## Income        1.253e-02  7.014e-03   1.786  0.08095 . 
## Price:Female -3.240e+00  1.018e+00  -3.183  0.00267 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.67 on 44 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4351, Adjusted R-squared:  0.3709 
## F-statistic: 6.778 on 5 and 44 DF,  p-value: 9.186e-05

2.5 Analice la presencia de multicolinealidad. Si la respuesta es afiramtiva, cómo lo solucionaría?

2.5.1 Análisis de multicolinialidad a través del VIF

Si se evalua el modelo 2, vemos que el indicador vif es muy alto en las variables Price y Female, asi como también la interacción entre ambas. Por lo que, por simplicidad del modelo podemos optar por el modelo sin interacción (modelo 1).

##          Age        Price       Female       Income Price:Female 
##     1.860602  3437.138657   146.989956     1.320652  3597.379621

##      Age   Income    Price 
## 1.116570 1.100580 1.094865

2.6 Analice los residuos e indique si aprecia la presencia de outliers o/y observaciones influyentes.

2.6.1 Análisis gráfico de normalidad

Según el análisis del QQplot a simple vista se observa falta de normalidad en los residuos por presencia de outliers, especificamente el registro 29 y 30 se marcan con valores atípicos.

## [1] 30 29

2.6.2 Análisis de Puntos influyentes

2.6.2.1 Gráfico de distancia Cook

2.6.2.2 Otros puntos influyentes

## Potentially influential observations of
##   lm(formula = Sales ~ Age + Income + Price, data = data) :
## 
##    dfb.1_ dfb.Age dfb.Incm dfb.Pric dffit   cov.r   cook.d hat    
## 2   0.15  -0.30    0.19     0.12     0.36    1.43_*  0.03   0.26_*
## 10 -0.03   0.03   -0.01     0.01     0.04    1.33_*  0.00   0.18  
## 29 -0.32  -0.18    0.48     0.52     0.89_*  0.61_*  0.17   0.09  
## 30  0.19   0.49    0.06    -1.03_*   1.42_*  0.08_*  0.27   0.05  
## 32  0.07  -0.10   -0.05     0.08     0.14    1.27_*  0.01   0.15  
## 38  0.00  -0.01    0.00     0.02    -0.02    1.30_*  0.00   0.16

La observación 9, 29 y 30 poseen un alto leverage y tienen influencia sobre la recta de regresión, por lo cual el modelo al no cumplir con lo supuestos no resulta el mas adecuado. Se podría armar un modelo robusto para restarle peso a estos los outliers.

2.7 Se justifica un método robusto? si la respuesta es afirmativa, constrúyalo.

modelo_robusto <- lmrob(Sales ~ Age + Income + Price, data = data)
resultados_robustos <- summary(modelo_robusto)
resultados_robustos

## 
## Call:
## lmrob(formula = Sales ~ Age + Income + Price, data = data)
##  \--> method = "MM"
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -42.2830  -5.4653  -0.8096  10.3100 142.5234 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 75.256351  30.357912   2.479 0.016904 *  
## Age          2.118118   1.913024   1.107 0.273960    
## Income       0.014023   0.003899   3.597 0.000783 ***
## Price       -1.870741   0.804314  -2.326 0.024487 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 10.86 
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4638, Adjusted R-squared:  0.4288 
## Convergence in 19 IRWLS iterations
## 
## Robustness weights: 
##  4 observations c(9,29,30,34) are outliers with |weight| <= 0.00031 ( < 0.002); 
##  6 weights are ~= 1. The remaining 40 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.09537 0.86780 0.95850 0.85710 0.98800 0.99780 
## Algorithmic parameters: 
##        tuning.chi                bb        tuning.psi        refine.tol 
##         1.548e+00         5.000e-01         4.685e+00         1.000e-07 
##           rel.tol         scale.tol         solve.tol       eps.outlier 
##         1.000e-07         1.000e-10         1.000e-07         2.000e-03 
##             eps.x warn.limit.reject warn.limit.meanrw 
##         9.239e-09         5.000e-01         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##            500             50              2              1            200 
##    maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
##            200              0           1000              0           2000 
##                   psi           subsampling                   cov 
##            "bisquare"         "nonsingular"         ".vcov.avar1" 
## compute.outlier.stats 
##                  "SM" 
## seed : int(0)

Si comparamos el error estandar del modelo robusto respecto al modelo 1, se nota como el modelo robusto disminuye el error a 10.86 vs 27.91 del mejor modelo lineal. Resulta la mejor alternativa en este caso, el modelo robusto.

Sin embargo vemos que el modelo no es adecuado ya que no cumple con los supuesto de normalidad, una alternativa sería intentar con modelos parametricos tipo GAMLSS.

2.7.1 Test de Normalidad

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_robusto$residuals
## W = 0.73152, p-value = 3.188e-08

3 Ejercicio 3.

A partir de los datos de la base insectos.xlsx con las variables disponibles: * grupo (Especie) * longitud_pata * circunf_abdomen * long_antena

Se pretende utilizar estas variables predecir la especie de los insectos.

3.1 Dataset

path <- "C:/Users/ander/Desktop/1. Ciencias de Datos/4. Regresión avanzada/9. Examen 11042023/"
insectos <- read_excel(paste(path,"insectos.xlsx", sep=""))
insectos$grupo <- as.factor(insectos$grupo)

head(insectos)

3.2 Ajustar un modelo logístico para predecir la probabilidad de pertenecer al grupo 1.

Se divide el dataset en test y train para validar las predicciones del modelo, en este caso se opto por incluir a todas las variables del excel y evaluar distintas metricas.

# Lanzamos una semilla para que salgan siempre los mismos datos
set.seed(12345)

# Creamos los dataframes

# Generamos una variable aleatoria con una distribución 70-30
insectos$random<-sample(0:1,size = nrow(insectos),replace = T,prob = c(0.3,0.7)) 

train<-filter(insectos,random==1)
test<-filter(insectos,random==0)

#Eliminamos ya la random para que no moleste
insectos$random <- NULL
train$random <- NULL
test$random <- NULL

rl <- glm(grupo ~ ., data = train, family = "binomial")
summary(rl)

## 
## Call:
## glm(formula = grupo ~ ., family = "binomial", data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0268  -0.1683   0.0000   0.0049   2.1294  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)     -58.3885    57.8054  -1.010    0.312
## longitud_pata     0.4925     0.4408   1.117    0.264
## circunf_abdomen  -1.2211     0.8600  -1.420    0.156
## long_antena       2.3795     1.6893   1.409    0.159
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 30.7891  on 22  degrees of freedom
## Residual deviance:  8.2266  on 19  degrees of freedom
## AIC: 16.227
## 
## Number of Fisher Scoring iterations: 10

pr2_rl <- 1 -(rl$deviance / rl$null.deviance)
pr2_rl

## [1] 0.7328072

3.3 Evaluar la calidad de ajuste del modelo con al menos dos criterios distintos.

Para evaluar el ajuste del modelo se opto por evaluar curva ROC, especificamente el AUC y la cantidad de aciertos del modelo (Accuracy).

Para ello, lo primero que hacemos es buscar el punto de corte optimo que maximice el F1 y el Accuracy, para tener la referencia el umbral para el pronostico del conjunto de test.

## [1] 0.05

Se observa que el umbral 0.05 es el que maximiza el F1 = 2* (precisión * cobertura)/ (precisión + cobertura) y la cantidad de aciertos

3.3.1 Matriz de confusión y curva ROC

confusion(test$grupo,rl_predict,umbral_final_rl2)

##     
## real FALSE TRUE
##    1     6    2
##    2     1    2

El modelo resulta poco adecuado para predecir el grupo de insecto sobre el conjunto de test, ya que solo logro un AUC de 0.58 pese a que el accuracy resulto ser del 72%

3.3.2 Métricas

##               [,1]
## umbral     0.05000
## acierto   72.72727
## precision 50.00000
## cobertura 66.66667
## F1        57.14286
## AUC       58.00000

3.4 Interpretar los coeficientes del modelo elegido.

La variable longitud_pata y 2.3795 cambia el logaritmo de las probabilidades de pertecer al grupo 1 a 0.492 en el caso de longitud _pata mientras que la circunf_abdomen disminuye la probabilidad de pertecer a este grupo de insecto.

exp(cbind(OR = coef(rl), confint(rl)))

##                           OR         2.5 %       97.5 %
## (Intercept)     4.387301e-26 3.980216e-111 3599.5237087
## longitud_pata   1.636415e+00  1.089527e+00    8.1556031
## circunf_abdomen 2.949046e-01  1.331519e-02    0.7622868
## long_antena     1.079967e+01  1.389075e+00 3397.7809766

En cuanto al ratio odds podemos decir que, al aumentar longitud_pata en una unidad, las probabilidades de pertenecer al grupo de insectos 1 aumenta un factor de 1.63 y en el caso de circunf_abdomen este ratio es de 2.94.

3.5 Evaluar la calidad de clasificación y compararlo con otro método de clasificasión

Se prueba con un arbol de clasificasión

library(tidyverse)
library(rpart)
library(rpart.plot)
library(caret)
library(rpart)
arbol <- rpart(grupo ~ ., train,
               xval=      0,
             cp=        0,   #esto significa no limitar la complejidad de los splits
             minsplit=  5,     #minima cantidad de registros para que se haga el split
             minbucket= 5,     #tamaño minimo de una hoja
             maxdepth=  5 )    #profundidad maxima del arbol


prediccion_1 <- predict(arbol, newdata = test, type = "prob")
prediccion_2 <- predict(arbol, newdata = test, type = "class")

library(tibble)
library(ggplot2)
library(rpart.plot)

#  Arbol generado

rpart.plot(arbol, 
                    # show fitted class, probs, percentages
           box.palette = "GnBu", # color scheme
           branch.lty = 3,       # dotted branch lines
           shadow.col = "grey",  # shadows under the node boxes
           nn = TRUE)            # display the node numbers 

3.5.1 Métricas

Comparando con el modelo de regresión logistica vemos como el AUC aumenta si usamos el arbol de clasificasión, sin embargo el Accuracy se mantiene igual.

confusionMatrix(prediccion_2, test[["grupo"]])

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction 1 2
##          1 6 1
##          2 2 2
##                                           
##                Accuracy : 0.7273          
##                  95% CI : (0.3903, 0.9398)
##     No Information Rate : 0.7273          
##     P-Value [Acc > NIR] : 0.6491          
##                                           
##                   Kappa : 0.3774          
##                                           
##  Mcnemar's Test P-Value : 1.0000          
##                                           
##             Sensitivity : 0.7500          
##             Specificity : 0.6667          
##          Pos Pred Value : 0.8571          
##          Neg Pred Value : 0.5000          
##              Prevalence : 0.7273          
##          Detection Rate : 0.5455          
##    Detection Prevalence : 0.6364          
##       Balanced Accuracy : 0.7083          
##                                           
##        'Positive' Class : 1               
## 

# loading the package
library(ROCR)


prediobj <-prediction(prediccion_1[,2],test$grupo)
perf <-  performance(prediobj, "tpr","fpr")

plot(perf,
     main = "Curva ROC",
     xlab="Tasa de falsos positivos", 
     ylab="Tasa de verdaderos positivos")
abline(a=0,b=1,col="blue",lty=2)
grid()
auc <- as.numeric(performance(prediobj,"auc")@y.values)
legend("bottomright",legend=paste(" AUC =",round(auc,4)))