Laboratorio N° 1

Alumnos:

• Marcos Irving Mera Sánchez.

• Robinson Elias Maqui Canaviri.

Caso: Chicago

#Invocación librería faraway
install.packages("faraway")

# Cargar la librería necesaria
library(faraway)

#Las primeras 5 filas del conjunto de datos
head(chicago, n=5)

#Dimensión del datoset (filas y columnas)
dim(chicago)

[1] 47  7

# Ver los nombres de las variables en el conjunto de datos
names(chicago)

[1] "race"     "fire"     "theft"    "age"      "volact"   "involact" "income"  

# O usar str() para ver la estructura del conjunto de datos
str(chicago)

'data.frame':   47 obs. of  7 variables:
 $ race    : num  10 22.2 19.6 17.3 24.5 54 4.9 7.1 5.3 21.5 ...
 $ fire    : num  6.2 9.5 10.5 7.7 8.6 34.1 11 6.9 7.3 15.1 ...
 $ theft   : num  29 44 36 37 53 68 75 18 31 25 ...
 $ age     : num  60.4 76.5 73.5 66.9 81.4 52.6 42.6 78.5 90.1 89.8 ...
 $ volact  : num  5.3 3.1 4.8 5.7 5.9 4 7.9 6.9 7.6 3.1 ...
 $ involact: num  0 0.1 1.2 0.5 0.7 0.3 0 0 0.4 1.1 ...
 $ income  : num  11744 9323 9948 10656 9730 ...

#Obtención informacion de dataset chicago
?faraway::chicago

# Cargar el conjunto de datos
data(chicago)

Preguntas parte 1: (30 puntos) Utilice el conjunto de datos “chicago” disponibles en la librerıa”faraway”. Considere Y = involact como variable respuesta, todas las demas seran variables explicativas.

1. Realizar una análisis descriptivo de las variables de la base de datos. Debe incluir indicadores y gráficas.

##Analisis Exploratorio
#Obtención de estadística descriptiva por cada variable
summary(chicago)

      race            fire           theft             age            volact         involact          income     
 Min.   : 1.00   Min.   : 2.00   Min.   :  3.00   Min.   : 2.00   Min.   : 0.50   Min.   :0.0000   Min.   : 5583  
 1st Qu.: 3.75   1st Qu.: 5.65   1st Qu.: 22.00   1st Qu.:48.60   1st Qu.: 3.10   1st Qu.:0.0000   1st Qu.: 8447  
 Median :24.50   Median :10.40   Median : 29.00   Median :65.00   Median : 5.90   Median :0.4000   Median :10694  
 Mean   :34.99   Mean   :12.28   Mean   : 32.36   Mean   :60.33   Mean   : 6.53   Mean   :0.6149   Mean   :10696  
 3rd Qu.:57.65   3rd Qu.:16.05   3rd Qu.: 38.00   3rd Qu.:77.30   3rd Qu.: 9.65   3rd Qu.:0.9000   3rd Qu.:11989  
 Max.   :99.70   Max.   :39.70   Max.   :147.00   Max.   :90.10   Max.   :14.30   Max.   :2.2000   Max.   :21480  

#Obtención de documentación de dataset chicago
data(chicago)

# Valores faltantes NA
sum(is.na(chicago))

[1] 0

# Ajustar los márgenes inferior, izquierdo, superior y derecho respectivamente
hist(chicago$involact, main = "Variable target involact", xlab = "Involact", ylab = "Valores", col = "lightblue")

#Distribución de los valores de la variable target del dataset
#hist(chicago$involact, main = "Variable target involact",
     #xlab="Involact", ylab = "Valores", col = "lightblue")

# Colores para los boxplots
boxplot_colors <- c( "#FF66CC", "#FFFF33", "#33FFCC", "#3366FF", "#FF3300", "#66FF33","#9933FF")

# Creamos un boxplot para cada variable con colores
par(mfrow = c(2, 4), mar = c(4, 4, 2, 1) + 0.1, cex.axis = 1.8, cex.lab = 1.8)

for (i in 1:7) {
  boxplot(chicago[, i], main = names(chicago)[i], ylab = "", col = boxplot_colors[i], cex.main = 2,  pch = 16, cex = 0.8)
}

## Gráfica de relación de campo fire con variable target involact
library(ggplot2)
ggplot(chicago, aes(x=fire, y=involact)) +  geom_point(size = 2, color = "#8800FF")

## Gráfica de relación de campo fire con variable target involact
library(ggplot2)
ggplot(chicago, aes(x=fire, y=involact)) +  geom_point(size = 2, color = "#0000FF")

## Gráfica de relación de campo race con variable target involact
library(ggplot2)
ggplot(chicago, aes(x=race, y=involact)) +  geom_point(size = 2, color = "#22FF88")

## Gráfica de relación de campo theft con variable target involact
library(ggplot2)
ggplot(chicago, aes(x=theft, y=involact)) +  geom_point(size = 2, color = "#9900FF")

## Gráfica de relación de campo age con variable target involact
library(ggplot2)
ggplot(chicago, aes(x=age, y=involact)) +  geom_point(size = 2, color = "#5500FF")

## Gráfica de relación de campo income con variable target involact
library(ggplot2)
ggplot(chicago, aes(x=age, y=involact)) +  geom_point(size = 2, color = "#0000FF")

# Diagrama de dispersión para dos variables numéricas
plot(chicago$age, chicago$involact)

#Calculo de valores outliers por colunmna
outliers_count <- sapply(chicago, function(x) {
  q <- quantile(x, probs = c(0.25, 0.75), na.rm = TRUE)
  iqr <- q[2] - q[1]
  lower_limit <- q[1] - 1.5 * iqr
  upper_limit <- q[2] + 1.5 * iqr
  sum(x < lower_limit | x > upper_limit, na.rm = TRUE)
})

print(outliers_count)

    race     fire    theft      age   volact involact   income 
       0        3        3        1        0        0        1 

#Diagrama de dispersión de las variables de dataset chicago
pairs(chicago, col = "#FF8800", cex = 1.5, pch = 16)

race: La variable “race” representa el porcentaje de composición racial de minorías en cada código postal.

fire: La variable “fire” representa la tasa de incendios por cada 100 unidades de vivienda.

theft: La variable “theft” representa la tasa de robos por cada 1000 habitantes.

age: La variable “age” representa el porcentaje de unidades de vivienda construidas antes de 1939.

volact: Los valores oscilan entre 0.50 y 14.30, con una mediana de 5.90 y una media de 6.53.

involact: La variable “involact” representa el número de nuevas pólizas y renovaciones del plan FAIR por cada 100 unidades de vivienda.

income: La variable “income” representa el ingreso familiar mediano en miles de dólares.

De acuerdo a lo observado las relaciones lineales con race y fire, con las variabless theft, income y age, no se ve claramente las relaciones.

2. Utilizando alguno de los criterios de selección de variables, determine el modelo lineal que mejor ajusta a la variable respuesta. Indique el criterio utilizado

```r
install.packages("car")
```

#Analisamos la Multicolienalidad entre las variables, menos la variable objeto
faraway::vif(chicago[,-6])

    race     fire    theft      age   volact   income 
3.491088 2.798840 1.684571 2.266203 4.851903 3.153110 

#Aplicando Matriz de correlación entre las variables
cor(chicago[,-6])

             race       fire      theft        age     volact     income
race    1.0000000  0.5927956  0.2550647  0.2505118 -0.7594196 -0.7037328
fire    0.5927956  1.0000000  0.5562105  0.4122225 -0.6864766 -0.6104481
theft   0.2550647  0.5562105  1.0000000  0.3176308 -0.3116183 -0.1729226
age     0.2505118  0.4122225  0.3176308  1.0000000 -0.6057428 -0.5286695
volact -0.7594196 -0.6864766 -0.3116183 -0.6057428  1.0000000  0.7509780
income -0.7037328 -0.6104481 -0.1729226 -0.5286695  0.7509780  1.0000000

Al obtener el valor de VIF de las columnas que resultaron bajos, se consideran aceptables y se sugiere que no hay una alta correlación entre las variables predictoras y no se presenta multicolinealidad que afecte la estabilidad del modelo dado.

#Construcción de Modelo Lineal
reg1 = lm(involact ~ ., data = chicago)

summary(reg1)

Call:
lm(formula = involact ~ ., data = chicago)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.84296 -0.14613 -0.01007  0.18386  0.81235 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4.862e-01  6.020e-01  -0.808 0.424109    
race         8.527e-03  2.863e-03   2.978 0.004911 ** 
fire         3.778e-02  8.982e-03   4.206 0.000142 ***
theft       -1.016e-02  2.908e-03  -3.494 0.001178 ** 
age          7.615e-03  3.330e-03   2.287 0.027582 *  
volact      -1.018e-02  2.773e-02  -0.367 0.715519    
income       2.568e-05  3.220e-05   0.798 0.429759    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3387 on 40 degrees of freedom
Multiple R-squared:  0.7517,    Adjusted R-squared:  0.7144 
F-statistic: 20.18 on 6 and 40 DF,  p-value: 1.072e-10

#Backward por AIC

El valor p es de 0.424109, nos dice que no se presenta suficiente evidencia para rechazar la hipótesis nula dado que el coeficiente del intercepto es igual a cero.

Asimismo, se indican que los coeficientes para la variable “fire” son significativos al nivel de 0.001, mientras que el coeficiente para la variable “theft” y “race” es significativo al nivel de 0.01 y “age” a 0.05 Los coeficientes para las variables “volact” e “income” no son significativos a estos niveles.

El valor de Multiple R-squared es de 0.7517 que explica el 75% las variables analizadas.

El estadístico F de 20.18 y un valor de p de 1.072e-10 indica que el modelo de regresión es estadísticamente significativo, hay suficiente evidencia para rechazar la hipótesis nula.

#Aplicación de Backward por AIC
mod1.step0=step(reg1)

Start:  AIC=-95.34
involact ~ race + fire + theft + age + volact + income

         Df Sum of Sq    RSS     AIC
- volact  1   0.01546 4.6047 -97.184
- income  1   0.07300 4.6622 -96.601
<none>                4.5892 -95.342
- age     1   0.59993 5.1892 -91.568
- race    1   1.01743 5.6067 -87.931
- theft   1   1.40048 5.9897 -84.825
- fire    1   2.02990 6.6191 -80.129

Step:  AIC=-97.18
involact ~ race + fire + theft + age + income

         Df Sum of Sq    RSS     AIC
- income  1   0.06710 4.6718 -98.504
<none>                4.6047 -97.184
- age     1   0.99296 5.5977 -90.007
- theft   1   1.46328 6.0680 -86.215
- race    1   1.74657 6.3513 -84.070
- fire    1   2.37807 6.9828 -79.615

Step:  AIC=-98.5
involact ~ race + fire + theft + age

        Df Sum of Sq    RSS     AIC
<none>               4.6718 -98.504
- age    1   0.99734 5.6691 -91.410
- theft  1   1.41436 6.0862 -88.074
- race   1   2.05375 6.7256 -83.379
- fire   1   2.38365 7.0554 -81.128

#Aplicación de Backward por BIC
mod1.step1=step(reg1, k=log(length(chicago$involact)))

Start:  AIC=-82.39
involact ~ race + fire + theft + age + volact + income

         Df Sum of Sq    RSS     AIC
- volact  1   0.01546 4.6047 -86.083
- income  1   0.07300 4.6622 -85.500
<none>                4.5892 -82.391
- age     1   0.59993 5.1892 -80.467
- race    1   1.01743 5.6067 -76.830
- theft   1   1.40048 5.9897 -73.724
- fire    1   2.02990 6.6191 -69.028

Step:  AIC=-86.08
involact ~ race + fire + theft + age + income

         Df Sum of Sq    RSS     AIC
- income  1   0.06710 4.6718 -89.254
<none>                4.6047 -86.083
- age     1   0.99296 5.5977 -80.756
- theft   1   1.46328 6.0680 -76.964
- race    1   1.74657 6.3513 -74.819
- fire    1   2.37807 6.9828 -70.364

Step:  AIC=-89.25
involact ~ race + fire + theft + age

        Df Sum of Sq    RSS     AIC
<none>               4.6718 -89.254
- age    1   0.99734 5.6691 -84.010
- theft  1   1.41436 6.0862 -80.674
- race   1   2.05375 6.7256 -75.978
- fire   1   2.38365 7.0554 -73.728

#Eliminacion y adicion de variables
reg2 = lm(involact ~race+fire+theft, data = chicago)
mod1.step2=step(reg2, scope=list(lower=reg2, upper=reg1), direction="both")

Start:  AIC=-91.41
involact ~ race + fire + theft

         Df Sum of Sq    RSS     AIC
+ age     1   0.99734 4.6718 -98.504
+ volact  1   0.47987 5.1893 -93.567
<none>                5.6691 -91.410
+ income  1   0.07148 5.5977 -90.007

Step:  AIC=-98.5
involact ~ race + fire + theft + age

         Df Sum of Sq    RSS     AIC
<none>                4.6718 -98.504
+ income  1   0.06710 4.6047 -97.184
+ volact  1   0.00955 4.6622 -96.601
- age     1   0.99734 5.6691 -91.410

#Evaluación de las metricas de los 03 nuevos modelos obtenidos
summary(mod1.step0)

Call:
lm(formula = involact ~ race + fire + theft + age, data = chicago)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.87108 -0.14830 -0.01961  0.19968  0.81638 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.243118   0.145054  -1.676 0.101158    
race         0.008104   0.001886   4.297 0.000100 ***
fire         0.036646   0.007916   4.629 3.51e-05 ***
theft       -0.009592   0.002690  -3.566 0.000921 ***
age          0.007210   0.002408   2.994 0.004595 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3335 on 42 degrees of freedom
Multiple R-squared:  0.7472,    Adjusted R-squared:  0.7231 
F-statistic: 31.03 on 4 and 42 DF,  p-value: 4.799e-12

summary(mod1.step1)

Call:
lm(formula = involact ~ race + fire + theft + age, data = chicago)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.87108 -0.14830 -0.01961  0.19968  0.81638 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.243118   0.145054  -1.676 0.101158    
race         0.008104   0.001886   4.297 0.000100 ***
fire         0.036646   0.007916   4.629 3.51e-05 ***
theft       -0.009592   0.002690  -3.566 0.000921 ***
age          0.007210   0.002408   2.994 0.004595 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3335 on 42 degrees of freedom
Multiple R-squared:  0.7472,    Adjusted R-squared:  0.7231 
F-statistic: 31.03 on 4 and 42 DF,  p-value: 4.799e-12

summary(mod1.step2)

Call:
lm(formula = involact ~ race + fire + theft + age, data = chicago)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.87108 -0.14830 -0.01961  0.19968  0.81638 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.243118   0.145054  -1.676 0.101158    
race         0.008104   0.001886   4.297 0.000100 ***
fire         0.036646   0.007916   4.629 3.51e-05 ***
theft       -0.009592   0.002690  -3.566 0.000921 ***
age          0.007210   0.002408   2.994 0.004595 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3335 on 42 degrees of freedom
Multiple R-squared:  0.7472,    Adjusted R-squared:  0.7231 
F-statistic: 31.03 on 4 and 42 DF,  p-value: 4.799e-12

Luego de haber obtenido la información de los 3 modelos, seleccionamos el modelo mod1.step2, se usará para comparar con el modelo reg1, y ver cual se ajusta mejor a los datos

#Comparación de modelos REG1 y mod1.step2
library(knitr)
df <- data.frame(Modelo = c("reg1",  "mod1.step2"),
                 AIC = c(AIC(reg1), AIC(mod1.step2)),
                 BIC = c(BIC(reg1), BIC(mod1.step2)))

tabla <- knitr::kable(df, caption = "Comparison of AIC and BIC", align = c("l", "r", "r"))
print(tabla)

Comparison of AIC and BIC

Modelo	AIC	BIC
reg1	40.03788	54.83906
mod1.step2	36.87586	47.97675

anova(reg1, mod1.step2)

Analysis of Variance Table

Model 1: involact ~ race + fire + theft + age + volact + income
Model 2: involact ~ race + fire + theft + age
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     40 4.5892                           
2     42 4.6718 -2 -0.082558 0.3598 0.7001

Luego de haber ejecutado la función ANOVA, el valor p obtenido es 0.7001, el cual es mayor que el umbral de significancia comúnmente utilizado 0.05). Esto nos indica que no hay una mejora significativa al eliminar variables predictoras.

install.packages("car")

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/robin/AppData/Local/R/win-library/4.3’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/car_3.1-2.zip'
Content type 'application/zip' length 1707562 bytes (1.6 MB)
downloaded 1.6 MB

package ‘car’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\robin\AppData\Local\Temp\RtmpCmWd9p\downloaded_packages

library(car)

Loading required package: carData

Attaching package: ‘car’

The following objects are masked from ‘package:faraway’:

    logit, vif

car::Anova(reg1,type=2)

Anova Table (Type II tests)

Response: involact
          Sum Sq Df F value    Pr(>F)    
race      1.0174  1  8.8679 0.0049112 ** 
fire      2.0299  1 17.6927 0.0001421 ***
theft     1.4005  1 12.2066 0.0011784 ** 
age       0.5999  1  5.2290 0.0275823 *  
volact    0.0155  1  0.1347 0.7155185    
income    0.0730  1  0.6363 0.4297589    
Residuals 4.5892 40                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Respuesta: Según lo observado en las pruebas de comparación el modelo mod1.step2 no se presenta una mejora significativa según la prueba ANOVA, pero si se observa que los estimadores AIC y BIC este modelo tiene menores valores entonces es mejor, considerar que se han omitido las covariables de volact e income que no son significativas para la variable predictoria. Si comparamos el valor del R-cuadrado (entre reg1 y mod1.step2) se tiene una disminución del 0.0045 en el modelo mod1.step2, lo cual no es muy pequeña, pero en el modelo mod1.step2 se tiene valores mas altos de R-cuadrado-ajustado y F-stadistico, y un p_valor mas pequeño. Por lo comentado, se considera al mejor modelo es mod1.step2, con el cual se realizará las predicciones.

#Diagnóstico
shapiro.test(reg1$residuals)

    Shapiro-Wilk normality test

data:  reg1$residuals
W = 0.98095, p-value = 0.6317

shapiro.test(mod1.step2$residuals)

    Shapiro-Wilk normality test

data:  mod1.step2$residuals
W = 0.98065, p-value = 0.6191

qqnorm(reg1$residuals, pch=20); qqline(reg1$residuals, col=2, lwd=2)

qqnorm(mod1.step2$residuals, pch=20); qqline(mod1.step2$residuals, col=2, lwd=2)

qqnorm(chicago$involact)

3. Analice la significancia del modelo obtenido luego del proceso de selección, y responda si:

• ¿Es el modelo obtenido significativo?

• ¿Existe alguna covariable no significativa?

• ¿En caso de existir alguna covariable no significativa, la quitar ́ıa del modelo?. Fundamente.

  RESPUESTAS:

  • ¿Es el modelo obtenido significativo?

    Respuesta:El modelo mod1.step2 es el que resulta como significativo. El estadístico F de 31.03 y un valor de p de 4.799e-12 indica que el modelo de regresión es estadísticamente significativo, por lo que hay suficiente evidencia para rechazar la hipótesis nula de que no hay relación entre las variables.

  • ¿Existe alguna covariable no significativa?

    Respuesta:De acuerdo al análisis en cada modelo y con la función ANOVA, las covariables volact (p_valor 0.715519) e income (p_valor 0.429759) no son significativas

  • ¿En caso de existir alguna covariable no significativa, la quitaría del modelo?. Fundamente.

    Respuesta:Si, se retirará del modelo, porque:

    • Facilita la interpretación del modelo y evita el riesgo de sobreajuste.

    • Tener un modelo más simple con la menor cantidad de covariables necesarias para explicar el fenómeno de interés.

    • Reducir los grados de libertad y se mejorará la estimación de los coeficientes de las covariables restantes.

4. Fijar al menos 5 valores para las covariables del modelo y con ellas realizar la predicción de la
media y la predicci ́on individual de la variable objetivo (incluir los intervalos de confianza).

#Seleccionar 5 valores del dataset
sample_data <- chicago[sample.int(47,5),]

#Visualizar valores seleccionados
sample_data

# Fijar los valores de las covariables
new_data <- data.frame(race = sample_data$race,
                       fire = sample_data$fire,
                       theft = sample_data$theft,
                       age = sample_data$age)
new_data

#Aplicando predicción de la media
# Predicción de la media
mean_prediction <- predict(mod1.step2, newdata = new_data, interval = "confidence")
mod1.step2

Call:
lm(formula = involact ~ race + fire + theft + age, data = chicago)

Coefficients:
(Intercept)         race         fire        theft          age  
  -0.243118     0.008104     0.036646    -0.009592     0.007210  

print(mean_prediction)

          fit         lwr        upr
1  0.46060120  0.27746564 0.64373676
2  0.75884612  0.57239299 0.94529926
3 -0.21258903 -0.50165730 0.07647923
4  0.07940122 -0.08662042 0.24542286
5  1.44268328  1.17969764 1.70566892

#Aplicando predicción individual
# Predicción individual
indiv_prediction <- predict(mod1.step2, newdata = new_data, interval = "prediction")

# Imprimir los resultados
print(indiv_prediction)

          fit         lwr       upr
1  0.46060120 -0.23693259 1.1581350
2  0.75884612  0.06043398 1.4572583
3 -0.21258903 -0.94510198 0.5199239
4  0.07940122 -0.61383606 0.7726385
5  1.44268328  0.72006568 2.1653009

Con ambas predicciones obtenemos resultados iguales, validamos con el tipo “response”:

reg1.pred1 <- predict(mod1.step2, new_data, type="response")
reg1.pred1

          1           2           3           4           5 
 0.46060120  0.75884612 -0.21258903  0.07940122  1.44268328 

library(ggplot2)

mydata <- cbind(sample_data, indiv_prediction )
p <- ggplot(mydata, aes(fit, involact)) +   geom_point(size = 3) +   stat_smooth(method = lm)
p + geom_line(aes(y = lwr), color = "red", linetype = "dashed") +   geom_line(aes(y = upr), color = "red", linetype = "dashed")

Resultados de la observaciones de error de datos aleatorios

err_obs=sample_data$involact-mydata$fit
var_err_obs=sqrt(sum(err_obs^2)/5)
score=1-sum(abs(err_obs))/5

print("Datos a predecir:")

[1] "Datos a predecir:"

print(sample_data$involact) #datos a predecir

[1] 0.0 0.4 0.0 0.0 1.9

print("Prediccion")

[1] "Prediccion"

print(mydata$fit) #prediccion

[1]  0.46060120  0.75884612 -0.21258903  0.07940122  1.44268328

print("Error observable")

[1] "Error observable"

print(err_obs) #error observable

[1] -0.46060120 -0.35884612  0.21258903 -0.07940122  0.45731672

print("Desviación estándar del error")

[1] "Desviación estándar del error"

print(var_err_obs) #desviacion estandar del error

[1] 0.3468606

print(score)

[1] 0.6862491

Podemos concluir que el modelo parece ajustarse razonablemente bien a los datos observados, ya que los errores son pequeños en comparación con los valores observados.

Caso: Wbca

# Cargar la librería necesaria
library(faraway)

#Las primeras 5 filas del conjunto de datos
head(wbca, n=5)

#Dimensión del datoset (filas y columnas)
dim(wbca)

[1] 681  10

# Ver los nombres de las variables en el conjunto de datos
names(wbca)

 [1] "Class" "Adhes" "BNucl" "Chrom" "Epith" "Mitos" "NNucl" "Thick" "UShap" "USize"

# O usar str() para ver la estructura del conjunto de datos
str(wbca)

'data.frame':   681 obs. of  10 variables:
 $ Class: int  1 1 1 1 1 0 1 1 1 1 ...
 $ Adhes: int  1 5 1 1 3 8 1 1 1 1 ...
 $ BNucl: int  1 10 2 4 1 10 10 1 1 1 ...
 $ Chrom: int  3 3 3 3 3 9 3 3 1 2 ...
 $ Epith: int  2 7 2 3 2 7 2 2 2 2 ...
 $ Mitos: int  1 1 1 1 1 1 1 1 5 1 ...
 $ NNucl: int  1 2 1 7 1 7 1 1 1 1 ...
 $ Thick: int  5 5 3 6 4 8 1 2 2 4 ...
 $ UShap: int  1 4 1 8 1 10 1 2 1 1 ...
 $ USize: int  1 4 1 8 1 10 1 1 1 2 ...

#Obtención informacion de dataset chicago
?faraway::wbca

# Cargar el conjunto de datos
data(wbca)

Preguntas parte 2: (30 puntos) Utilice el conjunto de datos “chicago” disponibles en la librerıa”faraway”. Considere Y = involact como variable respuesta, todas las demas ser an variables explica-tivas.

1. Realizar una análisis descriptivo de las variables de la base de datos. Debe incluir indicadores y gráficas.

# Obtener un resumen estadístico de las variables
summary(wbca)

     Class            Adhes            BNucl            Chrom            Epith            Mitos            NNucl            Thick            UShap       
 Min.   :0.0000   Min.   : 1.000   Min.   : 1.000   Min.   : 1.000   Min.   : 1.000   Min.   : 1.000   Min.   : 1.000   Min.   : 1.000   Min.   : 1.000  
 1st Qu.:0.0000   1st Qu.: 1.000   1st Qu.: 1.000   1st Qu.: 2.000   1st Qu.: 2.000   1st Qu.: 1.000   1st Qu.: 1.000   1st Qu.: 2.000   1st Qu.: 1.000  
 Median :1.0000   Median : 1.000   Median : 1.000   Median : 3.000   Median : 2.000   Median : 1.000   Median : 1.000   Median : 4.000   Median : 1.000  
 Mean   :0.6505   Mean   : 2.816   Mean   : 3.542   Mean   : 3.433   Mean   : 3.231   Mean   : 1.604   Mean   : 2.859   Mean   : 4.436   Mean   : 3.204  
 3rd Qu.:1.0000   3rd Qu.: 4.000   3rd Qu.: 6.000   3rd Qu.: 5.000   3rd Qu.: 4.000   3rd Qu.: 1.000   3rd Qu.: 4.000   3rd Qu.: 6.000   3rd Qu.: 5.000  
 Max.   :1.0000   Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :10.000  
     USize      
 Min.   : 1.00  
 1st Qu.: 1.00  
 Median : 1.00  
 Mean   : 3.14  
 3rd Qu.: 5.00  
 Max.   :10.00  

# Crear gráficas
# Histograma para la variable target del dataset
hist(wbca$Class)

# Diagrama de dispersión para todas las variables
plot(wbca)

# Diagrama de dispersión para dos variables numéricas
plot(wbca$Thick, wbca$Class)

#Generamos gráfico de dispersión para la variable involact y fire.
ggplot(wbca, aes(x=Thick, y=Class)) + geom_point()

# crea un histograma que visualiza la distribución de los valores de la variable "Thick" del conjunto de datos "wbca"
ggplot(wbca, aes(x=Thick)) + geom_histogram(binwidth=1)

Adhes, BNucl, Chrom, Epith, Mitos, NNucl, Thick, UShap, USize: Estas variables representan atributos y toman valores de 1 a 10.

Class: Esta variable representa una clasificación y toma valores 0 o 1.

#Validamos valores faltantes
sum(is.na(wbca))

[1] 0

#Revisamos valores outliers
# Colores para los boxplots
boxplot_colors <- c("#1366FF", "#AF3300", "#D6FF33", "#AE66CC", "#FFFF33", "#B3FFCC", "#6933FF", "#FF9933", "#CC33FF", "#A3FFFF")

# Se crea un color para cada variable
par(mfrow = c(2, 5), mar = c(4, 4, 2, 1) + 0.1, cex.axis = 1.8, cex.lab = 1.8)
for (i in 1:10) {
  boxplot(wbca[, i], main = names(wbca)[i], ylab = "", col = boxplot_colors[i], cex.main = 2, pch = 16, cex = 0.8)
}

# Contamos los outliers por columna
count_outliers <- sapply(wbca, function(x) {
  q <- quantile(x, probs = c(0.25, 0.75), na.rm = TRUE)
  iqr <- q[2] - q[1]
  lower_limit <- q[1] - 1.5 * iqr
  upper_limit <- q[2] + 1.5 * iqr
  sum(x < lower_limit | x > upper_limit, na.rm = TRUE)
})
print(count_outliers)

Class Adhes BNucl Chrom Epith Mitos NNucl Thick UShap USize 
    0    58     0    20    54   119    75     0     0     0 

Se tiene pocos valores out-liers a diferencia de la covariable Mitos, el cual tiene una mayor cantidad

# Ajustamos el tamaño de la imagen
par(mfrow = c(7, 7), mar = c(2, 2, 2, 2))
#  genera un gráfico de matriz que muestra la relación entre pares de variables en el dataframe "wbca"
pairs(wbca, col = "#A366FF", cex = 1.5, pch = 16)

#Cargamos libreria
library(gridExtra)

#crea un gráfico de dispersión utilizando ggplot2, donde los valores de la variable Class se mapean al eje x y los valores de la variable Adhes se mapean al eje y. 
p1 <- ggplot(wbca, aes(Class, Adhes)) +
  geom_point(size = 2, color = "#FF0000") +
  theme(plot.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12))

#crea un gráfico de dispersión utilizando ggplot2, donde los valores de la variable Class se mapean al eje x y los valores de la variable BNucl se mapean al eje y. 
p2 <- ggplot(wbca, aes(Class, BNucl)) +
  geom_point(size = 2, color = "#FF0000") +
  theme(plot.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12))

#crea un gráfico de dispersión utilizando ggplot2, donde los valores de la variable Class se mapean al eje x y los valores de la variable Chrom se mapean al eje y. 
p3 <- ggplot(wbca, aes(Class, Chrom)) +
  geom_point(size = 2, color = "#FF0000") +
  theme(plot.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12))

#crea un gráfico de dispersión utilizando ggplot2, donde los valores de la variable Class se mapean al eje x y los valores de la variable Epith se mapean al eje y. 
p4 <- ggplot(wbca, aes(Class, Epith)) +
  geom_point(size = 2, color = "#FF0000") +
  theme(plot.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12))

#crea un gráfico de dispersión utilizando ggplot2, donde los valores de la variable Class se mapean al eje x y los valores de la variable Adhes se mapean al eje y. 
p5 <- ggplot(wbca, aes(Class, Mitos)) +
  geom_point(size = 2, color = "#FF0000") +
  theme(plot.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12))

#crea un gráfico de dispersión utilizando ggplot2, donde los valores de la variable Class se mapean al eje x y los valores de la variable NNucl se mapean al eje y. 
p6 <- ggplot(wbca, aes(Class, NNucl)) +
  geom_point(size = 2, color = "#FF0000") +
  theme(plot.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12))

#crea un gráfico de dispersión utilizando ggplot2, donde los valores de la variable Class se mapean al eje x y los valores de la variable Thick se mapean al eje y. 
p7 <- ggplot(wbca, aes(Class, Thick)) +
  geom_point(size = 2, color = "#FF0000") +
  theme(plot.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12))

#crea un gráfico de dispersión utilizando ggplot2, donde los valores de la variable Class se mapean al eje x y los valores de la variable UShap se mapean al eje y. 
p8 <- ggplot(wbca, aes(Class, UShap)) +
  geom_point(size = 2, color = "#FF0000") +
  theme(plot.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12))

#crea un gráfico de dispersión utilizando ggplot2, donde los valores de la variable Class se mapean al eje x y los valores de la variable USize se mapean al eje y. 
p9 <- ggplot(wbca, aes(Class, USize)) +
  geom_point(size = 2, color = "#FF0000") +
  theme(plot.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12))

grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, p9, ncol = 3)

2. Utilizando alguno de los criterios de selección de variables, determine el modelo de regresión logística que mejor ajusta a la variable respuesta. Indique el criterio utilizado.

#Revisamos la correlación
cor_matrix <- cor(wbca[, -1])
cor_matrix

          Adhes     BNucl     Chrom     Epith     Mitos     NNucl     Thick     UShap     USize
Adhes 1.0000000 0.6750797 0.6661809 0.5945634 0.4221024 0.6031777 0.4887365 0.6833626 0.7046159
BNucl 0.6750797 1.0000000 0.6830874 0.5865168 0.3389872 0.5842596 0.5926375 0.7182476 0.6958544
Chrom 0.6661809 0.6830874 1.0000000 0.6180940 0.3479439 0.6644418 0.5534378 0.7343984 0.7547574
Epith 0.5945634 0.5865168 0.6180940 1.0000000 0.4811998 0.6291379 0.5236406 0.7228546 0.7540225
Mitos 0.4221024 0.3389872 0.3479439 0.4811998 1.0000000 0.4346642 0.3511003 0.4440921 0.4635913
NNucl 0.6031777 0.5842596 0.6644418 0.6291379 0.4346642 1.0000000 0.5326826 0.7194922 0.7208030
Thick 0.4887365 0.5926375 0.5534378 0.5236406 0.3511003 0.5326826 1.0000000 0.6556203 0.6444869
UShap 0.6833626 0.7182476 0.7343984 0.7228546 0.4440921 0.7194922 0.6556203 1.0000000 0.9065227
USize 0.7046159 0.6958544 0.7547574 0.7540225 0.4635913 0.7208030 0.6444869 0.9065227 1.0000000

# Calcular el VIF para cada variable predictora
library(car)
vif_values <- vif(glm(Class ~ ., family = binomial, data = wbca))

# Mostrar los valores del VIF
vif_values

   Adhes    BNucl    Chrom    Epith    Mitos    NNucl    Thick    UShap    USize 
1.215767 1.157928 1.180893 1.351148 1.050105 1.229996 1.228965 3.088274 3.110742 

#Analisis de multicolinealidad
# Ajustar un modelo de regresión
model <- lm(Class ~ ., data = wbca)

# Calcular los factores de inflación de la varianza (VIF)
vif_values <- car::vif(model)
cor(wbca[, -1])

          Adhes     BNucl     Chrom     Epith     Mitos     NNucl     Thick     UShap     USize
Adhes 1.0000000 0.6750797 0.6661809 0.5945634 0.4221024 0.6031777 0.4887365 0.6833626 0.7046159
BNucl 0.6750797 1.0000000 0.6830874 0.5865168 0.3389872 0.5842596 0.5926375 0.7182476 0.6958544
Chrom 0.6661809 0.6830874 1.0000000 0.6180940 0.3479439 0.6644418 0.5534378 0.7343984 0.7547574
Epith 0.5945634 0.5865168 0.6180940 1.0000000 0.4811998 0.6291379 0.5236406 0.7228546 0.7540225
Mitos 0.4221024 0.3389872 0.3479439 0.4811998 1.0000000 0.4346642 0.3511003 0.4440921 0.4635913
NNucl 0.6031777 0.5842596 0.6644418 0.6291379 0.4346642 1.0000000 0.5326826 0.7194922 0.7208030
Thick 0.4887365 0.5926375 0.5534378 0.5236406 0.3511003 0.5326826 1.0000000 0.6556203 0.6444869
UShap 0.6833626 0.7182476 0.7343984 0.7228546 0.4440921 0.7194922 0.6556203 1.0000000 0.9065227
USize 0.7046159 0.6958544 0.7547574 0.7540225 0.4635913 0.7208030 0.6444869 0.9065227 1.0000000

cor(wbca[,-10])

           Class      Adhes      BNucl      Chrom      Epith      Mitos      NNucl      Thick      UShap
Class  1.0000000 -0.7069017 -0.8266660 -0.7605170 -0.6913350 -0.4251902 -0.7221115 -0.7178680 -0.8222345
Adhes -0.7069017  1.0000000  0.6750797  0.6661809  0.5945634  0.4221024  0.6031777  0.4887365  0.6833626
BNucl -0.8266660  0.6750797  1.0000000  0.6830874  0.5865168  0.3389872  0.5842596  0.5926375  0.7182476
Chrom -0.7605170  0.6661809  0.6830874  1.0000000  0.6180940  0.3479439  0.6644418  0.5534378  0.7343984
Epith -0.6913350  0.5945634  0.5865168  0.6180940  1.0000000  0.4811998  0.6291379  0.5236406  0.7228546
Mitos -0.4251902  0.4221024  0.3389872  0.3479439  0.4811998  1.0000000  0.4346642  0.3511003  0.4440921
NNucl -0.7221115  0.6031777  0.5842596  0.6644418  0.6291379  0.4346642  1.0000000  0.5326826  0.7194922
Thick -0.7178680  0.4887365  0.5926375  0.5534378  0.5236406  0.3511003  0.5326826  1.0000000  0.6556203
UShap -0.8222345  0.6833626  0.7182476  0.7343984  0.7228546  0.4440921  0.7194922  0.6556203  1.0000000

De acuerdo a la correlación, la variable USize tiene una alta colinealidad con la variable Ushap, por lo tanto se retira de la prueba y se vuelve a realizar.

El valor F tiene valores bajos y se considera aceptable, esto sugiere a que no hay una alta correlación entre las variables predictoras. Como la multicolinealidad puede afectar la estabilidad del modelo, la interpretación de resultados y comprensión de como contribuye cada variable al modelo.

#Hacemos la regresion logistica con todas las variables
log1=glm(Class ~., data=wbca, family=binomial )
summary(log1)

Call:
glm(formula = Class ~ ., family = binomial, data = wbca)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 11.16678    1.41491   7.892 2.97e-15 ***
Adhes       -0.39681    0.13384  -2.965  0.00303 ** 
BNucl       -0.41478    0.10230  -4.055 5.02e-05 ***
Chrom       -0.56456    0.18728  -3.014  0.00257 ** 
Epith       -0.06440    0.16595  -0.388  0.69795    
Mitos       -0.65713    0.36764  -1.787  0.07387 .  
NNucl       -0.28659    0.12620  -2.271  0.02315 *  
Thick       -0.62675    0.15890  -3.944 8.01e-05 ***
UShap       -0.28011    0.25235  -1.110  0.26699    
USize        0.05718    0.23271   0.246  0.80589    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 881.388  on 680  degrees of freedom
Residual deviance:  89.464  on 671  degrees of freedom
AIC: 109.46

Number of Fisher Scoring iterations: 8

#Hacemos la regresion logistica sin la variable USize (que tiene alta colinealidad)
log2=glm(Class ~ Adhes + BNucl + Chrom + Epith + Mitos + NNucl + Thick + UShap, data= wbca, family=binomial)
summary(log2)

Call:
glm(formula = Class ~ Adhes + BNucl + Chrom + Epith + Mitos + 
    NNucl + Thick + UShap, family = binomial, data = wbca)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 11.11149    1.39160   7.985 1.41e-15 ***
Adhes       -0.39034    0.13101  -2.980  0.00289 ** 
BNucl       -0.41493    0.10238  -4.053 5.06e-05 ***
Chrom       -0.55806    0.18618  -2.998  0.00272 ** 
Epith       -0.06192    0.16544  -0.374  0.70818    
Mitos       -0.65047    0.36496  -1.782  0.07470 .  
NNucl       -0.28216    0.12484  -2.260  0.02381 *  
Thick       -0.62298    0.15847  -3.931 8.45e-05 ***
UShap       -0.23654    0.18264  -1.295  0.19529    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 881.388  on 680  degrees of freedom
Residual deviance:  89.523  on 672  degrees of freedom
AIC: 107.52

Number of Fisher Scoring iterations: 8

#Hacemos el Backward por AIC
log3=step(log1)

Start:  AIC=109.46
Class ~ Adhes + BNucl + Chrom + Epith + Mitos + NNucl + Thick + 
    UShap + USize

        Df Deviance    AIC
- USize  1   89.523 107.52
- Epith  1   89.613 107.61
- UShap  1   90.627 108.63
<none>       89.464 109.46
- Mitos  1   93.551 111.55
- NNucl  1   95.204 113.20
- Adhes  1   98.844 116.84
- Chrom  1   99.841 117.84
- BNucl  1  109.000 127.00
- Thick  1  110.239 128.24

Step:  AIC=107.52
Class ~ Adhes + BNucl + Chrom + Epith + Mitos + NNucl + Thick + 
    UShap

        Df Deviance    AIC
- Epith  1   89.662 105.66
- UShap  1   91.355 107.36
<none>       89.523 107.52
- Mitos  1   93.552 109.55
- NNucl  1   95.231 111.23
- Adhes  1   99.042 115.04
- Chrom  1  100.153 116.15
- BNucl  1  109.064 125.06
- Thick  1  110.465 126.47

Step:  AIC=105.66
Class ~ Adhes + BNucl + Chrom + Mitos + NNucl + Thick + UShap

        Df Deviance    AIC
<none>       89.662 105.66
- UShap  1   91.884 105.88
- Mitos  1   93.714 107.71
- NNucl  1   95.853 109.85
- Adhes  1  100.126 114.13
- Chrom  1  100.844 114.84
- BNucl  1  109.762 123.76
- Thick  1  110.632 124.63

summary(log3)

Call:
glm(formula = Class ~ Adhes + BNucl + Chrom + Mitos + NNucl + 
    Thick + UShap, family = binomial, data = wbca)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  11.0333     1.3632   8.094 5.79e-16 ***
Adhes        -0.3984     0.1294  -3.080  0.00207 ** 
BNucl        -0.4192     0.1020  -4.111 3.93e-05 ***
Chrom        -0.5679     0.1840  -3.085  0.00203 ** 
Mitos        -0.6456     0.3634  -1.777  0.07561 .  
NNucl        -0.2915     0.1236  -2.358  0.01837 *  
Thick        -0.6216     0.1579  -3.937 8.27e-05 ***
UShap        -0.2541     0.1785  -1.423  0.15461    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 881.388  on 680  degrees of freedom
Residual deviance:  89.662  on 673  degrees of freedom
AIC: 105.66

Number of Fisher Scoring iterations: 8

#Hacemos la Comparación de modelos
anova(log2, log1, test="Chisq")

Analysis of Deviance Table

Model 1: Class ~ Adhes + BNucl + Chrom + Epith + Mitos + NNucl + Thick + 
    UShap
Model 2: Class ~ Adhes + BNucl + Chrom + Epith + Mitos + NNucl + Thick + 
    UShap + USize
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       672     89.523                     
2       671     89.464  1 0.059177   0.8078

anova(log3, log1, test="Chisq")

Analysis of Deviance Table

Model 1: Class ~ Adhes + BNucl + Chrom + Mitos + NNucl + Thick + UShap
Model 2: Class ~ Adhes + BNucl + Chrom + Epith + Mitos + NNucl + Thick + 
    UShap + USize
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       673     89.662                     
2       671     89.464  2  0.19756   0.9059

anova(log2, log3, test="Chisq")

Analysis of Deviance Table

Model 1: Class ~ Adhes + BNucl + Chrom + Epith + Mitos + NNucl + Thick + 
    UShap
Model 2: Class ~ Adhes + BNucl + Chrom + Mitos + NNucl + Thick + UShap
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       672     89.523                     
2       673     89.662 -1 -0.13839   0.7099

anova(log3, log2, log1,  test="Chisq")

Analysis of Deviance Table

Model 1: Class ~ Adhes + BNucl + Chrom + Mitos + NNucl + Thick + UShap
Model 2: Class ~ Adhes + BNucl + Chrom + Epith + Mitos + NNucl + Thick + 
    UShap
Model 3: Class ~ Adhes + BNucl + Chrom + Epith + Mitos + NNucl + Thick + 
    UShap + USize
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       673     89.662                     
2       672     89.523  1 0.138389   0.7099
3       671     89.464  1 0.059177   0.8078

Esta información nos muestra que no hay diferencias significativas en la devianza residual(medida que se utiliza para evaluar el ajuste del modelo) entre los modelos log1, log2 y log3. Esto sugiere que estos modelos no son muy diferentes de manera significativa en su ajuste a los datos.

library(car)
car::Anova(log1,type=2)

Analysis of Deviance Table (Type II tests)

Response: Class
      LR Chisq Df Pr(>Chisq)    
Adhes   9.3798  1   0.002194 ** 
BNucl  19.5363  1  9.871e-06 ***
Chrom  10.3767  1   0.001276 ** 
Epith   0.1487  1   0.699763    
Mitos   4.0868  1   0.043220 *  
NNucl   5.7401  1   0.016582 *  
Thick  20.7744  1  5.167e-06 ***
UShap   1.1628  1   0.280887    
USize   0.0592  1   0.807802    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library(knitr)
# Creamos un data frame con los valores de AIC y BIC
df <- data.frame(Modelo = c("log1", "log2", "log3"),
                 AIC = c(AIC(log1), AIC(log2), AIC(log3)),
                 BIC = c(BIC(log1), BIC(log2), BIC(log3)))

# Generamos la tabla con una columna adicional
tabla <- knitr::kable(df, caption = "Comparison of AIC and BIC", align = c("l", "r", "r"))

# Se muestra la tabla
print(tabla)

Comparison of AIC and BIC

Modelo	AIC	BIC
log1	109.4642	154.6998
log2	107.5234	148.2354
log3	105.6618	141.8503

AIC y BIC tiene valores bajos, lo que indica que estos modelos se ajustan mejor a los datos.
De acuerdo a la tabla, el modelo “log3” tiene los valores más bajos de AIC y BIC en comparación con los otros modelo(log1 y log2), esto sugiere que el modelo “log3” se ajusta mejor en términos de complejidad del modelo y ajuste a los datos.

Las pruebas de comparación, nos indican que el modelo log3 no prensenta una mejora significativa según el test de annova, lo cual nos lleva a observar los estimadores AIC y BIC y de acuerdo a esto, este modelo tiene los menores valores, por lo tanto es el mejor modelo.

3. Analice la significancia del modelo obtenido luego del proceso de selección, y responda si:

• ¿Es el modelo obtenido significativo?

• ¿Existe alguna covariable no significativa?

• ¿En caso de existir alguna covariable no significativa, la quitar ́ıa del modelo?. Fundamente.

  Si, las retiraría del modelo, para:

  • Facilitad de la interpretación del modelo y evita el riesgo de sobreajuste.

  • Tener un modelo más simple con la menor cantidad de covariables necesarias para explicar el fenómeno de interés.

  • Reducir los grados de libertad y se mejorar la estimación de los coeficientes de las covariables restantes

• RESPUESTAS:

  • ¿Es el modelo obtenido significativo?

    Respuesta:  Si, ya que el modelo log3 es significativo en función del valor de bondad de ajuste (residual deviance) y los estadísticos AIC y BIC.

  • ¿Existe alguna covariable no significativa?

    Respuesta: De acuerdo al análisis realizado a cada modelo con la función ANOVA, obtenemos que las covariables Usize (p_valor 0,807802) y Ephith (p_valor 0,6997625) no son significativas.

  • ¿En caso de existir alguna covariable no significativa, la quitaría del modelo?. Fundamente.

    Respuesta: Si, las retiraríamos del modelo, porque:

    • Nos Facilita la interpretación del modelo y evita el riesgo de sobreajuste.

    • Nos permite tener un modelo más simple, con la menor cantidad de covariables necesarias para explicar el fenómeno de interés.

    • Permite reducir los grados de libertad y mejora la estimación de los coeficientes de las covariables restantes.

4. Fijar al menos 5 valores para las covariables del modelo y con ellas realizar la predicción de la
media y la predicción individual de la variable objetivo (incluir los intervalos de confianza).

# Seleccionamos 5 observaciones del dataset
sample_data_log <- wbca[sample.int(681, 5),]
sample_data_log

# Fijamos los valores de las covariables
new_data_log <- data.frame(Adhes = sample_data_log$Adhes,
                       BNucl = sample_data_log$BNucl,
                       Chrom  = sample_data_log$Chrom ,
                       Mitos  = sample_data_log$Mitos ,
                       NNucl  = sample_data_log$NNucl ,
                       Thick  = sample_data_log$Thick,
                       UShap  = sample_data_log$UShap   )
new_data_log

#realiza la predicción de valores utilizando un modelo de regresión logística
pred <- predict(log3, new_data_log, type="response", se.fit= TRUE)
pred$fit

        1         2         3         4         5 
0.9987698 0.2325307 0.9840732 0.9975856 0.9854111 

pred_end <- ifelse(pred$fit>0.8, "1", "0")
pred_end

  1   2   3   4   5 
"1" "0" "1" "1" "1" 

De acuerdo a este analisis los resultados en la predicción estan correctos.

#Se realizan los calculos de los intervalos de confianza
#Metodo 01
x =data.frame(puntual=pred$fit, LI= pred$fit - 2*pred$se.fit, LS= pred$fit + 2*pred$se.fit)
x

#Metodo 2 - empledo del exponente
pred2 <- predict(log3, new_data_log, type="link", se.fit= TRUE)
exp(data.frame(puntual=pred2$fit, LI= pred2$fit - 2*pred2$se.fit, LS= pred2$fit + 2*pred2$se.fit))/(1+exp(data.frame(puntual=pred2$fit, LI= pred2$fit - 2*pred2$se.fit, LS= pred2$fit + 2*pred2$se.fit)))

Los cálculos se llevan a cabo mediante dos enfoques. En ambos se utiliza el valor estimado (valor de la predicción); en el segundo enfoque (más comúnmente utilizado), se utiliza el exponente para restringir que los valores de los límites superior e inferior no sean inferiores a 0 ni superiores a 1.