Incertidumbre De Medición Estimada Por Modelos RSM Para Cumplir ISO-17025 Y Útil En Validación De Métodos Analíticos
Víctor Soto Del Toro, Químico Industrial, Ced. Prof. 4405730
2024-04-05
¿Cómo funciona la relación entre la incertidumbre de medición y las regresiones multivariables?
¿Y cómo nos ayuda en protocolos de aseguramiento de la calidad en el laboratorio?
Este reporte es extenso porque contiene detalles exactos.
Está dirigido a mis colegas de habla hispana que ejercen la química y trabajan con datos, pero que además requieren explorar o incluso, profundizar en la estadística multivariable aplicable para propósitos de control interno en laboratorios analíticos y de pruebas.
En este artículo veremos de dónde salen los datos, cómo se procesan y cómo se validan, para producir gráficos tridimensionales con propósitos predictivos mediante modelos RSM.
Una vez que elijamos un modelo RSM válido, la incertidumbre de medición asociada con el modelo se obtiene automáticamente para incluirlo en el presupuesto de incertidumbre. Y como es un asunto técnico que exige explicar de dónde sale cada cosa, la exposición toma tiempo. Pero la lectura es ágil.
Sin embargo, estos son conceptos avanzados y por eso esta presentación se ha organizado en tres partes:
Antecedentes
Caso tradicional RSM: regresión lineal múltiple
Caso avanzado RSM: regresión por componentes principales
Si el lector tiene interés en el tema, pero sus aptitudes corresponden a etapas iniciales, le conviene explorar la primera parte.
Si el lector ya tiene conocimientos de regresión simple por mínimos cuadrados, pero el tema de regresión múltiple le resulta nuevo, vaya a la segunda parte.
Si el lector está interesado en explorar uno de los métodos multivariables más utilizados en algoritmos de Machine Learning, PCA, vaya a la tercera parte.
Y si el lector prefiere una visión integral del proceso, lo invitamos a leer el artículo completo.
Parte 1, Antecedentes: Orígenes y Naturaleza de Datos Químicos Multivariables
Antes de comenzar, me gustaría primero ofrecer un poco de contexto.
Aún recuerdo cuando mi profesor de análisis instrumental, el Dr. Ernesto Vélez nos dijo en clase:
“Como químicos, ustedes primero deben aprender a medir. Ya después, hagan lo que quieran.”
En el trabajo científico cuantitativo, la exactitud y precisión en las mediciones instrumentales constituyen la norma. ¿Quién lo dice?
CODATA es el acrónimo de la Comisión de Datos para Ciencia y Tecnología y es la institución internacional que se encarga de recabar y verificar todos los resultados experimentales y de cómputo, de los valores aceptados para las constantes físicas universales y sus incertidumbres.
https://codata.org/about-codata/our-mission/
Ellos son la autoridad en materia de mediciones y precisión científica y estas son algunas de las constantes fundamentales:
Número de Avogadro, constante de Boltzmann, masa del eletrón, masa del neutrón, masa del protón, carga del electrón, constante de Newton, velocidad de la luz en el vacío y muchas otras constantes que describen las leyes de la Física.
¿Y lo que CODATA haga por qué nos debería interesar?
Por que es la institución encargada de promover la colaboración para el progreso de la Ciencia Abierta y mejorar la disponibilidad y usabilidad de datos para todas las áreas de investigación.
Por ejemplo, según CODATA y NIST, las constantes fundamentales terminan en una incertidumbre de dos dígitos, para todas las mediciones.
Eso es lo que los profesionales hacen y tú también deberías adoptar sus reglas.
¿Y cómo es que instituciones como CODATA, NIST y otras formalizan sus criterios para la expresión de incertidumbres?
Por un parte, la respuesta está es una disciplina científica: Análisis del Error.
Por otra parte, decimos que “el error se propaga en un sistema”.
¿Cuál sistema? En Análisis de Error nos referimos al Sistema de Medición.
¿Y qué es un sistema de medición?
En general, un sistema de medición es un arreglo múltiple conformado por instrumentos, sensores, estándares, herramientas de procesamiento de datos y elementos humanos.
¿Y qué es lo que hace un sistema de medición? ¿Para qué sirve?
El sistema de medición se diseña para que, de manera exacta y confiable, nos permita cuantificar cantidades físicas o atributos. Si trabajas en un laboratorio de pruebas, o estás a cargo de él, tú eres parte importante del sistema de medición.
Cada componente en el sistema de medición desempeña un papel crítico, que asegura la integridad y exactitud del proceso y actividades de medición.
¿Con qué beneficio?
La respuesta a esta pregunta contiene varias partes:
Primero y de manera muy general, podemos decir que el principal problema que existe dentro de un sistema de medición es que, mientras las mediciones se realicen y los datos sean procesados dentro del sistema, con cada operación individual el error se producirá, se acumulará y se propagará, impactando a los resultados finales de las mediciones.
Segundo, esto se debe simplemente al efecto natural de la entropía contenida en la información. En otras palabras, nuestros datos de medición exhiben variación porque contienen entropía. Esta es una ley universal y debemos atenderla.
Tercero, entender cómo es que estos errores se propagan y contribuyen individualmente a la incertidumbre final total, es esencial para evaluar la confiabilidad y exactitud de resultados de medición.
Cuarto, identificar las fuentes de error dentro de un sistema de medición y cuantificar sus niveles de impacto en la medición final, les permite a ingenieros y científicos tomar medidas correctivas para minimizar o mitigar estos errores.
Quinto, así es como se mejoran la calidad y confiabilidad de las mediciones obtenidas.
En resumen:
El estudio de la Propagación del Error es una parte integral del Análisis del Error, porque nos ofrece hallazgos valiosos asociados con nuestras mediciones, cálculos, modelos, estimaciones y predicciones.
¿Con qué objetivo?
El propósito es descubrir las magnitudes y orígenes de las fuentes individuales de error, porque así es como nos informanos lo más posible, antes de tomar decisiones.
¿Por dónde comenzar?
Pienso que actualmente, existe la necesidad de una mayor discusión acerca del Análisis del Error con el propósito de satisfacer los requerimientos de la Metrología.
Naturalmente, los vehículos que lo permiten son la Probabilidad y la Estadística.
¿Por qué?
Porque vistos desde la óptica estocástica, para un profesional de la estadística los valores numéricos de nuestras mediciones, son datos.
Pero esta óptica, aunque es válida, no está completa. ¿Por qué?
Porque vistos desde la óptica determinística, para un profesional de la ciencia los valores numéricos de nuestras mediciones, son observaciones.
Pero esta óptica, aunque es también es válida, tampoco está completa. ¿Por qué?
La respuesta a esta última pregunta está en la naturaleza misma del acto de medición.
Decimos “el acto de la medición”, porque una medición no es un número, sino una serie de actividades que se traducen en un número. ¿Y cuál es el efecto que tales actividades ejercen sobre esta cifra?
El efecto es la introducción de FLUCTUACIONES:
En un momento practicamos una medición sobre un objeto y obtenemos un valor. En otro momento practicamos otra medición sobre otro objeto equivalente, que en teoría, debería arrojar el mismo valor, pero eso no es lo que sucede, sino que obtenemos otro valor, un poco diferente, quizá mayor o menor. Y así sucesivamente. Son justamente esas diferencias las que introducen fluctuaciones en las cifras de medición.
¿Por qué sucede esto?
A estas fluctuaciones numéricas la Estadística les llama “variación”.
Por su parte, la Metrología les llama “incertidumbre”.
Pero la Física les llama “error”.
Mientras la Estadística y la Metrología son disciplinas altamente sofisticadas, perfectamente capaces de representar y cuantificar fluctuaciones, no explican su origen.
La única ciencia que sí lo puede hacer, es la Física. ¿Por qué? ¿Quién lo dice?
El mismísimo George. E.P. Box, el gran pionero de la Estadística aplicada a la industria, en la sección 15.1, página 600 de su ilustre libro Análisis de Series Temporales: Pronóstico y Control, él escribe y cito textualmente:
“El control de procesos es no menos que un intento por cancelar el efecto de una ley física fundamental: la Segunda Ley de la Termodinámica, lo cual implica que si lo dejamos por sí solo, la entropía o desorganización de cualquier sistema, nunca puede disminuir y usualmente, se incrementará.”
Ahí lo tenemos.
Así y con base en la definición del iluste Maestro Box (y otros más como el Dr. Harald Martens y Dr. Svante Wold, fundadores de la Quimiometría, que es la ciencia que estudia datos químicos por análisis multivariable), el autor de este artículo comparte la misma firme convicción:
Para hacer un buen trabajo de medición instrumental que refleje el mayor nivel de profesionalismo y produzca información analítica de la mayor calidad posible, el especialista requiere apoyarse sobre una estructura científica conformada por cuatro pilares cuantitativos que, en conjunto y de manera integrada, analizan y describen las fluctuaciones numéricas observadas durante los procesos de medición.
Estos cuatro pilares ya los conocemos y son:
Física, Matemáticas, Estadística y Metrología.
Eso suena bien, pero la pregunta es: ¿cómo las vamos a utilizar? ¿cómo funciona?
Lo que estamos buscando es ENTENDER la fluctuación numérica durantre el proceso de medición. Los datos cambian y nuestra misión es doble: descubrir y explicar por qué cambian.
Para satisfacer este propósito tecnológico tan específico, cada una de estas cuatro disciplinas nos va a informar una cosa por separado, acerca de la fluctuación durante la medición:
Su origen Determinístico lo explica la Física por medio de la Entropía.
Su conducta Matemática la explica la disciplina de Análisis del Error.
Su naturaleza Estocástica la explica la Estadísica.
Su aplicabilidad para propósitos mercantiles, legales y en toma de decisiones la define y se gobierna por medio de la Metrología.
Entonces, para hacer un buen trabajo cuantitativo que rinda frutos tangibles para atender objetivos metrológicos, todo lo anterior le exige al tecnólogo adaptarse intelectualmente a dos paradigmas, en simultáneo:
Determinismo y estocasticidad.
Científicos e ingenieros hemos sido educados bajo el paradigma determinista.
Estadísticos y actuarios han sido educados bajo el paradigma estocástico.
Sin embargo, un profesional que trabaje en la toma de decisiones basadas en datos obtenidos a partir de mediciones instrumentales sujetas a fluctuación, debe obligatoriamente poseer suficiente dominio, no sólo de los procesos determinísticos (físicos), sino de los procesos estocásticos (estadísticos).
¿Por qué?
Porque el resultado numérico que surge de la medición y que lo leemos en la pantalla de nuestro instrumento calibrado, está hecho de dos componentes:
El componente físico, porque la materia y la energía son los objetos sometidos a medición para obtener observaciones (lecturas). Este es el proceso determinístico.
El componente estadístico, porque es la entropía del sistema de medición lo que produce la variación en los datos (desviaciones). Este es el proceso estocástico.
¿Y estos conceptos son en verdad importantes? ¿O son sólo curiosidades numéricas de tipo opcional? ¿Existirá alguna implicación seria?
En México, el referente es la Ley Federal Sobre Metrología y Normalización, que en su artículo 3o fracción II a la letra dice:
“Calibración: el conjunto de operaciones que tiene por finalidad determinar los errores de un instrumento para medir y, de ser necesario, otras características metrológicas.”
Y ese es tan sólo el contexto nacional. En el internacional existen muchas normas.
Como vemos, la ciencia de las mediciones rige las economías de las naciones y regiones enteras. Su incumplimiento o mal uso se traduce en pérdidas cuantiosas.
Las personas se quedan sin empleo. Los negocios se van a la ruina. Especialistas pieden su prestigio y reputación. Juicios en los tribunales por escándalos de fraude.
Así que no, la variación, el error y la incertidumbre en mediciones no son una simple curiosidad matemática opcional.
Ahora que ya conocemos un poco del origen, relevancia e impacto de la fluctuación en las mediciones, sólo nos queda preguntarnos:
¿Qué es lo que vamos a hacer?
Existen cuatro grandes tareas y cada una de ellas es inmensa:
Predicción
Clasificación
Control
Pronóstico
En toda forma de actividad humana oganizada siempre existe la necesidad de predecir, clasificar, controlar y pronosticar lo que va a ocurrir en el futuro.
Proyectos, servicios, negocios e industrias enteras giran en torno a ellas.
¿Por qué?
Porque las formas más inteligentes e informadas de tomar decisiones dependen directamente de estos cuatro mecanismos y que por fortuna, todos ellos son cuantificables físicamente, analizables estadísticamente y programables computacionalmente.
Con la llegada de la relativamente nueva disciplina conocida como Ciencia de Datos, la pesada tarea de tomar decisiones en medio de procesos caracterizados por una alta complejidad debida a las numerosas variables involucradas, la toma de decisiones ahora es posible facilitarla, agilizarla y hasta cierto punto, automatizarla.
La ciencia de datos atiende procesos multivariables de manera eficaz y ofrece muchos métodos que reciben el nombre popular de Algoritmos de Aprendizaje Automático, conocidos en inglés como Machine Learning.
En el caso específico de mediciones de datos químicos multivariables en la industria de transformación de productos procesados y en tareas de análisis de muestras en el laboratorio químico, la Ciencia de Datos recibe el nombre de Quimiometría.
Para atender los intereses de la industria de productos procesados, la Quimiometría opera mediante dos mecanismos con dos objetivos:
Control estadístico de procesos multivariables, mediante la T2 Hotelling y Q.
Calibración multivariable, mediante regresiones múltiples y métodos RSM.
Y es aquí en donde las personas con orientación determinística comenzamos a hallarles más sentido y significado a diversos términos de la Ciencia de Datos.
Por ejemplo, cuando los científicos de datos dicen que van a “entrenar un algoritmo”, de lo que están hablando es sencillamente de una calibración. A final de cuentas, si lo que en la industria queremos es tomar decisiones informadas a partir de datos de alta calidad, necesariamente nuestros sistemas de medición deben calibrarse.
Pero los partidarios de la ciencia de datos van a decir que “entrenear un algoritmo” no es una calibración y así por el estilo. Pero científicos de las ciencias naturales e ingenieros en tecnologías rápidamente se dan cuenta de las similitudes y así, esto es materia de debate porque muchos especialistas de gran experiencia en la industria de proceso, saben desde hace décadas que los procesos multivariables comenzaron a analizarse hace más de 100 años con Pearson y Hotelling con el método clásico de Análisis por Componentes Principales (que hoy se considera uno de los principales algoritmos de machine learning).
Es cierto que estos y otros métodos de análisis multivariable ya eran accesibles desde los 80’s cuando las computadoras se volvieron accesibles, pero solamente especialistas conocían acerca de ellos y podían aplicarlos.
No fue hasta mediados de la década los 2000’s, con el advenimiento del software y código de fuente abierta, que las técnicas de análisis multivariable comenzaron a popularizarse y democratizarse y circular con gran velocidad entre las comunidades científicas.
Y ahora, 20 años más tarde, este autor ve con mucho gusto y celebra que el proceso de democratización de las metodologías multivariables no sólo continúa, sino que está ganando gran tracción entre los más jóvenes, aunque ahora bajo el estandarte de la popular Ciencia de Datos: ahora resulta que las matemáticas aplicadas están de moda. ¿Quién lo diría? Es un triunfo.
Fue justo en el año 2005 cuando este autor estaba actualizándose en estos asuntos que, al explorar la obra literaria del Profesor George E.P Box, “Statistics For Experimenters”, en el pie de página del prefacio de la segunda edición, textualmente dice:
“Todos los cómputos en este libro pueden hacerse con el lenguaje estadístico R”.
En esa época el lenguaje R aunque era incipiente, lo cierto es que ya era de los más competitivos contra paqueterías estadísticas comerciales. Esto fue así porque R era un producto creado y alimentado y mantenido por una comunidad de voluntarios internacionales. R es la “lingua franca” de los científicos que practican la estocástica.
En ese tiempo, este autor trabajaba para el Gobierno Mexicano con un equipo de compañeros en proyectos de diseño y desarrollo de programas de protección al medio ambiente y prevención de la contaminación, con tecnologías actuales en ese entonces.
Una de esas tecnologías fue el análisis estadístico multivariable a partir de datos obtenidos de muestras analizadas en laboratorio. Y todo el cómputo estadístico fue con el lenguaje de programación R.
En principio, la misión científica era simple:
Explorar y entender las relaciones cuentitativas empíricas entre diversos factores y contaminantes y sus impactos adversos sobre la calidad del entorno ambiental.
Expresado en términos determinísticos, el objetivo del proyecto era encontrar ecuaciones matemáticas que nos permitieran entender las correlaciones entre factores considerados como causas y contaminantes considerados como efectos, en diversos puntos de la ciudad.
En concreto, lo que buscábamos era estimar la conducta de los niveles de ozono troposférico y comparar contra diversos factores como monóxido de carbono, dióxido de nitrógeno, temperatura, humedad relativa, radiación solar y otros.
En la literatura existen numerosos estudios y modelos químico-matemáticos que describen los fenómenos de contaminación ambiental, pero debido a que los datos de medio ambiente exhiben elevados niveles de variabilidad, las estimaciones predictivas son de lo más desafiantes.
¿Por qué?
Porque a diferencia de los diseños experimentales multivariables en industria, donde las condiciones son controladas y en la mayoría de los casos los niveles de cada variable pueden definirse y modificarse, con los monitoreos de control de medio ambiente eso no sucede.
Los estudios de estimación estadística con base en datos de medio ambiente se parecen más a los estudios de carácter observacional, donde las variables exhiben sus conductas de acuerdo a las condiciones prevalecientes en el entorno.
En otras palabras: estamos observando a la naturaleza de la atmósfera respirable, mientrás está siendo impactada por agentes extraños, en plena acción.
Bajo esas condiciones, la investigación experimental es imposible.
Sólo nos queda tomar muestras, obtener lecturas de la mayor calidad posible y analizar los datos con el rigor que la investigación exige.
Y aunque todo esto suena muy interesante, la pregunta medular es: ¿para qué?
¿Cuáles son los propósitos y beneficios de estimar regresiones, correlaciones, erorres estándar, coeficientes de determinación, significancias, validación de datos y seguir al pie de la letra toda una serie de protocolos cuantitativos?
La toma decisiones, naturalmente.
En el contexto de vigilancia y control de medio ambiente, las decisiones a considerar son la gestión y evaluación de riesgos, desarrollo de politícas públicas, predicción de tendencias y la eventual ejecución de mecanismos de intervención.
Desde luego, la toma de decisiones es un acto humano. Pero tomar decisiones de interés público a partir de intuición y corazonadas, sin contar con datos ni tomando en serio toda evidencia cuantitativa, no es solo anti-profesional sino que conducirá al desperdicio, el despilfarro y hasta el desastre.
Y así, las decisiones más importantes deberían gozar de un firme apalancamiento al apoyarse en datos validados, porque de lo contrario, las decisiones sin datos no son decisiones: son adivinanzas.
Para ilustrar cómo se procesan datos de mediciones reales por métodos RSM para toma de decisiones en control de medio ambiente en el sector público, vamos a trabajar con los datos de un estudio realizado para monitoreo de calidad del aire.
El ozono es un producto de la descomposición fotoquímica del dióxido de nitrógeno. Es decir: cuando el dióxido de nitrógeno generado por combustibles fósiles entra en contacto con la luz solar, ocurre una reacción que produce ozono.
Este es un fenómeno sumamente dañino para la salud de las personas: el ozono debería estar en las zonas más elevadas de la atmósfera, filtrando la fracción ultavioleta, pero nunca en la biósfera, porque inhalar ozono con sus radicales libres produce cáncer.
Eso así como este asunto, además de ser un fenómeno de medio ambiente, al mismo tiempo es un tema de salud pública. Por eso en tareas de predicción y control es de crucial importancia explorar y conocer las conductas de las variables predictoras, con el mayor grado de detalle cuantitativo posible.
De otro modo, estamos a ciegas. Ignorantes de lo que pasa a nuestro alrededor. Literalmente.
Estudiaremos las relaciones que guardan la radiación solar y los niveles de dióxido de nitrógeno, en la generación de gas ozono a nivel del suelo.
Los objetivos de este ejercicio son los siguientes:
Exponer las formalidades científicas, reglas matemáticas y protocolos estadísticos aplicados durante el tratamiento de datos multivariables sujetos a regresión múltiple.
Detectar inconsistencias que producen resultados espurios.
Explorar opciones y mecanismos que producen resultados válidos.
Conocer cómo se construye un modelo RSM.
Identificar cómo obtener la incertidumbre de medición a partir de un modelo RSM.
Identificar el mayor problema en datos químicos: la MULTICOLINEALIDAD.
Hago una aclaración: este artículo no pretende en ninguna forma sustituir el aprendizaje formal de los diversos métodos analíticos de regresión múltiple y de análisis multivariable. Desarrollarse en esas áreas requiere años de preparación.
Esta presentación es una muestra de un caso real donde se podrán apreciar las ventajas y también las dificultades a las que el científico se enfrentará, cuando su tarea consista en desarrollar modelos empíricos útiles para toma de decisiones.
Es el deber del científico explicarles a los involucrados, en los términos más claros y suscintos posibles, el significado tangible de sus hallazgos numéricos.
Si el lector desea profundizar y obtener mayor dominio del tema, se le invita a explorar las numerosas obras literarias publicadas o matricularse en cursos y talleres especializados.
Vayamos a la segunda parte de este artículo.
PARTE 2, Caso Tradicional RSM: Regresión Lineal Múltiple
En el ejercicio del método científico, después de las etapas de observación, hipótesis y experimentación/mediciones, viene la etapa de análisis de datos.
Supongamos que toda la labor de las primeras tres ya se realizó y estamos listos para entrar en el análisis.
Por motivos de discreción, sólo analizaremos 10 observaciones para tres factores:
Ozono, radiación solar y dióxido de nitrógeno.
Las cifras se obtivieron por mecanismos instrumentales colocados en sitio. Los niveles de ozono se miden en unidades Dobson, la radiación solar en Watts/m2 y el dióxido de nitrógeno en partes por billón.
Una típica consulta en R luce así.
contaminacion.atmosfera[, c(2, 3, 4)]## Ozono (Y) Radiacion Solar (X1) NO2 (X2)
## [1,] 7 75 8
## [2,] 12 81 14
## [3,] 4 45 6
## [4,] 21 82 11
## [5,] 7 72 8
## [6,] 12 77 12
## [7,] 11 65 7
## [8,] 9 61 8
## [9,] 3 47 6
## [10,] 8 72 11Atendiendo a los requisitos del método RSM, vamos primero a analizarlos de manera formal con un poco de álgebra lineal y enseguida, con comandos y código.
Para empezar, el ozono es nuestra variable de respuesta y la representaremos con un vector Y, mientras que la radiación y el dóxido de nitrógeno son nuestras variables predictoras, X1 y X1.
¿Por qué las definimos así? ¿No suena esto arbitrario?
De ningún modo. El motivo por el cual el ozono es la variable de respuesta y la radiación y dióxido de nitrógeno las variables predictoras es porque, físico-químicamente, el ozono es un producto de la reacción foto-química entre la luz solar y el dióxido de nitrógeno.
En otras palabras:
Químicamente, la luz solar y el dóxido de nitrógeno son nuestros reactivos y el ozono es el producto.
Luego entonces
Matemática y estadísticamente, la luz solar y el dióxido de nitrógeno son nuestras variables predictoras y el ozono nuestra variable de respuesta.
Este arreglo descriptivo es válido, tanto físico-química como matemático-estadísticamente.
\[\quad Ozono \quad\quad\quad\quad\quad\quad\quad X_1 \quad X_2\]
\[ Y = \begin{bmatrix} 7 \\ 12 \\ 4 \\ 21 \\ 7 \\ 12 \\ 11 \\ 9 \\ 3 \\ 8 \\ \end{bmatrix} \quad\quad\quad\quad X = \begin{bmatrix} 75 & 8 \\ 81 & 14 \\ 44 & 6 \\ 82 & 11 \\ 72 & 8 \\ 77 & 12 \\ 65 & 7 \\ 61 & 8 \\ 47 & 6 \\ 72 & 11 \\ \end{bmatrix} \]
Con estas dos matrices es suficiente para obtener un modelo lineal de predicción.
Sin embargo, muchos científicos e ingenieros experimentados tienen la idea de que los modelos lineales no son los suficientemente descriptivos y no capturan la realidad.
Por eso se van directo a los modelos RSM, incluyendo los términos cuadráticos y de interacción en la ecuación, desde el inicio de su análisis, saltándose la etapa de orden lineal. ¿Por qué lo hacen así?
El razonamiento proviene de una regla experimental no escrita, que dice que los procesos de la naturaleza siguen patrones curvilíneos. Matemáticamente son cuadráticos. Y tienen razón: muchísimos procesos naturales son cuadráticos. Tal pareciera que las conductas de los fenómenos de la Naturaleza, el Universo y el Microcosmos pueden describirse matemáticamente con expresiones curvilíneas de segundo orden. (En lo personal, yo nunca he visto ni encontrado datos de fenómenos físico-químicos que sigan patrones cúbicos o de mayor orden. En todo caso, caerían en sobre-ajustes.)
Desde luego, los modelos cuadráticos no logran capturarlo todo y a veces, producen predicciones espurias. De hecho, esto es justo lo que sucede con nuestros datos atmosféricos. Y para demostrarlo, vamos a procesarlos obedeciendo a esta regla empírica cuadrática y veamos lo que pasa.
Estimación de Coeficientes en un Modelo Clásico de Regresión Múltiple Cuadrática
En esta primera iteración, nuestro objetivo es encontrar la ecuación cuadrática que describa lo mejor posible a nuestro conjunto de datos.
¿A qué nos referimos exactamente cuando decimos “lo mejor posible”?
Significa hallar aquellos coeficientes que, al multiplicarlos con datos nuevos de niveles de radiación solar y dióxido de nitrógeno, produzcan estimaciones aproximadas de los niveles reales de ozono. A esto le llamamos predicción.
Una de las primeras figuras de mérito que usamos, es la suma de distancias cuadradas.
Me estoy adelantando un poco, pero una de las cosas que veremos, es que los datos de nuestro conjunto original (que para efectos prácticos en este caso, podemos considerarlos como datos de calibración), gráficamente los conoceremos como distancias verticales. Si nuestro modelo final es válido, estas distancias deberían ser lo más cortas posible, estimadas desde los puntos de nuestros datos originales, hacia la línea de regresión.
Estas distancias verticales se abrevian como SSE y reciben el nombre de suma de errores cuadrados. La ecuación es:
\[SSE \; = \;\Sigma_i(Y_i \;- \hat{Y_i})^2\]
Como estas distancias cuadradas deben ser lo más pequeñas posible, es por eso que el procedimiento recibe el nombre de mínimos cuadrados.
Esa minimización de distancias se refleja en los valores SSE y proviene de los coeficientes de la ecuación de regresión. Esa ecuación es la siguiente:
\[\hat{B} \;= \;[ \;(X^TX)^{-1} \;] \;Y\]
Esta es una ecuación universal. ¿Qué es lo que vemos en ella?
Primero, la \(Y\): este es el vector de nuestras mediciones de ozono.
Segundo, la \(X\): es la matriz de datos de radiación solar y dióxido de nitrógeno.
Tercero, la \(\hat{B}\): es la matriz de coeficientes que vamos a estimar.
Además vemos un par de exponentes; estos describen operaciones especiales en álgebra lineal.
La ecuación anterior es simple, pero muy poderosa. Formalmente, es ecuación matricial que define las estimaciones mínimas cuadradas de los coeficientes de regresión.
Para ver cómo se aplica en un modelo RSM, primero es necesario definir la estructura de la ecuacuón cuadrática que queremos construir. A este proceso se le llama “ajuste”.
Lo que haremos será ajustar una ecuación cuadrática completa, a partir de nuestras dos variables predictoras, Radiación Solar y Dióxido de Nitrógeno.
La forma general de este tipo de ecuaciones en aplicaciones RSM es la siguiente:
\[\hat{Y} \;= \;b_0 \;+ \;b_1X_1 \;+ \;b_2X_2 \;+ \;b_{11}X_1^2 \;+ \;b_{22}X_2^2 \;+ \;b_{12}X_1X_2\]
La solución para juastar nestro modelo empírico de regresión curvilínea, está definida en esa ecuación. ¿Qué es lo que nos dice?
La \(\hat{Y}\) indica el vector que obtendremos de valores de predicción de ozono.
Las \(X_1\) y \(X_2\) son los vectores de valores en unidades originales de radiación solar y dióxido de nitrógeno. Esto indica las magnitudes del espacio euclidiano bidimensional donde navegaremos con nuestros datos.
Las \(X_1^2\) y \(X_2^2\) indican los mismos vectores de valores de radiación solar y dióxido de nitrógeno, pero elevados al cuadrado. Esto captura cualquier curvatura que pudiera estar presente en nuetro sistema.
La expresión \(X_1X_2\) indica algo que en estadística se conoce como “interacción” y lo que procura expresar es la naturaleza de la relación entre ambas variables predictoras (si es que acaso existe alguna), que puede ser sinergística o antagonística. No lo sabemos, pero es justo lo que vamos aaveriguar.
Lo que definitivamente sí sabemos, desde una perspectiva físicoquímica, es que las interacciones de carácter cuántico entre materia y energía sí existen y son totalmente reales. Por eso nos queda perfectamente claro que la radiación solar provoca cambios en la conducta del dióxido de nitrógeno, descomponiéndolo para entonces dar lugar al ozono. Pues bien, esta interacción físicoquímica entre la radiación solar y el dióxido de nitrógeno es justamente lo que la interacción estadística \(X_1X_2\) intenta retratar. Por eso se incluye en el modelo de regresión.
Finalmente, las literales \(b_n\) son los coeficientes que vamos a estimar.
Toda esta información la vamos a generar a partir de nuestros datos originales y lo que hacemos, es procesar cada uno para satisfacer el formato que la fórmula RSM exige.
Lo único que tienes que hacer es elevarlas al cuadrado y multiplicarlas. Con eso capturamos las posibles curvaturas e interacciones.
Los objetos matemáticos que producimos son las matrices siguientes:
\[\quad\quad\quad\quad I \quad X_1 \quad X_2 \quad X_1^2 \quad X_2^2 \quad X_1X_2 \quad\quad\quad\quad\quad Y \;(Ozono)\]
\[ X = \begin{bmatrix} 1 & 75 & 8 & 5625 & 64 & 600 \\ 1 & 81 & 14 & 6561 & 196 & 1134 \\ 1 & 44 & 6 & 2025 & 36 & 270 \\ 1 & 82 & 11 & 6724 & 121 & 902 \\ 1 & 72 & 8 & 5184 & 64 & 576 \\ 1 & 77 & 12 & 5929 & 144 & 924 \\ 1 & 65 & 7 & 4225 & 49 & 455 \\ 1 & 61 & 8 & 3721 & 64 & 488 \\ 1 & 47 & 6 & 2209 & 36 & 282 \\ 1 & 72 & 11 & 5184 & 121 & 792 \\ \end{bmatrix} \quad\quad\quad Y = \begin{bmatrix} 7 \\ 12 \\ 4 \\ 21 \\ 7 \\ 12 \\ 11 \\ 9 \\ 3 \\ 8 \\ \end{bmatrix} \]
Si observas con detalle, verás que la columna \(I\) representa el vector unitario para estimar el intercepto (si es que fuera estadísticamente significativo).
Ya lo tenemos todo. ¿Ahora qué sigue?
Para no tener que hacerlo manualmente, ni con lápiz ni papel ni calculadora, pero tampoco con software hojas de cálculo (porque ocurren muchos errores), yo prefiero apoyarme con el lenguaje de programación R.
Te sugiero que hagas lo mismo; desde luego que puedes usar Python, C, C++.
Además, están los clásicos de licencia comercial como Mathlab, Kaleidagraph o Graphpad Prism. Octave es gratuito y no veo diferencia. La decisión es tuya.
Yo uso R porque es con el que me siento más cómodo y además, es gratis. De hecho, todo este documento que estás leyendo fue creado dentro del entorno R, algo que también es posible con Python. Con todos los demás es imposible.
Enseguida los datos de las mediciones originales las organizamos en dos variables, X y Y para grabarlas en dos objetos del mismo nombre con el siguiente código R:
X <- matrix(c(rep(1, 10), contaminacion.atmosfera[ , c(3, 4)], contaminacion.atmosfera[ , c(3, 4)]^2,
contaminacion.atmosfera[ , 3] * contaminacion.atmosfera[ , 4]), ncol = 6, nrow = 10)
Y <- matrix(contaminacion.atmosfera[ , 2])Ya estamos listos para comenzar las maniobras.
Lo primero que hacemos es obtener el “inverso del producto de la matriz transpuesta”:
\[(X^TX)^{-1}\]
En lenguaje R se escribe de la siguiente manera y se convierte en esta matriz:
solve(t(X)%*%X)## [,1] [,2] [,3] [,4] [,5]
## [1,] 112.20365389 -3.929442353 1.82025265 0.0387225246 0.252406513
## [2,] -3.92944235 0.524977384 -3.03574239 -0.0080585144 -0.094263515
## [3,] 1.82025265 -3.035742392 22.90496680 0.0521942167 0.669251836
## [4,] 0.03872252 -0.008058514 0.05219422 0.0001365118 0.001947004
## [5,] 0.25240651 -0.094263515 0.66925184 0.0019470037 0.040180187
## [6,] -0.10700532 0.066382058 -0.48570211 -0.0012288244 -0.019902462
## [,6]
## [1,] -0.107005318
## [2,] 0.066382058
## [3,] -0.485702109
## [4,] -0.001228824
## [5,] -0.019902462
## [6,] 0.011890285El siguiente paso es ahora obtener el producto transposición de variables predictoras y de respuesta. La definición es la siguiente:
\[X^TY\]
De nuevo, para no hacerlo manualmente, nos apoyamos en la tecnología y en lenguaje R, el equivalente de la expresión algebráica anterior lo escribimos así:
t(X)%*%Y## [,1]
## [1,] 94
## [2,] 6808
## [3,] 934
## [4,] 502910
## [5,] 9852
## [6,] 69529Ya tenemos las dfos piezas. ¿Por qué estmaos haciendo todo esto?
Porque lo que estamos buscando son los coeficientes \(\hat{B}\) del modelo RSM:
\[\hat{B} = \begin{bmatrix} \hat{b_0} \\ \hat{b_1} \\ \hat{b_2} \\ \hat{b_{11}} \\ \hat{b_{22}} \\ \hat{b_{12}} \\ \end{bmatrix}\]
Para obtener los coeficientes, multiplicamos ambas matrices. Lo que ahora tenemos es la ecuación normal de mínimos cuadrados.
Para nuestro caso queda como sigue:
\[\hat{B}\quad = \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad (X^TX)^{-1} \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad X^TY\]
\[ \begin{bmatrix} \hat{b_0} \\ \hat{b_1} \\ \hat{b_2} \\ \hat{b_{11}} \\ \hat{b_{22}} \\ \hat{b_{12}} \\ \end{bmatrix} \quad \begin{bmatrix} 112.20365389 & -3.929442353 & 1.82025265 & 0.0387225246 & 0.252406513 & -0.107005318 \\ -3.92944235 & 0.524977384 & -3.03574239 & -0.0080585144 & -0.094263515 & 0.066382058 \\ 1.82025265 & -3.035742392 & 22.90496680 & 0.0521942167 & 0.669251836 & -0.485702109 \\ 0.03872252 & -0.008058514 & 0.05219422 & 0.0001365118 & 0.001947004 & -0.001228824 \\ 0.25240651 & -0.094263515 & 0.66925184 & 0.0019470037 & 0.040180187 & -0.019902462 \\ -0.10700532 & 0.066382058 & -0.48570211 & -0.0012288244 & -0.019902462 & 0.011890285 \\ \end{bmatrix} \quad \begin{bmatrix} 94 \\ 6808 \\ 934 \\ 502910 \\ 9852 \\ 69529 \\ \end{bmatrix} \]
Esa es la representacdión.
Lo que ahora sigue, es realizar la multiplicación matricial y así obtendremos los coeficientes. Nuevamente, dejamos que la tecnología lo haga por nosotros.
El código R es el siguiente:
(B <- (solve(t(X)%*%X))%*%(t(X)%*%Y))## [,1]
## [1,] 16.29693789
## [2,] 3.38155486
## [3,] -30.91078994
## [4,] -0.07693716
## [5,] -1.71403873
## [6,] 0.87723582Y ahí los tenemos.
Seamos formales y escribámoslos apropiadamente en su forma vectorial, con notación matricial:
\[ \hat{B} = \begin{bmatrix} 16.29693789 \\ 3.38155486 \\ -30.910789940 \\ -0.07693716 \\ -1.71403873 \\ 0.87723582 \\ \end{bmatrix} \]
Ahora lo que hacemos es tomar nuestros recién estimados coeficientes y los colocamos con sus respectivas variables, en la ecuación cuadrática.
Como vemos, tienen muchos dígitos así que conviene redondear
Y así es como construimos un modelo de regresión RSM. Aquí lo tenemos.
\[\hat{Y} \;= \;16.30 \;+ \;3.38X_1 \;- \;30.91X_2 \;- \;0.08X_1^2 \;- \;1.71X_2^2 \;+ \;0.88X_1X_2 \]
Por lo general, la primera pregunta de interés es:
“¿Qué tan bien se ajusta la ecuación a nuestros datos?”
El criterio más utilizado es \(R^2\), el famoso Coeficiente de Determinación y que es la fracción de la variabilidad total en los datos, que el modelo puede explicar.
Para calcular \(R^2\) necesitamos dos cosas:
- La Suma de Errores Totales
- La Suma de Errores Cuadrados
De nuevo, calcularlos manualmente es laborioso y propenso a errores, pero R lo hace mucho más rápido, prácticamente al instante:
(SST <- round(sum((Y-mean(Y))^2),2))## [1] 234.4(SSE <- round(sum((Y-(X%*%B))^2),2))## [1] 20.99Con estos resultados calculamos \(R^2\), que como ya dijimos, propiamente recibe el nombe de Coeficiente de Determinación. Esta es su definición:
\[R^2 = \frac{SST - SSE}{SST}\]
Y en lenguaje R sencillamente escribimos:
(SST-SSE)/SST## [1] 0.9104522Muy bien. Ya lo calculamos. ¿Pero qué significa? ¿Cuál es el mensaje que nos comunica?
En resumen, esta cifra indica que el modelo de regresión múltiple cuadrático es capaz de explicar el 91% de la variación de los datos de contaminación por ozono.
Esas son buenas noticias: con muy pocos datos fue posible explicar arriba del 90% del fenómeno de contaminación atmosférica por ozono, radiación solar y dióxido de nitrógeno.
¿Pero es esto realmente así? A mí me parece demasiado bueno para ser verdad.
Si observas con cuidado, te darás cuenta de que el Coeficiente de Determinación es simplemente una “regla de tres”:
A un total lo multiplicas por 100 y ahora es considerado como un 100%. Enseguida le restas una fracción; y como esta fracción es una parte de ese total, entonces naturalmente también puede verse como un porcentaje.
Y así es como el Coeficiente de Determinación en realidad resulta de una regla de tres.
De acuerdo. ¿Y qué más nos dice?
Lo que queremos que nos diga, es si nuestros datos reales de niveles de ozono son lo más cercanos posible a las predicciones estimadas de niveles de ozono.
Recordemos:
- Los niveles reales de ozono provienen de las mediciones practicadas.
- Los niveles estimados de ozono provienen de las predicciones del modelo.
Nuestra intención es que la diferencia entre ambos sea la menor posible y al mismo tiempo, queremos que esa diferencia se aproxime a la realidad. Pero tampoco queremos que la diferencia sea cero, porque eso significaría un “sobreajuste”, es decir, que nuestra ecuación únicamente representaría a nuestro conjunto de datos, incapaz de representar nada más.
En otras palabras: no queremos que las diferencias entre niveles medidos y niveles estimados sean amplias; y tampoco queremos que las diferencias sean cero.
Debe existir una diferencia que debe ser mínima. Mínima, pero elevada al cuadrado.
Lo que buscamos es que los niveles de medición y los niveles de estimación se aproximen unos con otros y eso es suficiente.
¿Por qué?
Porque el método de los mínimos cuadrados nos enseña que, mientras las diferencias cuadradas entre nuestros datos de medición y los valores estimados sean lo mínimas posible, tendremos la mejor aproximación cuantitativa que las ciencias deterministas pueden lograr.
(Desde luego, mínimos cuadrados es solamente uno de los numerosos métodos disponibles que existen para estas tareas de investigación.)
La \(R^2\) de nuestro modelo de regresión cuadrática es del 91%.
Puesto que esta cifra dice que nuestros datos representan en buena medida el fenómeno observado, resulta tentador proceder a predecir los niveles de ozono con base en nuevas mediciones de radiación solar y niveles de dióxido de nitrógeno.
Pero debemos ser cuidadosos.
La ecuación no debe usarse, hasta que hayamos verificado que todos los supestos sean correctos.
Una de las primeras pruebas que podemos aplicar es la exploración visual de la Superficie de Respuesta en tercera dimensión.
Aquí la tenemos.
Y como sabemos, la ecuación correspondiente es:
\[\hat{Y} \;= \;16.30 \;+ \;3.38X_1 \;- \;30.91X_2 \;- \;0.08X_1^2 \;- \;1.71X_2^2 \;+ \;0.88X_1X_2 \]
Estas representaciones tridimensionales son muy populares entre investigadores, científicos e ingenieros cuando se trata de proyectos de optimización de diseños experimentales DoE como Box-Bhenken y Diseño Compuesto Central para proyectos de sistemas de calida como Six Sigma.
Desde luego, el sistema \(6\sigma\) no es la única aplicación del método RSM. Ciertamente es de las más publicadas, más exploradas y más estudiadas. Pero en la química aplicada si algo existe de sobra, son proyectos de investigación científica publicados en el sistema mundial DOI, que utilizan regresiones cuadráticas para modelar fenómenos naturales con el método RSM.
¿Y qué es lo que vemos en nuestro modelo RSM en 3D?
A simple vista, lo que resalta al instante es la escala vertical de resultados espurios: predicciones falsas carentes de todo significado.
La escala de unidades Dobson va desde 0 en la parte superior, hasta -200 en la parte inferior. Esto no tiene ningún sentido. ¿Por qué?
Veamos una vez más el rango de nuestros datos originales para ozono:
\[\quad\quad\quad\quad\quad\quad Ozono \quad\quad\quad\quad\]
\[ Y = \begin{bmatrix} 7 \\ 12 \\ 4 \\ 21 \\ 7 \\ 12 \\ 11 \\ 9 \\ 3 \\ 8 \\ \end{bmatrix} \]
Está claro que el rango de medición de ozono en unidades Dobson va desde un mínimo de 3 hasta un máximo de 21: todas son cifras positivas superiores a cero.
Pero de alguna manera, el modelo cuadrático RSM indica que el rango de predicción va desde un mínimo de -200 hasta un máximo de 0: todas son cifras negativas inferiores a cero.
Naturalmente, el modelo no funciona ni es capaz de predecir nada. Y aún así, el Coeficiente de Determinación indica que es capaz de describir la variación en un 91%.
Esta es una obvia discrepancia.
¿Qué es lo que debería suceder?
Lo que esperamos que ocurra es que, si los modelos gráfico y matemático en verdad describen nuestro fenómeno con la mayor aproximación posible, entonces el rango de la variable de respuesta debería ser del mismo rango y orden. No esperamos que sea equivalente (y ya explicamos por qué), pero sí esperamos que sea semejante.
Y vemos que, definitivamente, los coeficientes que estimamos para el modelo de regresión cuadrática no representan a nuestros datos, en ninguna forma.
¿Por qué sucedió esto?
Se dice que “los coeficientes son inestables”. Pero eso suena ambiguo.
Lo cierto es que existen dos causas:
Los órdenes de magnitud entre las variables predictoras.
La correlación entre las variables predictoras.
Respecto a la primera causa, cuando al menos una variable predictora es al menos entre 5 y 10 veces mayor que el resto, los modelos producirán predicciones espurias.
Respecto a la segunda causa, cuando las variables predictoras exhiben conductas de cambio conjunto, sea directa o inversamente proporcional, los modelos producirán predicciones espurias.
Las estimaciones espurias obedecen a un fenómeno conocido como MULTICOLINEALIDAD, que termina por invalidar todo esfuerzo predictivo.
Por eso existen diversas técnicas y transformaciones que intentan mitigar los efectos de la multicolinealidad en los datos. Pero todos ellos se enfocan en transformar los valores de las variables predictoras como logaritmos y Box-Cox.
El más popular es sin duda la codificación en rangos -1, 0, +1.
Aunque todas estas transformaciones sí ayudan un poco con los efectos adversos de la multicolinealidad, la realidad es que están severamente limitadas. El motivo es que cuando dos o más variables predictoras están correlacionadas entre sí, al menos en la experiencia de este autor, cualquier intento de predicción por regresión múltiple, es totalmente inútil.
Al menos en el ejercicio de las ciencias químicas, la correlación entre predictores es de esperarse: lo normal para nosotros los químicos, que estamos acostumbrados a trabajar con materiales y equipos e instrumentos y diversas forma de energía en nuestros proyectos, lo que estamos seguros que va a ocurrir, es que las variables predictoras guarden algún tipo de correlación.
Parece que la multicolinealidad forma parte natural de las transformaciones fisicoquímicas de la materia, al someterla a cambios y procesos.
Entonces por un lado tenemos que los modelos tradicionales de regresión múltiple se ven seriamente impactados por la multicolinealidad de las variables, las predicciones no parecen funcionar y así, muchos esfuerzos resultan inútiles.
Por el otro lado, si de algo podemos estar seguros en el trabajo fisicoquímico de transformación de materias primas, es que la multicolinealidad de las variables siempre estará presente en los datos.
¿Por qué?
Porque nuestras “variables” no son simples variables: son sustancias y reactivos en diferentes concentraciones y que están sometidos a formas de energía que los activarán para reaccionar y producir efectos. Es tal el nivel de actividad existente entre la materia y la energía, que sencillamente resulta ingenuo pensar que las “variables” no interactúan entre sí y que por lo tanto, no guarden correlación.
¡Por supuesto que guardan correlación!
Pero la realidad transformativa de la materia y la energía es algo que ni matemáticos ni estadísticos pueden entender: el primero porque su paradigma es purista y el segundo, porque su paradigma es estocástico. Y nosotros en ciencias aplicadas y otras formas de tecnología, nuestro paradigma es determinístico. Por eso este autor considera que, aunque el matemático y el estadístico tienen siempre buenas intenciones y lo que quieren es ayudar, no siempre pueden hacerlo.
Al menos no, cuando insisten en recomendarnos metodologías que a todas luces, no funcionan.
¿Y cómo se ve la multicolinealidad?
Analíticamente, primero estimamos la correlación. En R la obtenemos así:
cor(x = contaminacion.atmosfera[, 3], y = contaminacion.atmosfera[, 4])## [1] 0.8134779Al instante vemos que es superior al 80%.
Esta elevada correlación entre ambas variables predictoras nos dice que todo esfuerzo y maniobra por transformarlas a otra escala para mitigar la colinealidad, es inútil.
Gráficamente, lo único que tenemos que hacer es comparar las variables predictoras y establecer una curva de regresión lineal simple.
En R lo hacemos así:
contaminacion.lm <- lm(contaminacion.atmosfera[,4] ~ contaminacion.atmosfera[, 3])
plot(x = contaminacion.atmosfera[, 3], y = contaminacion.atmosfera[, 4],
pch = 16, cex = 0.75, xlab = "Radiación Solar (Watts/m2)", ylab = "Dióxido de Nitrógeno (ppb)",
main = "Correlación entre Radiación Solar y Dióxido de Nitrógeno")
abline(contaminacion.lm)
Este tipo de evidencia entre predictores es de lo más común en el trabajo de estimaciones fisicoquímicas multivariables. Por eso, en la experiencia y opinión de este autor, no tiene ningún caso insistir con regresiones múltiples que no conducen a ningún lado.
En estas situaciones, la investigación de este tipo de modelos debe detenerse y no debemos tampoco invertir más tiempo ni recursos en ella.
Existen más pruebas que se pueden aplicar, como el índice de condición y el factor de inflación de variables. No los vamos a ver porque ya está claro que este modelo no representa a nuestros datos, pero los menciono para que el lector investigue acerca de ellos si gusta conocer más.
¿Tenemos alguna opción que sí funcione contra la multicolinealidad y que al mismo timepo sí sea capaz de producir predicciones confiables?
Sí existe, pero antes de explorarla, démosles una última oportunidad a la metdología de regresión clásica. Simplifiquemos las cosas y en lugar de aplicar el sofisticado modelo cuadrático, hagamos una segunda iteración pero ahora con un modelo lineal y veamos que ocurre.
Regresión Lineal Múltiple de Primer Orden
Los objetos matriciales del ejercicio anterior sirvieron como ilustación del origen del método. Ya no lo haremos así porque toma demasiado tiempo y espacio. Ahora iremos directo al código y que la tecnología ses la que se encargue de las labores mecánicas, mientras nosotros evaluamos el significado de cada resultado.
Nuestros datos originales eran los siguientes:
\[\quad Ozono \quad\quad\quad\quad\quad\quad\quad X_1 \quad X_2\]
\[ Y = \begin{bmatrix} 7 \\ 12 \\ 4 \\ 21 \\ 7 \\ 12 \\ 11 \\ 9 \\ 3 \\ 8 \\ \end{bmatrix} \quad\quad\quad\quad X = \begin{bmatrix} 75 & 8 \\ 81 & 14 \\ 44 & 6 \\ 82 & 11 \\ 72 & 8 \\ 77 & 12 \\ 65 & 7 \\ 61 & 8 \\ 47 & 6 \\ 72 & 11 \\ \end{bmatrix} \]
La regresión lineal múltiple de primer grado en R la escribimos así,
O3.Y <- lm(ozono.Y ~ radiacion.X1 + NO2.X2 )
summary(O3.Y)##
## Call:
## lm(formula = ozono.Y ~ radiacion.X1 + NO2.X2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.113 -2.520 -0.409 1.375 7.595
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.5773 7.0849 -1.352 0.218
## radiacion.X1 0.2588 0.1704 1.519 0.173
## NO2.X2 0.1603 0.8213 0.195 0.851
##
## Residual standard error: 3.907 on 7 degrees of freedom
## Multiple R-squared: 0.5441, Adjusted R-squared: 0.4139
## F-statistic: 4.178 on 2 and 7 DF, p-value: 0.06396Y el análisis de varianza de los coeficientes es el siguiente:
anova(O3.Y)## Analysis of Variance Table
##
## Response: ozono.Y
## Df Sum Sq Mean Sq F value Pr(>F)
## radiacion.X1 1 126.963 126.963 8.3173 0.02352 *
## NO2.X2 1 0.582 0.582 0.0381 0.85077
## Residuals 7 106.855 15.265
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La ecuación lineal obtenida es la siguiente:
\[\hat{Y} \;= \;-9.58 \;+ \;0.26X_1 \;+ \;0.16X_2\]
Si observamos la columna de valores t de Student y los valores P, de inmediato vemos que NINGUNO de los coeficientes de la regresión es estadísticamente significativo.
Esto significa que, para propósitos analíticos, los coeficientes son indistinguibles de cero.
Además, el Coeficiente de Determinación Ajustado es de 41.4%
¿Y esto qué significa?
Significa que si acaso existe alguna relación cuantitativa entre las variables predictoras y la variable de respuesta, el modelo lineal múltiple de primer orden tampoco alcanza a describirla.
No hay nada más por hacer.
Podemos aplicar algunas reglas analíticas, pero no tendrán significado.
Por ejemplo, la prueba de residuales:
shapiro.test(O3.Y$residuals)##
## Shapiro-Wilk normality test
##
## data: O3.Y$residuals
## W = 0.92446, p-value = 0.3957Los residuales del modelo provienen de una distribución normal.
Si obtenemos la gráfica de cuantiles, los residuales deberían ajustarse a una línea recta y vemos que así sucede.
plot(O3.Y, pch = 16, cex = 0.7, which = 2)
Incluso si aplicamos la prueba de factor de inflación de variables:
## Warning: package 'car' was built under R version 4.2.3## Loading required package: carData## Warning: package 'mdatools' was built under R version 4.2.3##
## Attaching package: 'mdatools'## The following object is masked from 'package:car':
##
## ellipsevif(O3.Y)## radiacion.X1 NO2.X2
## 2.956361 2.956361Vemos que los valores están por debajo de 5, lo que indica que bajo este modelo de primer orden, los coeficientes no exhiben multicolinealidad.
Finalmente, podemos incluso graficar el modelo RSM tridimensional y producirá un gráfico atractivo. Observa:
Esta gráfica parece comunicarnos un mensaje útil. Inclusive, vemos que la escala vertical de la unidades Dobson de ozono ahora quedan dentro del rango original de las mediciones. Esto es una buena señal y no deberíamos ignorarla.
Aún así, nada de esto es útil.
¿Por qué?
Porque ninguno de los coeficientes de la regresión es estadísticamente significativo y por lo tanto, el modelo no puede ser validado. No lo debemos usar para predecir. Pero ya estamos acercándonos.
Incluso se siente como haber realizado mucho trabajo analítico, a cambio de nada.
Todo se debe a la multicolinealidad existente en los datos.
Sin importar si se trata de diseños experimentales o de estudios observacionales, puedes contar con la certeza de que si hay cambios fisicoquímicos ocurriendo en tu sistema, la multicolinealidad casi siempre es la causa del problema.
Tengamos presente:
La correlación de las variables predictoras con la variables de respuesta, siempre será algo bueno. Pero la correlación entre variables predictoras siempre será un problema.
Para decirlo de manera más concreta:
Cuando dos o más variables predictoras guardan correlación entre sí, los modelos de predicción no son más que artefactos despojados de significado tangible.
Vayamos por fin a conocer la única opción analítica que sí funciona.
Parte 3, Caso Avanzado RSM: Regresión por Componentes Principales
En la Ciencia de Datos existen muchos algoritmos de Machine Learning, pero uno que es de particular interés para los químicos que trabajan con datos multivariables, es uno muy popular, clásico, conocido como Análisis de Componentes Principales.
A diferencia de lo que se enseña en cursillos de algoritmos de Machine Learning para Cencia de Datos, categorizándolo como un “algoritmo de aprendizaje no supervisado”, lo cierto es que PCA (Principal Components Analysis, por sus siglas en inglés) es un vasto campo de estudio de la Estadística Inferencial. Se han escrito obras enteras solamente acerca de este algoritmo.
¿Y por qué es tan especial?
Porque más que haber sido creado, fue concebido justamente para atender el problema de la correlación entre variables que ocurre en prácticamente, todas las ciencias aplicadas donde se generan y analizan datos multivariables.
Y en Quimiometría, que es la disciplina aplicada formal que los químicos debemos dominar antes de hablar de la supuesta “Ciencia de Datos”, se comienza por el estudio de los Componentes Principales.
¿Y qué son los Componentes Principales?
No quiero extenderme demasiado ni tampoco plagiar las prestigiadas obras de sobra publicadas, pero diré lo siguiente:
Como disciplina, el Análisis de Componentes Principales es solamente una de las numerosas técnicas de análisis estadístico multivariable.
Como objeto de estudio, los Componentes Principales son la versión transformada de nuestros datos originales.
El Análisis de Componentes Principales, como su nombbre lo indica, es un análisis, un proceso de exploración de datos que comienza con la matriz de covarianza, \(S\).
Y lo que este proceso analiza son las versiones transformadas de nuestros datos, que antes estaban expresados en unidades originales, pero que después de transformarlos, ahoro quedan expresados en la forma de Componentes Principales.
¿Por qué tenemos que hacer esto?
Porque la teoría del álgebra lineal estipula que la ortogonalidad se alcanza al rotar nuestro vector de datos, en un ángulo tal que, la correlación entre ejes adyacentes se minimice.
Voy a repetirlo:
La teoría del álgebra lineal estipula que la ortogonalidad se alcanza al rotar nuestro vector de datos, en un ángulo tal que, la correlación entre ejes adyacentes se minimice.
Cuando eso sucede, nuestros datos entran automáticamente en otro espacio Euclideano donde la multicolinealidad disminuye de manera notable. Es entonces que estamos listos para comenzar nuestro análisis… sola que ahora ya no analizaremos nuestros datos en su estado original, sino ahora en su versión transformada… como componentes principales.
Pero, ¿componentes principales de qué?
Los componenentes principales de… la variación presente en nuestros datos, por supuesto.
La variación es una entidad matemática compuesta y como tal, se puede descomponer, o separar, en sus elementos más simples. Esto facilita demasiado nuestro trabajo analítico.
Resulta que para poder analizar datos multivariables, su variación debe ser independiente.
La metodología matemática que hace que la variación de nuestras variables se transforme a un formato independiente en donde la colinealidad es mínima, sólo es posible gracias al álgebra lineal, por medio de los valores y vectores eigen.
No voy a presentar una demostración (hay muchísimos textos en el tema), sino solamente un recordatorio:
Como bien sabemos, los valores eigen se obtienen a partir de la famosa y fundamental ecuación característica (o ecuación eigen).
Los matemáticos la escriben así:
\[Ax = \lambda x \quad\quad\quad\quad y \;su \;forma \;equivalente \quad\quad\quad\quad (A - \lambda I)x = 0\]
Los estadísticos la escriben así:
\[Su = lu \quad\quad\quad\quad y \;su \;forma \;equivalente \quad\quad\quad\quad |S - lI|u = 0\]
Y los químicos la escribimos así:
\[\hat{H} \psi = E \psi \quad\quad\quad\quad y \;su \;forma \;equivalente \quad\quad\quad\quad (\hat{H} - EI) \psi = 0\]
Esta última forma es, ni más ni menos, la Ecuación de Schrodinger para la química cuántica.
¿Qué tiene que ver todo esto con el análisis de datos de un problema químico multivariable?
El argumento que estoy intentando establecer, es que los químicos ya tenemos un amplio y sólido precedente de familiaridad, con la obtención de valores eigen. ¿Por qué digo esto?
Voy a presentarte unas ideas que ya conoces, que parecieran no tener ninguna relación con el problema de multicolinealidad, pero en unos momentos quedará perfectamente claro.
Este es un nuevo recordatorio:
En la teoría de orbitales moleculares, el cálculo de los niveles cuánticos permitidos, para cada estado de energía, están definidos precisamente por los valores eigen, representados como bien sabemos, por la letra \(E\) en la ecuación de Schrodinger.
Cuando resolvemos la ecuación de Schrodinger al encontrar los niveles permitidos (cuantizados) de energía de una molécula, lo que matemáticamente estamos haciendo es obtener los valores eigen. Porque cada valor eigen representa un nivel permitido de energía cuántica.
Pues bien… cuando los químicos transitamos de las química teórica a la química aplicada en análisis de datos multivariables libres de multicolinealidad, resulta que vamos a seguir trabajando con valores eigen… es solo que bajo estas nuevas circunstancias, el significado de los valores eigen ya no son niveles permitidos de energía en un orbital molecular… no no no… los valores eigen ahora significan otra cosa.
¿Qué significan?
En el trabajo de análisis estadístico multivariable aplicando el método de análisis de componentes principales, los valores eigen representan la cantidad de variación presente en nuestras variables, pero ahora en su nuevo estado transformado, donde los vectores eigen son ortogonales.
Si es la primera vez que te enteras de estos tecnicismos y desconocías esta terminología, despreocúpate: no es tan difícil de aprender cuando ya sabes qué es lo que estás buscando. ¿Y qué es lo que buscamos?
Transformar nuestras variables colineales, desde sus unidades originales a nuevas unidades transformadas donde la colinealidad es mínima. Con eso es suficiente para obtener una nueva ecuación de predicción que producirá estimaciones aproximadas.
Y me pareció apropiado presentar el recordatorio de la química cuántica porque así, como profesional de la química aplicada que trabaja con datos multivariables, tu curva de aprendizaje se va a acortar en gran manera. ¿Por qué lo afirmo?
Porque a final de cuentas, los profesionales de la estadística que obtienen y analizan la variación por medio de componentes principales, están buscando lo mismo que tú:
Que las predicciones con modelos multivariables de regresión, sean lo más confiables posible, liberándolos de colinealidad. Y ellos lo hacen con vectores eigen.
Así que, si asististe a tus clases de química cuántica, el cálculo de valores y vectores eigen para ti no debería ser nada extraño ni sofisticado.
En resumen: el método PCA transforma nuestras variables a un formato donde la colinealidad se minimiza.
Y la expresión matemática que define la transformación es:
\[z_i \;= \;U^T \;[x_i \;- \;\bar{x}]\] Es con esta útima expresión que nuestros datos colineales se transforman y pasan de un espacio Euclideano, a otro. (Hay muchos libros que explican como funciona esto. Si desea conocer más, se le invita a lector a explorar la amplia literatura.)
En esta expresión, \(z_i\) es el componente principal, \(U^T\) es la matriz transpuesta de vectores eigen, \(x_i\) es nuestro vector de datos en unidades originales y \(\bar{x}\) es el vector promedio, también en unidades originales.
El producto entre la matriz de vectores eigen y el vector en datos originales, se conoce como “proyección”: básicamente lo que estamos haciendo al multiplicar los vectores eigen por nuestros datos originales, es proyectar nuestros datos originales desde el espacio donde casi todo es colineal, hacia el nuevo espacio donde casi todo es ortogonal.
Así obtenemos vectores que reciben el nombre de componentes principales. En la gran mayoría de los casos, los vectores de componentes principales son libres de multicolinealidad.
Suficiente de antecedentes y recordatorios.
¿Cómo producir modelos de regresión múltiple, libres de multicolinealidad?
Ya aprendimos que el algoritmo de Análisis de Componentes Principales es el camino que la Quimiometría nos ofrece para obtener modelos predictivos con coeficientes estables.
La pregunta ahora es, ¿cómo lo hacemos posible?
Existen tres maneras:
La instantánea, pero incomprensible.
La improductiva, pero con lujo de detalles.
La pragmática, con suficiente rigor técnico pero sin sacrificar productividad.
Todas son válidas y seleccionaremos aquella que mejor atienda nuestros intereses.
La primera es cuando compramos software de licencia y lo operamos mediante una interfase gráfica de usuario (GUI): seguimos el instructivo, hacemos clic en las opciones y obtenemos resultados y estimaciones. No tenemos idea de lo que sucedió, el proceso de análisis ocurrió dentro del software como una caja negra y lo único que sabemos es que debemos confiar en lo que el software nos diga. No está mal, sino sencillamente esta es la opción más acertada cuando únicamente queremos los resultados del análisis de nuestros datos y no tenemos ni el interés ni el tiempo para conocer ni entender los protocolos aplicados para obtenerlos. Confiamos en que el fabricante del software ha hecho un buen trabajo.
La segunda opción es cuando usamos lenguajes de programación aplicando únicamente funciones matemáticas, sin recurrir a librerías. Esta forma también funciona y nos brinda el privilegio de analizar nuestros datos con total libertad y sin depender de ningún fabricante, porque somos nosotros mismos quienes escribimos nuestros propios algoritmos analíticos a partir de ecuaciones y modelos publicados en la literatura especializada, lo cual es muy gratificante. Las desventajas son que este mecanismo consume demasiado tiempo y está sujeto a errores. Por eso, antes de aplicarlo en casos reales, debemos validarlo y someterlo a escrutinio y evaluación con colegas capacitados. Este mecanismo es más apropiado en entornos académicos, a menos que nuestro propósito sea justamente desarrollar un software y ofrecerlo como un producto bajo licencia comercial a personas y organizaciones interesadas en aprovecharlo bajo la primera opción.
La tercera opción es la que actualmente predomina y continúa creciendo en esta segunda década del siglo XXI.
¿De qué manera?
Cuando empleamos lenguajes de programación de fuente abierta, pero apoyándonos en librerías especializadas, previamente validadas y preferentemente producidas y publicadas por autoridades reconocidas en la materia de nuestro interés.
Es así como obtenemos lo mejor de los dos mundos:
Utilizamos herramientas especializadas de programación que sabemos funcionan y además, son gratuitas.
Gozamos de gran libertad de edición: si no estamos de acuerdo con algún método matemático empleado por alguna librería, es totalmente posible adaptar y modificar nuestro algoritmo de manera tal, que funcionará justo como lo necesitamos.
Es cierto que desarrollar un algoritmo personalizado bajo esta tercera opción toma tiempo, pero sólo será la primera vez. De ahí en adelante ya es sólo cuestión de correrlo y aplicarlo con cada nuevo conjunto de datos.
Como puedes darte cuenta, programar tus propios algoritmos analíticos, aunque impone sus desafíos, tiene sus ventajas. Y ahora veremos un ejemplo.
Entonces, regresando a nuestro asunto: ¿cómo obtenemos nuestro modelo, libre de multicolinealidad?
Seamos pragmáticos y apliquemos la tercera opción.
Primero aplicamos formatos para identificar a nuestros datos.
X.pca <- X[,-1]
Y.pca <- Y
colnames(X.pca) <- c("Radiacion X1", "NO2 X2", "Radiacion X1^2", "NO2 X2^2", "Radiacion*NO2 Int")
colnames(Y.pca) <- c("Ozono Y")
Y.pca; X.pca[, c(1, 2)]## Ozono Y
## [1,] 7
## [2,] 12
## [3,] 4
## [4,] 21
## [5,] 7
## [6,] 12
## [7,] 11
## [8,] 9
## [9,] 3
## [10,] 8## Radiacion X1 NO2 X2
## [1,] 75 8
## [2,] 81 14
## [3,] 45 6
## [4,] 82 11
## [5,] 72 8
## [6,] 77 12
## [7,] 65 7
## [8,] 61 8
## [9,] 47 6
## [10,] 72 11Aplicamos una librería especializada que producirá el modelo y rendirá el informe.
modelo.pca <- pls(x = X.pca[,c(1,2)], y = Y.pca, cv = 1, ncomp.selcrit = "min",
scale = TRUE, info = "Contaminación Por Ozono")
plot(modelo.pca)
¿Qué es lo que vemos?
A diferencia de la clásica regresión lineal múltiple por mínimos cuadrados ordinarios, aquí lo que tenemos es una moderna regresión por componentes principales; (y existe otra más sofisticada conocida en Quimiometría como regresión por estructuras latentes).
Estos son los resultados de un análisis gráfico muy completo que nos indica cosas como las distancias ortogonales q/q0, las distancias de Mahalanobis h/h0, el número de componentes principales recomendados por el algoritmo, la relación de los niveles de ozono de referencia y los de predicción y finalmente, el más impoortante de todos:
RMSE: la raíz del error cuadrado medio; en este caso, el modelo nos muestra que el RMSE es mínimo, cuando nuestro modelo consta de un componente principal.
Esta cifra está expresada en las mismas unidades de nuestra variable de respuesta y refleja el error del modelo. Más abajo veremos por qué este estadístico es tan importante.
Ahora, como nuestro interés es obtener un modelo de regresión confiable cuyos coeficientes sean estables y las variables reflejen la variación a partir de la menor multicolinealidad posible, lo que hacemos es obtener los coeficientes acompañados de sus respectivas significancias estadísticas.
Aquí los tenemos:
show(getRegcoeffs(modelo.pca, full = TRUE))## Estimated Std. err. t-value p-value 2.5% 97.5%
## Intercept -7.2595348 NA NA NA NA NA
## Radiacion X1 0.1585026 0.01959526 8.125103 1.954989e-05 0.1141751 0.2028302
## NO2 X2 0.6515283 0.16067610 4.114930 2.617674e-03 0.2880537 1.0150029
## attr(,"name")
## [1] "Regression coefficients for Ozono Y"Como te puedes dar cuenta, la t de Student y los valores p indican que ambas variables predictoras son estadísticamente significativas y en gran manera. También vemos sus intervalos de confianza en sus respectivas columnas 2.5% y 97.5%. En estos casos el intercepto podemos aceptarlo sin mayor requisito.
Este es el modelo por componentes principales:
\[\hat{Y} \;= \;-7.26 \;+ \;0.16 X_1 \;+ \;0.65X_2\]
Hagamos un poco de memoria y comparemos contra el modelo clásico por regresión lineal múltiple:
\[\hat{Y} \;= \;-9.58 \;+ \;0.26X_1 \;+ \;0.16X_2\]
Aunque son diferentes, podemos notar que los órdenes de magnitud son semejantes.
¿Significa esto que podemos utilizar el modelo clásico?
No.
¿Por qué?
Porque estadísticamente, los coeficientes del modelo clásico son indistinguibles de cero, lo cual no sucede con el modelo PCA.
Podemos asumir que el modelo PCA no producirá resultados espurios: al provenir de una regresión de componentes principales y sus coeficientes siendo estadísticamente significativos, esperamos que las predicciones de niveles de ozono sean confiables y se aproximen a la realidad.
Aquí tenemos más detalles del modelo.
##
## PLS model (class pls) summary
## -------------------------------
## Info: Contaminación Por Ozono
## Number of selected components: 1
## Cross-validation: full (leave one out)
##
## Response variable: Ozono Y
## X cumexpvar Y cumexpvar R2 RMSE Slope Bias RPD
## Cal 90.62115 51.88438 0.519 3.358 0.519 0.0000 1.52
## Cv NA NA 0.247 4.202 0.419 -0.0416 1.21Como vemos, la \(R^2\) es del 52% lo cual nos dice que, aún y cuando podemos confiar en las predicciones, en realidad el modelo apenas alcanza a explicar un poco más de la mitad de la variación observada en los datos. ¿Esto qué significa?
Que hace falta reunir más datos, más observaciones, posiblemente, más variables y repetir el análisis con más iteraciones.
Lo anterior tiene sentido porque a final de cuentas, este modelo lo construimos a partir de tan sólo 10 observaciones y 2 variables predictoras: definitivamente, no sería responsable ni profesional asumir que con un conjunto de datos tan abreviado, nuestro modelo tenga plena capacidad de producir predicciones de gran precisión.
Seguramente estarán dentro del rango observado para niveles de ozono, pero la variabilidad será amplia, lo que lo vuelve impráctico.
Uno de los protocolos más importantes en regresión multivariable, es el análsis de valores residuales: las diferencias estandarizadas entre valores observados y valores estimados de predicción.
El software de paquetería comercial te lo reporta de manera automática, pero como nosotros lo que queremos es entender de dónde es que los residuales etandarizados salen exactamente, aplicaremos el rigor matemático y estadístico.
Este es el código en R:
# Los residuales: diferencias entre valores observados y valores estimados.
residuales <- c(2.8403880, 2.7005735, -0.2177472, -8.0955088, 2.3648801,
0.7635064, -3.3961665, -1.3786487, 1.0992581, 3.3194650)
# Calcular la desviación estándar de los residuales.
residual_desv.std <- sd(residuales)
# Calcular residuales estandarizados
residuales_estandarizados <- residuales / residual_desv.std
# Ordenar residuales estandarizados
orden_residuales_estandarizados <- sort(residuales_estandarizados)
# Calcular cuantiles teóricos
cuantiles_teoricos <- qnorm(ppoints(length(orden_residuales_estandarizados)))
# Gráficos de residuales estandarizados cuantil-cuantil
plot(cuantiles_teoricos, orden_residuales_estandarizados, pch = 16,
xlim = c(-2.5, 2), ylim = c(-2.5, 2),
xlab = "Cuantiles Teóricos",
ylab = "Residuales Estandarizados Ordenados",
main = "Gráfico de Cuantiles de Residuales Estandarizados")
abline(0, 1, col = "red", lwd = 2) # Agregar línea de referencia para normalidad perfecta
¿Y este tipo de gráfico para qué sirve? ¿Qué nos dice?
Examinar los gráficos de cuantiles de residuales estandarizados nos permite evaluar la normalidad de los residuales. ¿Para qué queremos hacer eso?
Para identificar cualquier desviación contra la conducta Normal, lo cual puede impactar la validez y confiabilidad del modelo de regresión.
En otras palabras: si los residuales no exhiben una conducta semejante a la curva de la distribución Normal, el modelo de regresión no tiene validez y por lo tanto, debemos abstenernos de utilizarlo para hacer predicciones.
Si lo observas, el gráfico exhibe algo de desviación. Para corroborarlo analíticamente, uno de los siguientes formalismos es aplicar una prueba de hipótesis:
shapiro.test(residuales_estandarizados)##
## Shapiro-Wilk normality test
##
## data: residuales_estandarizados
## W = 0.85251, p-value = 0.06222El valor p de la prueba de Shapiro-Wilk es algo superior al 6% y como el nivel de significancia elegido (α) es del 5%, (comúnmente utilizado en pruebas de hipótesis estadísticas), un valor p del 6% excede este umbral.
Esto implica que no hay suficiente evidencia para concluir que los datos se apartan significativamente de la normalidad.
Lo anterior indica que podemos seguir adelante con nuestro análisis, aunque debemos ser cautelosos: un 6% en realidad es un resultado marginal.
¿Cómo luce visualmente nuestro nuevo modelo predictivo?
Aquí lo tienes:
Comparémoslo visualmente y veamos de nuevo el modelo lineal clásico.
¿Notas las diferencias entre ambos?
Como te puedes percatar, las predicciones de unidades Dobson de ozono son diferentes según los cambios de dióxido de nitrógeno y radiación solar. ¿Por qué las predicciones no son iguales?
Ya vimos que los tres coeficientes de ambas ecuaciones son semejantes en orden, pero también vemos que sus prediccines definitivamente no lo son. ¿Qué está pasando?
Son varias cosas, pero una que tiene pleno sentido fisicoquímico es la estequiometría y actividad de los gases presentes en el aire, entre el ozono y el dióxido de nitrógeno.
Si observas ambos gráficos verás que:
El modelo de componentes principales muestra que la acumulación de dióxido de nitrógeno a nivel troposférico impacta en el aumento de ozono, a medida que la radiación solar también vaya en aumento.
El modelo clásico lineal múltiple muestra que la acumulación de dióxido de nitrógeno a nivel troposférico impacto muy poco, o casi nada, sobre el aumento de ozono; pero vemos que además, el aumento de la radiación solar aparenta ser la principal
causante de la contaminación.
Ahora preguntémonos:
¿Cuál de las dos hipótesis es no-válida?
¿La primera, donde el incremento de dióxido de nitrógeno tiene un alto impacto?
¿O la segunda, donde el incremento de dióxido de nitrógeno tiene un bajo impacto?
Por supuesto que la segunda no tiene ninguna validez. Ni estadística ni fisicoquímica.
¿Por qué?
El modelo clásico de regresión lineal múltiple por mínimos cuadrados ordinarios no goza de validez estadística, porque sus coeficientes son no-significativos y tamboco goza de validez fisicoquímica porque estequiométricamente, las cantidades no concuerdan. (No mostrado.)
Es por estos motivos que las calibraciones multivariables con modelos no validados producen resultados espurios. A este autor le llegó a ocurrir:
Los datos se ven bien, tu algoritmo se ve bien, les ecuacions se ven bien, nos parece que el modelo es confiable y comenzamos a medir muestras, pero cuando analizamos un material de referencia los resultados son un disparate. Pensamos que nuestros instrumentos se descalibraron así que volvemos a calibrar, pero ocurre lo mismo otra vez. Y comenzamos a dar vueltas en círculo.
Sucede en espectroscopía, cromatografía y procesos de transformación.
Nos queda claro entonces que la relación empíricametne cuantitativa entre ozono generado por dióxido de nitrógeno que a su vez es activado por radiación solar, para nuestro modesto conjunto de datos, rsulta mejor descrita por el modelo de regresión de componentes principales. Es el que mejor parece aproximarse a la realidad.
El antídoto para estos fenómenos en donde las matemáticas convencionales no parecen funcionar, es la “Ciencia de Datos al servicio de la Ciencia Química”: Quimiometría.
Aún no hemos terminado.
¿Qué otra información adicional podemos obtener?
Algo que podemos hacer es averiguar ANALITICAMENTE cuáles son los niveles máximos y mínimos de ozono que el modelo alcanzaría a predecir.
Este es el mismo concepto de máximos y mínimos del cálculo diferencial, pero en lugar de hacerlo manualmente, hagámoslo con algo de código aplicando un ALGORITMO DE OPTIMIZACION.
Algunas personas me han llegado a preguntar:
“¿Para qué buscar los niveles de manera analítica, si ya los conocemos de manera gráfica?”
En análisis de datos, los métodos gráficos nos ofrecen una idea del espacio donde estamos navegando, pero los métodos analíticos nos indican con precisión las coordenadas de interés. Los resultados se deben informar con la mayor precisión.
En R, primero definimos nuestra FUNCION OBJETIVO de maximización:
# Función objetivo MAX
O3.funcion.max <- function(variables){
X1 <- variables[1]
X2 <- variables[2]
ozono <- -1 * ( -7.26 + 0.16* X1 + 0.65*X2 )
return(ozono)
}Ahora definimos nuestra FUNCION OBJETIVO de minimización:
# Función objetivo MIN
O3.funcion.min <- function(variables){
X1 <- variables[1]
X2 <- variables[2]
ozono <- 1 * ( -7.26 + 0.16* X1 + 0.65*X2 )
return(ozono)
} Una vez que ya hemos definido nuestras funciones, lo que sigue es aplicar un algoritmo de optimización, apropiado al caso y capaz de identificar los niveles.
No voy a profundizar porque el tema de algoritmos de optimización es verdaderamente amplísimo y es motivo de mucha investigación, pero sí puedo decir que si nuestro modelo validado, ahora visto como una función objetivo, satisface los requisitos del algoritmo de optimización que vayamos a seleccionar, tendremos en nuestro poder información sumamente valiosa.
¿Valiosa, de qué manera?
Valiosa, porque al conocer los rangos de nuestras variables predictoras, también podremos conocer los valores máximos y mínimos de nuestras variables de respuesta.
Traducido en términos prácticos, tendremos ESPECIFICACIONES de “lo mejor y lo peor” de las variables de respuesta y las condiciones precisas de las variables predictoras que las producen, donde “lo mejor” es lo óptimo y “lo peor” es lo no-óptimo.
En diseño de experimentos DoE para proyectos de investigación y desarrollo (R&D), cuando tenemos dos o más variables de respuesta y requerimos encontrar niveles óptimos para todas, el método recibe el nombre de FUNCIONES DE DESEABILIDAD: así es como se diseñan prototipos y secretos industriales, en el menor tiempo y costo posibles.
Para seleccionar un algoritmo de optimización, apropiado en el caso de nuestro modelo, lo que hacemos es categorizar la función.
Sin explicar pormenores, diremos que requerimos aplicar una “optimización global determinista y nuestra función objetivo es diferenciable de primer orden y que usa información por gradiente”.
(La terminología para clasificar funciones es necesaria para seleccionar algoritmos de optimización apropiados para cada aplicación. Como te darás cuenta, existen muchos algoritmos de optimización que satisfacen diferentes requisitos. Hay muchos textos en el tema.)
Vamos a usar uno clásico, popular, bajo el nombre BFGS, cuasi-Newtoniano, también conocido como algoritmo métrico variable, publicado en 1970 por Broyden, Fletcher, Goldfarb y Shanno; pero hecho del dominio público en código Pascal en 1990 por el mismísimo Dr. John Nash.
Este es el resultado de la maximización de niveles de ozono:
## $par
## [1] 82 14
##
## $value
## [1] -14.96
##
## $counts
## function gradient
## 6 6
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"Y aquí tenemos el resultado de la minimización de niveles de ozono:
## $par
## [1] 45 6
##
## $value
## [1] 3.84
##
## $counts
## function gradient
## 10 10
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"¿Qué significan estos números?
En el caso de la maximizción, las leyendas $par y $value indican las condiciones bajos las cuales los niveles de ozono serán máximos:
14.96 unidades Dobson de ozono, a 82 Watts/m2 de radiación solar y 14 ppb de dióxido de nitrógeno.
En el caso de la minimización, las leyendas $par y $value indican las condiciones bajos las cuales los niveles de ozono serán mínimos:
3.84 unidades Dobson de ozono, a 45 Watts/m2 de radiación solar y 6 ppb de dióxido de nitrógeno.
Por su parte, la leyenda $count nos indica los números de intentos que el algoritmo requirió para explorar el espacio y la leyenda $convergence en nivel 0 indica que la búsqueda terminó con éxito.
Estos resultados son evidencia de que el modelo de regresión por componentes principales, alcanza a describir con suficiente éxito la conducta del ozono bajo diferentes condiciones de radiación solar y niveles de dióxido de nitrógeno.
¿Y qué hay de la INCERTIDUMBRE DE MEDICION?
La promesa del autor al inicio de este artículo, fue indicar cómo obtener la incertidumbre de medición a partir del modelo RSM, para satisfacer la norma ISO-17025 en el laboratorio químico.
¿Por qué tuvimos que hacer toda esta investigación?
¿Y este ejercicio cómo encaja exactamente con el requisito normativo de la ISO-17025?
Primeramente, para hallar la incertidumbe medición debida a modelos de regresión, no siempre es necesario hacer todo este trabajo. Muchos modelos son simples, algunos requieren mayor atención y otros se complican cuando los datos son multicolineales y el análisis se convierte en un problema serio.
Como bien sabemos, la incertidumbre de medición es un concepto metrológico y las calibraciones de modelos regresión, ya sean univariables o multivariables, aportarán su parte correspondiente al”presupuesto de incertidumbre”.
Sin embargo, como la calibración de instrumentos espectroscópicos y cromatográficos es el corazón del método analítico completo, su contribución de incertidumbre al presupuesto requiere particular atención. ¿Cómo sabemos esto? ¿Quién lo dice?
En la Guía Eurachem CITAC “Cuantificación de Incertidumbre en Medidas Analíticas”, en su apéndice E bajo el título “Procedimientos Estadísticos Utiles”, apartado E.4.8, página 117 a la letra dice y cito textualmente:
“Mientras los métodos de los cálculos indicados anteriormente son aceptables para la mayoría de casos comunes de regresiones por mínimos cuadrados, no aplican a métodos de modelización más generales que consideran las incertidumbres de x o CORRELACIONES entre x y/o y.”
¿De qué habla la guía? ¿A qué se refiere?
Claramente está alertando que cuando las variables x están correlacionadas entre sí, los métodos clásicos y comunes de regresiones por mínimos cuadrados, no aplican.
Y después de nuestro larguísimo estudio, ahora sabemos por qué.
La guía Eurachem no dice más y solo atina señalar que la norma ISO TS 28038 podría, tal vez, contener algo de ayuda. Pero no sugiere nada más.
Finalmente, entonces, ¿qué tenemos?
Lo que tenemos es que, cuando las variables x están correlacionadas sabemos que por un lado, la modelización será no-válida y por el otro, la Guía Eurachem no ofrece ninguna indicación adicional.
¿Qué hacemos?
Es en estos casos cuando nosotros mismos definimos nuesta propuesta de modelización y señalamos por qué desechamos un método de regresión, aunque sea popular, y adoptamos otro, aunque sea menos conocido.
Volviendo a citar la misma Guía Eurachem, página 116, apartado E.4.3 que por subtítulo dice: “La incertidumbre \(u(x_{pred}, y)\) de un valor predicho \(x_{pred}\) debido a la variabilidad de \(y\) puede estimarse de varias maneras.”
La manera que en nuestro caso nos interesa es la que en su subtítulo dice:
“Mediante información dada por el software utilizado para derivar las curvas de calibración”
El párrafo dice así:
“Algunos programas dan el valor \(S\) descrito de varias maneras, por ejemplo RMS error o el error estándar residual.”
Eso es lo que la guía Eurachem nos indica.
Volvamos a revisar el informe de nuestro modelo de regresión por componentes principales.
##
## PLS model (class pls) summary
## -------------------------------
## Info: Contaminación Por Ozono
## Number of selected components: 1
## Cross-validation: full (leave one out)
##
## Response variable: Ozono Y
## X cumexpvar Y cumexpvar R2 RMSE Slope Bias RPD
## Cal 90.62115 51.88438 0.519 3.358 0.519 0.0000 1.52
## Cv NA NA 0.247 4.202 0.419 -0.0416 1.21Si observas, el lenguaje R nos ofrece las columna RMSE en la fila de datos de calibración (la fila de validación cruzada en este ejemplo podemos ignorarla).
Ahí leemos que la Raíz del Error Cuadrado Medio (Root Mean Square Error) tiene un valor de 3.358, pero ¿qué significa?
Como este estadístico está expresado en las mismas unidades Dobson de niveles de ozono em el aire, por regla general podemos tomarlo como la incertumbre de medición del modelo de regresión.
Así, podemos escribir que:
\[u(Ozono) = 3.358\]
Este valor ahora se puede agregar al presupuesto de incertidumbre. Este es el número que reportamos en cumplimiento de la ISO 17025.
¿Y por qué es válido utilizar el RMSE de una regresión multivariable por componentes principales, como un estadístico quer expresa la incertidumbre medición?
En la amplia bibliografía en Quimiometría y en la estadística general, existe un concepto que se conoce como “figura de mérito” y en los ajustes de regresión una de las figuras de mérito más importante (si no es que la más importante) es la del error que el modelo produce. Su definición formal es:
\[RMSE \;= \;\sqrt{\frac{\;\Sigma_{i=1}^N \;(y_i \;- \;\hat{y})^2}{N \;- \;1}} \]
Un modelo con un ajuste perfecto tendrá un RMSE igual a cero.
Y como ya mencionamos, una ventaja del RMSE es que se encuentra expresado en las mismas unidades que la variable de respuesta y por lo tanto, ofrece una buena evaluación absoluta del error predictivo del modelo.
Entonces, para responder la pregunta: sí, en efecto, es válido utilizar el RMSE como expresión del error de nuestro modelo, porque mide las diferencias entre los valores reales observados en las mediciones y los valores de predicción del modelo.
Para más información aquí tienes una fuente: https://en.wikipedia.org/wiki/Root-mean-square_deviation
Algo que podemos notar es que, al compararlo contra los valores de nuestros datos de ozono, esta incertidumbre es bastante elevada.
Una maniobra que podemos aplicar para cuantificar la magnitud del RMSE, es contrastarlo contra el rango de los valores observados, como un porcentaje de dicho rango.
En lenguaje R esto se analiza así:
paste(round(3.358 / (range(contaminacion.atmosfera[,2])[2] - range(contaminacion.atmosfera[,2])[1]) * 100, 2), "%")## [1] "18.66 %"La cifra es altísima. Naturalmente, ello se debe a varias cosas. ¿Qué es lo que nos dice?
Que nos faltan más datos y que nuestra investigación no está terminada. Por ejemplo:
Se necesitan incluir más observaciones para obtener más datos y así, aumentar la precisión.
Se requieren más variables predictoras para robustecer el modelo y encontrar más fuentes de variación.
Los fenómenos de contaminación ambiental son conocidos por contener y exhibir mucha variación.
Aún así, con tan sólo 10 observaciones, 2 variables y una variación notable de casi el 19%, es posible retratar y/o explicar la conducta, los cambios y relación entre ozono, radiación solar y dióxido de nitrógeno en casi un 52%.
Para un estudio modesto con recursos limitados, no está tan mal.
Y con eso hemos terminado.
Conclusiones
¿Qué opinas?
¿Te parece que estos mecanismos avanzados puedan ayudarte en la validación de métodos analíticos?
Esta fue una presentación demasiado larga y seguramente con exceso de detalles, pero quise publicar este trabajo porque estoy seguro que en el mundo de habla hispana existen otras personas como yo, con el mismo interés en resolver y atender este intrincado proceso y que quizá, aún no han terminado con la búsqueda.
Este reporte no pretende, en ninguna forma, comunicar que la labor de análisis de datos químicos multivariables se reduce a procesarlos con el método aquí expuesto.
Faltan más consideraciones por exponer y explicar, pero me parece que para el lector inexperto que está en la búsqueda de más respuestas en este tema, este reporte habrá sido exitoso si cumple lo siguiente:
Ahorrarle tiempo al lector que no conocía cuáles preguntas hacer, ni por dónde comenzar a buscar, porque no conocía que estas enseñanzas existen y son de uso convencional.
Aunque sí es cierto que muchos de los aspectos más relevantes han sido aclarados, también es cierto que este campo es muy extenso y existen personas e investigadores mucho, pero mucho más talentosos que yo… pero de alguna manera, me parece que no han publicado algún documento claro y útil que vaya directo al problema, sin rodeos y sin exceso de rigor matemático.
Desde luego, sin las formalidades matemáticas que explican este proceso, no es posible comprender cómo funciona; pero saturar un texto introductorio que primero explore los desafíos antes de entrar al método, tampoco favorece al lector inexperto.
Lo único que he querido hacer con esta publicación es compartir mi visión acerca de este método cuantitativo que encuentro apasionante: es una máquina analítica que responde muchas preguntas, pero que necesita afinarse antes de usarla rutinariamente.
Mi labor está hecha. Espero que con mi información este tema se democratice aún más y que tú te beneficies.
Todos los errores, omisiones, contradicciones e inconsistencias que el lector encuentre, son culpa mía.
Este reporte continuará actualizándose.
Si requieres apoyo, asesoría u orientación con gusto recibiré tus mensajes en:
victor[a]apestoff.com
Estoy en linkeIn: https://www.linkedin.com/in/qivictorsotodeltoro/
Víctor Soto Del Toro, químico industrial de profesión.
Ced. Prof. 4405730
Egresado de Universidad Autónoma de Baja California
Tijuana, Baja California, México