Tarea Visualización de Datos en R

Errores Comunes

library(readr)

## Warning: package 'readr' was built under R version 4.1.3

library(mice)

## Warning: package 'mice' was built under R version 4.1.3

## 
## Attaching package: 'mice'

## The following object is masked from 'package:stats':
## 
##     filter

## The following objects are masked from 'package:base':
## 
##     cbind, rbind

library(outliers)

## Warning: package 'outliers' was built under R version 4.1.3

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.1.3

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.1.3

## Warning: package 'ggplot2' was built under R version 4.1.3

## Warning: package 'tibble' was built under R version 4.1.3

## Warning: package 'tidyr' was built under R version 4.1.3

## Warning: package 'purrr' was built under R version 4.1.3

## Warning: package 'stringr' was built under R version 4.1.3

## Warning: package 'forcats' was built under R version 4.1.3

## Warning: package 'lubridate' was built under R version 4.1.3

## -- Attaching core tidyverse packages ------------------------ tidyverse 2.0.0 --
## v forcats   1.0.0     v stringr   1.5.0
## v ggplot2   3.4.1     v tibble    3.2.1
## v lubridate 1.9.2     v tidyr     1.3.0
## v purrr     1.0.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks mice::filter(), stats::filter()
## x dplyr::lag()    masks stats::lag()
## i Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(ggplot2)
library(EnvStats)

## 
## Attaching package: 'EnvStats'
## 
## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm

library(GGally)

## Warning: package 'GGally' was built under R version 4.1.3

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

df <- read.csv("C:/Users/pc/Desktop/Universidad/Semestre_VI/VIZ/Tarea2/Errores_comunes/diabetes.csv")
head(df)

##   Pregnancies Glucose BloodPressure SkinThickness Insulin  BMI
## 1           6     148            72            35       0 33.6
## 2           1      85            66            29       0 26.6
## 3           8     183            64             0       0 23.3
## 4           1      89            66            23      94 28.1
## 5           0     137            40            35     168 43.1
## 6           5     116            74             0       0 25.6
##   DiabetesPedigreeFunction Age Outcome
## 1                    0.627  50       1
## 2                    0.351  31       0
## 3                    0.672  32       1
## 4                    0.167  21       0
## 5                    2.288  33       1
## 6                    0.201  30       0

df <- df %>%
  mutate(across(c(Glucose, BloodPressure, SkinThickness, Insulin, BMI), ~na_if(.x, 0)))

summarize_na <- sapply(df, function(x) sum(is.na(x)))
summary_stats <- summary(df)

print(summarize_na)

##              Pregnancies                  Glucose            BloodPressure 
##                        0                        5                       35 
##            SkinThickness                  Insulin                      BMI 
##                      227                      374                       11 
## DiabetesPedigreeFunction                      Age                  Outcome 
##                        0                        0                        0

print(summary_stats)

##   Pregnancies        Glucose      BloodPressure    SkinThickness  
##  Min.   : 0.000   Min.   : 44.0   Min.   : 24.00   Min.   : 7.00  
##  1st Qu.: 1.000   1st Qu.: 99.0   1st Qu.: 64.00   1st Qu.:22.00  
##  Median : 3.000   Median :117.0   Median : 72.00   Median :29.00  
##  Mean   : 3.845   Mean   :121.7   Mean   : 72.41   Mean   :29.15  
##  3rd Qu.: 6.000   3rd Qu.:141.0   3rd Qu.: 80.00   3rd Qu.:36.00  
##  Max.   :17.000   Max.   :199.0   Max.   :122.00   Max.   :99.00  
##                   NA's   :5       NA's   :35       NA's   :227    
##     Insulin            BMI        DiabetesPedigreeFunction      Age       
##  Min.   : 14.00   Min.   :18.20   Min.   :0.0780           Min.   :21.00  
##  1st Qu.: 76.25   1st Qu.:27.50   1st Qu.:0.2437           1st Qu.:24.00  
##  Median :125.00   Median :32.30   Median :0.3725           Median :29.00  
##  Mean   :155.55   Mean   :32.46   Mean   :0.4719           Mean   :33.24  
##  3rd Qu.:190.00   3rd Qu.:36.60   3rd Qu.:0.6262           3rd Qu.:41.00  
##  Max.   :846.00   Max.   :67.10   Max.   :2.4200           Max.   :81.00  
##  NA's   :374      NA's   :11                                              
##     Outcome     
##  Min.   :0.000  
##  1st Qu.:0.000  
##  Median :0.000  
##  Mean   :0.349  
##  3rd Qu.:1.000  
##  Max.   :1.000  
## 

hist(df$Glucose, main = "Histograma de Glucosa", xlab = "Glucosa", border = "darkblue", col = "lightblue")

boxplot(list(df$Glucose, df$BMI, df$BloodPressure, df$SkinThickness), main = "Boxplot de Variables", col = "lightgreen")

boxplot(df$Insulin, col = "lightblue")

barplot(table(df$Outcome), main = "Distribución de Outcome", xlab = "Outcome", ylab = "Frecuencia", col = "salmon")

Imputación de datos

correlation_matrix <- cor(df %>% select(-Outcome) %>% na.omit(), use = "complete.obs")
corrplot::corrplot(correlation_matrix, method = "circle")

imputed_data_pmm <- mice(df, method = 'pmm', m=1, maxit = 5)

## 
##  iter imp variable
##   1   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   2   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   3   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   4   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   5   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI

# Aplicar imputación por normal prediction
imputed_data_norm_predict <- mice(df, method = 'norm.predict', m=1, maxit = 5)

## 
##  iter imp variable
##   1   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   2   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   3   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   4   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   5   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI

# Aplicar imputación por normal con observaciones completas (no bayesiano)
imputed_data_norm_nob <- mice(df, method = 'norm.nob', m=1, maxit = 5)

## 
##  iter imp variable
##   1   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   2   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   3   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   4   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   5   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI

# Aplicar imputación por normal (bayesiano)
imputed_data_norm <- mice(df, method = 'norm', m=1, maxit = 5)

## 
##  iter imp variable
##   1   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   2   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   3   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   4   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI
##   5   1  Glucose  BloodPressure  SkinThickness  Insulin  BMI

completed_data_pmm <- complete(imputed_data_pmm)
completed_data_norm_predict <- complete(imputed_data_norm_predict)
completed_data_norm_nob <- complete(imputed_data_norm_nob)
completed_data_norm <- complete(imputed_data_norm)

En esta parte se hizo la imputación de datos al dataframe df con los métodos estudiados en clase (pmm, norm predict, norm nob y norm)

Aquí hay algunas interpretaciones clave que podemos inferir del gráfico:

Círculos más grandes y oscuros: Indican una correlación más fuerte entre las variables correspondientes. Un círculo grande y oscuro en la diagonal superior izquierda (como el que aparece en la intersección de Glucose con Insulin) sugeriría una correlación positiva significativa entre estas variables. Esto podría significar que los niveles más altos de glucosa están asociados con niveles más altos de insulina.

Círculos pequeños o ausentes: Indican una correlación débil o inexistente. Por ejemplo, la correlación entre SkinThickness y Pregnancies parece ser muy baja, ya que el círculo es muy pequeño y claro.

Círculos azules claros: En general, los tonos más claros representan valores más bajos de correlación, y los tonos más oscuros indican valores más altos. No parece haber círculos negativos significativos, lo que implicaría una correlación negativa fuerte entre dos variables.

Relación con la variable de resultado: Dado que esta es una matriz de correlación y la variable de resultado Outcome (presencia o ausencia de diabetes) no está incluida, no se puede inferir nada directamente sobre la relación entre estas variables y el riesgo de diabetes.

Diagonal: La diagonal de la matriz de correlación siempre tendrá círculos oscuros y de tamaño máximo, ya que representa la correlación de cada variable consigo misma, lo cual es siempre perfectamente positivo (correlación de 1).

Implicaciones prácticas: En el contexto de la diabetes, se espera que ciertas variables como la glucosa y el IMC tengan una correlación positiva con otras medidas metabólicas. Las fuertes correlaciones positivas podrían ser útiles para predecir el riesgo de diabetes, mientras que las correlaciones débiles indican que las variables proporcionan información más independiente entre sí.

summary(df, na.rm = TRUE) 

##   Pregnancies        Glucose      BloodPressure    SkinThickness  
##  Min.   : 0.000   Min.   : 44.0   Min.   : 24.00   Min.   : 7.00  
##  1st Qu.: 1.000   1st Qu.: 99.0   1st Qu.: 64.00   1st Qu.:22.00  
##  Median : 3.000   Median :117.0   Median : 72.00   Median :29.00  
##  Mean   : 3.845   Mean   :121.7   Mean   : 72.41   Mean   :29.15  
##  3rd Qu.: 6.000   3rd Qu.:141.0   3rd Qu.: 80.00   3rd Qu.:36.00  
##  Max.   :17.000   Max.   :199.0   Max.   :122.00   Max.   :99.00  
##                   NA's   :5       NA's   :35       NA's   :227    
##     Insulin            BMI        DiabetesPedigreeFunction      Age       
##  Min.   : 14.00   Min.   :18.20   Min.   :0.0780           Min.   :21.00  
##  1st Qu.: 76.25   1st Qu.:27.50   1st Qu.:0.2437           1st Qu.:24.00  
##  Median :125.00   Median :32.30   Median :0.3725           Median :29.00  
##  Mean   :155.55   Mean   :32.46   Mean   :0.4719           Mean   :33.24  
##  3rd Qu.:190.00   3rd Qu.:36.60   3rd Qu.:0.6262           3rd Qu.:41.00  
##  Max.   :846.00   Max.   :67.10   Max.   :2.4200           Max.   :81.00  
##  NA's   :374      NA's   :11                                              
##     Outcome     
##  Min.   :0.000  
##  1st Qu.:0.000  
##  Median :0.000  
##  Mean   :0.349  
##  3rd Qu.:1.000  
##  Max.   :1.000  
## 

Aquí hay algunas observaciones:

Pregnancies: El número de embarazos varía de 0 a 17, con una media de aproximadamente 3.85, lo que indica una distribución con un sesgo hacia las mujeres con menos embarazos.

Glucose: Los niveles de glucosa varían de 44 a 199, con una media alrededor de 121.7, reflejando la importancia de esta medida en la evaluación del riesgo de diabetes.

BloodPressure: La presión arterial tiene una media de 72.41, con valores que van de 24 a 122, lo que sugiere una variación normal en la población.

SkinThickness: El grosor de la piel muestra una amplia gama de valores, desde 7 hasta 99, pero más de 200 registros tienen valores faltantes (NA’s), indicando posibles problemas de datos.

Insulin: Hay un número significativo de valores faltantes (374), y los niveles varían ampliamente, de 14 a 846, lo que podría indicar variaciones en las condiciones de salud de los sujetos.

BMI: El Índice de Masa Corporal (BMI) varía de 18.2 a 67.1, con una mediana cercana a la clasificación de sobrepeso, lo que es relevante para el estudio de la diabetes.

DiabetesPedigreeFunction: Una medida de la predisposición genética a la diabetes, con un rango de -0.078 a 2.420, que parece contener un error dado el valor mínimo negativo.

Age: Las edades van de 21 a 81 años, con la mediana en los 29 años, indicando una población relativamente joven.

Outcome: Es una variable binaria que indica la presencia o ausencia de diabetes, con la mayoría de los sujetos sin diabetes (valor 0).

IMPUTACIÓN CON PMM

summary(completed_data_pmm) # Después de la imputación con PMM

##   Pregnancies        Glucose      BloodPressure    SkinThickness  
##  Min.   : 0.000   Min.   : 44.0   Min.   : 24.00   Min.   : 7.00  
##  1st Qu.: 1.000   1st Qu.: 99.0   1st Qu.: 64.00   1st Qu.:21.00  
##  Median : 3.000   Median :117.0   Median : 72.00   Median :29.00  
##  Mean   : 3.845   Mean   :121.7   Mean   : 72.13   Mean   :28.94  
##  3rd Qu.: 6.000   3rd Qu.:141.0   3rd Qu.: 80.00   3rd Qu.:36.00  
##  Max.   :17.000   Max.   :199.0   Max.   :122.00   Max.   :99.00  
##     Insulin           BMI        DiabetesPedigreeFunction      Age       
##  Min.   : 14.0   Min.   :18.20   Min.   :0.0780           Min.   :21.00  
##  1st Qu.: 75.0   1st Qu.:27.40   1st Qu.:0.2437           1st Qu.:24.00  
##  Median :119.0   Median :32.15   Median :0.3725           Median :29.00  
##  Mean   :147.5   Mean   :32.40   Mean   :0.4719           Mean   :33.24  
##  3rd Qu.:183.0   3rd Qu.:36.60   3rd Qu.:0.6262           3rd Qu.:41.00  
##  Max.   :846.0   Max.   :67.10   Max.   :2.4200           Max.   :81.00  
##     Outcome     
##  Min.   :0.000  
##  1st Qu.:0.000  
##  Median :0.000  
##  Mean   :0.349  
##  3rd Qu.:1.000  
##  Max.   :1.000

Pregnancies: El número de embarazos varía ampliamente (0 a 17), con una mediana de 3, lo que es típico para poblaciones generales.

Glucose: Los valores oscilan entre 44 y 199, con una mediana de 117, lo cual es esperado en un conjunto de datos que incluye individuos con y sin diabetes.

BloodPressure: La presión arterial tiene un mínimo de 24 mm Hg y un máximo de 122 mm Hg, y una mediana de 72 mm Hg, indicando una distribución razonable para un conjunto de datos de adultos.

SkinThickness: El grosor de la piel varía de 7 a 99 mm, aunque con una mediana de 28 mm. Los datos imputados parecen haber llenado los valores previamente ausentes, con una distribución plausible.

Insulin: La imputación de insulina muestra un rango amplio (14 a 846), que puede reflejar tanto variabilidad biológica como la imputación de valores que representan casos de hiperinsulinemia.

BMI: Con valores que van de 18.2 a 67.1, el conjunto de datos cubre desde pesos saludables hasta obesidad severa. La mediana de 32.5 sugiere sobrepeso.

DiabetesPedigreeFunction: Esta función tiene un rango interesante de -0.078 a 2.420. El valor mínimo negativo podría ser un error de datos, ya que se esperaría que este valor siempre fuera positivo.

Age: La edad de los sujetos va de 21 a 81 años, con una mediana de 29, indicando que el conjunto de datos tiende hacia una población más joven.

Outcome: Como variable binaria, muestra que la mayoría de los individuos no tienen diabetes (0). La media (0.349) sugiere que alrededor del 34.9% de los individuos en este conjunto de datos tienen diabetes, lo cual es consistente con la variable siendo binaria y posiblemente desbalanceada hacia el grupo no diabético.

Este resumen nos indica que las técnicas de imputación PMM han sido aplicadas y han rellenado los valores faltantes.

densityplot(imputed_data_pmm)

Glucosa: Las curvas azul y roja muestran una superposición considerable, lo que indica que la imputación PMM ha mantenido la distribución original de los datos con una ligera elevación en la densidad post-imputación. Este cambio leve sugiere un mínimo ajuste en la centralidad de los datos, posiblemente acercando los valores imputados a la media o mediana de los datos observados.

Presión Arterial: La similitud entre las dos curvas es notable, aunque la curva roja presenta un pico ligeramente más pronunciado. Esto podría interpretarse como una tendencia de los valores imputados a agruparse alrededor de una medida central, reduciendo potencialmente la variabilidad en los datos imputados comparados con los originales.

Grosor de la Piel: Las curvas roja y azul aquí también coinciden en gran medida, indicando una preservación efectiva de la distribución original tras la imputación. Esto es indicativo de una buena práctica de imputación, ya que no parece haber un impacto significativo en la forma general de la distribución de los datos.

Insulina: La imputación parece haber generado un leve incremento en la frecuencia de valores más bajos, como se evidencia por la altura aumentada de la curva roja en el extremo izquierdo. Esta observación puede requerir un análisis más detallado para entender cómo la imputación PMM está afectando los valores más bajos de insulina.

BMI: Ambas curvas mantienen una congruencia en su forma, aunque con una densidad ligeramente mayor en la curva roja en torno al centro de la distribución. Este efecto puede ser un reflejo del método PMM, que tiende a seleccionar valores cercanos a los originales, minimizando así la alteración de la distribución.

Conclusión

La técnica de PMM ha mostrado ser efectiva en la conservación de la forma general de las distribuciones de las variables médicas analizadas. A pesar de las leves diferencias en densidad y altura de las curvas, no hay evidencia de distorsiones significativas que pudieran comprometer los análisis subsecuentes.

IMPUTACIÓN CON NORM

summary(imputed_data_norm)

## Class: mids
## Number of multiple imputations:  1 
## Imputation methods:
##              Pregnancies                  Glucose            BloodPressure 
##                       ""                   "norm"                   "norm" 
##            SkinThickness                  Insulin                      BMI 
##                   "norm"                   "norm"                   "norm" 
## DiabetesPedigreeFunction                      Age                  Outcome 
##                       ""                       ""                       "" 
## PredictorMatrix:
##               Pregnancies Glucose BloodPressure SkinThickness Insulin BMI
## Pregnancies             0       1             1             1       1   1
## Glucose                 1       0             1             1       1   1
## BloodPressure           1       1             0             1       1   1
## SkinThickness           1       1             1             0       1   1
## Insulin                 1       1             1             1       0   1
## BMI                     1       1             1             1       1   0
##               DiabetesPedigreeFunction Age Outcome
## Pregnancies                          1   1       1
## Glucose                              1   1       1
## BloodPressure                        1   1       1
## SkinThickness                        1   1       1
## Insulin                              1   1       1
## BMI                                  1   1       1

Las variables SkinThickness, Insulin, y BMI han sido imputadas utilizando métodos basados en la distribución normal (“norm”), lo que sugiere que la imputación para estas variables se hizo asumiendo que seguían una distribución aproximadamente normal.

Las variables Pregnancies, Glucose, BloodPressure, DiabetesPedigreeFunction, Age, y Outcome no han sido imputadas, lo cual está indicado por “” en los métodos de imputación.

La PredictorMatrix muestra con qué variables se ha realizado la imputación para cada variable imputada. Un “1” indica que la variable de esa fila se utilizó como predictor en la imputación de la variable de esa columna. Por ejemplo, Pregnancies fue utilizada como predictor para SkinThickness, Insulin, y BMI, pero no para sí misma (indicado con un “0”).

densityplot(imputed_data_norm)

IMPUTACIÓN CON NORM NOB

Las distribuciones pre y post-imputación parecen similares, lo que indicaría que el método de imputación preserva la estructura general de los datos.

No hay cambios drásticos en las formas de las distribuciones, lo que sugiere que no se introducen sesgos significativos.

En la imputación con el método NORM, las curvas de densidad de los datos imputados se alinean estrechamente con los datos originales, lo que sugiere una buena conservación de las propiedades estadísticas de los datos.

summary(imputed_data_norm_nob)

## Class: mids
## Number of multiple imputations:  1 
## Imputation methods:
##              Pregnancies                  Glucose            BloodPressure 
##                       ""               "norm.nob"               "norm.nob" 
##            SkinThickness                  Insulin                      BMI 
##               "norm.nob"               "norm.nob"               "norm.nob" 
## DiabetesPedigreeFunction                      Age                  Outcome 
##                       ""                       ""                       "" 
## PredictorMatrix:
##               Pregnancies Glucose BloodPressure SkinThickness Insulin BMI
## Pregnancies             0       1             1             1       1   1
## Glucose                 1       0             1             1       1   1
## BloodPressure           1       1             0             1       1   1
## SkinThickness           1       1             1             0       1   1
## Insulin                 1       1             1             1       0   1
## BMI                     1       1             1             1       1   0
##               DiabetesPedigreeFunction Age Outcome
## Pregnancies                          1   1       1
## Glucose                              1   1       1
## BloodPressure                        1   1       1
## SkinThickness                        1   1       1
## Insulin                              1   1       1
## BMI                                  1   1       1

Indica que se ha realizado una única imputación y se han utilizado métodos de imputación normal no bayesiana (“norm.nob”) para las variables Glucose, BloodPressure, SkinThickness, y BMI. La matriz predictora revela qué variables se utilizaron como predictores en la imputación de cada una: un “1” indica que la variable de la fila se usó para predecir la variable de la columna, mientras que un “0” indica que no se utilizó.

densityplot(imputed_data_norm_nob)

Las curvas azules y rojas podrían representar los datos antes y después de la imputación, respectivamente. La similitud en la forma de las curvas sugiere que la imputación ha mantenido la distribución general de los datos originales. Además, no se observan cambios significativos en la distribución, lo que indica que el método de imputación ha sido efectivo y no ha introducido distorsiones notables en las variables imputadas.

IMPUTACIÓN CON NORM PREDICT

summary(imputed_data_norm_predict)

## Class: mids
## Number of multiple imputations:  1 
## Imputation methods:
##              Pregnancies                  Glucose            BloodPressure 
##                       ""           "norm.predict"           "norm.predict" 
##            SkinThickness                  Insulin                      BMI 
##           "norm.predict"           "norm.predict"           "norm.predict" 
## DiabetesPedigreeFunction                      Age                  Outcome 
##                       ""                       ""                       "" 
## PredictorMatrix:
##               Pregnancies Glucose BloodPressure SkinThickness Insulin BMI
## Pregnancies             0       1             1             1       1   1
## Glucose                 1       0             1             1       1   1
## BloodPressure           1       1             0             1       1   1
## SkinThickness           1       1             1             0       1   1
## Insulin                 1       1             1             1       0   1
## BMI                     1       1             1             1       1   0
##               DiabetesPedigreeFunction Age Outcome
## Pregnancies                          1   1       1
## Glucose                              1   1       1
## BloodPressure                        1   1       1
## SkinThickness                        1   1       1
## Insulin                              1   1       1
## BMI                                  1   1       1

Se emplearon métodos basados en la predicción normal (“norm.predict”) para las variables Glucose, BloodPressure, SkinThickness, Insulin, y BMI. No se aplicó imputación a Pregnancies, DiabetesPedigreeFunction, Age, y Outcome. La matriz de predictores muestra que todas las variables disponibles fueron utilizadas como predictores en la imputación de cada variable, lo cual es típico para un enfoque de imputación múltiple que intenta aprovechar toda la información disponible en el conjunto de datos.

densityplot(imputed_data_norm_predict)

Las líneas azules representan la distribución original de los datos antes de la imputación, y las rojas después de la imputación.

Para Glucose, BloodPressure, y SkinThickness, las curvas después de la imputación siguen estando muy alineadas con las curvas originales, lo que sugiere que la imputación no ha alterado significativamente la distribución subyacente de estos datos.

Para Insulin y BMI, hay una marcada diferencia entre las curvas pre y post-imputación, particularmente en BMI donde la curva roja post-imputación tiene picos múltiples. Esto podría indicar que la imputación ha introducido o exacerbado ciertas irregularidades en la distribución de estos datos. Sería importante investigar más a fondo por qué ocurre esto, ya que podría deberse a la naturaleza de los valores faltantes o a particularidades en la técnica de imputación utilizada.

PRUEBAS

Para las pruebas utilizaremos el dataframe completed data norm, para analizar los datos con uno de los métodos de imputación y su comportamiento, se debe tener en cuenta que puede que cambie dependiendo del método seleccionado, para este caso utilizaremos este como ejemplo.

PRUEBA GRUBBS

insulin_test <- grubbs.test(completed_data_norm$Insulin, opposite = F)
insulin_test

## 
##  Grubbs test for one outlier
## 
## data:  completed_data_norm$Insulin
## G = 5.74304, U = 0.95694, p-value = 2.48e-06
## alternative hypothesis: highest value 846 is an outlier

Al nivel de significancia del 5%, concluimos que el valor más alto 846 es un valor atípico

preg_test <- grubbs.test(completed_data_norm$Pregnancies, opposite = F)
preg_test

## 
##  Grubbs test for one outlier
## 
## data:  completed_data_norm$Pregnancies
## G = 3.9040, U = 0.9801, p-value = 0.03367
## alternative hypothesis: highest value 17 is an outlier

Al nivel de significancia del 5%, concluimos que el valor más alto 17 es un valor atípico

glu_test <- grubbs.test(completed_data_norm$Glucose, opposite = F)
glu_test

## 
##  Grubbs test for one outlier
## 
## data:  completed_data_norm$Glucose
## G = 2.54972, U = 0.99151, p-value = 1
## alternative hypothesis: lowest value 44 is an outlier

Al nivel de significancia del 5%, concluimos que no tiene un valor atípico

bp_test <- grubbs.test(completed_data_norm$BloodPressure, opposite = T)
bp_test_ <- grubbs.test(completed_data_norm$BloodPressure, opposite = F)

print(list(bp_test, bp_test_))

## [[1]]
## 
##  Grubbs test for one outlier
## 
## data:  completed_data_norm$BloodPressure
## G = 3.90490, U = 0.98009, p-value = 0.03355
## alternative hypothesis: lowest value 24 is an outlier
## 
## 
## [[2]]
## 
##  Grubbs test for one outlier
## 
## data:  completed_data_norm$BloodPressure
## G = 4.02776, U = 0.97882, p-value = 0.01984
## alternative hypothesis: highest value 122 is an outlier

Al nivel de significancia del 5%, concluimos que el valor más alto 122 es un valor atípico y el valor más bajo 24 es un valor atípico

st_test <- grubbs.test(completed_data_norm$SkinThickness, opposite = F)
st_test

## 
##  Grubbs test for one outlier
## 
## data:  completed_data_norm$SkinThickness
## G = 6.71235, U = 0.94118, p-value = 3.693e-09
## alternative hypothesis: highest value 99 is an outlier

Al nivel de significancia del 5%, concluimos que el valor más alto 99 es un valor atípico

bmi_test <- grubbs.test(completed_data_norm$BMI, opposite = F)
bmi_test

## 
##  Grubbs test for one outlier
## 
## data:  completed_data_norm$BMI
## G = 5.00540, U = 0.96729, p-value = 0.0001738
## alternative hypothesis: highest value 67.1 is an outlier

Al nivel de significancia del 5%, concluimos que el valor más alto 67.1 es un valor atípico

dpf_test <- grubbs.test(completed_data_norm$DiabetesPedigreeFunction, opposite = F)
dpf_test

## 
##  Grubbs test for one outlier
## 
## data:  completed_data_norm$DiabetesPedigreeFunction
## G = 5.87973, U = 0.95487, p-value = 1.056e-06
## alternative hypothesis: highest value 2.42 is an outlier

Al nivel de significancia del 5%, concluimos que el valor más alto 2.42 es un valor atípico

age_test <- grubbs.test(completed_data_norm$Age, opposite = F)
age_test

## 
##  Grubbs test for one outlier
## 
## data:  completed_data_norm$Age
## G = 4.06107, U = 0.97847, p-value = 0.01716
## alternative hypothesis: highest value 81 is an outlier

Al nivel de significancia del 5%, concluimos que el valor más alto 81 es un valor atípico

PRUEBA DIXON

subdata <- completed_data_norm[1:25,] # Reducir el dataset para hacer la prueba de dixon

preg_test <- dixon.test(subdata$Pregnancies, opposite = F)
preg_test

## 
##  Dixon test for outliers
## 
## data:  subdata$Pregnancies
## Q = 0.1, p-value = 0.3494
## alternative hypothesis: highest value 11 is an outlier

Con un nivel de significancia del 5% tenemos no tiene valores atípicos

glu_test <- dixon.test(subdata$Glucose, opposite = F)
glu_test

## 
##  Dixon test for outliers
## 
## data:  subdata$Glucose
## Q = 0.074074, p-value = 0.2051
## alternative hypothesis: highest value 197 is an outlier

Con un nivel de significancia del 5% tenemos no tiene valores atípicos

bp_test <- dixon.test(subdata$BloodPressure, opposite = F)
bp_test

## 
##  Dixon test for outliers
## 
## data:  subdata$BloodPressure
## Q = 0.32258, p-value = 0.3247
## alternative hypothesis: lowest value 30 is an outlier

Con un nivel de significancia del 5% tenemos no tiene valores atípicos

st_test <- dixon.test(subdata$SkinThickness, opposite = F)
st_test

## 
##  Dixon test for outliers
## 
## data:  subdata$SkinThickness
## Q = 0.29664, p-value = 0.4319
## alternative hypothesis: lowest value 8.03438793546078 is an outlier

Con un nivel de significancia del 5% tenemos no tiene valores atípicos

insulin_test <- dixon.test(subdata$Insulin, opposite = F)
insulin_test

## 
##  Dixon test for outliers
## 
## data:  subdata$Insulin
## Q = 0.42367, p-value = 0.07381
## alternative hypothesis: highest value 846 is an outlier

Con un nivel de significancia del 5% tenemos que el valor máximo 846 es un valor atípico

bmi_test <- dixon.test(subdata$BMI, opposite = F)
bmi_test

## 
##  Dixon test for outliers
## 
## data:  subdata$BMI
## Q = 0.135, p-value = 0.5796
## alternative hypothesis: highest value 45.8 is an outlier

Con un nivel de significancia del 5% tenemos no tiene valores atípicos

dpf_test <- dixon.test(subdata$DiabetesPedigreeFunction, opposite = F)
dpf_test

## 
##  Dixon test for outliers
## 
## data:  subdata$DiabetesPedigreeFunction
## Q = 0.74682, p-value < 2.2e-16
## alternative hypothesis: highest value 2.288 is an outlier

Con un nivel de significancia del 5% tenemos que el valor máximo 2.288 es un valor atípico

age_test <- dixon.test(subdata$Age, opposite = F)
age_test

## 
##  Dixon test for outliers
## 
## data:  subdata$Age
## Q = 0.15625, p-value = 0.7207
## alternative hypothesis: highest value 59 is an outlier

Con un nivel de significancia del 5% tenemos no tiene valores atípicos

PRUEBA ROSNER

all_stats_results <- list()

for (column_name in names(completed_data_norm)[-which(names(completed_data_norm) == "DiabetesPedigreeFunction")]) {
  if (is.numeric(completed_data_norm[[column_name]])) {
    test <- tryCatch({
      rosnerTest(completed_data_norm[[column_name]], k = 6)
    }, error = function(e) { 
      return(NULL) 
    })

    if (!is.null(test) && !is.null(test$all.stats)) {
      all_stats_results[[column_name]] <- test$all.stats
    } else {
      all_stats_results[[column_name]] <- NA
    }
  } else {
    all_stats_results[[column_name]] <- NA
  }
}

print(all_stats_results)

## $Pregnancies
##   i   Mean.i     SD.i Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 3.845052 3.369578    17     160 3.904034   3.974092   FALSE
## 2 1 3.827901 3.338063    15      89 3.346880   3.973762   FALSE
## 3 2 3.813316 3.315699    14     299 3.072258   3.973432   FALSE
## 4 3 3.800000 3.297310    14     456 3.093431   3.973102   FALSE
## 5 4 3.786649 3.278714    13      29 2.810050   3.972771   FALSE
## 6 5 3.774574 3.263821    13      73 2.826572   3.972440   FALSE
## 
## $Glucose
##   i   Mean.i     SD.i Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 121.8166 30.51965    44      63 2.549720   3.974092   FALSE
## 2 1 121.9180 30.40969   199     662 2.534784   3.973762   FALSE
## 3 2 121.8174 30.30150   198     562 2.514153   3.973432   FALSE
## 4 3 121.7178 30.19564   197       9 2.493149   3.973102   FALSE
## 5 4 121.6193 30.09209   197     229 2.505002   3.972771   FALSE
## 6 5 121.5205 29.98759   197     409 2.517026   3.972440   FALSE
## 
## $BloodPressure
##   i   Mean.i     SD.i Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 72.24111 12.35399   122     107 4.027758   3.974092    TRUE
## 2 1 72.17624 12.23045    24     598 3.939041   3.973762   FALSE
## 3 2 72.23913 12.11369    30      19 3.486892   3.973432   FALSE
## 4 3 72.29434 12.02478    30     126 3.517267   3.973102   FALSE
## 5 4 72.34970 11.93471   114     692 3.489847   3.972771   FALSE
## 6 5 72.29512 11.84671   110      44 3.182730   3.972440   FALSE
## 
## $SkinThickness
##   i   Mean.i      SD.i       Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 28.88758 10.445283 99.00000000     580 6.712352   3.974092    TRUE
## 2 1 28.79617 10.140048 63.00000000     446 3.373142   3.973762   FALSE
## 3 2 28.75152 10.070934 60.80159303     194 3.182433   3.973432   FALSE
## 4 3 28.70962 10.010504 60.00000000      58 3.125754   3.973102   FALSE
## 5 4 28.66867  9.952720 57.87803163     155 2.934812   3.972771   FALSE
## 6 5 28.63039  9.902803  0.03784939      34 2.887318   3.972440   FALSE
## 
## $Insulin
##   i   Mean.i     SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 150.1245 121.1684 846.0000      14 5.743045   3.974092    TRUE
## 2 1 149.2173 118.6084 744.0000     229 5.014677   3.973762    TRUE
## 3 2 148.4408 116.7189 680.0000     248 4.554184   3.973432    TRUE
## 4 3 147.7459 115.1990 600.0000     585 3.925851   3.973102   FALSE
## 5 4 147.1540 114.1043 588.5249     580 3.868137   3.972771   FALSE
## 6 5 146.5755 113.0526 579.0000     410 3.824986   3.972440   FALSE
## 
## $BMI
##   i   Mean.i     SD.i Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 32.46809 6.918916  67.1     178 5.005396   3.974092    TRUE
## 2 1 32.42294 6.809265  59.4     446 3.961817   3.973762   FALSE
## 3 2 32.38772 6.743451  57.3     674 3.694292   3.973432   FALSE
## 4 3 32.35515 6.687321  55.0     126 3.386236   3.973102   FALSE
## 5 4 32.32551 6.641229  53.2     121 3.143166   3.972771   FALSE
## 6 5 32.29816 6.602364  52.9     304 3.120374   3.972440   FALSE
## 
## $Age
##   i   Mean.i     SD.i Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 33.24089 11.76023    81     460 4.061069   3.974092    TRUE
## 2 1 33.17862 11.64053    72     454 3.335018   3.973762   FALSE
## 3 2 33.12794 11.56315    70     667 3.188755   3.973432   FALSE
## 4 3 33.07974 11.49346    69     124 3.125278   3.973102   FALSE
## 5 4 33.03272 11.42714    69     685 3.147531   3.972771   FALSE
## 6 5 32.98558 11.36006    68     675 3.082239   3.972440   FALSE
## 
## $Outcome
##   i    Mean.i      SD.i Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 0.3489583 0.4769514     1       1 1.365006   3.974092   FALSE
## 2 1 0.3481095 0.4766818     1       3 1.367559   3.973762   FALSE
## 3 2 0.3472585 0.4764098     1       5 1.370126   3.973432   FALSE
## 4 3 0.3464052 0.4761355     1       7 1.372708   3.973102   FALSE
## 5 4 0.3455497 0.4758587     1       9 1.375304   3.972771   FALSE
## 6 5 0.3446920 0.4755795     1      10 1.377915   3.972440   FALSE

Viendo los resultados de la prueba de Rosner podemos observar que las variables Outcome, Pregnancies y Glucose no tienen valores atípicos, el resto podemos ver que si tienene especialmente la variable Insulin con 5 valores atípicos.

rosner_insulin <- rosnerTest(completed_data_norm$DiabetesPedigreeFunction, k = 12)

## Warning in rosnerTest(completed_data_norm$DiabetesPedigreeFunction, k = 12): The true Type I error may be larger than assumed.
## Although the help file for 'rosnerTest' has a table with information
## on the estimated Type I error level,
## simulations were not run for k > 10 or k > floor(n/2).

rosner_insulin$all.stats

##     i    Mean.i      SD.i Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1   0 0.4718763 0.3313286 2.420     446 5.879733   3.974092    TRUE
## 2   1 0.4693364 0.3239768 2.329     229 5.740114   3.973762    TRUE
## 3   2 0.4669086 0.3171301 2.288       5 5.742410   3.973432    TRUE
## 4   3 0.4645281 0.3104137 2.137     371 5.387880   3.973102    TRUE
## 5   4 0.4623390 0.3046509 1.893      46 4.696067   3.972771    TRUE
## 6   5 0.4604640 0.3004070 1.781      59 4.395823   3.972440    TRUE
## 7   6 0.4587310 0.2967633 1.731     372 4.287150   3.972108    TRUE
## 8   7 0.4570591 0.2933457 1.699     594 4.233710   3.971776    TRUE
## 9   8 0.4554250 0.2900522 1.698     622 4.283971   3.971443    TRUE
## 10  9 0.4537879 0.2867083 1.600     396 3.997834   3.971110    TRUE
## 11 10 0.4522757 0.2838528 1.476     331 3.606533   3.970776   FALSE
## 12 11 0.4509234 0.2815864 1.461     623 3.587093   3.970442   FALSE

Y para la variables Insulin podemos ver que es la que tiene mayor datos atípicos con un total de 10.

PRUEBA HAMPEL

lower_bound <- median(df$Outcome) - 3 * mad(df$Outcome, constant = 1)
lower_bound

## [1] 0

upper_bound <- median(df$Outcome) + 3 * mad(df$Outcome, constant = 1)
upper_bound

## [1] 0

outlier_ind <- which(df$Outcome < lower_bound | df$Outcome > upper_bound)
outlier_ind

##   [1]   1   3   5   7   9  10  12  14  15  16  17  18  20  23  24  25  26  27
##  [19]  32  38  39  40  44  46  49  54  57  62  65  67  71  73  79  85  89  94
##  [37] 100 101 110 111 112 115 116 117 121 125 126 129 130 131 132 133 144 153
##  [55] 155 156 160 165 166 171 172 176 178 180 186 187 188 189 190 193 194 196
##  [73] 198 199 200 207 208 210 214 215 216 217 219 220 221 222 228 231 232 236
##  [91] 237 238 239 243 244 246 255 256 260 262 265 267 270 271 277 281 284 285
## [109] 288 292 293 294 297 299 301 302 304 307 309 310 313 315 318 320 322 323
## [127] 324 327 329 333 338 339 340 350 356 357 358 360 361 364 367 370 371 376
## [145] 379 387 388 389 392 395 398 400 401 403 405 407 409 410 415 416 418 420
## [163] 425 426 428 430 436 441 444 445 446 449 452 456 459 469 477 481 485 486
## [181] 494 499 503 507 511 516 517 524 536 540 541 542 543 546 547 561 562 570
## [199] 578 580 581 585 587 589 591 593 596 599 604 605 607 612 613 615 619 620
## [217] 631 636 639 643 647 648 649 656 660 662 663 664 665 667 668 676 677 679
## [235] 682 684 690 692 694 696 697 702 703 707 709 710 713 716 717 720 723 731
## [253] 732 733 740 741 744 747 749 750 751 754 755 756 758 760 762 767

boxplot(completed_data_norm$Outcome)

hampel_outliers <- function(x, k = 10) {
  mediana <- median(x, na.rm = TRUE)
  mad <- mad(x, constant = 1, na.rm = TRUE) 
  outlier_threshold <- k * mad
  
  lower_bound <- mediana - outlier_threshold
  upper_bound <- mediana + outlier_threshold

  outlier_ind <- which(x < lower_bound | x > upper_bound)
  
  return(list(lower_bound = lower_bound, upper_bound = upper_bound, outliers = outlier_ind))
}

outliers_hampel <- list()
exclude_columns <- c("Outcome")

for (column_name in setdiff(names(completed_data_norm), exclude_columns)) {
  if (is.numeric(completed_data_norm[[column_name]])) {
    outliers_hampel[[column_name]] <- hampel_outliers(completed_data_norm[[column_name]], k = 3)
  }
}

print(outliers_hampel)

## $Pregnancies
## $Pregnancies$lower_bound
## [1] -3
## 
## $Pregnancies$upper_bound
## [1] 9
## 
## $Pregnancies$outliers
##  [1]   8  12  13  25  26  29  35  37  73  87  89 144 160 194 216 247 255 260 271
## [20] 275 282 299 307 324 328 334 358 359 376 437 456 459 465 506 511 519 543 559
## [39] 560 579 583 591 615 635 636 649 659 661 668 673 692 707 713 718 741 745 746
## [58] 764
## 
## 
## $Glucose
## $Glucose$lower_bound
## [1] 57
## 
## $Glucose$upper_bound
## [1] 177
## 
## $Glucose$outliers
##  [1]   3   9  14  23  41  46  57  63 155 176 186 187 207 210 213 229 237 238 246
## [20] 259 261 318 320 328 333 340 360 361 400 409 426 428 441 446 490 499 507 546
## [39] 547 550 562 580 596 605 607 623 648 662 676 681 716 749 754 760
## 
## 
## $BloodPressure
## $BloodPressure$lower_bound
## [1] 48
## 
## $BloodPressure$upper_bound
## [1] 96
## 
## $BloodPressure$outliers
##  [1]   5  19  44  81  85 107 126 178 188 208 304 347 363 370 380 388 441 465 550
## [20] 576 577 598 600 659 663 673 674 692 698 708 742
## 
## 
## $SkinThickness
## $SkinThickness$lower_bound
## [1] 8
## 
## $SkinThickness$upper_bound
## [1] 50
## 
## $SkinThickness$outliers
##  [1]  34  58  68  87 100 103 117 121 155 194 197 212 255 276 304 379 446 520 533
## [20] 580 700
## 
## 
## $Insulin
## $Insulin$lower_bound
## [1] -63.49687
## 
## $Insulin$upper_bound
## [1] 316.1896
## 
## $Insulin$outliers
##  [1]   3   9  14  23  27  34  55  56  81 112 124 132 140 154 155 156 184 186 187
## [20] 200 211 221 222 229 232 236 238 248 249 259 263 269 275 281 284 287 297 320
## [39] 361 371 376 379 393 400 410 416 437 438 441 475 481 487 496 506 520 537 575
## [58] 580 585 601 613 623 635 636 646 656 659 676 685 692 696 708 711 713 716 738
## [77] 750 754 758 760 768
## 
## 
## $BMI
## $BMI$lower_bound
## [1] 18.4
## 
## $BMI$upper_bound
## [1] 46.3
## 
## $BMI$outliers
##  [1]  58  85  93 100 121 126 155 156 178 194 240 248 304 336 379 419 439 446 488
## [20] 527 591 674 682 747 748
## 
## 
## $DiabetesPedigreeFunction
## $DiabetesPedigreeFunction$lower_bound
## [1] -0.13
## 
## $DiabetesPedigreeFunction$upper_bound
## [1] 0.875
## 
## $DiabetesPedigreeFunction$outliers
##  [1]   5  13  36  40  46  58  59 101 105 132 148 153 180 188 199 219 221 229 244
## [20] 246 255 260 266 267 268 271 280 293 309 315 331 337 339 357 371 372 380 383
## [39] 384 390 396 397 409 413 417 435 444 446 459 465 488 494 509 533 535 540 547
## [58] 575 589 591 594 607 614 619 622 623 658 660 662 674 693 700 703 717 745 748
## [77] 751 756
## 
## 
## $Age
## $Age$lower_bound
## [1] 8
## 
## $Age$upper_bound
## [1] 50
## 
## $Age$outliers
##  [1]   9  10  13  14  15  25  29  31  40  44  54  68  94 116 124 130 141 149 187
## [20] 205 207 208 213 222 224 237 260 264 275 279 285 286 295 300 320 345 362 363
## [39] 364 376 389 402 454 457 459 460 476 480 488 490 496 499 510 513 517 518 520
## [58] 538 543 547 553 558 561 580 583 585 589 604 659 661 667 675 677 685 703 718
## [77] 720 735 758 760 764

Viendo los resultados de la prueba de Hampel, podemos observar que todas las variables seleccionadas tienen datos atípicos, BMI y Skin Thickness son las que menos tienen. Pero debemos tener en cuenta que esta prueba es muy sensible al basarse en la mediana y el MAD, es sensible a la forma de la distribución de los datos. Para variables con una distribución sesgada o con valores “legítimamente altos” (como podría ser el número de embarazos para algunas mujeres), el método puede identificar incorrectamente estos valores como atípicos.

Capping

boxplot(df$Glucose, main = "Boxplot de Glucosa antes del capping", col = "lightblue")

Podemos obeservar hay un valor atípico 0

x <- df$Glucose
qnt <- quantile(x, probs=c(.25, .75), na.rm = T)
caps <- quantile(x, probs=c(.05, .95), na.rm = T)

H <- 1.5 * IQR(x, na.rm = T)
x[x < (qnt[1] - H)] <- caps[1]
x[x > (qnt[2] + H)] <- caps[2]

boxplot(x, col = "salmon")

Con el capping quitamos en valor atípico y queda la grafíca sin estos valores

boxplot(df$Pregnancies, main = "Boxplot de Pregnancies antes del capping", col = "lightblue")

Podemos observar que para pregnancies tenemos 3 valores atípicos

x <- df$Pregnancies
qnt <- quantile(x, probs=c(.25, .75), na.rm = T)
caps <- quantile(x, probs=c(.05, .95), na.rm = T)

H <- 1.5 * IQR(x, na.rm = T)
x[x < (qnt[1] - H)] <- caps[1]
x[x > (qnt[2] + H)] <- caps[2]

boxplot(x, main = "Boxplot de Pregnancies después del capping", col = "salmon")

Con el capping quitamos los valores atípicos del boxplot como se puede ver en la grafica.

boxplot(df$BloodPressure, main = "Boxplot de Blood Pressure antes del capping", col = "lightblue")

En el boxplot podemos observar que la variable blood pressure tiene valores atípicos antes y despues de los valores extremos, 3 abajo del minimo y 4 arriba del maximo

x <- df$BloodPressure
qnt <- quantile(x, probs=c(.25, .75), na.rm = T)
caps <- quantile(x, probs=c(.05, .95), na.rm = T)

H <- 1.5 * IQR(x, na.rm = T)
x[x < (qnt[1] - H)] <- caps[1]
x[x > (qnt[2] + H)] <- caps[2]

boxplot(x, main = "Boxplot de Blood Pressure después del capping", col = "salmon")

Con el capping eliminamos estos valores atípicos, quedando el boxplot sin estos.

boxplot(df$SkinThickness, main = "Boxplot de Skin Thickness antes del capping", col = "lightblue")

para skin thickness podemos observar en el boxplot que tiene solo un valor atípico

x <- df$SkinThickness
qnt <- quantile(x, probs=c(.25, .75), na.rm = T)
caps <- quantile(x, probs=c(.05, .95), na.rm = T)

H <- 1.5 * IQR(x, na.rm = T)
x[x < (qnt[1] - H)] <- caps[1]
x[x > (qnt[2] + H)] <- caps[2]

boxplot(x, main = "Boxplot de Skin Thickness después del capping", col = "salmon")

Con capping podemos quitar ese unico valor atípico, como podemos ver en la gráfica

boxplot(df$Insulin, main = "Boxplot de Insluin antes del capping", col = "lightblue")

Podemos observar que la variable insulin es una de las que más tiene valores atípicos.

x <- df$Insulin
qnt <- quantile(x, probs=c(.25, .75), na.rm = T)
caps <- quantile(x, probs=c(.05, .95), na.rm = T)

H <- 1.5 * IQR(x, na.rm = T)
x[x < (qnt[1] - H)] <- caps[1]
x[x > (qnt[2] + H)] <- caps[2]

boxplot(x, main = "Boxplot de Insulin después del capping", col = "salmon")

Con capping podemos ver que todos los valores atípicos fueron eliminados de la gráfica, y queda un boxplot más limpio para leer

boxplot(df$BMI, main = "Boxplot de BMI antes del capping", col = "lightblue")

En el boxplot podemos ver que la variable BMI tiene valores atípicos debajo del mínimo y arriba del máximo

x <- df$BMI
qnt <- quantile(x, probs=c(.25, .75), na.rm = T)
caps <- quantile(x, probs=c(.05, .95), na.rm = T)

H <- 1.5 * IQR(x, na.rm = T)
x[x < (qnt[1] - H)] <- caps[1]
x[x > (qnt[2] + H)] <- caps[2]

boxplot(x, main = "Boxplot de BMI después del capping", col = "salmon")

Con el capping quitamos estos valores atípicos y queda un boxplot más facil de analizar

boxplot(df$DiabetesPedigreeFunction, main = "Boxplot de Diabetes Pedigree Function antes del capping", col = "lightblue")

Podemos ver en el boxplot que la variable diabetes pedigree function igual que insulin es la que tiene la gran mayor cantidad de números atípicos.

x <- df$DiabetesPedigreeFunction
qnt <- quantile(x, probs=c(.25, .75), na.rm = T)
caps <- quantile(x, probs=c(.05, .95), na.rm = T)

H <- 1.5 * IQR(x, na.rm = T)
x[x < (qnt[1] - H)] <- caps[1]
x[x > (qnt[2] + H)] <- caps[2]

boxplot(x, main = "Boxplot de Diabetes Pedigree Function después del capping", col = "salmon")

Con el metodo de capping podemos ver que la gráfica queda mejor para el analisis del boxplot

boxplot(df$Age, main = "Boxplot de Age antes del capping", col = "lightblue")

La variable age es otro que tiene una cantidad considerable de datos atípicos según el boxplot.

x <- df$Age
qnt <- quantile(x, probs=c(.25, .75), na.rm = T)
caps <- quantile(x, probs=c(.05, .95), na.rm = T)

H <- 1.5 * IQR(x, na.rm = T)
x[x < (qnt[1] - H)] <- caps[1]
x[x > (qnt[2] + H)] <- caps[2]

boxplot(x, main = "Boxplot de Age después del capping", col = "salmon")

Podemos observar con el método de capping que la gráfica queda sin los valores atípicos.

Visualizaciones Interactivas

# Cargar las librerías necesarias
library(plotly)

## Warning: package 'plotly' was built under R version 4.1.3

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(rjson)
library(readr)
library(dplyr)
library(ggplot2)
library(tidyr)
library(quantmod)

## Warning: package 'quantmod' was built under R version 4.1.3

## Loading required package: xts

## Warning: package 'xts' was built under R version 4.1.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.1.3

## 
## Attaching package: 'zoo'

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

## 
## ################################### WARNING ###################################
## # We noticed you have dplyr installed. The dplyr lag() function breaks how    #
## # base R's lag() function is supposed to work, which breaks lag(my_xts).      #
## #                                                                             #
## # Calls to lag(my_xts) that you enter or source() into this session won't     #
## # work correctly.                                                             #
## #                                                                             #
## # All package code is unaffected because it is protected by the R namespace   #
## # mechanism.                                                                  #
## #                                                                             #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## # You can use stats::lag() to make sure you're not using dplyr::lag(), or you #
## # can add conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop   #
## # dplyr from breaking base R's lag() function.                                #
## ################################### WARNING ###################################

## 
## Attaching package: 'xts'

## The following objects are masked from 'package:dplyr':
## 
##     first, last

## Loading required package: TTR

## Warning: package 'TTR' was built under R version 4.1.3

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(gapminder)

## Warning: package 'gapminder' was built under R version 4.1.3

df_stocks <- read.csv("C:/Users/pc/Desktop/Universidad/Semestre_VI/VIZ/Tarea2/Errores_comunes/stocks.csv")

p <- plot_ly(data = df_stocks, x = ~date, y = ~AMZN, type = 'scatter', mode = 'lines',
             line = list(color = 'salmon')) %>%
  layout(title = 'Amazon Stock Prices Over Time')
p

Esta gráfica muestra la evolución del precio de las acciones de Amazon a lo largo del tiempo. La gráfica ilustra cómo el precio de las acciones de Amazon ha aumentado significativamente desde sus inicios, lo que refleja el crecimiento y la expansión exitosa de la empresa en el comercio electrónico y otros sectores.

fig <- plot_ly(data = df_stocks) %>%
  add_trace(x = ~date, y = ~AAPL, mode = 'lines', name = 'Apple', type = 'scatter') %>%
  add_trace(x = ~date, y = ~AMZN, mode = 'lines+markers', name = 'Amazon', type = 'scatter') %>%
  add_trace(x = ~date, y = ~FB, mode = 'lines+markers', name = 'Facebook', type = 'scatter') %>%
  add_trace(x = ~date, y = ~GOOG, mode = 'lines+markers', name = 'Google', type = 'scatter',
            line = list(color = 'red', width = 2, dash = 'dashdot')) %>%
  layout(title = 'Stock Price data',
         xaxis = list(title = 'Date'),
         yaxis = list(title = 'Price'))

fig

Esta gráfica muestra la evolución de los precios de las acciones de cuatro grandes empresas tecnológicas: Apple, Amazon, Facebook y Google, a lo largo de un período de tiempo específico.

Se pueden observar varias tendencias interesantes:

Apple (línea azul) y Amazon (línea naranja) muestran una tendencia general al alza, con períodos de volatilidad. Apple parece haber experimentado un crecimiento más constante, mientras que Amazon exhibe mayores fluctuaciones.

Facebook (línea verde) tuvo un período de crecimiento significativo seguido por una caída sustancial, lo que sugiere posibles desafíos o eventos que afectaron su desempeño.

Google (línea roja) muestra una tendencia ascendente más suave, pero con un crecimiento notable en los últimos períodos de tiempo.

df_gapminder <- read.csv("https://raw.githubusercontent.com/plotly/datasets/master/gapminderDataFiveYear.csv")

df_pakistan <- subset(df_gapminder, country == "Pakistan")

# Crea un gráfico de barras usando plot_ly
plot_ly(df_pakistan, x = ~year, y = ~pop, type = 'bar', name = 'Population of Pakistan Over the Years', marker = list(color = '#0099ff')) %>%
  layout(title = 'Population of Pakistan Over the Years',
         xaxis = list(title = 'Year'),
         yaxis = list(title = 'Population'),
         showlegend = FALSE)

Este gráfico muestra el crecimiento de la población de Pakistán a lo largo de los años, desde 1950 hasta 2000. Las barras representan el tamaño de la población, que ha aumentado constantemente con el tiempo. El gráfico ilustra el rápido crecimiento de la población de Pakistán, que pasó de alrededor de 40 millones en 1950 a más de 160 millones en el año 2000.

df_tips <- read_csv("C:/Users/pc/Desktop/Universidad/Semestre_VI/VIZ/Tarea2/Errores_comunes/tips.csv")

## New names:
## Rows: 244 Columns: 8
## -- Column specification
## -------------------------------------------------------- Delimiter: "," chr
## (4): sex, smoker, day, time dbl (4): ...1, total_bill, tip, size
## i Use `spec()` to retrieve the full column specification for this data. i
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## * `` -> `...1`

fig <- plot_ly(data = df_tips, x = ~day, y = ~tip, color = ~sex, type = 'bar') %>%
  layout(title = 'Tips by Sex on Each Day',
         xaxis = list(title = 'Day of the Week'),
         yaxis = list(title = 'Tip Amount'))

# Muestra el gráfico
fig

Este gráfico de barras muestra las propinas recibidas por trabajadores de servicios masculinos y femeninos en diferentes días de la semana. Las barras representan los montos de las propinas, con barras azules que indican trabajadores varones y barras verdes para trabajadoras. Parece que los trabajadores tienden a recibir propinas más altas la mayoría de los días, particularmente los jueves, mientras que las trabajadoras reciben propinas más bajas en comparación.

fig <- plot_ly(data = df_tips, x = ~sex, y = ~total_bill, color = ~smoker, type = 'bar', barmode = 'group') %>%
  layout(title = 'Total Bill by Sex and Smoker Status',
         xaxis = list(title = 'Sex'),
         yaxis = list(title = 'Total Bill'))

# Muestra el gráfico
fig

Este gráfico compara los montos totales de las facturas de fumadores y no fumadores, clasificados por género. Las barras representan la factura total, las barras verdes indican no fumadores y las barras azules representan fumadores. Tanto para hombres como para mujeres, la factura total es mayor para los no fumadores que para los fumadores. Además, la factura total de los hombres no fumadores es mayor que la de las mujeres no fumadoras.

df_asia <- subset(df_gapminder, continent == "Asia" & year == 2007 & pop > 2e6)

# Crea el gráfico de barras
plot <- plot_ly(data = df_asia, y = ~pop, x = ~country, text = ~pop, color = ~country, type = 'bar') %>%
  layout(title = 'Population of Asian Countries in 2007 (Population > 2 million)',
         xaxis = list(title = 'Country', tickangle = -45),
         yaxis = list(title = 'Population'),
         uniformtext = list(minsize = 8))

# Muestra el gráfico
plot

Esta gráfica de barras muestra la población de varios países asiáticos en 2007, para aquellos con una población superior a 2 millones de personas. Se puede observar que China e India son los países más poblados, seguidos por Indonesia, Pakistán y Bangladesh.

data("iris")
df_iris <- iris

fig <- plot_ly(data = df_iris, x = ~Sepal.Width, y = ~Sepal.Length, color = ~Species,
                size = ~Petal.Length, text = ~paste("Petal Width: ", Petal.Width, "cm"),
                type = 'scatter', mode = 'markers') %>%
  layout(title = 'Sepal Width vs Sepal Length by Species',
         xaxis = list(title = 'Sepal Width (cm)'),
         yaxis = list(title = 'Sepal Length (cm)'))

fig

Este gráfico de dispersión muestra la relación entre el ancho y la longitud de los sépalos de diferentes especies de flores iris. Los puntos se agrupan por especie y se puede apreciar ciertos patrones y separaciones entre las diferentes especies según estas dos características.

fig <- plot_ly(data = df_iris, x = ~Sepal.Width, y = ~Sepal.Length,
                text = ~Species, type = 'scatter', mode = 'markers',
                marker = list(color = ~Sepal.Width, size = 10, line = list(width = 2), showscale = TRUE)) %>%
  layout(title = 'Sepal Width vs Sepal Length',
         xaxis = list(title = 'Sepal Width (cm)'),
         yaxis = list(title = 'Sepal Length (cm)'))

fig

Este es otro gráfico de dispersión que muestra la misma relación entre ancho y longitud de sépalos, pero sin distinguir por especie. Se observa una forma general de nube de puntos que sugiere una correlación positiva entre ambas variables.

set.seed(123) 
random_data <- data.frame(
  x = rnorm(100000),
  y = rnorm(100000),
  color = rnorm(100000)
)

fig <- plot_ly(data = random_data, x = ~x, y = ~y, type = 'scattergl', mode = 'markers',
                marker = list(color = ~color, colorscale = 'Viridis', line = list(width = 1))) %>%
  layout(title = 'Random Scatter Plot',
         xaxis = list(title = 'X'),
         yaxis = list(title = 'Y'))

fig

Esta gráfica muestra una nube de puntos aleatorios dispersos en un plano cartesiano. No parece haber un patrón específico, por lo que podría representar datos o muestras aleatorias.

df_samer <- subset(df_gapminder, year == 2007 & continent == 'Asia')

plot_ly(df_samer, labels = ~country, values = ~pop, type = 'pie', 
        text = ~paste(country, "<br>", pop), 
        colors = c('blue', 'red', 'green', 'yellow', 'orange', 'purple')) %>%
  layout(title = 'Population of Asian continent')

Esta gráfica circular o pie chart muestra la población relativa de los principales países asiáticos. China e India son los más poblados, seguidos por Indonesia, Pakistán, Bangladesh y otros países de la región.

colors <- c('blue', 'green', 'black', 'purple', 'red', 'brown')

data <- list(
  labels = c('Water', 'Grass', 'Normal', 'Psychic', 'Fire', 'Ground'),
  values = c(110, 90, 80, 80, 70, 60)
)

fig <- plot_ly(labels = data$labels, values = data$values, type = 'pie') %>%
  layout(
    hoverinfo = 'label+percent',
    textinfo = 'label+percent',
    pull = c(0.1, 0, 0.2, 0, 0, 0),
    marker = list(colors = colors, line = list(color = '#FFFFFF', width = 2))
  )

fig

Este gráfico de barras circular muestra la distribución de diferentes tipos de elementos. El agua representa el 22,4% de la distribución, mientras que la hierba y el tipo normal juntos constituyen alrededor del 35% de los elementos. El psíquico y el fuego tienen porcentajes más pequeños del 14,3% y el 12,2%, respectivamente.

set.seed(123) 

dice_1 <- sample(1:6, 5000, replace = TRUE)
dice_2 <- sample(1:6, 5000, replace = TRUE)
dice_sum <- dice_1 + dice_2

# Histograma
histogram <- plot_ly(x = dice_sum, type = 'histogram', nbinsx = 11) %>%
  layout(
    xaxis = list(title = 'Dice Roll'),
    yaxis = list(title = 'Frequency'),
    title = '5000 Dice Roll Histogram',
    bargap = 0.2,
    showlegend = FALSE
  )

violin <- plot_ly(x = rep("Dice Sum", length(dice_sum)), y = dice_sum, type = 'violin') %>%
  layout(
    xaxis = list(title = ''),
    yaxis = list(title = 'Dice Roll'),
    showlegend = FALSE
  )

subplot <- subplot(
  histogram, violin,
  nrows = 2,
  margin = 0.05
)

subplot

Esta es un histograma que muestra los resultados de 5000 lanzamientos de dados. La distribución sigue una forma de campana aproximada, lo que indica que los resultados se distribuyen de manera cercana a la distribución normal esperada. Los valores centrales de 6, 7 y 8 tienen las frecuencias más altas, mientras que los valores extremos de 2 y 12 tienen las frecuencias más bajas.

fig <- plot_ly(data = df_tips, x = ~total_bill, color = ~sex, type = 'histogram') %>%
  layout(title = 'Total Bills by Sex',
         barmode = 'stack',
         barnorm = '') 
fig

Este gráfico de barras agrupado muestra el recuento total de facturas para hombres y mujeres. El eje horizontal representa los montos de las facturas, y el eje vertical muestra las frecuencias. La distribución parece sesgada hacia la izquierda, con más facturas de montos más bajos tanto para hombres como para mujeres.

fig <- plot_ly(df_tips, x = ~sex, y = ~tip, type = 'box',
               boxpoints = 'all') %>%
  layout(title = 'Tip Distribution by Sex',
         xaxis = list(title = 'Sex'),
         yaxis = list(title = 'Tip'))

fig

Este diagrama de caja muestra la distribución de propinas para hombres y mujeres. Los hombres tienden a dar propinas ligeramente más altas en promedio, con una mediana y un rango intercuartílico más altos. Sin embargo, ambos grupos tienen valores atípicos tanto altos como bajos.

ggplot(df_tips, aes(x = day, y = tip, color = sex)) +
  geom_boxplot() +
  labs(x = "Day", y = "Tip", color = "Sex") +
  theme_minimal()

Este gráfico de cajas y bigotes muestra la distribución de propinas por día de la semana y género. Tanto para hombres como para mujeres, las propinas tienden a ser más altas los fines de semana (sábado y domingo) en comparación con los días de semana. Se observan varios valores atípicos, particularmente en los días de semana.

fig <- plot_ly(df_tips, x = ~sex, y = ~tip, type = 'box', marker = list(color = 'blue'), boxmean = 'sd')

fig

Este gráfico muestra una representación de árbol de decisión para clasificar el sexo (femenino o masculino) en función de ciertas características o variables. Los rectángulos representan nodos de decisión y las ramas que salen muestran los posibles resultados.

fig <- plot_ly() %>%
  add_boxplot(data = df_stocks, y = ~GOOG, name = 'Google', boxpoints = 'all', fillcolor = 'blue', jitter = 0.5, whiskerwidth = 0.2) %>%
  add_boxplot(data = df_stocks, y = ~AAPL, name = 'Apple', boxpoints = 'all', fillcolor = 'red', jitter = 0.5, whiskerwidth = 0.2)

# Mostrar el diagrama
fig

Es un diagrama de caja que compara la distribución de un conjunto de datos (posiblemente ingresos o valoraciones) para Google y Apple. Se puede ver que la mediana de Google es menor que la de Apple, pero Google tiene más valores atípicos superiores.

fig <- plot_ly(df_tips, y = ~total_bill, type = 'violin', box = list(visible = TRUE), points = 'all')

fig

Es un gráfico de dispersión que muestra una gran cantidad de puntos concentrados formando una distribución similar a una curva normal o gaussiana. Parece representar la distribución de alguna variable o métrica.

fig <- subplot(
  plot_ly(df_tips[df_tips$sex == "Male",], y = ~tip, x = ~smoker, type = 'violin', color = I("blue"), box = list(visible = TRUE), points = 'all', name = "Masculino"),
  plot_ly(df_tips[df_tips$sex == "Female",], y = ~tip, x = ~smoker, type = 'violin', color = I("red"), box = list(visible = TRUE), points = 'all', name = "Femenino"),
  nrows = 2
)

fig

Estas son distribuciones de densidad de kernels que muestran cómo se distribuyen ciertas características o variables para diferentes grupos. Las curvas azules y rojas representan dos grupos distintos, posiblemente sexo masculino y femenino. Los paneles muestran la distribución condicionada a una respuesta “Sí” o “No”.

df_smokers <- df_tips[df_tips$smoker == "Yes", ]

fig <- plot_ly(df_smokers, x = ~day, y = ~total_bill, type = 'violin', legendgroup = 'Yes', scalegroup = 'Yes', name = 'Yes', side = 'negative', line = list(color = 'blue'))

fig

Este gráfico muestra cuatro distribuciones de una variable (posiblemente una métrica o conteo) a lo largo de diferentes días de la semana (viernes, sábado, domingo, jueves). Las distribuciones tienen forma de curva normal con un pico pronunciado, lo que sugiere que la variable tiene valores concentrados alrededor de una media específica cada día.

df_non_smokers <- df_tips[df_tips$smoker == "No", ]

fig <- fig %>%
  add_trace(data = df_non_smokers, x = ~day, y = ~total_bill, type = 'violin', legendgroup = 'No', scalegroup = 'Yes', name = 'No', side = 'positive', line = list(color = 'red'))

fig

Este gráfico muestra la distribución del total de facturas en diferentes días de la semana, divididas por una categoría binaria (Sí/No). El ancho de cada ‘violín’ muestra la densidad de puntos de datos en diferentes importes de billetes. Parece que la mediana de la factura total es generalmente más alta el domingo para ambas categorías, con una distribución más amplia el sábado, lo que indica una mayor variabilidad en los montos totales de la factura ese día.

flights <- read.csv("C:/Users/pc/Desktop/Universidad/Semestre_VI/VIZ/Tarea2/Errores_comunes/flights.csv")

fig <- plot_ly(flights, x = ~year, y = ~month, z = ~passengers, type = "heatmap", colors = colorRamp(c("blue", "yellow", "red")))

# Mostrar el gráfico
fig

El mapa de calor muestra el número de pasajeros por mes y año, con colores que indican el volumen de pasajeros. Los colores más cálidos representan un mayor número de pasajeros, mientras que los colores más fríos representan un número menor. Este mapa de calor sugiere una tendencia de aumento del número de pasajeros a lo largo de los años y posiblemente tendencias estacionales con ciertos meses que muestran mayores volúmenes de pasajeros de manera constante.

fig <- plot_ly(flights, x = ~year, y = ~month, z = ~passengers, color = ~year,
               type = "scatter3d", mode = "markers", opacity = 0.7)

# Definir el tamaño del gráfico
fig <- layout(fig, width = 800, height = 400)

# Mostrar el gráfico
fig

Este gráfico es una representación 3D del número de pasajeros por mes y año. Cada punto representa el número de pasajeros en un mes y año específico, y el degradado de color muestra el paso del tiempo de 1950 a 1960. La tendencia sugiere un aumento general de pasajeros a lo largo de los años, y también se pueden observar patrones estacionales.

fig <- plot_ly(flights, x = ~year, y = ~month, z = ~passengers, color = ~year, type = "scatter3d", mode = "lines")

# Mostrar el gráfico
fig

El gráfico es similar al gráfico de dispersión 3D pero muestra los datos como una serie de líneas que conectan los puntos del mismo año. Esto ayuda a rastrear el patrón del número de pasajeros a lo largo de cada año, destacando más claramente las tendencias estacionales, así como la tendencia al alza en el número de pasajeros a lo largo del tiempo.

fig <- plot_ly(data = flights, color = ~month) %>%
  add_trace(
    type = 'splom',
    dimensions = list(
      list(label = 'year', values = ~year),
      list(label = 'month', values = ~month),
      list(label = 'passengers', values = ~passengers)
    )
  )

fig

Año vs. Año (panel superior izquierdo): Dado que las variables son las mismas en ambos ejes, los puntos forman una línea perfecta ascendente, indicando simplemente la progresión temporal de un año a otro.

Mes vs. Año (paneles superiores e inferiores del medio): Estos gráficos muestran la distribución de los vuelos a lo largo de cada mes en diferentes años. No parece haber una tendencia o patrón claro en la variación anual, pero sí se observa la recurrencia estacional por el color de los puntos.

Pasajeros vs. Año (panel superior derecho): Existe una tendencia creciente en el número de pasajeros a lo largo de los años, lo que indica un crecimiento en la industria de vuelos o un aumento en la demanda de viajes aéreos durante este período.

Mes vs. Mes (diagonal central): Al igual que con el año vs. año, estos puntos forman una línea diagonal perfecta, ya que es la misma variable en ambos ejes.

Pasajeros vs. Mes (paneles inferiores a la derecha y superior a la izquierda): Estos paneles muestran variaciones estacionales en el número de pasajeros, con algunos meses teniendo consistentemente más pasajeros que otros, como se puede ver por la concentración de colores. Por ejemplo, hay meses con picos de pasajeros (posiblemente los meses de vacaciones).

Pasajeros vs. Pasajeros (panel inferior derecho): De nuevo, como es la misma variable en ambos ejes, los datos forman una línea perfectamente ascendente.

library(countrycode)
df <- gapminder %>% filter(year == 2007)

df$iso_alpha <- countrycode(df$country, origin = "country.name", destination = "iso3c")


fig <- plot_geo(df, locations = ~iso_alpha,
                color = ~continent,
                hoverinfo = "text",
                text = ~paste("Country: ", country, "<br>Population: ", pop),
                size = ~pop,
                marker = list(line = list(width = 0.5, color = "black")),
                projection = "orthographic")

fig

## No scattergeo mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Este mapa muestra burbujas en varios lugares del mundo, con diferentes colores que representan diferentes continentes. El tamaño de las burbujas probablemente represente una cantidad, como la población, el PIB o alguna otra estadística. La distribución de los tamaños de las burbujas podría indicar que Asia y Europa tienen valores más altos de la cantidad medida, aunque sin una leyenda para los tamaños, es difícil interpretar los valores exactos.

url <- 'https://raw.githubusercontent.com/plotly/datasets/master/geojson-counties-fips.json'
counties <- rjson::fromJSON(file=url)
url2<- "https://raw.githubusercontent.com/plotly/datasets/master/fips-unemp-16.csv"
df <- read.csv(url2, colClasses=c(fips="character"))

g <- list(
  scope = 'usa',
  projection = list(type = 'albers usa'),
  showlakes = TRUE,
  lakecolor = toRGB('white')
)
fig <- plot_ly()
fig <- fig %>% add_trace(
    type="choropleth",
    geojson=counties,
    locations=df$fips,
    z=df$unemp,
    colorscale="Viridis",
    zmin=0,
    zmax=12,
    marker=list(line=list(
      width=0)
    )
  )
fig <- fig %>% colorbar(title = "Unemployment Rate (%)")
fig <- fig %>% layout(
    title = "US Unemployment"
)

fig <- fig %>% layout(
    geo = g
  )

fig

El mapa de coropletas muestra los Estados Unidos con diferentes áreas coloreadas para representar la tasa de desempleo. La escala de colores indica que el verde más oscuro representa tasas de desempleo más bajas, mientras que el amarillo indica tasas más altas. Este mapa permite una comparación rápida del desempleo en diferentes regiones del país.

wind <- read_csv("C:/Users/pc/Desktop/Universidad/Semestre_VI/VIZ/Tarea2/Errores_comunes/wind.csv")

## New names:
## Rows: 128 Columns: 4
## -- Column specification
## -------------------------------------------------------- Delimiter: "," chr
## (2): direction, strength dbl (2): ...1, frequency
## i Use `spec()` to retrieve the full column specification for this data. i
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## * `` -> `...1`

fig <- plot_ly(data = wind, type = 'scatterpolar', mode = 'markers',
               r = ~frequency, theta = ~direction,
               color = ~strength, size = ~frequency, symbol = ~strength)

fig <- fig %>% layout(
  polar = list(
    radialaxis = list(visible = TRUE, range = c(min(wind$frequency), max(wind$frequency))),
    angularaxis = list(direction = "clockwise")
  )
)

fig

Este gráfico de rosa de los vientos proporciona una descripción general de la dirección del viento y las frecuencias de velocidad. La trama se divide en direcciones cardinales e intercardinales, con diferentes símbolos y colores que representan categorías de velocidad. Vemos que la mayoría de las velocidades del viento caen en las categorías inferiores, con relativamente menos casos de velocidades más altas.

colors <- c('blue', 'cyan', 'green', 'yellow', 'orange', 'red', 'darkred')

p <- plot_ly(wind, type = 'scatterpolar', mode = 'lines',
             r = ~frequency, theta = ~direction,
             line = list(color = ~as.numeric(strength), width = 2),
             color = ~strength, colors = colors,
             opacity = 0.7)

p <- p %>% layout(
  polar = list(
    radialaxis = list(
      visible = TRUE,
      range = c(0, max(wind$frequency))
    ),
    bgcolor = "black" 
  ),
  paper_bgcolor = "black", 
  plot_bgcolor = "black", 
  showlegend = TRUE,
  width = 800,
  height = 400
)
p

El gráfico de radar es una forma de mostrar múltiples puntos de datos y la variación entre ellos. Aquí, parece representar los mismos datos que la rosa de los vientos: dirección del viento y frecuencia de la velocidad. Cada línea conecta puntos de la misma categoría de velocidad en diferentes direcciones, lo que facilita comparar la distribución de las velocidades del viento en cada dirección.

df_exp <- read_csv("C:/Users/pc/Desktop/Universidad/Semestre_VI/VIZ/Tarea2/Errores_comunes/experiment.csv")

## New names:
## Rows: 100 Columns: 6
## -- Column specification
## -------------------------------------------------------- Delimiter: "," chr
## (2): gender, group dbl (4): ...1, experiment_1, experiment_2, experiment_3
## i Use `spec()` to retrieve the full column specification for this data. i
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## * `` -> `...1`

df_exp <- df_exp %>%
  mutate(gender = as.factor(gender))

colors <- c('female' = '#FF69B4', 'male' = '#1E90FF') # Rosa para 'female', Azul para 'male'

fig <- plot_ly(data = df_exp, type = 'scatterternary',
               mode = 'markers',
               a = ~experiment_1,
               b = ~experiment_2,
               c = ~experiment_3,
               hoverinfo = 'text',
               text = ~paste("Group:", group, "<br>Gender:", gender),
               marker = list(size = 10,
                             color = ~I(colors[as.character(gender)]))) 

fig

Un gráfico ternario muestra las proporciones de tres variables que suman una constante. En este gráfico, el Componente A, el Componente B y el Componente C deben sumar 1 para cualquier punto dado. La agrupación de puntos alrededor del centro sugiere que los valores de estos componentes están relativamente equilibrados para la mayoría de las observaciones.

ggplot(df_tips, aes(x = total_bill, y = tip, color = smoker)) +
  geom_point() +
  facet_wrap(~sex) +
  labs(title = "Scatter plot of Total Bill vs. Tip", x = "Total Bill", y = "Tip") +
  theme_minimal()

Este diagrama de dispersión compara la factura total con la propina dada, dividida por género (mujer a la izquierda y hombre a la derecha), con codificación de colores para fumadores y no fumadores. Esto puede proporcionar información sobre el comportamiento de dejar propinas por género y preferencia de fumar, aunque no muestra un patrón claro de correlación a primera vista.

df_att <-read_csv("C:/Users/pc/Desktop/Universidad/Semestre_VI/VIZ/Tarea2/Errores_comunes/att.csv")

## New names:
## Rows: 60 Columns: 6
## -- Column specification
## -------------------------------------------------------- Delimiter: "," chr
## (1): attention dbl (5): ...1, Unnamed: 0, subject, solutions, score
## i Use `spec()` to retrieve the full column specification for this data. i
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## * `` -> `...1`

p <- ggplot(df_att, aes(x = solutions, y = score)) + 
  geom_line() + 
  facet_wrap(~subject, ncol = 5) + 
  labs(title = 'Scores Based on Attention', x = 'Solutions', y = 'Score')

print(p)

Se utilizan múltiplos pequeños para comparar muchos gráficos que muestran el mismo tipo de datos. Cada minigráfico tiene un título con un número, que sugiere un individuo o una categoría, y traza la ‘Puntuación’ frente a las ‘Soluciones’. Este diseño permite una comparación rápida entre diferentes temas para ver cómo varían sus puntuaciones con la cantidad de soluciones.

df_cnt <-read_csv("C:/Users/pc/Desktop/Universidad/Semestre_VI/VIZ/Tarea2/Errores_comunes/cnt.csv")

## New names:
## Rows: 1704 Columns: 9
## -- Column specification
## -------------------------------------------------------- Delimiter: "," chr
## (3): country, continent, iso_alpha dbl (6): ...1, year, lifeExp, pop,
## gdpPercap, iso_num
## i Use `spec()` to retrieve the full column specification for this data. i
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## * `` -> `...1`

fig <- df_cnt %>%
  plot_ly(
    x = ~gdpPercap, 
    y = ~lifeExp, 
    size = ~pop, 
    color = ~continent, 
    frame = ~year, 
    text = ~country, 
    hoverinfo = "text",
    type = 'scatter',
    mode = 'markers'
  )
fig <- fig %>% layout(
    xaxis = list(
      type = "log"
    )
  )

fig

El gráfico de burbujas interactivo rastrea la esperanza de vida en comparación con el PIB per cápita, donde el tamaño de la burbuja posiblemente represente el tamaño de la población y los colores representen diferentes continentes. Hay una clara tendencia positiva, que indica que a medida que aumenta el PIB per cápita, la esperanza de vida también tiende a aumentar. El control deslizante en la parte inferior sugiere que puede pasar los años para ver cómo estas variables han cambiado con el tiempo.

fig <- df_cnt %>%
  plot_ly(
    x = ~continent,        
    y = ~pop,             
    color = ~continent,    
    frame = ~year,         
    text = ~country,       
    hoverinfo = "text",
    type = 'bar',          
    marker = list(size = ~pop)
  )

fig

Este gráfico de barras animado muestra la población (población) por continente a lo largo del tiempo. El botón “Reproducir” en la parte inferior sugiere que puedes observar cómo cambian las poblaciones a lo largo de los años. En 1952, podemos ver que Asia tenía una población significativamente mayor que otros continentes. Este tipo de gráfico sería útil para observar las tendencias de crecimiento de la población.

Visualizaciones Geo Espaciales

library(sf) # para trabajar con datos espaciales

## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE

library(ggplot2) 
library(dplyr)
library(plotly)
library(leaflet)

## 
## Attaching package: 'leaflet'

## The following object is masked from 'package:xts':
## 
##     addLegend

library(leaflet.extras)

datos_sf <- st_read("C:/Users/pc/Desktop/factoriesindustrial-sites.geojson")

## Reading layer `factoriesindustrial-sites' from data source 
##   `C:\Users\pc\Desktop\factoriesindustrial-sites.geojson' using driver `GeoJSON'
## Simple feature collection with 732 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 2.730691 ymin: 6.382559 xmax: 4.190527 ymax: 6.687934
## Geodetic CRS:  WGS 84

datos_sf <- datos_sf %>%
  mutate(lon = st_coordinates(geometry)[,1], 
         lat = st_coordinates(geometry)[,2])
names(datos_sf)

##  [1] "id"         "state_code" "source"     "name"       "status"    
##  [6] "ward_code"  "global_id"  "geometry"   "lon"        "lat"

plot_ly(datos_sf, x = ~lon, y = ~lat, type = "scatter", mode = "markers", 
        marker = list(color = "blue")) %>%
  layout(title = "A Scatter Plot of the Latitude and Longitude of Factories in Lagos",
         xaxis = list(title = "Longitude", gridcolor = "lightgrey", showgrid = TRUE),
         yaxis = list(title = "Latitude", gridcolor = "lightgrey", showgrid = TRUE),
         plot_bgcolor = 'rgba(0,0,0,0)') 

Este gráfico muestra la latitud y longitud de las fábricas en Lagos. En el mismo gráfico se superponen dos trazas, que probablemente representan diferentes tipos de fábricas o quizás diferentes épocas. Los grupos pueden indicar áreas industriales dentro de la ciudad.

mark <- leaflet() %>%
  addTiles(urlTemplate = "https://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png") %>%
  setView(lng = 3.3792, lat = 6.5244, zoom = 12)

for(i in 1:nrow(datos_sf)) {
  mark <- mark %>%
    addMarkers(lng = datos_sf$lon[i], lat = datos_sf$lat[i], popup = datos_sf$name[i])
}

mark

Aquí, las ubicaciones individuales de las fábricas en Lagos están marcadas en un mapa con íconos puntuales, lo que permite una visualización geográfica de la distribución de las fábricas. El mapa está ampliado y muestra una alta densidad de fábricas en ciertas áreas, lo que sugiere que se trata de zonas industriales clave.

bubble <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lng = 3.3792, lat = 6.5244, zoom = 12)

mc <- markerClusterOptions()
bubble <- bubble %>%
  addMarkers(lng = datos_sf$lon, lat = datos_sf$lat, clusterOptions = mc)

bubble

Este mapa representa diferentes áreas de Lagos con burbujas, donde el tamaño de la burbuja probablemente corresponde al número de fábricas en esa área. Las burbujas más grandes representan una mayor concentración de fábricas, lo que proporciona una comparación visual rápida de la actividad industrial en diferentes regiones.

heat <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lng = 3.3792, lat = 6.5244, zoom = 10)

heat <- heat %>%
  addHeatmap(lng = ~lon, lat = ~lat, data = datos_sf, radius = 10, blur = 20, max = 0.05)

heat

Un mapa de calor indica la densidad con color, mostrando cuán concentradas están las fábricas en ciertos lugares. Los colores más cálidos (como el rojo y el amarillo) sugieren mayores concentraciones de fábricas, mientras que los colores más fríos (como el azul) representan densidades más bajas. Este mapa destaca las zonas más industrializadas.

status <- aggregate(cbind(count = status) ~ status, data = datos_sf, FUN = length)
colnames(status)[2] <- "count"

ggplot(status, aes(x = status, y = count, fill = status)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("In Use" = "blue", "Unused Good" = "yellow", "Unused Poor" = "purple")) +
  labs(title = "Status of Factories and Industrial Sites in Lagos", x = "Status", y = "Count")

El gráfico de barras muestra el estado de las fábricas y sitios industriales en Lagos, clasificados como “En uso”, “Bien no utilizado” y “Pobres no utilizados”. La mayoría de las fábricas están en uso, y un pequeño número se clasifica como sin uso, ya sea en buenas o malas condiciones.

circle <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lng = 3.3792, lat = 6.5244, zoom = 10)

color_producer <- function(val) {
  if (val == 'In Use') {
    return('blue')
  } else if (val == 'Unused Good') {
    return('yellow')
  } else {
    return('purple')
  }
}

for(i in 1:nrow(datos_sf)) {
  circle <- circle %>%
    addCircles(lng = datos_sf$lon[i], lat = datos_sf$lat[i], radius = 20,
               color = color_producer(datos_sf$status[i]), fill = TRUE)
}

circle

fig_scatter <- plot_ly(data = datos_sf, type = 'scattermapbox', mode = 'markers',
               lat = ~lat, lon = ~lon,
               color = ~status, colors = c('In Use' = 'blue', 'Unused Good' = 'yellow', 'Unused Poor' = 'purple'),
               text = ~name, hoverinfo = 'text') %>%
  layout(mapbox = list(
    style = 'carto-positron',
    zoom = 10,
    center = list(lat = 6.5244, lon = 3.3792) 
    )
  )

fig_scatter

Al igual que la segunda imagen, este mapa también tiene íconos señalados pero con códigos de colores para diferenciar el estado de las fábricas: “En uso”, “Bueno no usado” o “Pobre no usado”. Esto permite una evaluación del estado operativo de las fábricas de toda la ciudad.

m <- leaflet(datos_sf) %>%
  addProviderTiles("CartoDB.Positron") %>% 
  setView(lng = 3.3792, lat = 6.5244, zoom = 9) %>%
  addHeatmap(lng = ~lon, lat = ~lat, blur = 15, max = 0.05, radius = 15) 

m

Este mapa de densidad, a menudo denominado “mapa de puntos críticos”, utiliza colores para representar la densidad de las fábricas. Las zonas con mayores concentraciones están en rojo, mientras que las que tienen menos están en azul. Esta visualización ayuda a identificar qué áreas tienen la mayor agregación de sitios industriales.

Series de Tiempo

Load data

end_date <- as.Date("2023-12-31")
start_date <- as.Date("2023-01-01")
data <- getSymbols("TGLS", src = "yahoo", from = start_date, to = end_date)
tgls <- as.data.frame(TGLS)

summary(tgls$TGLS.Close)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   29.97   34.24   38.09   39.32   44.17   53.88

frequency(tgls$TGLS.Close)

## [1] 1

Frequency equals 1 so we don’t have to eliminate a tendancy.

Autocorrelation function

plot(tgls$TGLS.Close,
     type = "l",
     main = "TGLS Closing Price",
     ylab = "Price ($)",
     xlab = "Date")

acf(tgls$TGLS.Close, lag.max = 300, ylim = c(-1, 1), main = "ACF of the original sery")

We can assume that the time sery follow a tendancy bceause the values go out of the lines.

tgls_d1 <- diff(tgls$TGLS.Close, lag=1)[-1]
plot(tgls_d1,
     type = "l",
     main = "TGLS - First Differencing",
     ylab = "Price ($)",
     xlab = "Date")

acf(tgls_d1, lag.max = 300, ylim = c(-0.5, 0.5), main = "ACF of 1st Differencing")

We can assume that the time sery doesn’t follow a tendancy because the values mainly stay between the lines.

pacf(tgls$TGLS.Close, lag.max = 300, ylim = c(-0.5, 0.5), main = "Partial ACF of the original sery")

With the partial autocorrelation function, it seems we don’t have a tendancy.

ARIMA models

Best orders

best_order_arima_aic <- function(input, rge_p, rge_d, rge_q) {
  best_aic <- Inf
  best_pdq <- NULL
  fit <- NULL
  for(p in rge_p) {
    for(d in rge_d) {
      for (q in rge_q) {
        tryCatch({
          fit <- arima(input, order=c(p, d, q))
          aic <- AIC(fit)
          if (aic < best_aic) {
            best_aic <- aic
            best_pdq <- c(p, d, q)
          }
        }, error=function(e){})
      }
    }
  }
  return(best_pdq)
}

best_order_arima_bic <- function(input, rge_p, rge_d, rge_q) {
  best_bic <- Inf
  best_pdq <- NULL
  fit <- NULL
  for(p in rge_p) {
    for(d in rge_d) {
      for (q in rge_q) {
        tryCatch({
          fit <- arima(input, order=c(p, d, q))
          bic <- BIC(fit)
          if (bic < best_bic) {
            best_bic <- bic
            best_pdq <- c(p, d, q)
          }
        }, error=function(e){})
      }
    }
  }
  return(best_pdq)
}

best_order_arima_hqic <- function(input, rge_p, rge_d, rge_q) {
  best_hqic <- Inf
  best_pdq <- NULL
  fit <- NULL
  for(p in rge_p) {
    for(d in rge_d) {
      for (q in rge_q) {
        tryCatch({
          fit <- arima(input, order=c(p, d, q))
          hqic <- AIC(fit)
          if (hqic < best_hqic) {
            best_hqic <- hqic
            best_pdq <- c(p, d, q)
          }
        }, error=function(e){})
      }
    }
  }
  return(best_pdq)
}

MAPE <- function(actual, forecast) {
  mean(abs((actual - forecast) / actual)) * 100
}

MAE <- function(actual, forecast) {
  mean(abs(actual - forecast))
}

RMSE <- function(actual, forecast) {
  sqrt(mean((actual - forecast)^2))
}

MSE <- function(actual, forecast) {
  mean((actual - forecast)^2)
}

R2 <- function(actual, forecast) {
  cor(actual, forecast)^2
}

rge_p = c(0:7)
rge_d = c(0:1)
rge_q = c(0:7)

best_order_aic = best_order_arima_aic(tgls$TGLS.Close, rge_p, rge_d, rge_q)
best_order_bic = best_order_arima_bic(tgls$TGLS.Close, rge_p, rge_d, rge_q)
best_order_hqic = best_order_arima_hqic(tgls$TGLS.Close, rge_p, rge_d, rge_q)

Best order AIC:

best_order_aic

## [1] 2 1 2

Best order BIC:

best_order_bic

## [1] 0 1 0

Best order HQIC:

best_order_hqic

## [1] 2 1 2

Rolling Forecast

n = 14

nb_data <- nrow(tgls)
nb_day_predicted = c(14, 21, 28)
train_size = nb_data - max(nb_day_predicted)
train <- tgls$TGLS.Close[1:train_size]

test <- tgls$TGLS.Close[(train_size + 1):(train_size + 14)]
test <- test[complete.cases(test)]
test_dates <- index(tgls)[(train_size + 1):(train_size + 14)]
test_dates <- test_dates[complete.cases(test_dates)]

forecast_aic <- c()
forecast_bic <- c()
forecast_hqic <- c()

for (i in 1:14) {
  fit_aic <- arima(train, order = best_order_aic)
  fit_bic <- arima(train, order = best_order_bic)
  fit_hqic <- arima(train, order = best_order_hqic)
  
  forecast_aic <- c(forecast_aic, predict(fit_aic, n.ahead = 1)$pred)
  forecast_bic <- c(forecast_bic, predict(fit_bic, n.ahead = 1)$pred)
  forecast_hqic <- c(forecast_hqic, predict(fit_hqic, n.ahead = 1)$pred)

  train <- c(train, test[i])
}

n_last = 50
last_dates <- tail(index(tgls), n_last)
last_tgls <- tail(tgls$TGLS.Close, n_last)
last_forecast_aic <- tail(forecast_aic, n_last)
last_forecast_bic <- tail(forecast_bic, n_last)
last_forecast_hqic <- tail(forecast_hqic, n_last)

plot(last_dates, last_tgls, type = "l", col = "black", ylim = range(c(last_tgls, last_forecast_aic, last_forecast_bic, last_forecast_hqic)), 
     main = "TGLS Closing Price", ylab = "Price ($)", xlab = "Date")
lines(test_dates, forecast_aic, col = "red")
lines(test_dates, forecast_bic, col = "green")
lines(test_dates, forecast_hqic, col = "purple")

legend("topleft", legend = c("Actual", "Forecast AIC", "Forecast BIC", "Forecast HQIC"), 
       col = c("black", "red", "green", "purple"), lty = c(1, 1, 1, 1))

mape_aic <- MAPE(test, forecast_aic)
mape_bic <- MAPE(test, forecast_bic)
mape_hqic <- MAPE(test, forecast_hqic)

mae_aic <- MAE(test, forecast_aic)
mae_bic <- MAE(test, forecast_bic)
mae_hqic <- MAE(test, forecast_hqic)

rmse_aic <- RMSE(test, forecast_aic)
rmse_bic <- RMSE(test, forecast_bic)
rmse_hqic <- RMSE(test, forecast_hqic)

mse_aic <- MSE(test, forecast_aic)
mse_bic <- MSE(test, forecast_bic)
mse_hqic <- MSE(test, forecast_hqic)

r2_aic <- R2(test, forecast_aic)
r2_bic <- R2(test, forecast_bic)
r2_hqic <- R2(test, forecast_hqic)

metrics_table <- data.frame(
  Method = c("AIC", "BIC", "HQIC"),
  MAPE = c(mape_aic, mape_bic, mape_hqic),
  MAE = c(mae_aic, mae_bic, mae_hqic),
  RMSE = c(rmse_aic, rmse_bic, rmse_hqic),
  MSE = c(mse_aic, mse_bic, mse_hqic),
  R2 = c(r2_aic, r2_bic, r2_hqic)
)

print(metrics_table)

##   Method     MAPE       MAE      RMSE       MSE        R2
## 1    AIC 1.875319 0.6783717 0.8655771 0.7492237 0.7948870
## 2    BIC 1.965542 0.7114285 0.8975364 0.8055715 0.7771502
## 3   HQIC 1.875319 0.6783717 0.8655771 0.7492237 0.7948870

plot(test, forecast_aic, main="AIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_bic, main="BIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_hqic, main="HQIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

residuals_aic <- test - forecast_aic
residuals_bic <- test - forecast_bic
residuals_hqic <- test - forecast_hqic

plot(residuals_aic, main="AIC Residuals", ylab="Residuals", col="blue")
abline(0,0)

plot(residuals_bic, main="BIC Residuals", ylab="Residuals", col="blue")
abline(0,0)

plot(residuals_hqic, main="HQIC Residuals", ylab="Residuals", col="blue")
abline(0,0)

shapiro.test(residuals_aic)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_aic
## W = 0.96072, p-value = 0.7349

shapiro.test(residuals_bic)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_bic
## W = 0.96032, p-value = 0.7283

shapiro.test(residuals_hqic)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_hqic
## W = 0.96072, p-value = 0.7349

n = 21

nb_data <- nrow(tgls)
nb_day_predicted = c(14, 21, 28)
train_size = nb_data - max(nb_day_predicted)
train <- tgls$TGLS.Close[1:train_size]

test <- tgls$TGLS.Close[(train_size + 1):(train_size + 21)]
test <- test[complete.cases(test)]
test_dates <- index(tgls)[(train_size + 1):(train_size + 21)]
test_dates <- test_dates[complete.cases(test_dates)]

forecast_aic <- c()
forecast_bic <- c()
forecast_hqic <- c()

for (i in 1:21) {
  fit_aic <- arima(train, order = best_order_aic)
  fit_bic <- arima(train, order = best_order_bic)
  fit_hqic <- arima(train, order = best_order_hqic)
  
  forecast_aic <- c(forecast_aic, predict(fit_aic, n.ahead = 1)$pred)
  forecast_bic <- c(forecast_bic, predict(fit_bic, n.ahead = 1)$pred)
  forecast_hqic <- c(forecast_hqic, predict(fit_hqic, n.ahead = 1)$pred)

  train <- c(train, test[i])
}

n_last = 50
last_dates <- tail(index(tgls), n_last)
last_tgls <- tail(tgls$TGLS.Close, n_last)
last_forecast_aic <- tail(forecast_aic, n_last)
last_forecast_bic <- tail(forecast_bic, n_last)
last_forecast_hqic <- tail(forecast_hqic, n_last)

plot(last_dates, last_tgls, type = "l", col = "black", ylim = range(c(last_tgls, last_forecast_aic, last_forecast_bic, last_forecast_hqic)), 
     main = "TGLS Closing Price", ylab = "Price ($)", xlab = "Date")
lines(test_dates, forecast_aic, col = "red")
lines(test_dates, forecast_bic, col = "green")
lines(test_dates, forecast_hqic, col = "purple")

legend("topleft", legend = c("Actual", "Forecast AIC", "Forecast BIC", "Forecast HQIC"), 
       col = c("black", "red", "green", "purple"), lty = c(1, 1, 1, 1))

mape_aic <- MAPE(test, forecast_aic)
mape_bic <- MAPE(test, forecast_bic)
mape_hqic <- MAPE(test, forecast_hqic)

mae_aic <- MAE(test, forecast_aic)
mae_bic <- MAE(test, forecast_bic)
mae_hqic <- MAE(test, forecast_hqic)

rmse_aic <- RMSE(test, forecast_aic)
rmse_bic <- RMSE(test, forecast_bic)
rmse_hqic <- RMSE(test, forecast_hqic)

mse_aic <- MSE(test, forecast_aic)
mse_bic <- MSE(test, forecast_bic)
mse_hqic <- MSE(test, forecast_hqic)

r2_aic <- R2(test, forecast_aic)
r2_bic <- R2(test, forecast_bic)
r2_hqic <- R2(test, forecast_hqic)

metrics_table <- data.frame(
  Method = c("AIC", "BIC", "HQIC"),
  MAPE = c(mape_aic, mape_bic, mape_hqic),
  MAE = c(mae_aic, mae_bic, mae_hqic),
  RMSE = c(rmse_aic, rmse_bic, rmse_hqic),
  MSE = c(mse_aic, mse_bic, mse_hqic),
  R2 = c(r2_aic, r2_bic, r2_hqic)
)

print(metrics_table)

##   Method     MAPE       MAE     RMSE      MSE        R2
## 1    AIC 1.947686 0.7535679 1.078031 1.162151 0.9249987
## 2    BIC 2.016224 0.7780951 1.085067 1.177371 0.9254211
## 3   HQIC 1.947686 0.7535679 1.078031 1.162151 0.9249987

plot(test, forecast_aic, main="AIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_bic, main="BIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_hqic, main="HQIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

residuals_aic <- test - forecast_aic
residuals_bic <- test - forecast_bic
residuals_hqic <- test - forecast_hqic

plot(residuals_aic, main="AIC Residuals", ylab="Residuals", col="blue")
abline(0,0)

plot(residuals_bic, main="BIC Residuals", ylab="Residuals", col="blue")
abline(0,0)

plot(residuals_hqic, main="HQIC Residuals", ylab="Residuals", col="blue")
abline(0,0)

shapiro.test(residuals_aic)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_aic
## W = 0.89981, p-value = 0.0347

shapiro.test(residuals_bic)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_bic
## W = 0.91446, p-value = 0.06728

shapiro.test(residuals_hqic)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_hqic
## W = 0.89981, p-value = 0.0347

n = 28

nb_data <- nrow(tgls)
nb_day_predicted = c(14, 21, 28)
train_size = nb_data - max(nb_day_predicted)
train <- tgls$TGLS.Close[1:train_size]

test <- tgls$TGLS.Close[(train_size + 1):(train_size + 28)]
test <- test[complete.cases(test)]
test_dates <- index(tgls)[(train_size + 1):(train_size + 28)]
test_dates <- test_dates[complete.cases(test_dates)]

forecast_aic <- c()
forecast_bic <- c()
forecast_hqic <- c()

for (i in 1:28) {
  fit_aic <- arima(train, order = best_order_aic)
  fit_bic <- arima(train, order = best_order_bic)
  fit_hqic <- arima(train, order = best_order_hqic)
  
  forecast_aic <- c(forecast_aic, predict(fit_aic, n.ahead = 1)$pred)
  forecast_bic <- c(forecast_bic, predict(fit_bic, n.ahead = 1)$pred)
  forecast_hqic <- c(forecast_hqic, predict(fit_hqic, n.ahead = 1)$pred)

  train <- c(train, test[i])
}

n_last = 50
last_dates <- tail(index(tgls), n_last)
last_tgls <- tail(tgls$TGLS.Close, n_last)
last_forecast_aic <- tail(forecast_aic, n_last)
last_forecast_bic <- tail(forecast_bic, n_last)
last_forecast_hqic <- tail(forecast_hqic, n_last)

plot(last_dates, last_tgls, type = "l", col = "black", ylim = range(c(last_tgls, last_forecast_aic, last_forecast_bic, last_forecast_hqic)), 
     main = "TGLS Closing Price", ylab = "Price ($)", xlab = "Date")
lines(test_dates, forecast_aic, col = "red")
lines(test_dates, forecast_bic, col = "green")
lines(test_dates, forecast_hqic, col = "purple")

legend("topleft", legend = c("Actual", "Forecast AIC", "Forecast BIC", "Forecast HQIC"), 
       col = c("black", "red", "green", "purple"), lty = c(1, 1, 1, 1))

mape_aic <- MAPE(test, forecast_aic)
mape_bic <- MAPE(test, forecast_bic)
mape_hqic <- MAPE(test, forecast_hqic)

mae_aic <- MAE(test, forecast_aic)
mae_bic <- MAE(test, forecast_bic)
mae_hqic <- MAE(test, forecast_hqic)

rmse_aic <- RMSE(test, forecast_aic)
rmse_bic <- RMSE(test, forecast_bic)
rmse_hqic <- RMSE(test, forecast_hqic)

mse_aic <- MSE(test, forecast_aic)
mse_bic <- MSE(test, forecast_bic)
mse_hqic <- MSE(test, forecast_hqic)

r2_aic <- R2(test, forecast_aic)
r2_bic <- R2(test, forecast_bic)
r2_hqic <- R2(test, forecast_hqic)

metrics_table <- data.frame(
  Method = c("AIC", "BIC", "HQIC"),
  MAPE = c(mape_aic, mape_bic, mape_hqic),
  MAE = c(mae_aic, mae_bic, mae_hqic),
  RMSE = c(rmse_aic, rmse_bic, rmse_hqic),
  MSE = c(mse_aic, mse_bic, mse_hqic),
  R2 = c(r2_aic, r2_bic, r2_hqic)
)

print(metrics_table)

##   Method     MAPE       MAE     RMSE      MSE        R2
## 1    AIC 1.906612 0.7641704 1.038412 1.078300 0.9495750
## 2    BIC 1.957943 0.7821428 1.057291 1.117864 0.9478544
## 3   HQIC 1.906612 0.7641704 1.038412 1.078300 0.9495750

plot(test, forecast_aic, main="AIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_bic, main="BIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_hqic, main="HQIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

residuals_aic <- test - forecast_aic
residuals_bic <- test - forecast_bic
residuals_hqic <- test - forecast_hqic

plot(residuals_aic, main="AIC Residuals", ylab="Residuals", col="blue")
abline(0,0)

plot(residuals_bic, main="BIC Residuals", ylab="Residuals", col="blue")
abline(0,0)

plot(residuals_hqic, main="HQIC Residuals", ylab="Residuals", col="blue")
abline(0,0)

shapiro.test(residuals_aic)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_aic
## W = 0.92826, p-value = 0.05574

shapiro.test(residuals_bic)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_bic
## W = 0.94246, p-value = 0.1276

shapiro.test(residuals_hqic)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_hqic
## W = 0.92826, p-value = 0.05574

With more values predicted, we obtain better r2 and a low p-value for normality test so we can assume with a good confidence thqt the residuals follow a normal distribution.

Without Rolling Forecast

n = 14

nb_data <- nrow(tgls)
nb_day_predicted = c(14, 21, 28)
train_size = nb_data - max(nb_day_predicted)
train <- tgls$TGLS.Close[1:train_size]

test <- tgls$TGLS.Close[(train_size + 1):(train_size + 14)]
test <- test[complete.cases(test)]
test_dates <- index(tgls)[(train_size + 1):(train_size + 14)]
test_dates <- test_dates[complete.cases(test_dates)]

fit_aic <- arima(train, order = best_order_aic)
fit_bic <- arima(train, order = best_order_bic)
fit_hqic <- arima(train, order = best_order_hqic)

forecast_aic <- predict(fit_aic, n.ahead = 14)$pred
forecast_bic <- predict(fit_bic, n.ahead = 14)$pred
forecast_hqic <- predict(fit_hqic, n.ahead = 14)$pred

n_last = 50
last_dates <- tail(index(tgls), n_last)
last_tgls <- tail(tgls$TGLS.Close, n_last)
last_forecast_aic <- tail(forecast_aic, n_last)
last_forecast_bic <- tail(forecast_bic, n_last)
last_forecast_hqic <- tail(forecast_hqic, n_last)

plot(last_dates, last_tgls, type = "l", col = "black", ylim = range(c(last_tgls, last_forecast_aic, last_forecast_bic, last_forecast_hqic)), 
     main = "TGLS Closing Price", ylab = "Price ($)", xlab = "Date")
lines(test_dates, forecast_aic, col = "red")
lines(test_dates, forecast_bic, col = "green")
lines(test_dates, forecast_hqic, col = "purple")

legend("topleft", legend = c("Actual", "Forecast AIC", "Forecast BIC", "Forecast HQIC"), 
       col = c("black", "red", "green", "purple"), lty = c(1, 1, 1, 1))

mape_aic <- MAPE(test, forecast_aic)
mape_bic <- MAPE(test, forecast_bic)
mape_hqic <- MAPE(test, forecast_hqic)

mae_aic <- MAE(test, forecast_aic)
mae_bic <- MAE(test, forecast_bic)
mae_hqic <- MAE(test, forecast_hqic)

rmse_aic <- RMSE(test, forecast_aic)
rmse_bic <- RMSE(test, forecast_bic)
rmse_hqic <- RMSE(test, forecast_hqic)

mse_aic <- MSE(test, forecast_aic)
mse_bic <- MSE(test, forecast_bic)
mse_hqic <- MSE(test, forecast_hqic)

r2_aic <- R2(test, forecast_aic)
r2_bic <- R2(test, forecast_bic)
r2_hqic <- R2(test, forecast_hqic)

metrics_table <- data.frame(
  Method = c("AIC", "BIC", "HQIC"),
  MAPE = c(mape_aic, mape_bic, mape_hqic),
  MAE = c(mae_aic, mae_bic, mae_hqic),
  RMSE = c(rmse_aic, rmse_bic, rmse_hqic),
  MSE = c(mse_aic, mse_bic, mse_hqic),
  R2 = c(r2_aic, r2_bic, r2_hqic)
)

print(metrics_table)

##   Method     MAPE      MAE     RMSE      MSE         R2
## 1    AIC 4.377584 1.606444 1.838467 3.379960 0.01850756
## 2    BIC 4.407885 1.615714 1.834709 3.366157         NA
## 3   HQIC 4.377584 1.606444 1.838467 3.379960 0.01850756

plot(test, forecast_aic, main="AIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_bic, main="BIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_hqic, main="HQIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

n = 21

nb_data <- nrow(tgls)
nb_day_predicted = c(14, 21, 28)
train_size = nb_data - max(nb_day_predicted)
train <- tgls$TGLS.Close[1:train_size]

test <- tgls$TGLS.Close[(train_size + 1):(train_size + 21)]
test <- test[complete.cases(test)]
test_dates <- index(tgls)[(train_size + 1):(train_size + 21)]
test_dates <- test_dates[complete.cases(test_dates)]

fit_aic <- arima(train, order = best_order_aic)
fit_bic <- arima(train, order = best_order_bic)
fit_hqic <- arima(train, order = best_order_hqic)

forecast_aic <- predict(fit_aic, n.ahead = 21)$pred
forecast_bic <- predict(fit_bic, n.ahead = 21)$pred
forecast_hqic <- predict(fit_hqic, n.ahead = 21)$pred

n_last = 50
last_dates <- tail(index(tgls), n_last)
last_tgls <- tail(tgls$TGLS.Close, n_last)
last_forecast_aic <- tail(forecast_aic, n_last)
last_forecast_bic <- tail(forecast_bic, n_last)
last_forecast_hqic <- tail(forecast_hqic, n_last)

plot(last_dates, last_tgls, type = "l", col = "black", ylim = range(c(last_tgls, last_forecast_aic, last_forecast_bic, last_forecast_hqic)), 
     main = "TGLS Closing Price", ylab = "Price ($)", xlab = "Date")
lines(test_dates, forecast_aic, col = "red")
lines(test_dates, forecast_bic, col = "green")
lines(test_dates, forecast_hqic, col = "purple")

legend("topleft", legend = c("Actual", "Forecast AIC", "Forecast BIC", "Forecast HQIC"), 
       col = c("black", "red", "green", "purple"), lty = c(1, 1, 1, 1))

mape_aic <- MAPE(test, forecast_aic)
mape_bic <- MAPE(test, forecast_bic)
mape_hqic <- MAPE(test, forecast_hqic)

mae_aic <- MAE(test, forecast_aic)
mae_bic <- MAE(test, forecast_bic)
mae_hqic <- MAE(test, forecast_hqic)

rmse_aic <- RMSE(test, forecast_aic)
rmse_bic <- RMSE(test, forecast_bic)
rmse_hqic <- RMSE(test, forecast_hqic)

mse_aic <- MSE(test, forecast_aic)
mse_bic <- MSE(test, forecast_bic)
mse_hqic <- MSE(test, forecast_hqic)

r2_aic <- R2(test, forecast_aic)
r2_bic <- R2(test, forecast_bic)
r2_hqic <- R2(test, forecast_hqic)

metrics_table <- data.frame(
  Method = c("AIC", "BIC", "HQIC"),
  MAPE = c(mape_aic, mape_bic, mape_hqic),
  MAE = c(mae_aic, mae_bic, mae_hqic),
  RMSE = c(rmse_aic, rmse_bic, rmse_hqic),
  MSE = c(mse_aic, mse_bic, mse_hqic),
  R2 = c(r2_aic, r2_bic, r2_hqic)
)

print(metrics_table)

##   Method     MAPE      MAE     RMSE      MSE          R2
## 1    AIC 8.098980 3.297485 4.333084 18.77562 0.008232079
## 2    BIC 8.091165 3.291905 4.314561 18.61544          NA
## 3   HQIC 8.098980 3.297485 4.333084 18.77562 0.008232079

plot(test, forecast_aic, main="AIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_bic, main="BIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_hqic, main="HQIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

n = 28

nb_data <- nrow(tgls)
nb_day_predicted = c(14, 21, 28)
train_size = nb_data - max(nb_day_predicted)
train <- tgls$TGLS.Close[1:train_size]

test <- tgls$TGLS.Close[(train_size + 1):(train_size + 28)]
test <- test[complete.cases(test)]
test_dates <- index(tgls)[(train_size + 1):(train_size + 28)]
test_dates <- test_dates[complete.cases(test_dates)]

fit_aic <- arima(train, order = best_order_aic)
fit_bic <- arima(train, order = best_order_bic)
fit_hqic <- arima(train, order = best_order_hqic)

forecast_aic <- predict(fit_aic, n.ahead = 28)$pred
forecast_bic <- predict(fit_bic, n.ahead = 28)$pred
forecast_hqic <- predict(fit_hqic, n.ahead = 28)$pred

n_last = 50
last_dates <- tail(index(tgls), n_last)
last_tgls <- tail(tgls$TGLS.Close, n_last)
last_forecast_aic <- tail(forecast_aic, n_last)
last_forecast_bic <- tail(forecast_bic, n_last)
last_forecast_hqic <- tail(forecast_hqic, n_last)

plot(last_dates, last_tgls, type = "l", col = "black", ylim = range(c(last_tgls, last_forecast_aic, last_forecast_bic, last_forecast_hqic)), 
     main = "TGLS Closing Price", ylab = "Price ($)", xlab = "Date")
lines(test_dates, forecast_aic, col = "red")
lines(test_dates, forecast_bic, col = "green")
lines(test_dates, forecast_hqic, col = "purple")

legend("topleft", legend = c("Actual", "Forecast AIC", "Forecast BIC", "Forecast HQIC"), 
       col = c("black", "red", "green", "purple"), lty = c(1, 1, 1, 1))

mape_aic <- MAPE(test, forecast_aic)
mape_bic <- MAPE(test, forecast_bic)
mape_hqic <- MAPE(test, forecast_hqic)

mae_aic <- MAE(test, forecast_aic)
mae_bic <- MAE(test, forecast_bic)
mae_hqic <- MAE(test, forecast_hqic)

rmse_aic <- RMSE(test, forecast_aic)
rmse_bic <- RMSE(test, forecast_bic)
rmse_hqic <- RMSE(test, forecast_hqic)

mse_aic <- MSE(test, forecast_aic)
mse_bic <- MSE(test, forecast_bic)
mse_hqic <- MSE(test, forecast_hqic)

r2_aic <- R2(test, forecast_aic)
r2_bic <- R2(test, forecast_bic)
r2_hqic <- R2(test, forecast_hqic)

metrics_table <- data.frame(
  Method = c("AIC", "BIC", "HQIC"),
  MAPE = c(mape_aic, mape_bic, mape_hqic),
  MAE = c(mae_aic, mae_bic, mae_hqic),
  RMSE = c(rmse_aic, rmse_bic, rmse_hqic),
  MSE = c(mse_aic, mse_bic, mse_hqic),
  R2 = c(r2_aic, r2_bic, r2_hqic)
)

print(metrics_table)

##   Method     MAPE      MAE     RMSE      MSE          R2
## 1    AIC 11.30916 4.840822 6.067108 36.80981 0.004342446
## 2    BIC 11.28383 4.827857 6.043484 36.52369          NA
## 3   HQIC 11.30916 4.840822 6.067108 36.80981 0.004342446

plot(test, forecast_aic, main="AIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_bic, main="BIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

plot(test, forecast_hqic, main="HQIC", xlab="Observed", ylab="Predicted", col="blue")
abline(0,1)

Without rolling forecast, the model doesn’t manage to predict the price as we can see with the graphics. Consequently we won’t do the analysis of the residuals.