Imputación de valores perdidos y detección de valores atípicos en datos funcionales: una aplicación con datos de PM₁₀

1. Introducción

El manejo y análisis de los datos de calidad de aire con idoneidad facilitan comprobar el cumplimiento de las normas de calidad y evaluar la exposición de la población a los contaminantes del aire [1, 2]. Por ejemplo, la información en el monitorio de particulado generado en las zonas minería es de gran relevancia ambiental y para la salud de los habitantes cercanos a estas zonas por la fuerte evidencia de efectos adversos como la disminución de la función pulmonar, dificultad respiratoria e hipoxemia, cambios adversos en la función autonómica cardíaca [3], [4] y [5].

La mina de carbón del Cerrejón, ubicada en el departamento de La Guajira al norte de Colombia, representa una de las minas más grandes a cielo abierto del mundo, con aporte del 38% en la producción y más del 50% en la exportación, convirtiéndose en la principal compañía de obtención de carbón térmico en el país [6]. Esta capacidad de producción hace que las operaciones mineras a cielo abierto realizadas generen material particulado en cantidades significativas. La masa de partículas menor de 10 μm (PM₁₀) es uno de los contaminantes de referencias más importantes que se utilizan para evaluar la calidad del aire en Colombia.

En Colombia a través de la Resolución 650 de 2010, por la cual se adopta el Protocolo para el Monitoreo y Seguimiento de la Calidad del Aire, establece Sistemas Especiales de Vigilancia de Calidad de Aire (SEVCA) a cualquier población con problemáticas específicas de calidad del aire como minería y alto nivel de industrialización. La Corporación Ambiental de La Guajira (Corpoguajira) cuenta con una red de diez estaciones (10) ubicada en la zona de influencia de la Mina el Cerrejón. Estas estaciones registran mediciones de partícula sedimentable, TSP, PM₁₀ y PM_2.5.

En la practica el análisis de datos generados por particulado en zonas carboníferas de cielo abierto es realizado con métodos estadísticos tradicionales univariado y multivariados, sin tener en cuenta el comportamiento en el tiempo y sin considérala como una trayectoria. La importancia de darle el enfoque funcional o de que las observaciones sean curvas, es que el objeto de estudio (la estación de monitoreo) es considerada como una unidad total en un periodo determinado de tiempo (año). A partir de considerar cada estación como curvas, se realiza el analice y tratamiento que dé a lugar. El monitoreo de la calidad del aire se realiza con el fin de detectar cualquier concentración significativa de contaminante que puede tener posibles efectos adversos para la salud humana. Sin embargo, dicho análisis se complica si existen gran cantidad de observaciones faltantes, que pueden ser originado por diferentes motivos [7]. Los conjuntos de datos incompletos pueden causar sesgo, por lo tanto, es importante encontrar la mejor manera de estimar valores faltantes que garanticen que los datos analizados presenten información de alta calidad [8].

En la literatura del ADF se ha propuesto diversas técnicas para imputar datos, por ejemplo, [8] compara el método de interpolación lineal y el método de remplazar el dato por la media en un estudio de particulados PM₁₀.

Además, en el caso de datos de precipitación y modelos hidrológicos, generalmente, se aplican métodos de regresión o un método de distancia ponderado que son los más comunes para estimar los datos faltantes en esto temas [9]. Mientras en problemas de particulado de aire se aplican métodos como, [10] y [11], donde se emplearon técnicas simples, como la “mean top bottom” para remplazar los valores faltantes en un conjunto de concentraciones de PM₁₀ [11]. Sin embargo, en [11] se encontró que este método presenta buenos resultados sólo cuando el número de los valores faltantes es pequeño. Mientras tanto, [12] aplica el método de imputación del vecino más cercano para estos casos.

Los diferentes métodos de imputación presentan sus ventajas y desventajas en cada escenario de trabajo. Es por esto que distintos procedimientos podrían dar lugar a resultados desiguales.

Una de las maneras de hacer control sobre la generación de particulados de aire en el monitoreo ambiental es a través de técnicas de identificación de valores atípicos, que hasta hace poco se basaban en herramientas de detección univariados o bivariados [13], [14] y [15].

La detección de valores atípicos es un campo muy amplio, que se ha estudiado en diferentes contextos. [16], [17] y [18], además de proporcionar una amplia descripción de las diferentes técnicas, se ha examinado en distintos dominios de datos, incluidos los de alta dimensión [19] y en series de tiempo [ 20], [21].

La detección de datos irregulares es muy popular en estudios de la industria y monitoreo ambiental y, por lo tanto, se han creado algunas herramientas en software para la detección eficiente de estos valores, como el software R (paquetes "outliers’1 y" outlierD " [22]).

En el ADF se ha implementado algunas herramientas para detectar valores extremos, como por ejemplo [23] en estudios sobre calidad del agua y [24] sobre valores atípicos en la calidad del aire.

Adicionalmente, la detección de valores atípicos es otro tema importante en la investigación y análisis de calidad del aire, sobre todo en el contexto funcional. Esto incluyen la detección de atípicos temporales en términos de magnitud e identificación de patrones inusuales (es decir, curvas atípicas de forma), que también proporciona información muy útil a través de sus métodos gráficos [25]. Una visión más global en el caso funcional sobre temas de los datos faltantes se puede encontrar en [26] y [27].

En el ADF existen diverso enfoque para abordar el problema de datos faltantes, entre las técnicas implementadas en la última década están los autores [26], [27]. Además, existe una gran cantidad de literatura que aborda algunos métodos para imputar datos, en el caso multivariado [28] y para datos longitudinales esta [29]. Recientemente, [30] propuso abordar estos temas utilizando el análisis probabilístico de componentes principales (ACPP) y el análisis de componentes principales bayesiano (ACPB) para los algoritmos de imputación.

En cuanto a los métodos de detección de valores atípicos, un punto esencial son las observaciones periféricas. Mientras la detección de datos multivariados atípicos se ha desarrollado a lo largo de varias décadas, pero en observaciones funcionales sólo ha sido discutida en los últimos años [30] y [31]. Es de aclarar que observaciones atípicas pueden conducir de manera adversa a errores de modelación, estimación de parámetros y resultados incorrectos.

Los mecanismos de detección de atípicos en el caso funcional incluyen el uso de análisis de componentes principales robusta [32], parcelas de descomposición de valores singulares [33] y métodos gráficos como el rainbow plot, bagplot y boxplot para datos funcionales [25]. Hace dos décadas el ADF introdujo algunos métodos estadísticos para estudiar este tipo de datos y ha sido ampliamente desarrollados por [25] y [26].

En este trabajo se proponen el análisis de datos funcionales al considerar las concentraciones semanales de PM₁₀ como una trayectoria en el tiempo, este es su primer aporte. El enfoque de imputar datos funcionales es muy resiente y de gran relevancia porque es muy común encontrar observaciones faltantes en estudios de monitoreo ambiental. La metodología propuesta para imputar curvas se basa en determinar la esperanza condicional del ACPF, ver (4) y (5), para obtener los escores generados por este mismo y así, ser utilizado como elementos para imputar datos faltantes en el análisis de concentraciones de PM₁₀ de las diferentes estaciones de muestreó. Para esto se ha utilizado la metodología expuesta por [34] y obtener la muestra completa.

2. Metodología

Es muy común encontrar datos faltantes en procesos de observación, sobre todos cuando se trata de monitoreo ambiental. Existen muchos motivos entre ellas fallas de equipos, problemas con la electricidad, errores humanos entre otros. Esto nos coloca en un primer escenario de dificultades. Incluso, el ADF tampoco escapa a esta situación, sin olvidar que si faltan algunas observaciones no es un problema grave, solamente cuando se encuentra intervalos muy grandes se dificultan recuperar la estructura completa de la curva.

La imputación de datos funcionales ha sido muy poca tratada en la literatura, es por esto, que hemos querido adoptar la metodología expuesta por [34] al caso de imputación da datos de particulado cuando la observación es una trayectoria.

El enfoque se basa en recuperar la estructura de las curvas utilizado la herramienta del ACPF y la expansión de Karhunen–Loève KL, esta última expresión (KL) no es más que la media del proceso más un término de error, la cual es estimada a partir de los escores del ACPF [34].

Inicialmente se describe la ubicación de las estaciones de muestreo y los equipos utilizados para este fin, siguiendo con la metodología de la reconstrucción de las observaciones faltantes para su imputación. Finalmente, se obtiene las muestras completas y a partir de ahí se procede con la búsqueda de curvas atípicas.

2.1. Área de estudio

El Cerrejón es una mina de explotación de carbón a cielo abierto localizada al norte de Colombia en el departamento de La Guajira, frontera con Venezuela. Es una de las minas a cielo abierto más grandes de Latinoamérica y decimo segunda en el mundo. Está ubicada en las coordenadas 11°5’2’’ Norte y 72°40’31’’ Oeste (Figura 1).

Los depósitos de carbón se encuentran en una extensión aproximada de 70.000 hectáreas. La topografía del terreno es representada por una llanura aluvial entre los sistemas montañosos de la Serranía del Perijá por el oriente y la Sierra Nevada de Santa Marta por el occidente.

2.2. Datos de estudios PM₁₀

Se tomaron para el estudio los datos de material particulado PM₁₀, determinados en 7 estaciones que conforman la Red de Calidad de Aire de el Cerrejón. Las concentraciones fueron recolectadas por la empresa entre el 4 de enero y el 29 de diciembre de 2012. Las partículas PM₁₀ se determinaron por el método gravimétrico, utilizando muestreadores de alto de volumen (Hi-Vol) con cabezal de separación inercial para las PM₁₀, según lo establecido en el método de referencia para PM₁₀ de la Agencia Ambiental de los Estados Unidos EPA. Los equipos Hi-Vol operaron cada seis días con tiempo de muestreo de 24 horas, desde las 12:00 h hasta las 12:00 h del día siguiente. Se tomaron 61 muestras en cada estación. La Tabla 1 presenta la ubicación de las estaciones, los equipos utilizados en el monitoreo y los parámetros medidos en cada estación.

2.3. Modelos de valores perdidos

En general, los datos faltantes pueden ser de naturaleza aleatoria y algunas veces son causados por fallas del detector o las herramientas de medición. Hay tres clases de patrones típicos de datos faltantes [35], [36] para los datos multivariados y para datos longitudinales [36]. Primero cuando los datos faltantes son completamente al azar, segundo, cuando los faltantes presentan un patrón determinado y por últimos una combinación de los dos, pero en el caso funcional este enfoque es diferente.

En el ADF se adoptan los métodos anteriormente expuesto, y se describe las siguientes situaciones que podrían ocurrir [34].

1. Los puntos faltantes PF están aislados o dispersos aleatoriamente en la curva observada.

2. Se presentan intervalos faltantes IF, o pequeños grupos de puntos perdidos aleatoriamente a través de la curva.

3. Se presentan tanto puntos faltantes PF como intervalos faltantes IF en la curva observada.

En nuestro caso se presentaron intervalos grandes de observaciones faltantes que dificultan la obtención de la curva completa.

A continuación, se describe la técnica para reconstruir las curvas a través del enfoque de ACPF.

2.3.1. Análisis de Componente principales funcional ACPF

El análisis de componentes principales (ACP) es un enfoque estándar para la exploración de la variabilidad en datos multivariados. El ACP utiliza una descomposición de valor propio de la matriz de covarianza de los datos para encontrar direcciones en el espacio de observaciones a lo largo de las cuales los datos tienen la mayor variabilidad. Para cada componente principal, el análisis produce un vector de peso que proporciona la dirección de la variabilidad correspondiente a ese componente [37].

En el contexto funcional, cada componente principal está especificado por una función de peso o función propia 𝜙(𝑡) definida en el mismo rango de 𝜏=[0,𝑇] que los datos funcionales. De forma similar el ACPF resuelve el problema obteniendo las funciones propias 𝜙_𝑖(𝑡),𝑖=1,…,𝑚 , valores propios 𝜆_𝑘 y la función de covarianza 𝐶𝑜𝑣_𝑋(𝑠,𝑡)=𝑐𝑜𝑣(𝑋(𝑠),𝑋(𝑡))=𝐺(𝑠,𝑡),𝑠,𝑡∈𝜏, que en la práctica es estimada a partir de la muestra de datos funcionales. De igual manera esta técnica permite encontrar los escores de cada curva, que se proyectan como un vector de datos en el espacio reducido [38].

Primero, se adopta la idea de que cada trayectoria de particulado PM₁₀ es una realización de una función aleatoria 𝑋(𝑡). Donde se asume que el proceso estocástico 𝑋(∙) es la función aleatoria de particulado que tiene una curva media suavizada desconocida 𝐸𝑋(𝑡)=𝜇(𝑡) con función de covarianza definida, así 𝑐𝑜𝑣(𝑋(𝑠),𝑋(𝑡))=𝐺(𝑠,𝑡),𝑠,𝑡∈𝜏, donde 𝜏=[0,𝑇] se encuentra en el espacio cuadrado integrable L^². Aquí G tiene una expansión ortogonal en L^². La solución del problema resulta de resolver la siguiente ecuación [39].

así que los {𝜆_𝑘} son el conjunto de valores propios y {𝜙_𝑘} el conjunto de funciones propias que forman un conjunto base ortonormal con norma uno en L^². Entonces cada trayectoria aleatoria de particulado PM₁₀ tiene la siguiente expansión Karhunen - Loéve.

(1)

En la práctica, la función aleatoria 𝑋_𝑖 está contaminada por errores de medida y la i-ésima curva puede expresarse como:

(2)

(3)

Donde los errores de medición aleatorios 𝜀_𝑖𝑗 se asumen no correlacionados entre sí e independientes de 𝜉_𝑖𝑘 con 𝐸(𝜉_𝑖𝑘)=0 y 𝑣𝑎𝑟(𝜉_𝑖𝑘)=𝜎^². Para obtener la función estimada (1), se utiliza la función del modelo de componentes 𝜇 y 𝜙_𝑘 aplicando localmente el método suavizado mínimos cuadrados ponderado con la combinación de todas las trayectorias para estimar la función media. Para obtener la estimación de G, se suaviza las covarianzas empíricas. Seguidamente se obtiene los 𝜆_𝑘 valores propios estimados y los 𝜙_𝑘 aplicando el procedimiento de la descomposición espectral para la estimación de la función de covarianza en el ACPF, similarmente como se procede en el ACP tradicional.

La estimación de coeficiente aleatorio de 𝜉_𝑖𝑘 no se puede obtener fácilmente a través de esta aproximación

Primero, cuando faltan muchas observaciones esta integral de aproximación se dificultad. Segundo, 𝑋_𝑖(𝑡) no se puede observar directamente, sino solo a través 𝑌_𝑖 que están afectadas por los errores de medición, y estimar 𝜉_𝑖𝑘 sustituyendo 𝑌_𝑖 por 𝑋_𝑖 puede conducir a escores sesgados en el ACPF. Para superar estas dificultades, se adopta el enfoque de [38]. Sea Σ_𝑌𝑖 la matriz de covarianza de 𝑌_𝑖 , donde los escores 𝜉_𝑖𝑘 y los errores 𝜀_𝑖𝑘 cumplen conjuntamente el supuesto Gaussiano.

(4)

La estimación condicional de los escores de la ecuación (33) son obtenidos de la siguiente manera:

(5)

Donde, la matriz de covarianza de 𝑌_𝑖 es estimada por y 𝛿_𝑖𝑗 es el delta de kronocker. Para la estimación de valores perdidos se utilizaría la metodología anteriormente detallada.

A continuación, se aplica el criterio de la fracción de la varianza explicada con la regla de los primeros L componentes que explican al menos 𝜏_𝜆 ×100% de la varianza total.

(6)

donde M es el número mayor de componentes con 0 y 𝜏_𝜆 es un valor predeterminado.

Basado en la estimación de un modelo de componente 𝜇̂ ,𝜙̂ y {𝜉_𝑖𝑘}_{𝑘=1,…,𝐿} para todo i, la predicción de la función de la i-ésima curva, está dada por

(7)

Los datos ajustados del modelo (4) se pueden utilizar para imputar una observación faltante por el valorObserve que las trayectorias pronosticadas incluyen las componentes de la función media suavizada y una combinación lineal de las funciones propias, que recuperar las trayectorias individuales de las mediciones del ruido, mientras que el error depende del número de componentes seleccionados en el modelo ACPF.

A partir de la muestra de observaciones completas, se continua con la búsqueda de valores atípicos que podría entregar algún indicio que las concentraciones de PM₁₀ sobrepasan los estándares establecidos por la ley ambiental o presentar una alarma al respecto.

A continuación, se describe como se implementó estas herramientas gráficas para la búsqueda de atípicos funcionales en nuestro estudio.

2.3.2. Herramienta gráficas para detectar atípicos funcionales

Las herramientas gráficas para obtener valores atípicos en el caso funcional son de mucha utilidad. Entre los más utilizados están las gráficas bagplot funcional y el HDR boxplot funcional, propuesta por [25]. Ambos métodos se basan en tomar los dos primeros escores de los componentes principales del ACPF. Para este propósito, [36] aplica un algoritmo robusto que permite intuir que las proyecciones desde la descomposición de componentes principales deben ser sensibles a valores atípicos. Sin embargo, este método no considera valores perdidos y se requiere un conjunto de datos completos para proyectar los datos en un espacio de dimensión menor de tal manera que una medida robusta de la varianza de los datos proyectados sea máxima.

Como primer paso se imputaron las trayectorias de PM₁₀ que presentaban intervalos de faltantes muy grandes, implementando la metodología desarrollada por [34]. Con la muestra completa de datos se aplicó los métodos gráficos HDR boxplot y bagplot funcional con la idea de encontrar algunos valores atípicos.

La herramienta grafica consiste en aplicar ACPF a las curvas reconstruidas, que en este caso corresponde a las trayectorias de las diferentes estaciones de monitoreo de PM₁₀, luego con los dos primeros escores del ACPF se construyen regiones bivariadas. Los métodos pretenden construir una bolsa que contenga el 95% de los escores y los puntos que este por fuera se ponderan como valores atípicos.

2.3.3. Detección de valores atípicos funcionales

Para, llevar a cabo el procedimiento en la detección de valores atípicos en los gráficos funcionales, se necesita pre especificar una probabilidad de cobertura para determinar la región periférica. La probabilidad de cobertura generalmente suele ser de 99%, 95% y 90%, correspondiente a un nivel de significancia de 𝛼=0.01,0.05 𝑦 0.10.

3. Resultados

Los resultados de la Figura 2 no evidencia diferencia significativa entre los promedios de las siete (7) estaciones. Se observa mayores concentraciones en la estación RD7 (Casitas) y la menor concentración en la estación RD1 (Sol y Sombra), la cual está ubicada viento arriba de la zona de explotación. El rango de las concentraciones de PM₁₀ para el periodo de estudio varió entre 2,00 μg/m³ y 87 μg/m³. La estación RD7 reportó el mayor promedio en las concentraciones de PM₁₀ con un nivel de 41 μg/m³ y RD1 presentó el menor promedio con un nivel de 17 μg/m³.

Los resultados muestran que el material particulado PM₁₀ son inferiores a los estándares nacionales y de la EPA, a pesar de estar en una zona de influencia minera. Un reciente proyecto de la EPA sugiere nuevo estándar de 24 horas para PM₁₀ en el rango de 65 - 85 μg/m³. Aun para los niveles proyectados se estarían cumpliendo estos preceptos.

En la Figura 3 y Figura 4 ilustra el HDR boxplot y bagplot funcional, donde la estrella negra corresponde a la curva mediana, mientras la región gris oscura pertenece al 50% de las curvas y la región gris claro corresponde al 95% de ellas. Mientras que, el punto 10 (Figura 3) representa a la estación de monitoreo (Provicial) que se muestra como una trayectoria atípica a un nivel de significancia de 𝛼=0.05. Esta curva presenta los niveles más bajos de concentración a lo largo del tiempo de estudio, lo que en este caso es muy positivo. Mientras que la Figura 4 muestra el bagplot funcional, donde la estrella representa la curva mediana y la región gris clara corresponde al 95% de las curvas. Este grafico muestra como atípica la estación de monitoreo Barranca, debido a que presenta los niveles de concentración más alta a lo largo del tiempo de estudio.

La estación Barrancas se presenta como atípica dado que las variaciones de concentraciones son las más altas, con respecto a la curva mediana. A pesar de estar en el rango permitidos por la norma ambiental. Aunque esta curva es atípica es de forma muestra sus picos más altos cerca de los meses de septiembre y octubre.

Los resultados muestran que la estación (Barrancas) es una curva atípica, lo que indica que se debe proponer estudios más profundos y monitoreo más seguido en el tiempo sobre esta zona, que permita una evaluación constante.

En diferentes sitios de muestreo los resultados de las concentraciones de PM₁₀ no excedieron los estándares diarios de la Norma Colombiana y la National Ambient Air Quality Standards (NAAQS) de los Estados Unidos.

La región gris oscuro y gris claro de la Figura 3 muestran el 50% de HDR y el HDR externo, respectivamente. La línea negra es la curva modal. Las curvas fuera de la región exterior son valores atípicos (Provincial).

Las regiones oscuras y gris claro en la figura 4 muestran la bolsa de regiones. El asterisco es la mediana de profundidad de Tukey. A la izquierda, la línea negra es la curva mediana, rodeada por intervalos de confianza puntuales del 95%. Las curvas fuera de la región exterior se muestran como valores atípicos. (estación Barrancas).

4. Conclusiones

Este artículo presenta un estudio de caso donde se implementa la imputación de curvas en el caso funcional y de detección de atípico en el mismo contexto.

Observaciones faltantes en estudios de monitoreo ambiental son más comunes de lo que se cree, en el contexto funcional no escapa a este problema, por esto se presentó una metodología que permite reconstruir trayectorias y para ilustrar esta situación, consideramos el caso de trayectorias de concentración PM₁₀ a la periferia de una zona de mina de carbón a cielo abierto.

Una de las maneras para evaluar y controlar la calidad del aire es a través de detección de valores atípicos. En el contexto funcional se ha implementado últimamente técnicas graficas que permite obtener curvas atípicas tanto de escala como de forma. Es así que estas curvas pueden suministrar valiosa información para su intervención preventiva en estudios de monitoreo ambiental.

El resultado evidencia que la estación Barrancas como una curva atípica, mostrando niveles de concentración PM₁₀ más altas que el promedio, esto es un indicativo para tomar acciones sobre un monitoreo más seguidos y más exhaustivos sobre la calidad del aire. Aunque la curva 10 (Provincial) se presenta como una curva atípica de forma, no debe ser tenida en cuenta para su intervención en un monitoreo más detallado porque presenta bajas concentración de particulado, como es lo deseado.

Nuestros resultados numéricos indican que el enfoque de datos funcionales para la imputación de valores perdidos y la detección de valores atípicos puede aplicarse y ajustarse muy bien a este tipo de estudios.

Agradecimientos

Los autores expresan sus agradecimientos a la Corporación Autónoma Regional de La Guajira por el suministro de datos de la zona de estudio. A Colciencias por el financiamiento y apoyo del proyecto "Aplicación de los modelos de receptor en el aporte de fuentes a la contaminación del aire por material particulado" Código: 111552430465l.

Referencias

Información adicional

Cómo citar: R. Meléndez, S. Bolívar, R. Rojano, “Imputación de valores perdidos y detección de valores atípicos en datos funcionales: una aplicación con datos de PM₁₀,” Rev. UIS Ing., vol. 19, no. 2, pp. 1-10, 2020. doi: 10.18273/revuin.v19n2-2020001

Revista UIS Ingenierías
ISSN: 1657-4583
Vol. 19
Num. 2
Año. 2020

Imputación de valores perdidos y detección de valores atípicos en datos funcionales: una aplicación con datos de PM10

RafaelStevensonRoberto MeléndezBolívarRojano
Universidad de La GuajiraPontificia Universidad JaverianaUniversidad de La Guajira,ColombiaColombiaColombia