ImputacionR

Code

library(mice)

Warning: package 'mice' was built under R version 4.2.3

Attaching package: 'mice'

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

    filter

The following objects are masked from 'package:base':

    cbind, rbind

Code

library(dplyr)

Warning: package 'dplyr' was built under R version 4.2.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

Code

library(VIM)

Warning: package 'VIM' was built under R version 4.2.3

Loading required package: colorspace

Loading required package: grid

The legacy packages maptools, rgdal, and rgeos, underpinning this package
will retire shortly. Please refer to R-spatial evolution reports on
https://r-spatial.org/r/2023/05/15/evolution4.html for details.
This package is now running under evolution status 0 

VIM is ready to use.

Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues

Attaching package: 'VIM'

The following object is masked from 'package:datasets':

    sleep

Code

library(foreign)

Warning: package 'foreign' was built under R version 4.2.2

Code

library(ggplot2)

Warning: package 'ggplot2' was built under R version 4.2.3

Code

library(lattice)

Warning: package 'lattice' was built under R version 4.2.3

Code

library(outliers)
library(EnvStats)

Warning: package 'EnvStats' was built under R version 4.2.3

Attaching package: 'EnvStats'

The following objects are masked from 'package:stats':

    predict, predict.lm

Imputacion

Se hace lectura de los datos

Code

Diabetes<-read.csv("https://raw.githubusercontent.com/lihkir/Data/main/diabetes.csv")

Se reemplaza con NA los valores de las mediciones medicas que no tienen sentido que sean 0, lo cual representaria la falta de informacion de ese valor.

Code

columnas<-colnames(Diabetes[,-c(ncol(Diabetes),1)])

for (i in columnas){
  Diabetes[[i]][Diabetes[[i]]==0]<-NA
  }

Observamos las primeras filas del conjunto de datos con NA (Pregnancies y Outcome si pueden tomar valores de 0)

Code

head(Diabetes)

  Pregnancies Glucose BloodPressure SkinThickness Insulin  BMI
1           6     148            72            35      NA 33.6
2           1      85            66            29      NA 26.6
3           8     183            64            NA      NA 23.3
4           1      89            66            23      94 28.1
5           0     137            40            35     168 43.1
6           5     116            74            NA      NA 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

Code

colnames(Diabetes)<-c("Prg","Glu","BlP","SkT",
                      "Ins","BMI","DPF","Age","Out")

md.pattern(Diabetes,plot=TRUE,rotate.names=TRUE)

    Prg DPF Age Out Glu BMI BlP SkT Ins    
392   1   1   1   1   1   1   1   1   1   0
140   1   1   1   1   1   1   1   1   0   1
192   1   1   1   1   1   1   1   0   0   2
2     1   1   1   1   1   1   0   1   0   2
26    1   1   1   1   1   1   0   0   0   3
1     1   1   1   1   1   0   1   1   1   1
1     1   1   1   1   1   0   1   1   0   2
2     1   1   1   1   1   0   1   0   0   3
7     1   1   1   1   1   0   0   0   0   4
1     1   1   1   1   0   1   1   1   1   1
4     1   1   1   1   0   1   1   1   0   2
      0   0   0   0   5  11  35 227 374 652

se puede observar que hay 392 filas sin datos faltantes, 140 filas que solo tienen faltantes en la caracteristica de insulina, 26 filas que contienen faltantes en BLB, SkinThickness e Insuline, y observando la matriz dada se evidencia que esas 3 caracteristicas son las que representan casi todos los datos faltantes, donde hay unos 5 y 11 para glucose y BMI los cuales son pocos.

Code

aggr(Diabetes, numbers=TRUE)

Como se evidencio anteriormente, las caracteristicas que mas valores faltantes presentan son: Ins,SKT,MIP. Donde SkT prsenta casi un 30% de datos faltantes y en Ins podemos observar mas del 40% de datos faltantes. El 25% de los datos les hace falta los datos de Ins y Skt y el 18% solo el dato de Ins y el 51% de filas no tienen datos faltantes, lo cual es muy poco.

Se haran imputacion de datos por emparejamiento predictivo medio para los valores faltantes y se evaluara el resultado en comparacion de los datos originales, con el fin de encontrar la mejor manera de manejar estos valores faltantes.

Code

#Metodo pmm
imp_pmm <- mice(Diabetes, m=5, maxit=50, method ='pmm', seed=500, printFlag = FALSE)
imp_pmm_df<-complete(imp_pmm)
#Metodo norm.predict
imp_norm.pred <- mice(Diabetes, m=5, maxit=50, method ='norm.predict', seed=500, printFlag = FALSE)
imp_norm.pred_df<-complete(imp_norm.pred)
#Metodo norm.nob
imp_norm.nob <- mice(Diabetes, m=5, maxit=50, method ='norm.nob', seed=500, printFlag = FALSE)
imp_norm.nob_df<-complete(imp_norm.nob)
#Metodo norm
imp_norm <- mice(Diabetes, m=5, maxit=50, method ='norm', seed=500, printFlag = FALSE)
imp_norm_df<-complete(imp_norm)

Code

head(imp_pmm_df)

  Prg Glu BlP SkT Ins  BMI   DPF Age Out
1   6 148  72  35  83 33.6 0.627  50   1
2   1  85  66  29  55 26.6 0.351  31   0
3   8 183  64  20 175 23.3 0.672  32   1
4   1  89  66  23  94 28.1 0.167  21   0
5   0 137  40  35 168 43.1 2.288  33   1
6   5 116  74  24 175 25.6 0.201  30   0

Code

dat<-bind_rows(Diabetes %>% mutate(df= "Diabetes"),
               imp_pmm_df %>% mutate(df= "imp_pmm"),
               imp_norm.pred_df %>% mutate(df= "imp_norm.pred"),
               imp_norm.nob_df %>% mutate(df= "imp_norm.nob"),
               imp_norm_df %>% mutate(df= "imp_norm"))

Code

variable<-"Ins"
# Crear el gráfico de frecuencias de líneas
ggplot(dat, aes(x = Ins, color = df)) +
  geom_freqpoly(linewidth=1) +
  theme_minimal() +
  labs(x = "Ins", y = "Frecuencia", title = "Frecuencia de lineas Ins")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 374 rows containing non-finite values (`stat_bin()`).

Code

datos_lista<-list(Diabetes$Ins,imp_pmm_df$Ins,imp_norm_df$Ins,imp_norm.nob_df$Ins,imp_norm.pred_df$Ins)

minimos <- numeric(length(datos_lista))

for (i in seq_along(datos_lista)) {
  minimos[i] <- min(datos_lista[[i]], na.rm = TRUE)
}

print(minimos)

[1]   14.00000   14.00000 -211.49713 -156.09009  -22.18447

Debido a que los metodos “norm”,“norm.nob”,“norm.pred” resultan en datos negativos no es util para imputar dichos datos, debido a que no es posible valores negativos para estas caracteristicas, lo cual se usaria la imputacion usando el metodo “pmm”

Code

xyplot(imp_pmm, Ins ~ Prg + Glu + BlP + SkT + BMI + DPF + Age + Out, pch=18, cex=1)

realizando unos graficos de dispersion podemos observar que los datos originales y los imputados siguen la mismas tendencias.

Code

densityplot(imp_pmm)

Como se dijo previamente, los datos inputados siguen distribuciones muy similares a los datos originales.

Datos Atipicos

Histograma

Code

par(mfrow=c(2,4))

hist(imp_pmm_df$Prg,xlab = "Prg",main = "Histogram of Prg ",breaks = sqrt(nrow(imp_pmm_df)))
hist(imp_pmm_df$Glu,xlab = "Glu",main = "Histogram of Glu",breaks = sqrt(nrow(imp_pmm_df)))
hist(imp_pmm_df$BlP,xlab = "BlP",main = "Histogram of BlP",breaks = sqrt(nrow(imp_pmm_df)))
hist(imp_pmm_df$SkT,xlab = "SkT",main = "Histogram of SkT",breaks = sqrt(nrow(imp_pmm_df)))
hist(imp_pmm_df$Ins,xlab = "Ins",main = "Histogram of Ins",breaks = sqrt(nrow(imp_pmm_df)))
hist(imp_pmm_df$BMI,xlab = "BMI",main = "Histogram of BMI",breaks = sqrt(nrow(imp_pmm_df)))
hist(imp_pmm_df$DPF,xlab = "DPF",main = "Histogram of DPF",breaks = sqrt(nrow(imp_pmm_df)))
hist(imp_pmm_df$Age,xlab = "Age",main = "Histogram of Age",breaks = sqrt(nrow(imp_pmm_df)))

Se puede observar que la mayoria de caracteristicas son sesgadas a la izquierda, en la cual se pueden observar claros datos atipicos en SkT, en Age con un valor de 80, BMI con valores aproximadamente 70 y en INS valores por encima de los 600, lo cual teniendo en cuenta que en los datos originales casi el 50% eran faltantes, por medio de practicidad academica se dejara la variable, la cual es candidata a eliminar por la cantidad de datos faltantes.

Boxplot

Code

par(mfrow=c(2,4))

boxplot(imp_pmm_df$Prg,ylab = "Prg")
boxplot(imp_pmm_df$Glu,ylab = "Glu")
boxplot(imp_pmm_df$BlP,ylab = "BlP")
boxplot(imp_pmm_df$SkT,ylab = "SkT")
boxplot(imp_pmm_df$Ins,ylab = "Ins")
boxplot(imp_pmm_df$BMI,ylab = "BMI")
boxplot(imp_pmm_df$DPF,ylab = "DPF")
boxplot(imp_pmm_df$Age,ylab = "Age")

Se observan Muchos valores atipicos para DPF INS y BMi, unos cuantos para Age, SkT, y Prg.

Code

out1<-boxplot.stats(imp_pmm_df$Prg)$out
out3<-boxplot.stats(imp_pmm_df$BlP)$out
out4<-boxplot.stats(imp_pmm_df$SkT)$out
out5<-boxplot.stats(imp_pmm_df$Ins)$out
out6<-boxplot.stats(imp_pmm_df$BMI)$out
out7<-boxplot.stats(imp_pmm_df$DPF)$out
out8<-boxplot.stats(imp_pmm_df$Age)$out

out1_ind <- which(imp_pmm_df$Prg %in% c(out1))
out3_ind <- which(imp_pmm_df$BlP %in% c(out3))
out4_ind <- which(imp_pmm_df$SkT %in% c(out4))
out5_ind <- which(imp_pmm_df$Ins %in% c(out5))
out6_ind <- which(imp_pmm_df$BMI %in% c(out6))
out7_ind <- which(imp_pmm_df$DPF %in% c(out7))
out8_ind <- which(imp_pmm_df$Age %in% c(out8))

Estos serian los datos atipicos que corresponden a cada caracteristica.

Code

imp_pmm_df_copy<-imp_pmm_df

imp_pmm_df_copy[out1_ind,which(colnames(imp_pmm_df)=="Prg")]<-NA
imp_pmm_df_copy[out3_ind,which(colnames(imp_pmm_df)=="BlP")]<-NA
imp_pmm_df_copy[out4_ind,which(colnames(imp_pmm_df)=="SkT")]<-NA
imp_pmm_df_copy[out5_ind,which(colnames(imp_pmm_df)=="Ins")]<-NA
imp_pmm_df_copy[out6_ind,which(colnames(imp_pmm_df)=="BMI")]<-NA
imp_pmm_df_copy[out7_ind,which(colnames(imp_pmm_df)=="DPF")]<-NA
imp_pmm_df_copy[out8_ind,which(colnames(imp_pmm_df)=="Age")]<-NA

Se eliminan los datos Atipicos encontrados y se asignan NA, para luego hacer imputacion de datos mediante el metodo pmm y tener unos datos sin valores atipicos.

Code

imp_pmm_2 <- mice(imp_pmm_df_copy, m=5, maxit=50, method ='pmm', seed=500, printFlag = FALSE)
imp_pmm_2_df<-complete(imp_pmm_2)
head(imp_pmm_2_df)

  Prg Glu BlP SkT Ins  BMI   DPF Age Out
1   6 148  72  35  83 33.6 0.627  50   1
2   1  85  66  29  55 26.6 0.351  31   0
3   8 183  64  20 175 23.3 0.672  32   1
4   1  89  66  23  94 28.1 0.167  21   0
5   0 137  40  35 168 43.1 0.227  33   1
6   5 116  74  24 175 25.6 0.201  30   0

A modo de practica academica se procedera a realizar las otras maneras de identificacion de datos atipicos de la columna “BIP” sin imputacion, solo extrayendo los valores atipicos.

percentiles

Code

lower_bound <- quantile(imp_pmm_2_df$BMI, 0.025)
upper_bound <- quantile(imp_pmm_2_df$BMI, 0.975)

outlier_ind <- which(imp_pmm_2_df$BMI < lower_bound | imp_pmm_2_df$BMI > upper_bound)
outlier_ind

 [1]  10  34  51  58  69  85  91  93  98 100 155 156 194 204 232 240 289 317 336
[20] 379 382 419 439 446 454 461 470 488 520 527 559 591 608 618 640 682 690 747
[39] 748

Datos atipicos mediante el metodo de percentiles

Filtro de Hampel

Code

lower_bound <- median(imp_pmm_2_df$BMI) - 3 * mad(imp_pmm_2_df$BMI, constant = 1)
upper_bound <- median(imp_pmm_2_df$BMI) + 3 * mad(imp_pmm_2_df$BMI, constant = 1)

outlier_ind <- which(imp_pmm_2_df$BMI < lower_bound | imp_pmm_2_df$BMI > upper_bound)
outlier_ind

 [1]  10  58  85  93 100 155 156 194 232 240 336 379 419 439 446 470 488 527 559
[20] 591 682 690 747 748

Datos atipicos mediante el metodo de Hampel

Prueba de Grubbs

Code

test <- grubbs.test(imp_pmm_2_df$BMI)
test

    Grubbs test for one outlier

data:  imp_pmm_2_df$BMI
G = 2.70371, U = 0.99046, p-value = 1
alternative hypothesis: highest value 50 is an outlier

El valor p es de 1. Al nivel de significación del 5%, no rechazamos la hipótesis de que el valor más alto 50 no es un valor atípico. No tenemos evidencia suficiente para decir que 50 es un valor atípico.

Code

test <- grubbs.test(imp_pmm_2_df$BMI,opposite = TRUE)
test

    Grubbs test for one outlier

data:  imp_pmm_2_df$BMI
G = 2.15842, U = 0.99392, p-value = 1
alternative hypothesis: lowest value 18.2 is an outlier

El valor p es de 1. Al nivel de significación del 5%, no rechazamos la hipótesis de que el valor más bajo 18.2 no es un valor atípico. No tenemos evidencia suficiente para decir que 18.2 es un valor atípico.

Prueba de Dixon

Code

#dado que la prueba dixon solo acepta conjuntos de datos de 3 a 30 obs se hace un sample en modo de demostracion
test <- dixon.test(sample(imp_pmm_2_df$BMI,30))
test

    Dixon test for outliers

data:  sample(imp_pmm_2_df$BMI, 30)
Q = 0.096899, p-value = 0.3924
alternative hypothesis: highest value 47.9 is an outlier

El valor p es de 0.3536. Al nivel de significación del 5%, no rechazamos la hipótesis de que el valor más alto 45.6 no es un valor atípico. No tenemos evidencia suficiente para decir que 45.6 es un valor atípico.

Code

test <- dixon.test(sample(imp_pmm_2_df$BMI,30),opposite = TRUE)
test

    Dixon test for outliers

data:  sample(imp_pmm_2_df$BMI, 30)
Q = 0.048913, p-value = 0.1211
alternative hypothesis: lowest value 21.8 is an outlier

El valor p es de 1. Al nivel de significación del 5%, no rechazamos la hipótesis de que el valor más bajo 21.8 no es un valor atípico. No tenemos evidencia suficiente para decir que 21.8 es un valor atípico. ### Prueba de Rosner

Code

test <- rosnerTest(imp_pmm_2_df$BMI, k = 10)
test$all.stats

   i   Mean.i     SD.i Value Obs.Num    R.i+1 lambda.i+1 Outlier
1  0 32.31680 6.540350  50.0     156 2.703709   3.974092   FALSE
2  1 32.29374 6.513315  49.7     100 2.672411   3.973762   FALSE
3  2 32.27102 6.487077  49.6     682 2.671308   3.973432   FALSE
4  3 32.24837 6.460935  49.3     747 2.639190   3.973102   FALSE
5  4 32.22605 6.435590  48.8      85 2.575359   3.972771   FALSE
6  5 32.20433 6.411724  48.3     379 2.510351   3.972440   FALSE
7  6 32.18320 6.389315  47.9     155 2.459857   3.972108   FALSE
8  7 32.16255 6.368015  47.9     336 2.471328   3.971776   FALSE
9  8 32.14184 6.346519  46.8      10 2.309638   3.971443   FALSE
10 9 32.12253 6.328318  46.8      58 2.319332   3.971110   FALSE

La prueba de Rosner arroja que no existen Outliers para la columna observada

Tratamiento de datos atipicos con Capping

Code

x <- imp_pmm_2_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]
head(x)

[1] 33.6 26.6 23.3 28.1 43.1 25.6