¿Como tratar los valores atipicos?

El método, conocido como filtro de Hampel, se basa en la identificación de valores atípicos dentro de un conjunto de datos mediante la delimitación de un intervalo específico.

Intervalo = [ mediana(X) - 3 * MAD(X) ; mediana(X) + 3 * MAD(X) ]

hayamos el intervalo utilizando nuestra variable dependiente

limite_inferior <-  median(datos$`HH En la Maquina`) - 3 * 
  mad(datos$`HH En la Maquina`, constant = 1)

limite_inferior

[1] -8

limite_superior <-  median(datos$`HH En la Maquina`) + 3 * 
  mad(datos$`HH En la Maquina`, constant = 1)

limite_superior

[1] 22

Una vez definimos las margenes para el filtro extraemos las filas consideradas atípicas:

atipicos <- which(datos$`HH En la Maquina` <
                    limite_inferior | 
                    datos$`HH En la Maquina` > 
                    limite_superior)


### vemos la cantidad de valores atípicos 
atipicos |> length()

[1] 5251

De acuerdo a los resultados del filtro de hampel existen 5251 observaciones en donde las horas utilizadas para reparar la maquina son atípicas y muy diferentes a los datos esperados.

## filtramos los valores atípicos
datos <- datos[-atipicos,  ]

par(mfrow = c(1, 4))


boxplot(datos$`HH Actv Gnerales`, 
        main = "HH Actv Gnerales")


boxplot(datos$`HH En la Maquina`, 
        main = "HH En la Maquina")


boxplot(datos$`Otras HH`,
        main = "Otras HH")

boxplot(datos$dias_para_diferencia,
        main = "Dias de diferencia")

Al aplicar el filtro de hampel, el gráfico de caja de la variable que representa las horas reales utilizadas en la maquina refleja una menor cantidad de valores atípicos.

Ahora aplicamos el mismo procedimiento para las otras 3 variables.

## creamos una función para extraer el valor atipico

atipicos <- function(x){
  
## encontramos los valores atipicos para la variable que representa horas generales

limite_inferior <-  median(x) - 3 * 
  mad(x, constant = 1)

limite_superior <-  median(x) + 3 * 
  mad(x, constant = 1)


## encontramos los valores atipicos
atipicos <- which(x < limite_inferior | 
                   x  > limite_superior)

  return(atipicos)
}


## filtramos las horas generales
datos <- datos[-atipicos(datos$`HH Actv Gnerales`), ]

## filtramos por otras horas h

datos <- datos[-atipicos(datos$`Otras HH`), ]

## filtramos por dias generales

datos <- datos[-atipicos(datos$dias_para_diferencia), ]

Vemos los resultados luego de aplicar los filtros:

par(mfrow = c(1, 4))


boxplot(datos$`HH Actv Gnerales`, 
        main = "HH Actv Gnerales")


boxplot(datos$`HH En la Maquina`, 
        main = "HH En la Maquina")


boxplot(datos$`Otras HH`,
        main = "Otras HH")

boxplot(datos$dias_para_diferencia,
        main = "Dias de diferencia")

Se observa una menor dispersión y una menor cantidad de valores atípicos.

Entendiendo el modelo

• (Intercept): El coeficiente sugiere que, en promedio, el número de horas en la máquina aumenta en 3.271 cuando todas las demás variables son cero.

• DisciplinaInstrumentos: Un coeficiente positivo de 0.384 sugiere que la presencia de la disciplina de Instrumentos se asocia con un aumento en el número de horas en la máquina.

• DisciplinaLinea de Mtto: Un coeficiente negativo de -1.566 sugiere que la presencia de la disciplina de Línea de Mantenimiento está asociada con una disminución en el número de horas en la máquina.

• DisciplinaMecanica: Un coeficiente positivo de 1.037 sugiere que la presencia de la disciplina Mecánica está asociada con un aumento en el número de horas en la máquina.

• DisciplinaMecanica CBM: Un coeficiente positivo de 0.862 sugiere que la presencia de la disciplina Mecánica CBM está asociada con un aumento en el número de horas en la máquina.

• DisciplinaMecanica VAL: Un coeficiente positivo de 2.803 sugiere que la presencia de la disciplina Mecánica VAL está asociada con un aumento significativo en el número de horas en la máquina.

• DisciplinaPozoz: Un coeficiente negativo de -1.933 sugiere que la presencia de la disciplina Pozoz está asociada con una disminución en el número de horas en la máquina.

• EpocaAñoVerano: Un coeficiente negativo de -0.242 sugiere que la temporada de verano está asociada con una ligera disminución en el número de horas en la máquina.

• HH_Actv_Gnerales: Un coeficiente positivo de 2.040 sugiere que un aumento en las horas de trabajo generales está asociado con un aumento en el número de horas en la máquina.

• Otras_HH: Un coeficiente positivo de 1.693 sugiere que un aumento en las horas de otros trabajos está asociado con un aumento significativo en el número de horas en la máquina.

Comprobando supuestos

par(mfrow = c(1, 4))

plot(modelo_lineal)

Los graficos del modelo parecen mostrar que no se cumplen los supuestos de linealidad entre las variables, el grafico de q-q residuales no parece mostrar evidencia de normalidad, ademas parece existir heterocedasticidad en los residuales.

Prueba de normalidad

Prueba de Kolmogorov-Smirnov: Esta prueba es una opción para evaluar la normalidad de los datos. La función ks.test() en R puede ser utilizada para esto.

Planteamiento de Hipótesis:

• \(H_0\): Los residuales se distribuyen de manera normal.
• \(H_1\): Los residuales NO se distribuyen de maneran normal.

ks.test(modelo_lineal$residuals, "pnorm")

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  modelo_lineal$residuals
D = 0.12, p-value = 0.0000008
alternative hypothesis: two-sided

Se observa que el p-value es inferior al 0.05 lo que indica evidencia a favor de la hipotesis alternativa, lo que indica que en nuestro modelo los residuales no se distribuyen de manera normal.

Prueba de homocedasticidad

La prueba de Breusch-Pagan es una prueba de hipótesis utilizada para evaluar la homocedasticidad en un modelo de regresión lineal. Esta prueba se utiliza para determinar si la varianza de los errores en el modelo es constante en toda la gama de valores de las variables independientes.

Planteamiento de Hipótesis:

• \(H_0\): La varianza de los errores es constante en toda la gama de valores de las variables independientes (homocedasticidad).

• \(H_1\): La varianza de los errores no es constante en toda la gama de valores de las variables independientes (heterocedasticidad).

bgtest(modelo_lineal)

    Breusch-Godfrey test for serial correlation of order up to 1

data:  modelo_lineal
LM test = 2.2, df = 1, p-value = 0.1

El P-value es superir al 0.05 lo que indica evidencia a favor de la homosedasticidad, lo que indica que la varianza de los errores es constante.

Prueba de multicolinealidad

El factor de inflación de la varianza para cada variable independiente en el modelo de regresión. Un valor de VIF mayor que 10 es a menudo considerado indicativo de multicolinealidad problemática.

vif(modelo_lineal)

                  GVIF Df GVIF^(1/(2*Df))
Disciplina       1.726  6           1.047
EpocaAño         1.010  1           1.005
HH_Actv_Gnerales 1.831  1           1.353
Otras_HH         1.898  1           1.378

Ninguno de los valores de VIF parece indicar un problema grave de multicolinealidad, ya que todos son menores que 10.

Conclusiones

El modelo de regresión lineal da buenos resultados, sin embargo no cumple el supuesto de normalidad, lo que podria resultar en que el valor de los coeficientes no es real ni confiable.

Planteamiento del modelo

# Ajustar el modelo lineal mixto usando el grupo de disciplina
modeloLMM <- lmer(HH_En_la_Maquina ~ EpocaAño + 
                    HH_Actv_Gnerales + 
                    Otras_HH  + (1 | Disciplina), 
                  
                  data = diario)


tidy(modeloLMM)

# A tibble: 6 × 6
  effect   group      term             estimate std.error statistic
  <chr>    <chr>      <chr>               <dbl>     <dbl>     <dbl>
1 fixed    <NA>       (Intercept)         3.22      0.488      6.60
2 fixed    <NA>       EpocaAñoVerano     -0.351     0.143     -2.45
3 fixed    <NA>       HH_Actv_Gnerales    2.70      0.388      6.95
4 fixed    <NA>       Otras_HH            1.78      0.172     10.4 
5 ran_pars Disciplina sd__(Intercept)     1.18     NA         NA   
6 ran_pars Residual   sd__Observation     1.66     NA         NA   

Al analizar los coeficientes vemos que el sentido es el mismo al aplicar el modelo de regresión lineal. Lo cual puede darnos una idea del verdadero valor de los coeficientes.

Comprobando supuestos

Normalidad

Evaluamos la normalidad de la misma manera en que lo hicimos anteriormente.

Prueba de Kolmogorov-Smirnov: Esta prueba es una opción para evaluar la normalidad de los datos. La función ks.test() en R puede ser utilizada para esto.

Planteamiento de Hipótesis:

• \(H_0\): Los residuales se distribuyen de manera normal.
• \(H_1\): Los residuales NO se distribuyen de maneran normal.

resid(modeloLMM) |> ks.test("pnorm")

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  resid(modeloLMM)
D = 0.12, p-value = 0.0000009
alternative hypothesis: two-sided

Se observa que el p-value es inferior al 0.05 lo que indica evidencia a favor de la hipotesis alternativa, lo que indica que en nuestro modelo los residuales no se distribuyen de manera normal.

Homocedasticidad

plot(modeloLMM)

Al analizar el grafico de los valores predichos frente a los residuales del modelo parecer haber evidencia de homosedasticidad, lo que llevaria a aceptar el supuesto de homocedasticidad. Es decir la varianza de los residuales de nuestro modelo es homogénea.

Multicolinealidad

vif(modeloLMM)

        EpocaAño HH_Actv_Gnerales         Otras_HH 
           1.001            1.297            1.297 

El factor de inflación de la varianza para cada variable independiente en el modelo de regresión. Un valor de VIF mayor que 10 es a menudo considerado indicativo de multicolinealidad problemática.

Planteamiento del modelo

# Ajustar el modelo glmm usando el grupo de disciplina 

## con la distribución poisson y una aproximación de 15
modeloGLMM <- glmer(HH_En_la_Maquina ~ EpocaAño + 
                    HH_Actv_Gnerales + 
                    Otras_HH  + (1 | Disciplina), 
                    data = diario,
                    family = poisson(link = "log"), nAGQ = 15)

Comprobando supuestos

plot(modeloGLMM)

Parece haber evidencia a favor de la hipótesis homocedastica

Metodologia de box jenkins

En esta etapa, se identifican los modelos potenciales para la serie de tiempo. Esto implica identificar si la serie de tiempo es estacionaria o no estacionaria, y seleccionar un modelo ARIMA (Autoregressive Integrated Moving Average) apropiado. La identificación también puede incluir la detección de estacionalidad y tendencias, y la determinación del orden de la diferencia y de los términos autorregresivos y de media móvil.

Identificando si la serie es estacionaria

Para determinar la estacionariedad de una serie temporal, recurrimos a la Prueba de Dickey-Fuller Aumentada (ADF). Esta prueba se encarga de evaluar la presencia de una raíz unitaria en la serie, lo que constituye la hipótesis nula de no estacionariedad. En caso de que se rechace esta hipótesis nula, se establece evidencia de que la serie es estacionaria.

adf.test(diario$HH_En_la_Maquina)

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag   ADF p.value
[1,]   0 -6.85  0.0100
[2,]   1 -4.32  0.0100
[3,]   2 -3.18  0.0100
[4,]   3 -2.43  0.0164
[5,]   4 -1.99  0.0460
[6,]   5 -1.65  0.0959
Type 2: with drift no trend 
     lag    ADF p.value
[1,]   0 -22.29    0.01
[2,]   1 -16.79    0.01
[3,]   2 -13.98    0.01
[4,]   3 -12.01    0.01
[5,]   4 -10.69    0.01
[6,]   5  -9.55    0.01
Type 3: with drift and trend 
     lag    ADF p.value
[1,]   0 -22.27    0.01
[2,]   1 -16.78    0.01
[3,]   2 -13.97    0.01
[4,]   3 -12.00    0.01
[5,]   4 -10.68    0.01
[6,]   5  -9.54    0.01
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

El p-value es inferior al 0.05 lo que muestra evidencia a favor de la hipotesis alternativa, es decir la serie es estacionaria lo que implica que no es necesario aplicar diferencias para nuestro modelo arima. q = 0.

Test de autocorrelación

La función acf calcula los coeficientes de autocorrelación para una serie temporal y los representa en un gráfico de barras o puntos contra el rezago. Esto nos permite identificar patrones de autocorrelación en la serie, lo que puede ser útil para identificar estacionalidad, tendencias y patrones cíclicos en los datos.

acf(diario$HH_En_la_Maquina, lag.max = 20)

Parece mostrar un numero de medias moviles igual a 1, de acuerdo al numero de rezagos significativos. por lo tanto q = 1.

Test de autocorrelación parcial

Mientras que la función de autocorrelación (ACF) muestra la correlación entre una serie temporal y sus valores rezagados, la PACF controla los efectos de los rezagos intermedios y muestra la correlación entre dos puntos separados por un número específico de intervalos de tiempo, pero sin considerar los efectos de los rezagos intermedios.

pacf(diario$HH_En_la_Maquina, lag.max = 20)

Al analizar el grafico de autocorrelación parcial el numero de autoregresivos parece ser igual a 1, p = 1.

Planteando el modelo

De acuerdo a los resultados de aplicar la metodologia de box jenkins definimos los valores para plantear el modelo arima

modelo_arima <- arima(diario$HH_En_la_Maquina, order = c(1, 0, 1))


summary(modelo_arima)

          Length Class  Mode     
coef        3    -none- numeric  
sigma2      1    -none- numeric  
var.coef    9    -none- numeric  
mask        3    -none- logical  
loglik      1    -none- numeric  
aic         1    -none- numeric  
arma        7    -none- numeric  
residuals 546    ts     numeric  
call        3    -none- call     
series      1    -none- character
code        1    -none- numeric  
n.cond      1    -none- numeric  
nobs        1    -none- numeric  
model      10    -none- list