“ LA CALIBRACIÓN CALIBRACIÓN EN EN QUÍ QUÍMI MICA CA ANALÍT ANAL ÍTIC ICA” A” AUTORES: Dr. Héctor C. Goicoechea Investigador Adjunto de CONICET. Profesor Adjunto Cátedra de Química Analítica I, Director del Laboratorio de Desarrollo Analítico y Quimiometría (LADAQ), Facultad de Bioquímica y Ciencias Biológicas, Universidad Nacional del Litoral, Santa Fe, Argentina
Dr. Alejandro C. Oliv Olivieri ieri Investigador Principal CONICET. Profesor Titular Cátedra de Química Analítica, Facultad de Ciencias Bioquímicas y Farmacéuticas, Universidad Nacional de Rosario, Rosario, Argentina
“ LA CALIBRACIÓN CALIBRACIÓN EN EN QUÍ QUÍMI MICA CA ANALÍT ANAL ÍTIC ICA” A” AUTORES: Dr. Héctor C. Goicoechea Investigador Adjunto de CONICET. Profesor Adjunto Cátedra de Química Analítica I, Director del Laboratorio de Desarrollo Analítico y Quimiometría (LADAQ), Facultad de Bioquímica y Ciencias Biológicas, Universidad Nacional del Litoral, Santa Fe, Argentina
Dr. Alejandro C. Oliv Olivieri ieri Investigador Principal CONICET. Profesor Titular Cátedra de Química Analítica, Facultad de Ciencias Bioquímicas y Farmacéuticas, Universidad Nacional de Rosario, Rosario, Argentina
LA CALIBRACIÓN CAL IBRACIÓN EN QUÍM QUÍMICA ICA ANALÍTICA ANA LÍTICA RESUMEN En el libro se presentan los aspectos más importantes de la calibración metodológica, una de las herramientas más importantes de la química analítica. En primer lugar se describe la técnica de calibración univariada, basada en las regresiones por cuadrados mínimos clásica y ponderada. Luego se profundiza en la aplicación de las técnicas de regresión aplicadas en comparación de métodos analíticos. A continuación se introduce al alumno en el manejo de calibración bivariada, para desarrollar los conceptos más significativos que sirven de base en calibración multivariada. De esta última se presentan las calibraciones: a) directa (método CLS), b) indirecta (ILS, PCR y PLS) y c) para datos no-lineales (ANN). Finalmente se realiza una presentación de los avances recientes en calibración multivariada, incluyendo la de orden superior a uno con utilización de la denominada “ventaja de segundo orden”. En cada uno de los capítulos se presentan ejercicios resueltos. Se anexa un capítulo de introducción al álgebra matricial y un CD conteniendo programas escritos en Basic y entorno MATLAB para la aplicación de las técnicas presentadas, como también lecturas adicionales correspondientes a producción científica de los autores.
2
ÍNDICE CAPÍTULO 1. INTRODUCCIÓN Definición de calibración metodológica y sus diferencias con calibración instrumental. Distintas aplicaciones analíticas de la calibración metodológica. Resolución de problemas originados por interferencias, efectos de matriz y variaciones instrumentales.
CAPÍTULO 2. REGRESIÓN LINEAL PARTE I: CALIBRACIÓN UNIVARIADA. Determinación del extremo superior del rango lineal. Preparación de patrones. Medición de la respuesta de los patrones. Estimación de los parámetros de la regresión. Predicción de muestras incógnita. Cifras de mérito del método. Sensibilidad de calibración. Sensibilidad analítica. Límite de detección. Límite de cuantificación. Rango dinámico. Rango lineal. Programas de computación. Ejercicios. Referencias.
CAPÍTULO 3. REGRESIÓN LINEAL PARTE II: EXACTITUD Y COMPARACIÓN DE MÉTODOS ANALÍTICOS Exactitud de un método analítico. Región de confianza en el caso homoscedástico. Regresión ponderada. Región de confianza en el caso heteroscedástico. Comparación de métodos analíticos. Programas de computación. Ejercicios. Referencias.
CAPÍTULO 4. CALIBRACIÓN BIVARIADA Determinación de dos analitos usando dos sensores. La etapa de calibración. La calibración en notación matricial. Etapa de predicción. Coeficientes de regresión. Colinealidad. Cifras de mérito. Ejercicios. Referencias.
CAPÍTULO 5. CALIBRACIÓN CALIBRACIÓN DIRECTA
MULTIVARIADA
PARTE
I:
Determinación de multianalitos usando múltiples sensores. El modelo cuadrados mínimos clásicos (CLS) en notación matricial: etapa de calibración. Etapa de predicción y coeficientes de regresión. Cifras de mérito. Colinealidad espectral. Interferentes no modelados. Ventajas y desventajas de CLS. Comparación de métodos. Ejercicios. Referencias.
CAPÍTULO 6. CALIBRACIÓN MULTIVARIADA CALIBRACIÓN INVERSA (ILS Y PCR)
PARTE
II:
Regresión por cuadrados mínimos inversos (ILS). Calibración. Predicción. Ventajas y desventajas de ILS. Regresión por componentes principales (PCR). Compresión de la información. Componentes principales y fuentes de variación espectral. Calibración. Predicción. Validación cruzada. Residuos espectrales. Cifras de mérito. Ventajas y desventajas de PCR. Más allá de PCR. Ejercicios. Referencias.
3
CAPÍTULO 7. CALIBRACIÓN MULTIVARIADA PARTE III: REGRESIÓN EN CUADRADOS MÍNIMOS PARCIALES (PLS) Regresión por cuadrados mínimos parciales (PLS). Un algoritmo para PLS. Calibración. Predicción. Residuos espectrales y cifras de mérito. Ventajas y desventajas de PLS. Más allá de PLS. Ejercicios. Referencias.
CAPÍTULO 8. CALIBRACIÓN RESPUESTAS NO LINEALES CONCENTRACIÓN DEL ANALITO
MULTIVARIADA PARA RESPECTO DE LA
Redes Neuronales Artificiales. Elementos simples de cálculo. Estructura de red óptima. Redes neuronales artificiales en calibración multivariada. Desarrollo de modelos de calibración: 1. Detección de no linealidad, 2. Muestras de entrenamiento, monitoreo y validación. 3. Arquitectura. Entrenamiento de una red. Modelos basados en métodos lineales. NL-PLS y SPL-PLS. Referencias
CAPÍTULO 9. RECIENTES DESARROLLOS EN CALIBRACIÓN MULTIVARIADA Nuevas tendencias en calibración multivariada. Calibración multivariada de orden superior a uno. Ventaja de segundo orden. Diferentes modelos utilizados. Referencias
ANEXO I. REPASO DE ÁLGEBRA MATRICIAL Operaciones con matrices: suma, resta y multiplicación Matrices especiales. Norma o longitud de un vector. Normalización de un vector. Determinante. Inversa. Inversa Generalizada. Seudoinversa. Autovectores y autovalores. Descomposición en valores singulares. Referencias
ANEXO II Programas de computación. Programas ejecutables escritos en QB. Programas ejecutables escritos en Matlab.
4
CAPÍTULO 1 INTRODUCCIÓN La palabra calibración es utilizada varias veces al día en el laboratorio químico, y en especial en el laboratorio analítico. La calibración se encuentra dentro de las actividades desarrolladas por una disciplina científica que ha tenido un crecimiento sorprendente en los últimos años: la quimiometría. Esta disciplina científica tuvo sus inicios por el año 1969, cuando se publican los primeros trabajos científicos. Con posterioridad aparecen las primeras revistas especializadas de investigación en quimiometría (Chemometrics and Intelligent Laboratory Systems y Journal of Chemometrics).1 La quimiometría trata, específicamente, de todos aquellos procesos que transforman señales analíticas y datos más o menos complejos en información, utilizando métodos de origen matemático, estadístico y otros procedentes del campo de la lógica formal para conseguir sus fines. 2 Algunos autores proponen una división de la quimiometría en tres áreas generales: optimización de los procedimientos experimentales y mediciones químicas, extracción de la máxima información química de los datos analíticos y calibración, validación y significancia de las mediciones analíticas. 3 En el contexto de la tercer área definida dentro de la quimiometría, es necesario distinguir entre dos tipos de calibraciones: calibración instrumental y calibración metodológica. En la práctica, en muchas oportunidades ambas palabras se emplean como si significaran lo mismo, aunque constituyen dos aspectos diferentes y complementarios que son claves de la trazabilidad química. 4 La calibración instrumental tiene que ver con la seguridad que el analista tiene de que un instrumento y/o aparato funciona correctamente. Por lo tanto, el objetivo perseguido al realizar tal operación es el de corregir la respuesta instrumental hasta que la misma alcance el valor considerado como verdadero correspondiente al patrón de comparación utilizado. Ejemplos de este tipo de calibración se tienen cuando se corrige la velocidad una centrífuga, la temperatura de un termómetro o la longitud de onda suministrada por una red de difracción. Por otro lado, la calibración metodológica analítica consiste en caracterizar la señal que proporciona un instrumento en función de las propiedades de un analito o de un grupo de ellos. 5,6 Concretamente se busca establecer una relación inequívoca entre la señal instrumental y la concentración del analito. Generalmente, la calibración se relaciona con el análisis cuantitativo. Si se analizan la técnicas analíticas utilizadas para el desarrollo del método analítico, la calibración puede ser dividida en diferentes grupos. Cuando se aplica un método absoluto como la gravimetría, solo se calibra la balanza. Por otra parte, cuando se aplica un método estequiométrico (por ejemplo, una volumetría), se calibra la bureta (instrumentos) y se emplea un patrón primario o secundario previa valoración. Finalmente, se deben considerar los métodos relativos, involucrados en la mayoría de las técnicas analíticas. En primer lugar se debe calibrar el instrumento, y realizar las acciones correctivas que correspondan. Luego se establece la relación señal-concentración para cada proceso de medida químico. El siguiente cuadro muestra de manera esquemática los tipos de calibración aplicados a los diferentes métodos, según lo discutido previamente:
5
Métodos Primarios
Calibración Instrumental Metodológica
Absolutos (Gravimetría, culombimetría, dilución isotópica seguida de espectrometría de masas)
Estequiométricos (Volumetrías)
Secundarios Relativos (La mayoría de las técnicas instrumentales)
En lo que sigue, utilizaremos la palabra calibración para referirnos a calibración metodológica y centraremos la atención en diferentes tipos de calibración, originados según las características de los datos que proporciona el instrumento utilizado para recoger la señal relacionada con la concentración del analito de interés. Resolución de problemas originados por variaciones ins tru mentales, interf erencias y efectos de matriz Cuando el sistema analítico presenta fluctuaciones instrumentales que llevan a la pérdida de repetibilidad entre medición y medición, el empleo de un estándar interno aparece como una solución inteligente. Este problema está presente en algunas técnicas instrumentales separativas tales como la cromatografía líquida de alta resolución (CLAR o HPLC, del inglés High Performance Liquid Chromatography ) o la electroforesis capilar (CE, del inglés Capillary Electrophoresis). Se puede decir que este tipo de calibración tiene implicancias tanto instrumentales como metodológicas, ya que tiene el doble fin, por un lado de minimizar las fluctuaciones instrumentales y por otro, de establecer referencias con fines cualitativos y cuantitativos. 4 El método consiste en utilizar un estándar interno diferente del estándar empleado en la calibración convencional, que es agregado tanto a las muestras como a los patrones usados para la calibración. Una vez obtenidas las respuestas instrumentales, éstas se normalizan respecto de la medida del estándar interno. Por ejemplo, se pueden normalizar los tiempos de retención (fines cualitativos) o las áreas debajo de un cromatograma (fines cuantitativos). En el último caso, la relación concentración-respuesta buscada en la calibración utiliza la áreas normalizadas como respuesta instrumental. El analista se encuentra muy a menudo con otro tipo de problema, la existencia de interacciones entre el analito de interés con el entorno en el que se encuentra en la solución problema. Así, se presentan exaltaciones o inhibiciones de las señales, lo que se denomina efecto matriz. Esta exaltación o inhibición es proporcional a la concentración del analito presente en la muestra, por lo que una buena solución al problema consiste en realizar la calibración sobre la misma muestra. Este procedimiento es conocido como adición estándar o calibración interna. Si bien este método permite solucionar la modificación de la pendiente de la curva de calibrado por el efecto matríz, presenta algunos inconvenientes que vale la pena citar: a) implica un mayor trabajo experimental, b) necesidad de grandes cantidades de muestra, lo que es prohibitivo por ejemplo cuando se trata de algunas muestras biológicas, c) no se puede aplicar a métodos no lineales, y d) está sujeto a un mayor nivel de incertidumbre en la medición por no realizarse las mediciones cerca del centroide de la curva de calibrado, detalle que se verá más adelante. 6
Finalmente, el problema más común que debe enfrentar el analista es la aparición de interferencias. Por muy sofisticado que que sea el procedimiento procedimiento utilizado para medir un analito, una parte de la señal obtenida se debe a causas diferentes de la concentración o la masa de ese analito en cuestión. Las interferencias originan la presencia de un error sistemático constante en la medición. Ese error puede ser positivo o negativo dependiendo del tipo de efecto que produzca la interferencia. Para poder tener señales separadas tanto de analito como de interferencia, se puede recurrir a diferentes estrategias, una de las cuales consiste en aplicar una técnica separativa acoplada a un procedimiento de detección de la señal generada por el analíto (principio en el que se basan las técnicas separativas como HPLC y CE). Fuera del ámbito de las técnicas separativas, se puede recurrir a la preparación de un blanco, es decir una solución que contenga todas las especies presentes en la muestra problema, menos el analito de interés. Esa solución ideal no está siempre al alcance del analista debido a la complejidad que pueden presentar determinadas matrices. Como veremos a lo largo del libro, la calibración multivariada aparece como una solución interesante, al alcance del laboratorio de mediana complejidad, para resolver este problema realizando una separación matemática de las señales. Referencias 1. G. Ramis Ramos, M. C. García García Álvarez Coque, Quimiometría, Editorial Síntesis, Madrid, 2001, Capítulo 1. 2. D. L. Massart, B. G. M. Vandeginste, L. M. C. Buydens, Buydens, S. De Jong, P. J. Lewi y J. Smeyers-Verbeke, Handbook of Chemometrics and Qualimetrics, Elsevier, Amsterdam, 1997, Capítulo 1. 3. R. G. Brereton, Chemometrics. Applications of Mathematics Mathematics and and Satistics Satistics to Laboratory Systems, Ellis Horwood, New York, 1990. Capítulo 1. 4. M. Valcárcel, Principios de química analítica, Springer-Verlag Ibérica, Barcelona, 1999, p. 150. 5. K. Danzer y L. L. A. Currie, Guidelines for calibration calibration in analytical chemistry. Part 1. 1. Fundamentals and single component calibration, Pure & Appl. Chem . 1998, 70, 993. 6. 3. J. N. Miller y J. C. Miller, Estadística Estadística y quimiometría quimiometría para para química química analítica, analítica, 4ta. Edición, Prentice Hall, Madrid, 2002.
7
CAPÍTULO 2 REGRESIÓN LINEAL PARTE I: CALIBRACIÓN UNIVARIADA. En este capítulo estudiaremos una de las más populares aplicaciones de la regresión lineal en química analítica: la recta de calibración univariada. La teoría se expone a continuación, pero se recomienda consultar paralelamente el ejemplo concreto que se analiza al final del capítulo. El análisis mediante recta de calibración puede hacerse cuando sólo el analito de interés presenta señal analítica o respuesta (absorbancia, fluorescencia, potencial eléctrico, corriente, etc.), o cuando la señal del blanco es constante. Las etapas que deben seguirse en un análisis mediante recta de calibración son: • Determinación del extremo superior del rango lineal li neal • Preparación de patrones • Medición de la respuesta de los patrones • Estimación de los parámetros de la regresión • Cálculo de las cifras de mérito del método • Predicción en muestras incógnita Las expresiones matemáticas que se presentarán a continuación y su empleo en el análisis univariado están tomadas, en general, del trabajo de referencia clásico de Danzer y Currie, preparado para la Unión Internacional de Química Pura y Aplicada (IUPAC).1 De la amplia literatura que existe en este campo, recomendamos también los libros de Gardiner 2 y Miller y Miller.3 Determinació Determinació n del extr extr emo superior del rango lineal Esta etapa es fundamental, ya que la regresión lineal está basada en la suposición de que los datos de respuesta analítica están linealmente relacionados con la concentración del analito. Si se sospecha que existen desvíos de la linealidad, se recomienda realizar un análisis exploratorio previo cuyo objeto es extender el rango de aplicabilidad de la técnica analítica a la máxima concentración posible. En dicho análisis, se incluyen patrones de concentración conocida del analito desde cero hasta valores que se desvíen visiblemente de la linealidad. Una prueba estadística apropiada permitirá luego decidir hasta qué concentración se cumple la relación lineal respuesta-concentración. Sin embargo, dado que los parámetros a emplear en esta prueba se obtienen del análisis matemático-estadístico de la regresión, diferiremos el cálculo detallado para más adelante. Preparación Preparación d e patro patro nes Una vez estimado el extremo superior del rango lineal de la técnica, deben prepararse patrones de concentración conocida dentro de dicho rango, e incluyendo el valor cero de concentración del analito (blanco). Usualmente, se preparan varios patrones (como mínimo cinco) con concentraciones concentraciones igualmente espaciadas entre cero y el extremo superior del rango lineal, y cada patrón se analiza por triplicado. Debe ponerse especial cuidado en la preparación de los patrones del analito para la calibración, de manera que las concentraciones de calibrado se conozcan con la máxima precisión posible. Este requisito se relaciona con con el hecho de que la recta de regresión regresión se ajusta mediante ecuaciones que suponen que los valores del eje x (concentraciones) tienen una incertidumbre considerablemente menor que los del eje y (respuestas). Sólo a modo de ejemplo, si se realizan mediciones de absorbancia como respuesta, podemos suponer que el nivel de incertidumbre en la respuesta puede ser de alrededor 8
de 0,005 unidades de absorbancia. Si los valores de las respuestas son, en promedio, de 1 unidad de absorbancia, esto implica un nivel relativo de incertidumbre de aproximadamente 0,5% en la respuesta. Por lo tanto, se deben preparar patrones de calibrado cuyas concentraciones se conozcan con un error menor al 0,5%. Preparar soluciones de calibrado, por ejemplo, con incertidumbres del orden del 0,1% en promedio, requiere pesar más de 200 mg de patrón, preparar soluciones en matraces calibrados de al menos 100 mL, tomar alícuotas con pipetas aforadas calibradas, etc. Medic Medic ión d e la respuesta de los patro nes Una vez preparados los patrones de concentración conocida, se miden sus respuestas analíticas, incluyendo réplicas de cada medición. Usualmente cada patrón se mide por triplicado. Es importante establecer la siguiente nomenclatura: si se emplean 6 patrones, cada uno por triplicado, entonces el número de niveles diferentes de concentración ( p) es 6, y el número total de puntos de la recta de calibrado ( m) es 18. Estimación d e los p arámetros arámetros d e la regresió regresió n El análisis de los datos de calibrado mediante regresión lineal implica el cálculo de la pendiente ( A) y ordenada al origen ( B) de la recta ajustada a la ecuación y = A x + B. Los valores estimados de A y B se calculan mediante las siguientes ecuaciones: ecuaciones: m
A =
∑=1 ( x − x )( y − y ) i
Q xy
=
Q xx
i
i
m
(2.1)
∑=1 ( x − x ) 2 i
i
B = y − A x
(2.2)
donde xi es la concentración de cada uno de los m patrones de calibrado, x es el promedio de las concentraciones concentraciones de calibrado, yi es la respuesta en cada punto e y es el promedio de las respuestas de los los patrones de calibrado. calibrado. Además de los valores individuales de A y B, es importante tener una idea de su incertidumbre asociada, ya que los datos instrumentales llevan asociados un error que depende del ruido instrumental, y el ajuste por cuadrados mínimos sólo provee estimaciones de la pendiente y ordenada al origen. Los desvíos estándar en los parámetros A y B se calculan con las l as siguientes ecuaciones: s A =
s y / x Q xx
s B = s y / x
(2.3) 1
m
+
x 2 Q xx
(2.4)
En las ecuaciones precedentes, el parámetro s y/x es el desvío estándar de los l os residuos de la regresión y está dado por: m
∑=1 ( y − yˆ ) 2 i
s y/x =
i
m−2
i
(2.5)
9
donde yi es la respuesta experimental de cada patrón de calibrado e yˆ i representa la respuesta estimada en cada punto, esto es, yˆ i = A xi + B. En la ecuación (2.5) se emplean m – 2 grados de libertad, ya que hay m datos disponibles, y 2 parámetros estimados en la regresión ( A y B). Estos parámetros estadísticos dan también una idea de la bondad de la regresión. Es deseable que s y/x sea lo más pequeña posible, no obstante su valor está limitado por el ruido instrumental. La distribución de los residuos, es decir, el modo en que los valores de ( yi – yˆ i ) varían con la respuesta, cumple también un papel importante en el análisis de la adecuación de los datos al modelo lineal, como veremos más adelante. Predicción en muestras in cógnita Los valores de A y B se requieren para realizar predicciones en muestras incógnitas, a través de la ecuación yinc = A xinc + B, de donde puede obtenerse la concentración estimada del analito en la muestra: xinc = ( yinc – B) / A
(2.6)
donde yinc es, en general, un promedio de las respuestas obtenidas para un determinado número de réplicas de la incógnita (habitualmente tres). Un resultado no es tal, sin embargo, si no está acompañado por su correspondiente nivel de incertidumbre. Para informar xinc con su incertidumbre asociada, y establecer su número correcto de cifras significativas, es necesario calcular el error estándar en la concentración predicha s( xinc), lo cual se lleva a cabo mediante la siguiente expresión: s( xinc) =
s y / x
1
A
n
+
1 m
+
( yinc − y ) 2 A2 Q xx
=
s y / x
1
A
n
+
1 m
+
( xinc − x ) 2 Q xx
(2.7)
donde s y/x es el desvío estándar de los residuos de la regresión dado por la ecuación (2.5), A es la pendiente de la recta de regresión, n es el número de réplicas de la muestra incógnita, m es el número total de patrones de calibrado, yinc es el promedio de las respuestas de las réplicas de la incógnita, y es el promedio de las respuestas de los patrones de calibrado, y Q xx fue definido en la ecuación (2.1). La ecuación (2.7) es responsable de que la incertidumbre en la predicción dependa de cada muestra y no de la calibración en forma global, ya que para cada muestra incógnita hay un valor predicho de la concentración ( xinc) y por lo tanto un valor asociado del desvío estándar s( xinc). La forma de la ecuación (2.7) proviene de un análisis de la propagación de las distintas fuentes de error a la concentración predicha. Puede demostrarse que hay dos fuentes principales de incertidumbre: 1) la señal medida para la muestra incógnita y 2) las señales medidas para las muestras de calibrado. La primera contribuye con el término (1/ n) dentro de la raíz cuadrada de la ecuación (2.7), ⎛ 1 ( x − x ) 2 ⎞ ⎟⎟ , que colectivamente reciben el nombre y la segunda con los términos ⎜⎜ + inc m Q xx ⎝ ⎠ de leva (del inglés leverage ). La leva mide, de algún modo, la "distancia" de la muestra incógnita al centro de la calibración. Dado que la leva es mínima cuando la concentración de la incógnita es igual al promedio de las concentraciones de calibrado (esto es, cuando xinc = x ), se concluye que el método posee su máxima precisión en este último caso. De ahí que se recomiende analizar muestras cuya concentración de analito sea cercana al centro de las concentraciones de calibrado. La extrapolación a
10
concentraciones mucho mayores o menores que el promedio de la calibración aumenta la leva y con ello el error en la predicción. Otra conclusión que puede extraerse de la ecuación (2.7) es que el efecto de la calibración sobre el error de predicción será también menor si m > n, es decir, cuando el número de patrones de calibrado es superior al de réplicas empleadas para predecir. En todo caso, el análisis de la ecuación (2.7) muestra que, para muestras no demasiado alejadas del centro de la calibración, y dado que en general se cumple que m > n, el error estándar en la concentración se puede aproximar por s( xinc) = s y/x / ( A n1/2). Debe notarse finalmente que el intervalo de confianza para la concentración predicha puede calcularse multiplicando el valor del desvío estándar dado por la ecuación (2.7) por el correspondiente coeficiente de student para un dado nivel de confianza (usualmente 95%) y un número de grados de libertad igual a ( m – 2). Cifras de mérito d el métod o Las cifras de mérito de un método analítico se utilizan regularmente con el propósito de calificar un determinado método y comparar sus propiedades analíticas con las provistas por otras técnicas. Incluyen, entre otras, las siguientes: • Sensibilidad de calibración • Sensibilidad analítica • Límite de detección • Límite de cuantificación • Rango dinámico • Rango lineal Debe notarse que la expresión "cifras de mérito" es la traducción correcta del inglés figures of merit , la cual no debe traducirse como "figuras de mérito". Sensibilidad de calibración La sensibilidad de calibración es igual a la pendiente de la recta de calibrado: SEN = A
(2.8)
Indica la variación de respuesta producida por una unidad de variación de concentración del analito, y sus unidades son de señal × concentración –1. Sensib ilid ad analítica La sensibilidad de calibración no es adecuada para comparar dos métodos analíticos cuando estos están basados en respuestas de diferente naturaleza (por ejemplo, absorbancia y fluorescencia, o absorbancia y medidas electroquímicas, etc.). Para ello es preferible utilizar la llamada sensibilidad analítica γ, definida por la relación entre la sensibilidad y el ruido instrumental: γ = SEN / s y (2.9) donde s y es una medida conveniente del nivel de ruido en la respuesta. Para estimar el nivel de ruido pueden usarse dos procedimientos, que en teoría deberían coincidir. En el primero, se estima el ruido instrumental (s y) a través de los desvíos de las réplicas de las mediciones de calibrado respecto de sus promedios: p
s y =
r
( y ∑∑ =1 =1 i
ij
− yi ) 2
j
m − p
(2.10)
11
donde p es el número de niveles de concentración estudiados en la recta, r es el número de réplicas de cada punto, yij es el valor de la respuesta correspondiente a cada nivel y réplica, e yi es el promedio de las respuestas de las réplicas para cada nivel de concentración. En la ecuación (2.10), el número de grados de libertad es m – p, ya que de los m datos disponibles, p grados de libertad se reservan para el cálculo de las p medias y i . Este cálculo se ilustra en forma detallada en el ejercicio resuelto que acompaña al presente libro. En el segundo método de estimación del nivel de ruido, se lo calcula como el desvío estándar de los residuos de la regresión lineal, el parámetro ya definido s y/x [véase la ecuación (2.5)]. Si los datos estudiados cumplen la relación lineal entre respuesta y concentración, los dos métodos anteriormente descritos deben proveer resultados similares en cuanto a la estimación del ruido instrumental. Límite de detección Es la mínima concentración detectable de manera confiable por la técnica. En la definición moderna, el límite de detección (LOD) se calcula en función del desvío estándar de la concentración predicha para una muestra blanco ( s0).4 Para estimar s0 se recurre a la ecuación (2.7), escrita del modo siguiente: s( xinc) =
s y / x
1
A
n
+
1 m
+
( xinc − x ) 2 Q xx
(2.11)
Si suponemos que se analiza una muestra por triplicado (lo más usual es n = 3) en la que el analito no está presente ( xinc = 0), la ecuación (2.11) se reduce a: s0 =
s y / x A
1 1 x 2 + + 3 m Q xx
(2.12)
aunque s0 será diferente si se emplea un número diferente de réplicas. En todo caso, es importante informar qué valor de n se considera en el cálculo de s0 y por lo tanto del LOD. Como se muestra en la Figura 2.1, el LOD se calcula mediante una prueba de hipótesis estadística. En primer lugar se fija una concentración llamada nivel crítico (L C en la Figura 2.1), a partir de la cual se toman decisiones respecto de la detección del analito. Para concentraciones superiores a L C, existe una probabilidad α de cometer el llamado error de tipo I o falso positivo. Este último consiste en aceptar erróneamente la hipótesis alternativa, admitiendo que el analito está presente cuando en realidad está ausente. Como se aprecia en la Figura 2.1, la probabilidad de cometer este error de tipo I está dada por la zona sombreada a la derecha de L C (área α), siendo la "distancia" de LC al cero de la escala igual al producto de s0 por el coeficiente t α,ν. Si α se toma igual a 0,05, entonces una concentración superior a L C tendrá sólo un 5% de probabilidad de constituir un falso positivo. Del mismo modo, existe una probabilidad β de cometer un error de tipo II o falso negativo, en el que se acepta erróneamente la hipótesis nula, admitiendo que el analito está ausente cuando en realidad está presente (zona sombreada a la izquierda de L C en la Figura 2.1, con probabilidad igual a β). Si β se toma también como 0,05, la probabilidad de obtener un falso negativo será del 5%. En este caso la distancia de L C a la concentración correspondiente a dicho valor de β es el producto del coeficiente t β,ν por 12
s0,
considerando que este último parámetro es muy cercano al desvío estándar en la concentración de una muestra blanco. Puede notarse entonces que el valor de LOD depende de α y β, y de los desvíos estándar de las dos curvas gaussianas de la Figura 2.1. En general, ambas probabilidades se toman como iguales 0,05, mientras que los desvíos estándar se suponen ambos iguales a s0. De este modo, el LOD está dado por: 5 LOD = 2 × t 0,05,m –2 × s0 (2.13) 6 7 definición que ha sido adoptada también por IUPAC e ISO. En la práctica, dado que m es un número relativamente grande, el valor de (2 ×t 0,05,m –2) tiende a 3,3, por lo que una ecuación aproximada para el límite de detección es LOD = 3,3 s0. Nótese que antiguamente se definía el LOD contemplando únicamente errores de tipo I, como la concentración correspondiente a una relación señal/ruido igual a 3, lo que equivale a fijar el límite de detección como LOD = 3 s bl / A, donde s bl es el desvío estándar en la señal del blanco. En esta aproximación, la probabilidad de cometer errores de tipo I era de 0,1%, que corresponde a t 0,001,ν = 3 (para un número muy grande de grados de libertad). Esta definición, ya abandonada por la IUPAC, no contempla los errores de tipo II. (tα,ν + tβ,ν) s0
Hipótesis nula: analito ausente
Hipótesis alternativa: analito presente a este nivel
α
β 0
LC
LOD
Predicción
Figura 2.1. Prueba de significación empleada para estimar el límite de
detección. L C es el nivel crítico, LOD el límite de detección, α y β las probabilidades correspondientes a errores de tipo I y II respectivamente, s0 el desvío estándar del blanco (en unidades de concentración) y t α,ν y t β,ν los coeficientes de student para ν grados de libertad. Límite de cuantificación Es la mínima concentración cuantificable en forma confiable. Este parámetro (LOQ) se toma como la concentración correspondiente a 10 veces el desvío estándar (en unidades de concentración) del blanco, con lo cual: LOQ = 10 s0
(2.14)
De este modo, el desvío estándar relativo (DSR) para una concentración igual al LOQ es del 10%, nivel que se toma convencionalmente como el máximo DSR aceptable para cuantificar el analito en una muestra. Rango dinámico Se considera que va desde la menor concentración detectable (el LOD) hasta la pérdida de relación entre respuesta y concentración; véase la Figura 2.2, adaptada de la 13
excelente obra de Valcárcel. 8 El rango dinámico es también el rango de aplicabilidad de la técnica. En la zona de pérdida de la linealidad, podría aplicarse, en principio, un método de regresión polinómica para la calibración (o algún otro de naturaleza no lineal), de modo que nada impide que dicha zona sea utilizada con propósitos predictivos.
Rango dinámico
a t s e u p s e R
Rango lineal Pérdida de la relación respuesta-concentración Extremo superior del rango lineal
Concentración LOD
LOQ
Figura 2.2. Rangos dinámico y lineal de un método analítico.
Rango lineal Se considera que el rango lineal comprende desde la menor concentración que puede medirse (el LOQ) hasta la pérdida de la linealidad (Figura 2.2). Una manera conveniente de medir el cumplimiento de la linealidad es a través de la relación que existe entre la varianza de la regresión, medida por ( s y/x)2 [ecuación (2.5)], y la del ruido instrumental, medida por ( s y)2 [ecuación (2.10)]. Si la primera es significativamente mayor que la segunda, se supone que hay causas de desvío de la ley lineal que son estadísticamente superiores al ruido en la respuesta. Para emplear esta prueba es esencial que se cumpla el supuesto bajo el cual se realiza el ajuste lineal, esto es, que los errores en concentración de calibrado sean menores que en respuesta. De lo contrario, se acumularían en (s y/x)2 incertidumbres derivadas de la imprecisión en las concentraciones de los patrones, que nada tienen que ver con el ruido instrumental o las pérdidas de la linealidad. La prueba estadística que se utiliza para determinar si los datos se ajustan a la ley lineal es la F : en primer lugar se calcula un valor "experimental" de F , dado por: F exp =
(s y / x )2 2 (s y )
(2.15)
14
Luego se compara este valor con el crítico que se encuentra en tablas de F (de una cola) para m – 2 y m – p grados de libertad, y un determinado nivel de confianza, por ejemplo 95%. Si F exp < F , se acepta que los datos se comportan linealmente. Alternativamente, se calcula la probabilidad pF asociada a este valor de F exp, y se considera que la prueba de linealidad es aceptada si pF > 0,05. Esta prueba se describe en detalle en el trabajo de Danzer y Currie. 1
A s o u d i 0 s e R
B s o u d i 0 s e R
C s o u d i 0 s e R
Concentración Residuos de la regresión. A) Comportamiento lineal. B) Comportamiento no lineal. C) Comportamiento lineal con alta incertidumbre en la concentración de los patrones. Figura
2.3.
15
También es útil, como en todo ajuste por cuadrados mínimos, examinar visualmente la distribución de los residuos de la regresión. Un gráfico de residuos ( yi – A xi + B) en función de xi puede ser muy informativo respecto de la presencia de no linealidades, ya que el valor de F exp puede resultar significativo no solamente porque la relación entre las variables no sea lineal, sino por incertidumbres en la preparación de los patrones. La Figura 2.3 ilustra casos representativos al respecto. En el caso A), el comportamiento es lineal: se espera que la distribución de los residuos sea al azar, y que la variabilidad interna de las réplicas a cada nivel de concentración sea comparable a la variabilidad global (precisamente este es el sentido de la prueba estadística F antes comentada). En el caso B) se aprecia visualmente que los residuos poseen un comportamiento parabólico, caso típico de desvíos de la ley lineal. Finalmente, en el caso C), los residuos muestran una variabilidad global significativamente mayor que la que presentan las réplicas a cada nivel. Esta situación es típica de la presencia de mayor incertidumbre en las concentraciones nominales de los patrones de calibrado que en la señal instrumental, aunque el sistema se comporte linealmente. De ahí que se haya puesto hincapié en la necesidad de contar con patrones cuya concentración se conozca con mayor precisión que el ruido instrumental. En general, sin embargo, la distribución de los residuos no es tan clara como los casos presentados en la Figura 2.3, por lo que es importante aplicar el criterio estadístico F . Debe notarse que no hemos empleado, en todo este documento, al parámetro r , el coeficiente de correlación, aún cuando popularmente se recurre a él como prueba de linealidad o de bondad del ajuste. En este sentido, vale la pena repetir textualmente el siguiente pasaje del trabajo de Danzer y Currie: "el coeficiente de correlación, que es una medida de la relación de dos variables azarosas, no tiene ningún significado en la calibración analítica, debido a que los valores de x no están distribuidos al azar". 1 El coeficiente de correlación se emplea para responder preguntas tales como: ¿está correlacionada la concentración de antimonio con la de plomo en muestras de agua de una zona productora de metales?. En este caso se trata de analizar si existe correlación entre variables sobre las que el operador tiene muy poco control. Programas de comp utación Los métodos descritos en este capítulo pueden aplicarse con cualquier programa comercial que sea capaz de efectuar una regresión por cuadrados mínimos. Los parámetros faltantes pueden calcularse luego "a mano" con las ecuaciones provistas en este documento. En este sentido, la obra de Gardiner. 2 hace una excelente descripción del uso de la planilla de cálculo EXCEL para propósitos analíticos en general, y para estudios mediante regresión univariada en particular. Similar tratamiento se puede encontrar en la obra de Harris. 9 Para quienes deseen introducirse al mundo del entorno matricial MATLAB, esencial para cálculos avanzados en quimiometía, se proveen dos rutinas que calculan todos los parámetros aquí descritos, y permiten calibrar y predecir a partir de datos univariados. Confiamos que la discusión del ejercicio resuelto que se acompaña, el contenido del documento 'COMO OPERAR CON MATLAB.DOC', así como las rutinas 'LR_CAL.M' y 'LR_PRED.M', proveerán la información requerida para organizar los datos e implementar las rutinas. También se proveen programas independientes ejecutables en QB, como alternativa para quienes no puedan acceder a MATLAB: 'LR_CAL.EXE' y 'LR_PRED.EXE'. Para operarlos puede consultarse el documento 'COMO OPERAR CON QB.DOC'. Material suministrado Para este capítulo se provee el siguiente material suplementario: Archivos de texto (*.TXT) conteniendo datos típicos.
16
Rutinas (*.M) para el entorno de programación MATLAB. COMO OPERAR CON MATLAB.PDF, documento de Word que explica el empleo del entorno MATLAB. Programas ejecutables en QB (*.EXE). COMO OPERAR CON QB.PDF, documento de Word que explica el uso de los programas en QB.
Ejercicio 1) La Tabla 2.1 proporciona un ejemplo de datos de respuesta-concentración para su análisis, incluyendo respuestas medidas por triplicado. Grafique los datos de respuesta en función de la concentración y compruebe en forma visual que se desvían de la linealidad. Establezca un límite superior del rango lineal en forma cualitativa, para luego compararlo con el calculado mediante una prueba estadística apropiada. Tabla 2.1. Concentraciones y respuestas para un rango en el que se sospecha que existen desvíos de la linealidad. Concentración Respuesta 1 Respuesta 2 Respuesta 3 del patrón 0,06 0,08 –0,06 0,00 1,41 1,00 1,44 1,56 2,82 2,76 2,90 2,00 4,08 3,00 4,15 4,20 5,46 5,52 4,00 5,29 6,69 5,00 6,61 6,54 7,70 7,69 6,00 7,79 8,83 7,00 8,89 8,97 9,88 9,77 8,00 10,03 10,91 10,65 9,00 10,84 11,90 11,87 11,81 10,00 Note que los valores de concentración están dados con una precisión de ±0,01, lo cual implica un error relativo porcentual promedio de 0,01 ×100/5 = 0,2% (tomamos 5 como el valor promedio de las concentraciones de calibrado). Los valores de respuesta también están informados con una incertidumbre de ±0,01 unidades, si bien un análisis cualitativo de la variabilidad de los replicados indica que la incertidumbre en esta medición es mayor que lo informado en la Tabla 1. Posteriormente haremos un análisis más detallado, pero en principio es importante verificar que la incertidumbre relativa es mayor en la respuesta que en la concentración. Usuarios de MATLAB: los datos de la Tabla 2.1 están contenidos, en el formato apropiado para ser estudiados por la rutina 'LR_CAL.M' de Matlab, en el archivo de texto 'DATOS_EJ_RES_COMPLETOS.TXT'. Usuarios de QB: los datos están en el archivo de texto 'D_E_R_C.TXT', para ser estudiados por el programa 'LR_CAL.EXE'. 2) La Tabla 2.2 muestra los mismos datos que la Tabla 2.1, restringidos hasta un límite superior de concentración para el cual se cumple la linealidad (más adelante se muestra cómo se llegó a esta conclusión).
17
Tabla 2.2. Concentraciones y respuestas para un rango en el que existe linealidad. Concentración Respuesta 1 Respuesta 2 Respuesta 3 del patrón 0,06 0,08 –0,06 0,00 1,41 1,00 1,44 1,56 2,82 2,76 2,90 2,00 3,00 4,15 4,20 4,08 5,29 5,46 5,52 4,00 6,54 6,69 5,00 6,61 Usuarios de MATLAB: los datos de la Tabla 2.2 están contenidos, en el formato apropiado para ser estudiados por la rutina 'LR_CAL.M' de Matlab, en el archivo de texto 'DATOS_EJ_RES_LINEAL.TXT'. Usuarios de QB: los datos están disponibles para ser estudiados por el programa 'LR_CAL.EXE' en el archivo de texto 'D_E_R_L.TXT'. Calcule los valores de la pendiente y ordenada al origen para la recta ajustada con los datos de la Tabla 2.2. 3) Estime los desvíos estándar en la pendiente y ordenada al origen, e informe los valores de A y B con el número correcto de cifras significativas. 4) La Tabla 2.3 muestra los valores de la respuesta para cuatro muestras incógnita, todos por triplicado. Tabla 2.3. Respuestas para cuatro muestras incógnita. Muestra Respuesta 1 Respuesta 2 1 0,69 0,65 2,20 2,13 2 3 3,55 3,41 4,71 4 4,82
Respuesta 3 0,75 2,05 3,52 4,70
Los datos de la Tabla 2.3 están contenidos, en el formato apropiado para ser estudiados por la rutina 'LR_PRED.M' de Matlab, en el archivo de texto 'DATOS_EJ_RES_TEST.TXT'. Estime la concentración del analito en las cuatro muestras de la Tabla 2.3, calcule sus desvíos estándar e informe el resultado con el número apropiado de cifras significativas. 5) Calcule las cifras de mérito del método.
Respuesta 1) El análisis de estos datos mediante los programas LR_CAL.M (Matlab) o LR_CAL.EXE (QB) indica que los datos no se comportan en forma lineal. En particular, se obtiene un valor de F exp de 8,88, con una probabilidad asociada pF de 0,001. La gráfica de los residuos es informativa al respecto:
18
Figura 2.4. Gráfica de los residuos del Ejercicio 1
2) Los valores estimados, dados por las ecuaciones (2.1) y (2.2) son, para el ejemplo de la Tabla 2.2, A = 1,3174 y B = 0,1237. Estos últimos números tienen, probablemente, más cifras significativas que lo permitido por sus desvíos estándar. Para acotarlos al número correcto de cifras es necesario estimar sus incertidumbres. 3) Los desvíos estándar calculados son s y/x = 0,1, s A = 0,01 y s B = 0,04. Lo correcto es informar la pendiente y ordenada al origen de la recta ajustada del modo que sigue: A = 1,32(1) B = 0,12(4)
En la Tabla 2.4 encontrará un resumen de todos los cálculos intermedios necesarios para estimar A, B y sus errores estándar.
19
Tabla 2.4. Parámetros necesarios para el cálculo de A, B, s A y s B. i xi xi – x yi ( xi – x )2 ( xi – x ) ( yi – y ) yi – y 8,39 –3,36 6,25 –2,50 0,06 0,00 1 2,97 2,25 1,44 –1,98 –1,50 2 1,00 0,30 –0,60 0,25 –0,50 2,82 2,00 3 0,37 0,25 4,15 0,73 0,50 4 3,00 2,81 1,87 2,25 5,29 4,00 1,50 5 6,25 7,98 6,61 3,19 2,50 6 5,00 8,34 –3,34 6,25 0,08 0,00 –2,50 7 2,25 2,79 1,56 –1,86 –1,50 8 1,00 0,33 –0,66 0,25 2,76 2,00 –0,50 9 0,39 0,78 0,25 4,20 3,00 0,50 10 2,25 3,06 2,04 1,50 5,46 4,00 11 7,81 3,12 6,25 6,54 5,00 2,50 12 6,25 8,69 –3,48 –2,50 –0,06 0,00 13 3,01 2,25 1,41 –2,01 1,00 –1,50 14 0,25 0,26 –0,52 –0,50 2,90 2,00 15 0,33 0,25 4,08 0,66 0,50 16 3,00 2,25 3,15 2,10 1,50 5,52 4,00 17 8,18 6,25 6,69 3,27 2,50 18 5,00 Total Q xx = Q xy = 69,17 52,5 Promedio x = 2,50 y = 3,42 4) Los valores de predicción se muestran en la Tabla 2.5. Tabla 2.5. Predicciones en muestras incógnita. Muestra Respuesta Concentración Desvío DSR = 100 a s( xinc) / xinc promedio predicha estándar ( yinc) ( xinc) (%) s( xinc) 1 0,70 0,44 0,05 12 2 2,13 1,52 0,05 3,3 3,49 2,56 0,05 1,9 3 4 4,74 3,51 0,05 1,4 a A partir de la ecuación (2.6), insertando s y/x = 0,1; A = 1,32; n = 3; m = 18; yinc de la columna 2 de la Tabla 2.5, y = 3,42 y Q xx = 52,5. Note que los valores pueden aproximarse por s( xinc) = s y/x / ( A n1/2), tal como se dijo en la parte teórica. Puede notarse que la concentración predicha se acotó a dos cifras decimales significativas, teniendo en cuenta que los desvíos estándar son todos aproximadamente de 0,05 unidades. Nótese que los valores de s( xinc) son iguales en la Tabla 2.5 porque se informan con una sola cifra significativa, aunque su cálculo detallado demuestra que difieren entre sí, de la manera prevista por el efecto de la leva. Es importante destacar también que el desvío estándar relativo (DSR) dado en la Tabla 2.5 es alto para la primera muestra, y razonablemente bajo para las otras. En el primer caso, la concentración predicha es también baja. Estas consideraciones se
20
relacionan con la mínima concentración detectable por la técnica, que se considerará a continuación. También pueden fijarse los intervalos de confianza alrededor de una predicción, empleando los coeficientes de student de dos colas para un 95% de confianza y ( m – 2) grados de libertad. Por ejemplo, para la muestra Nº 4 en la Tabla 2.5: xinc = 3,51 ± t ( p = 0,05; 16 GL) × s( xinc) = 3,51 ± 2,1 × 0,05 = 3,5 ± 0,1 5) Es importante analizar la gráfica de los residuos para este caso.
Figura 2.5. Gráfica de los
residuos del Ejercicio 1, para menos puntos.
Como puede verse en la figura anterior, la distribución de los residuos conserva aún rastros de la falta de linealidad de los datos, pero la prueba F dice que esta impresión no es estadísticamente relevante: F exp = 1,58, pF = 0,21. La Tabla 2.6 ilustra el cálculo detallado de s y para esta prueba. En el presente ejemplo, la sensibilidad está dada por SEN = 1,32 (Unidades de respuesta)×(Unidades de concentración) –1 Para el cálculo de la sensibilidad analítica se requiere una estimación del nivel de ruido instrumental. Para los datos de la Tabla 2.2, p = 6, r = 3, s y = 0,08 (véase la Tabla 2.6 para el detalle del cálculo).
21
Tabla 2.6. Parámetros requeridos para el cálculo de s y. i
1 2 3 4 5 6
j
yij
1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
0,06 0,08 –0,06 1,44 1,56 1,41 2,82 2,76 2,90 4,15 4,20 4,08 5,29 5,46 5,52 6,61 6,54 6,69
0,03 1,47 2,83 4,14 5,42 6,61 p
Total
( yij – yi )2 0,0009 0,0025 0,0081 0,0009 0,0081 0,0036 0,0001 0,0049 0,0049 0,0001 0,0036 0,0036 0,0169 0,0016 0,0100 0,0000 0,0049 0,0064
yi
r
( y ∑∑ =1 =1 i
ij
− yi ) 2 = 0,081
j
A partir de los resultados de la tabla anterior, se puede calcular un nivel de ruido instrumental de (0,081/12) 1/2 = 0,08. Dado que, para los mismos datos, s y/x = 0,1, puede notarse que ambos procedimientos para estimar el ruido producen resultados similares. Empleando 0,1 unidades de respuesta como nivel de ruido, podemos calcular la sensibilidad analítica para el ejemplo en estudio a partir de la ecuación (2.10), como γ = SEN / s y/x = 13 (Unidades de concentración) –1. El parámetro γ se interpreta mejor en términos de su inversa. El valor de γ –1 (0,08 unidades de concentración en nuestro caso) indica la menor diferencia de concentración que puede apreciarse a lo largo del intervalo de aplicación de la técnica analítica. Con respecto al límite de detección, puede estimarse como LOD = 2 ×t 0,05,16 × 0.06 = 0,2. Se interpreta este último resultado diciendo que la técnica es capaz de detectar al analito cuando está en concentraciones superiores a 0,2. Para el ejemplo de la Tabla 2.2 el LOQ se calcula como 0,6 unidades de concentración. Se interpreta como la menor concentración que se puede cuantificar, esto es, en el intervalo de concentración entre 0,2 y 0,6 la técnica puede detectar pero no cuantificar al analito. Con esto se comprueba que la concentración predicha para la muestra incógnita Nº 1 de la Tabla 2.5 está por debajo del LOQ, lo cual explica el alto valor de DSR. Con respecto al rango dinámico, la máxima concentración probada fue de 10,00 unidades (Tabla 2.1). Hasta esa concentración existe un cambio de respuesta al cambiar la concentración, por lo que, a falta de mayor información, supondremos que el rango dinámico está entre 0,3 y 10 unidades de concentración. Para estimar el rango lineal, se recurre a los datos de la Tabla 2.1, y se comprueba que para este caso, si se incluyen todos los datos, F exp = 8,88, pF = 0,001, con lo cual dichos datos se declaran no lineales. Si vamos quitando datos, comenzando con los de mayor concentración, y recalculamos los valores de F exp y sus pF asociadas, se obtienen los resultados informados en la Tabla 2.7. 22
Tabla 2.7. Rangos de concentración y estudio de la linealidad mediante la prueba F . Rango de concentración F exp pF 0-10 8,88 0,001 0-9 6,69 0,001 0-8 4,62 0,001 3,50 0,007 0-7 0-6 2,73 0,031 1,58 0,214 0-5 Estos resultados indican que a partir de una concentración de analito igual a 6 unidades se pierde la linealidad. En realidad, la no-linealidad se mantiene. Debería decirse que a partir de 6 unidades de concentración no es posible distinguir la incertidumbre por falta la linealidad de la incertidumbre intrínseca de la respuesta analítica. La Tabla 2.8 resume las cifras de mérito calculadas. Tabla 2.8. Cifras de mérito. Cifra de mérito Valor (unidades) Sensibilidad SEN = 1,32 (Unidades de respuesta) ×(Unidades de concentración) –1 Sensibilidad analítica γ = SEN / s y/x = 13 (Unidades de concentración) –1 Límite de detección LOD = 0,2 (Unidades de concentración) Límite de LOQ = 0,6 (Unidades de concentración) cuantificación Rango dinámico 0,2-10,0 (Unidades de concentración) Rango lineal 0,6-6,0 (Unidades de concentración) Ejercicios complementarios 1) Se analiza una serie de muestras patrones mediante dos métodos analíticos, uno basado en medidas de absorbancia y otro basado en medidas de fluorescencia. Los resultados se muestran en la siguiente tabla: Tabla 2.9. Concentraciones de patrones y respuestas obtenidas mediante dos métodos analíticos. Concentración Método A Método B del patrón Respuesta Respuesta Respuesta Respuesta Respuesta Respuesta 1 2 3 1 2 3 0,02 0,02 2,0 1,9 1,9 0,000 0,01 0,17 0,17 0,17 17,4 17,4 17,3 0,100 0,200 0,32 0,33 0,32 32,5 32,6 32,6 0,300 0,48 0,48 0,48 47,8 47,8 48,0 0,400 0,64 0,64 0,64 63,2 63,3 63,3 0,79 78,4 78,5 78,4 0,500 0,79 0,79
Calcule las cifras de mérito para cada método. ¿Cuál de estos métodos puede considerarse más sensible? ¿Qué parámetro(s) emplea para justificar la mayor sensibilidad de un método sobre el otro?. 2) Se mide por triplicado una muestra incógnita, usando ambos métodos descriptos en el problema anterior. Los resultados se presentan en la siguiente tabla: 23
Tabla 2.10. Medición de muestra incógnita Método A Respuesta Respuesta Respuesta Respuesta 1 2 3 1 0,25 0,26 0,25 25,2
Método B Respuesta 2 25,1
Respuesta 3 25,3
Calcular la concentración del analito por ambos métodos, y estimar su desvío estándar. ¿Qué comentarios pueden hacerse respecto de estos resultados? Se recomienda emplear las rutinas de MATLAB 'LR_CAL.M' y 'LR_PRED.M' (o sus versiones respectivas en QB), organizando los datos de los ejercicios propuestos de la manera que se presenta en los archivos de texto correspondientes al ejercicio resuelto. 3) En el análisis fluorimétrico de un compuesto, se realizan dos curvas de calibrado, empleando dos longitudes de onda diferentes para la excitación. En el caso A, la emisión del compuesto está superpuesta con la dispersión Raman del solvente, y el analista observa por lo tanto la presencia de un blanco constante de intensidad significativa. Decide modificar la longitud de onda de excitación, en este caso generando los datos del caso B, donde el blanco parece ser menor. En la tabla siguiente se informan los datos de calibración para cada caso, en sus respectivos rangos lineales. ¿Qué conclusiones pueden extraerse respecto de las cifras de mérito de estos dos casos? Tabla 2.11. Datos de calibración Muestra 1 2 3 4 5 6
Concentración 0,000 0,198 0,392 0,583 0,769 0,950
Muestra 1 2 3 4 5 6 7 8
Concentración 0,000 0,198 0,392 0,583 0,769 0,950 1,130 1,310
Caso A Respuesta 1 0,78 3,38 5,75 8,53 10,97 13,40 Caso B Respuesta 1 0,01 1,96 3,75 5,59 7,30 9,07 10,83 12,08
Respuesta 2 0,80 3,44 6,16 8,51 11,04 13,08
Respuesta 3 0,82 3,51 6,01 8,68 10,89 13,37
Respuesta 2 0,03 1,88 3,75 5,52 7,35 8,95 10,71 12,11
Respuesta 3 0,04 1,90 3,80 5,56 7,27 9,03 10,46 12,21
24
Respuestas a los ejercicios complementarios 1) Empleando las ecuaciones de regresión lineal y cálculo de cifras de mérito expuestas con anterioridad, se obtienen los siguientes resultados respecto de la sensibilidad: Tabla 2.12. Sensibilidad de los métodos Método Sensibilidad de calibración A 1,552 B 153.0
Sensibilidad analítica 4,3×102 1,8×103
Observe que la sensibilidad de calibración tiene cifras significativas compatibles con su desvío estándar. En cambio, la sensibilidad analítica se informa con un número de cifras significativas que depende del cociente sensibilidad/ruido. Dado que el ruido se conoce con una o a lo sumo dos cifras significativas, la sensibilidad analítica se informa con dos cifras como máximo. Estos resultados indican que tanto la sensibilidad de calibración como la sensibilidad analítica es significativamente mayor para el método B. Sin embargo, la sensibilidad de calibración es dos órdenes de magnitud mayor para B, mientras que la sensibilidad analítica es superior, pero en menos de un orden de magnitud. La sensibilidad analítica es un mejor parámetro para la comparación. 2) Las concentraciones predichas para la incógnita y sus desvíos estándar, usando ambos métodos, son: Tabla 2.13. Concentraciones predichas Método Concentración (desvío estándar) A 0,153(1) B 0,1517(4) Como puede apreciarse, el desvío estándar calculado mediante el método B es menor, debido a su mayor sensibilidad analítica. Como comentario, la sensibilidad analítica parece comportarse mejor, en cuanto cifra de mérito, para calificar el desempeño de estos dos métodos, ya que se correlaciona con la precisión de cada cálculo de concentración. 3) Cifras de mérito en cada caso: Tabla 2.14. Cifras de mérito Caso Sensibilidad Sensibilidad analítica A 13,2 114,1 B 9,3 91,2
1/γ
LOD
LOQ
0,009 0,011
0,02 0,03
0,06 0,08
Rango lineal 0,06-0,95 0,08-1,31
Debe notarse que el caso A posee efectivamente un blanco significativo, puesto que la ordenada al origen es significativamente distinta de cero. En cuanto a las cifras de mérito, son algo mejores en el caso A, aunque el rango lineal es también sensiblemente menor.
25
La elección entre estos dos casos es un ejemplo de que no se puede tener todo en la vida: habría que decidir qué es más importante para aplicaciones concretas, si el rango lineal extendido o la mayor sensibilidad.
Referencias 1. K. Danzer y L. A. Currie, Guidelines for calibration in analytical chemistry. Part 1. Fundamentals and single component calibration, Pure & Appl. Chem. 1998, 70, 993-1014. 2. W. P. Gardiner, Statistical analysis methods for chemists. A software-based approach, The Royal Society of Chemistry, Cambridge, 1997. 3. J. N. Miller y J. C. Miller, Estadística y quimiometría para química analítica, 4ta. Edición, Prentice Hall, Madrid, 2002. 4. C. A. Clayton, J. W. Hines y P. D. Elkins, Detection limits with specified assurance probabilities, Anal. Chem. 1987, 59, 2506-2514. 5. L. A. Currie, Detection and quantification limits: origins and historical perspective, Anal. Chim. Acta 1999, 391, 127-134. 6. L. A. Currie, Recommendations in Evaluation of Analytical Methods including Detection and Quantification Capabilities, Pure Appl. Chem. 1995, 67 , 1699-1723. 7. P. Wilrich, ISO/DIS 11843-1,2 (1995), Capability of Detection, ISO/TC69/SC6, ISO Standard, 11843-1, 1977. 8. M. Valcárcel, Principios de química analítica, Springer-Verlag Ibérica, Barcelona, 1999, p. 81. 9. D. Harris, Análisis Químico Cuantitativo, 2da Ed. Reverté S.A., 2001, p. 106.
26
CAPÍTULO 3 REGRESIÓN LINEAL PARTE II: EXACTITUD Y COMPARACIÓN DE MÉTODOS ANALÍTICOS En este capítulo sobre regresión lineal exploraremos su uso para el análisis de la exactitud de un método analítico y para la comparación de dos métodos analíticos diferentes. La teoría se expone a continuación, pero se recomienda consultar paralelamente el ejemplo concreto que se analiza a continuación. La discusión que sigue está basada en trabajos recientes acerca del empleo de ensayos de recuperación para la validación y comparación de métodos, 1–3 así como en la obra clásica de Massart y colaboradores. 4 Para el estudio de la exactitud de un método analítico, es usual preparar una serie de patrones con concentraciones conocidas del analito de interés, diferentes a las utilizadas en la etapa de calibración. Luego se determina la concentración del analito en cada uno de ellos por interpolación en la recta de calibrado, y se analiza la exactitud de la determinación a través de la recuperación de las concentraciones nominales del analito. Por otro lado, cuando se desean comparar dos métodos analíticos, se determina, por ambos métodos, el contenido de un analito en una serie de muestras en las que su concentración es variable (dentro del rango lineal de cada uno de ellos). En ambos casos se trata de comparar parejas de valores que idealmente serían iguales, y estudiar el posible desvío de esta situación ideal, en un contexto estadístico y con un cierto nivel de confianza. Es por esta razón que ambos procedimientos se incluyen en el presente capítulo. Exactit ud de un método analítico Si se dispone de una serie de patrones de concentración conocida para la validación de un método analítico, se procede del modo siguiente. En primer lugar se miden sus respuestas, incluyendo réplicas de cada medición (usualmente cada patrón se mide por triplicado). Se estima la concentración a partir de cada respuesta analítica, se promedian los valores para cada nivel y se calcula el desvío estándar asociado. Luego se realiza una regresión lineal de los promedios en función de las concentraciones nominales a cada nivel. El análisis difiere en ciertas sutilezas respecto del realizado en el caso del Capítulo 2. La nomenclatura empleada aquí se describe a continuación: x indica la variable concentración nominal de cada nivel, y la variable concentración promedio predicha para las réplicas de cada nivel, n el número de réplicas, q el número de niveles de validación estudiados, y s( yi) el desvío estándar en la señal para cada nivel de concentración ( xi). Hay q desvíos estándar, dados por: n
s( yi) =
∑=1 ( y
ij
− y i ) 2
j
n −1
(3.1)
En la ecuación (3.1), yij indica la concentración para el patrón i en la réplica j, e yi es el promedio de las n réplicas para el nivel i. Debemos notar que una de las premisas para realizar un estudio por regresión lineal simple es que la varianza de la variable y sea aproximadamente constante, u 27
homoscedástica.4
La Figura 3.1 muestra las diferencias entre una varianza homoscedástica y otra heteroscedástica.
Figura 3.1. Arriba,
varianza homoscedástica; abajo, varianza heteroscedástica.
En la calibración de datos analíticos se supone que la distribución del ruido instrumental es constante a lo largo del rango de calibración, o en otras palabras, que la respuesta analítica es homoscedástica. Esto no es necesariamente así, sin embargo, si la variable y es la concentración predicha para patrones de validación, y no la respuesta analítica. Como se estudió antes, el desvío estándar en la concentración predicha mediante una recta de calibrado no es constante para diferentes muestras, sino que varía con la concentración del analito. Es decir que, en principio, la variable y que estamos considerando ahora no es homoscedástica. En estos casos, se recomienda realizar una regresión lineal mediante cuadrados mínimos ponderados (WLS, por weighted leastsquares) y no una regresión ordinaria (OLS, por ordinary least-squares ) como la empleada en el Capítulo 2.
28
Dado que el método WLS es más complicado que el OLS, lo recomendable es previamente verificar si efectivamente la varianza no es constante, para utilizar el primero en los casos en los que es estrictamente necesario. Una prueba de constancia de la varianza (o prueba de la homoscedasticidad) puede realizarse mediante el uso del parámetro estadístico F , calculando el valor "experimental" F exp definido por el cociente entre el máximo y el mínimo valor de las varianzas en las réplicas de los patrones [se toma como medida de cada varianza el valor de s( yi)2]: max[s( yi ) 2 ] F exp = min[s( yi ) 2 ]
(3.2)
Este valor se compara luego con el valor crítico de tablas para n – 1 y n – 1 grados de libertad para numerador y denominador, respectivamente (usualmente con el 95% de confianza). Si F exp > F crit entonces se recomienda calcular los parámetros A y B de la regresión con el método WLS que se describe más adelante. Región de conf ianza en el caso ho mosc edástic o Si se ha podido aplicar el método OLS descrito antes, debido a que las varianzas son aproximadamente constantes, se dispone de los valores ajustados de A y B y de sus desvíos estándar. Estos parámetros han sido utilizados tradicionalmente para determinar si las concentraciones estimadas de los patrones de validación se diferencian estadísticamente (o no), de las nominales. El procedimiento consistía en verificar si los valores ideales de A y B (1 y 0 respectivamente) estaban contenidos dentro de los correspondientes intervalos de confianza para la pendiente y ordenada al origen ajustadas. Sin embargo, actualmente se considera que este procedimiento es incorrecto, puesto que no tiene en cuenta que A y B no son variables estadísticamente independientes, y que siempre existe un cierto grado de correlación entre ellas. El procedimiento correcto debe considerar el intervalo de confianza conjunto entre la pendiente y la ordenada al origen. Este intervalo es una región en el plano de las dos variables (pendiente y ordenada al origen) que tiene forma elíptica. Por este motivo, la prueba estadística correcta consiste en investigar si el punto (1,0) está contenido en la región elíptica de confianza conjunta de la pendiente y la ordenada al origen. La prueba se conoce como EJCR (por elliptical joint confidence region). Específicamente, la región elíptica está descripta por la siguiente ecuación: 1 q(β − B) 2
q
+ 2(α − A)(β − B )∑ xi i =1
+ (α − A) 2
q
∑=1 x 2 = 2s 2 / i
y x F 2, q − 2
(3.3)
i
En la ecuación precedente, α y β son las variables que corresponden a las dos dimensiones del plano en que se representa la región elíptica, y F 2,q –2 es el valor del parámetro estadístico F con 2 y q – 2 grados de libertad para un dado nivel de confianza (usualmente 95%). Por lo tanto, debe dibujarse en un gráfico bidimensional la región anterior y verificar si contiene al punto (1,0). Detalles de cómo se dibuja esta elipse en un caso particular se dan en el ejercicio resuelto del documento que se acompaña. La Figura 3.2 ilustra este tipo de región para un caso típico: si el punto (1,0) no está contenido dentro de la elipse, esto implica que el método no es exacto. Es importante remarcar que el tamaño de la elipse, que está controlado, entre otros parámetros, por el desvío estándar de la regresión s y/x, da una idea de la precisión del método analítico que se está probando. En este sentido, es importante utilizar un número significativo de niveles de concentración para la prueba de exactitud, de manera que s y/x 29
sea representativo de la regresión. De lo contrario, si se emplean sólo unos pocos niveles de concentración, se corre el riesgo de que la elipse abarque un área considerable, e incluya al punto ideal (1,0) sólo por azar. Véase la Figura 3.3 para aclarar este punto. Nótese que el valor de s y/x en este caso es similar al parámetro usualmente empleado en la comparación de concentraciones predichas y nominales, llamado RMSE (por root mean square error ):
∑ ( y predicho − ynominal ) 2 RMSE =
q
(3.4)
Se divide el numerador por q (y no por q – 1) debido a que RMSE no es un desvío estándar, sino la raíz cuadrada de una media de desvíos. 0.2
0.2
n e g i r o l a a 0.0 d a n e d r O
n e g i r o l a a 0.0 d a n e d r O
-0.2
-0.2 1.0
1.1
1.0
Pendiente
1.1
Pendiente
Figura 3.2. Dos regiones elípticas de confianza conjunta. Izquierda, método exacto.
Derecha, método no exacto. El cuadrado marca el punto ideal (1,0).
30
■
n e g i r o l a a d a n e d r O
Pendiente
Figura 3.3.
Distintos tipos de elipses, de acuerdo con la exactitud y precisión: exacta y precisa; exacta e imprecisa; inexacta e imprecisa; inexacta y precisa. El cuadrado negro marca el punto ideal (1,0). Regresión pond erada Si los datos no cumplen con la prueba de homoscedasticidad, el análisis de los datos de validación debe hacerse mediante regresión lineal ponderada. En este caso se calculan la pendiente ( A) y ordenada al origen ( B) de la recta ajustada a la ecuación y = A x + B, minimizando la siguiente suma ponderada de cuadrados (SC): q
SC =
∑=1 w ( y i
i
− yˆ i ) 2
(3.5)
i
donde wi es el "peso" o "ponderación" aplicado a cada punto de la regresión, q el número de puntos, yi el valor de la variable y en cada punto (los promedios yi de las réplicas) e y es el promedio de los valores de la variable y. En el método OLS utilizado en calibración, la suma de cuadrados no incluye peso o ponderación alguna. Cuando los datos son heteroscedásticos, el peso wi se define como inversamente proporcional a la varianza de la variable en el punto i: wi =
1 s( yi ) 2
(3.6)
El efecto concreto del pesado de los datos en forma inversamente proporcional a su varianza es dar mayor contribución, en la regresión, a los datos más precisos y comparativamente menor peso a los menos precisos. Los valores estimados de A y B de una regresión lineal ponderada se calculan mediante las siguientes ecuaciones:
31
q
∑=1 w ( x i
A =
− x w )( yi − y w )
i
i
q
∑=1 w ( x i
i
(3.7)
− x w ) 2
i
B = y w – A x w
(3.8)
donde xi es la concentración de cada uno de los q patrones de validación, y los parámetros xw e y w son las coordenadas del centro de gravedad pesado por donde pasa la recta ajustada, que están dadas por: q
∑=1 w x
i i
xw =
i
q
(3.9)
(3.10)
∑=1 w
i
i
q
∑=1 w y i
y w =
i
i
q
∑=1 w
i
i
En el método WLS el parámetro s y/x (el desvío estándar de los residuos de la regresión) está dado por: q
∑=1 w ( y i
s y/x =
i
− yˆ i ) 2
i
q−2
(3.11)
donde yi es la respuesta experimental, e yˆ i representa la respuesta estimada en cada punto, esto es, yˆ i = A xi + B. El lector podrá comprobar que si todos los wi son idénticos entre sí (homoscedasticidad perfecta), las ecuaciones anteriores se reducen al caso OLS tratado en el Capítulo 2. Región de conf ianza en el caso heteroscedástico Cuando se aplica el método WLS para determinar A y B, la prueba de exactitud del método analítico es idéntica a la descrita en el caso OLS, excepto que la ecuación que describe la elipse de confianza conjunta es: (β − B) 2
q
q
∑=1 w + 2(α − A)(β − B)∑=1 w x i
i
i i
i
+ (α − A) 2
q
∑=1 w x 2 = 2s 2 / i i
y x F 2, q − 2
(3.12)
i
32
Comparació n de métodos analíticos La comparación de dos métodos se lleva a cabo disponiendo de una serie de muestras para las que se ha determinado el contenido de un analito por dos métodos alternativos. Usualmente se mide cada muestra por triplicado por ambos métodos, y se aplica un modelo de regresión lineal para verificar si los resultados provistos por ambos métodos son comparables. Cada muestra estudiada proporciona entonces una concentración predicha por cada uno de los dos métodos, acompañadas por sus respectivas varianzas. Supongamos que los resultados determinados por el método 1 se consideran la variable x y los provistos por el método 2 la variable y (en la comparación de un método dado frente a otro considerado como referencia, este último se toma como método 1). Ambas variables, por lo tanto, tienen asociada una incertidumbre finita. La regresión lineal de y vs. x difiere, en este caso, tanto del método OLS como del WLS, ya que en estos dos últimos la suposición básica es que no hay error en la variable x, aunque en realidad debería decirse que en OLS y WLS la incertidumbre asociada a la variable x (concentración nominal de patrones) es significativamente menor que la asociada a la variable y (respuesta analítica de los patrones, o concentración predicha por un dado método). Este supuesto no se cumple en la comparación de métodos analíticos, y es necesario recurrir a un método de regresión que tenga en cuenta los errores en ambos ejes. Un método popular para estos casos es el de cuadrados mínimos bivariados o BLS (por bivariate 6 least-squares). En la técnica BLS la pendiente y la ordenada al origen de la recta ajustada se obtienen minimizando una función idéntica a la mostrada en la ecuación (3.5), excepto que los pesos son una función de las varianzas en ambas variables: wi
= [s( yi ) 2 + A 2 s( xi ) 2 ]−1
(3.13)
En otras palabras, los pesos de la regresión "doblemente ponderada" BLS se eligen como inversamente proporcionales a una combinación de las varianzas en x y en y. Lamentablemente no existen fórmulas explícitas para estimar la pendiente y la ordenada al origen cuando los pesos tienen la forma dada por la ecuación (3.13), y debe recurrirse a un algoritmo matemático iterativo que no está disponible en los programas comerciales de ajuste por cuadrados mínimos. Esto es así porque en la ecuación (3.13) interviene la pendiente estimada A, que a su vez depende de los pesos. Sin embargo, hay ocasiones en que no es imprescindible aplicar el método BLS: cuando la varianza en la variable x es significativamente menor que en la variable y, la comparación puede realizarse con éxito empleando el método WLS, considerando que no hay error en la variable x. De hecho, si s( xi)2 << s( yi)2, la ecuación (3.13) se reduce al caso WLS en que wi = s( yi) –2. Por este motivo se aconseja asignar, para la regresión lineal, la variable x a los valores hallados por el método más preciso, y la variable y al método menos preciso. Si puede hacerse esta última aproximación, la comparación de métodos consiste en el cálculo de la pendiente y ordenada al origen mediante WLS, y consideración de la región elíptica de confianza conjunta, tal como se describió para el estudio de exactitud. Si el punto ideal (1,0) está contenido dentro de la elipse, los métodos son comparables estadísticamente en cuanto a la predicción de la concentración del analito en las muestras de validación. Se recomienda consultar la Referencia número 3, en el que se ilustran los peligros de no emplear el método correcto de regresión para la comparación de métodos 33
analíticos. También se discute el hecho de que en ciertos casos los métodos WLS y BLS pueden producir resultados similares, pero muy diferentes a los provistos por OLS. Programas de comp utación Usuarios de MATLAB: se provee acceso a la rutina EJCR.M que puede usarse para aplicar los métodos OLS, WLS y BLS, y generar la elipse correspondiente. Usuarios de QB: se provee acceso al programa EJCR.EXE, que realiza las operaciones necesarias pero no grafica la elipse. Esta última puede obtenerse importando los datos generados por el programa en un entorno gráfico apropiado. Véase también el ejercicio resuelto detalladamente que se acompaña. Material suministrado Para este capítulo se provee la siguiente siguiente información complemetaria: LECTURA ADICIONAL CAP 3.PDF, documento de Adobe con un trabajo educativo para lectura adicional. Archivos de texto (*.TXT) conteniendo datos típicos para estudios de exactitud y comparación de métodos. Archivos (*.M) con rutinas para el entorno de programación MATLAB. Archivos (*.EXE) con programas ejecutables en QB.
Ejercicio 1) La Tabla 3.1 muestra datos para analizar la exactitud de un método analítico. Determine si el método es exacto mediante regresión lineal y estudio de la región elíptica de confianza conjunta para A y B.
Tabla 3.1. Concentraciones nominales de patrones, y valores hallados por un método analítico (con sus desvíos estándar). Muestra Nominal Hallada Desvío estándar (promedio de cinco réplicas) 1 0,05 0,06 0,06 2 5,16 5,02 0,05 3 9,91 10,00 0,04 4 14,90 15,20 0,02 5 19,80 19,90 0,03 6 24,90 25,00 0,04 7 30,00 30,00 0,06 2) La Tabla 3.2 muestra datos para la comparación de dos métodos analíticos (promedios de tres réplicas en cada caso), incluyendo los desvíos estándar de cada uno. Compare los resultados mediante regresión WLS y análisis de la región elíptica conjunta.
34
Tabla 3.2. Concentraciones halladas por dos métodos analíticos con sus desvíos estándar. Muestra Método 1 Desvío Método 2 Desvío estándar estándar 1 0,05 0,03 0,06 0,06 2 5,16 0,02 5,02 0,05 3 9,91 0,02 10,00 0,04 4 14,90 0,01 15,20 0,02 5 19,80 0,02 19,90 0,03 6 24,90 0,01 25,00 0,04 7 30,00 0,03 30,00 0,06
Respuesta 1) En primer lugar debemos determinar si los datos de la Tabla 1 son homoscedásticos. Para ello calculamos el cociente: max[s( yi ) 2 ] (0,06) 2 = =9 min[s( yi ) 2 ] (0,02) 2
F exp =
Dado que este último valor es mayor que el de tabla [ F crit (95%,4,4) = 6,5] concluimos que los datos son heteroscedásticos, y que debemos emplear el método WLS para el análisis por regresión lineal. Calculamos entonces los pesos wi de cada dato, los que se reúnen en la Tabla 3. El cálculo de cada peso se realiza mediante la ecuación: qs( yi ) −2
wi =
q
∑=1 s( y ) −2 i
i
De esta manera, se consigue que la suma de los pesos sea igual a q, lo que facilita los cálculos. Tabla 3.3. Datos xi, yi y pesos wi para exactitud de métodos. i
xi
yi
wi
1 2 3 4 5 6 7
0,05 5,16 9,91 14,90 19,80 24,90 30,00
0,06 5,02 10,00 15,20 19,90 25,00 30,00
0,33 0,48 0,75 3,00 1,33 0,75 0,33
Note que los pesos son mayores para datos con menor desvío estándar. Para la muestra número 1, por ejemplo, tendremos: 7 (0,06) 2 = 0,33 w1 = 1 1 1 1 1 1 1 + + + + + + (0,06) 2 (0,05) 2 (0,04) 2 (0,02) 2 (0,03) 2 (0,04) 2 (0,06) 2 35
Luego debemos calcular los valores de los diferentes productos de variables y pesos, que se muestran en la Tabla 3.4. Tabla 3.4. Cálculos parciales para el método WLS. i
wi xi
wi xi2
wi yi
wi xi yi
1 2 3 4 5 6 7 Total
0,0167 2,4839 7,4538 44,8281 26,4756 18,7285 10,0287 110,0153
0,0008 12,8169 73,8671 667,9384 524,2178 466,3399 300,8596 2.046,0405
0,0201 2,4165 7,5215 45,7307 26,6094 18,8037 10,0287 111,1304
0,0010 12,4692 74,5380 681,3868 526,8653 468,2128 300,8596 2.064,3327
Con los resultados anteriores, calculamos: xw = 110,0153 / 7 = 15,72 y w = 111,1304 / 7 = 15,88 q
∑=1 w ( x i
A =
− x w )( yi − y w )
i
=
i
q
∑=1 w ( x i
i
− x w ) 2
i
q
∑=1 w x y i i
=
i
− v x w y w =
i
q
∑=1 w x 2 − v x 2 i i
2.064,3327 − 7 × 15,72 × 15,88 = 1,0022 2.046,0405 − 7 × (15,72) 2
w
i
B = y w – A x w = 15,88 – 1,0022
× 15,72 = 0,12
Estos valores deben acotarse al número correcto de cifras significativas conociendo los desvíos estándar correspondientes. Los desvíos estándar en la pendiente y la ordenada al origen, estimadas por el método WLS de regresión lineal, están dados por ecuaciones análogas a las empleadas en el método OLS, pero con los valores de x e y pesados convenientemente: s A =
s y / x Q xx
s B = s y / x
1 m
+
x w2 Q xx
donde s y/x se determina mediante la ecuación apropiada para datos pesados (WLS), tal como se describió en la parte teórica:
36
q
∑=1 w ( y i
s y/x =
i
− yˆ i ) 2 = 0,16
i
q−2
Por su parte, Q xx está dado por: q
Q xx =
∑=1 w x 2 − q x 2 = 316,2 i i
w
i
A partir de estos parámetros, se obtiene (redondeando a una cifra significativa): s A = 0,01 s B = 0,2
Por lo tanto, la pendiente y la ordenada al origen se informan como A = 1,00(1) y B = 0,1(2). Para el estudio de la región elíptica, necesitamos los siguientes parámetros: q = 7 q
∑=1 w x
i i
= 110,0153
i
q
∑=1 w x 2 = 2.046,0405 i i
i
s y2 / x
= 0,026
F 2,q −2
= 8,6
Por lo tanto, la ecuación de la elipse estará dada por: 7(β − 0,1) 2 + 220,0306(α − 1)(β − 0,1) + 2.046,0405(α − 1) 2 = 0,44
La ecuación anterior tiene la siguiente forma: a1 (α − A) 2
+ a 2 (α − A)(β − B) + a3 (β − B) 2 = a 4
donde a1, a2, a3, a4, A y B son constantes y α y β son las variables. Los valores de las constantes son: a1 = 2,046×103 a2 = 220,03 a3 = 7
37
a4 = 0,44 A = 1 B = 0,1
La ecuación describe una elipse en el plano ( α,β). Para dibujar esta elipse es necesario conocer sus límites en el eje de las abscisas ( α). Estos límites se pueden calcular a partir de las siguientes consideraciones. En primer lugar re-escribimos la ecuación anterior como de segundo grado en ( β – B): a3 (β − B ) 2
+ a 2 (α − A)(β − B) + [a1 (α − A) 2 − a 4 ] = 0
Luego calculamos los valores de ( β – B) a partir de la resolvente de segundo grado: (β – B) =
− a 2 (α − A) ±
a2
2
(α − A) 2 − 4a3 [a1 (α − A) 2 − a 4 ] 2a3
Observamos que sólo se obtendrán valores reales de ( β – B) si se cumple que la expresión dentro de la raíz cuadrada es positiva; los límites se encuentran cuando esta expresión se iguala a cero: a2
2
(α − A) 2 − 4a3 [a1 (α − A) 2 − a 4 ] = 0
de donde se pueden calcular los límites superior e inferior de ( α – A) como: LIM(α – A) = ±
4a3 a 4 = ± 0,0373 2 − a 2 + 4a3 a1
Para construir una tabla de valores de α y β, y graficar la elipse se calculan los correspondientes valores de β dentro de estos límites de α mediante la ecuación:
β = B +
− a 2 (α − A) ±
a2
2
(α − A) 2 − 4a3 [a1 (α − A) 2 − a 4 )] 2a3
Ejemplos de pares de valores de α y β calculados con la ecuación anterior son: Tabla 3.5. Valores de α y β calculados. α – A α –0,0373 0,9627 –0,0273 0,9727 –0,0173 0,9827 –0,0073 0,9927 0,0027 1,0027 0,0127 1,0127 0,0227 1,0227 0,0327 1,0327
β 0,7110 0,6971 0,5903 0,4563 0,3027 0,1306 –0,0642 –0,3022
0,6520 0,3516 0,1441 –0,0362 –0,1970 –0,3393 –0,4587 –0,5350 38
La gráfica de la elipse correspondiente, construida con datos de la tabla anterior, es la siguiente (el cuadrado sólido marca el punto ideal de pendiente 1 y ordenada 0):
1 ) β ( n e g i r o l a 0 a d a n e d r O
-1 0.96
0.98
1.00
1.02
1.04
1.06
Pendiente (α) Figura 3.4. Elipse construida con los datos de la Tabla 3.5
Se aprecia claramente que el punto ideal (1,0) está contenido en la elipse, por lo que el método analizado es exacto. Usuarios de MATLAB: los datos de la tabla están incluidos en el archivo de texto 'DATOS_EXACT_WLS.TXT', y organizados de tal modo que pueden estudiarse mediante la rutina de MATLAB 'EJCR.M', de la manera descrita en el Capítulo 2. Esta rutina proporciona los valores ajustados de pendiente y ordenada al origen, produce una figura con la correspondiente elipse, y genera un archivo de texto que contiene los valores numéricos necesarios para graficar la región elíptica mediante programas gráficos: la primera columna de este archivo cuenta con los valores de pendiente y la segunda y tercera los valores de ordenada al origen que corresponden a las dos mitades de la elipse. Usuarios de QB: los datos están en el archivo 'D_E_WLS.TXT' para ser estudiados por EJCR.EXE. 2) En este caso se trata de comparar dos métodos analíticos. Los resultados del análisis mediante WLS son idénticos a los discutidos para la parte 1) (¿porqué?). Cuando se realiza un análisis BLS se calculan los siguientes valores de pendiente y ordenada al origen: A = 1.00(1) B = 0,1(2)
39
Obsérvese que son idénticos a los hallados mediante la técnica WLS. La explicación es que los valores de la variable x (las concentraciones estimadas mediante el método analítico 1) tienen desvíos estándar menores que los de y (las concentraciones estimadas mediante el método analítico 2). Como consecuencia, es prácticamente lo mismo realizar el análisis mediante WLS o mediante BLS. Usuarios de MATLAB: los datos de la tabla están contenidos en el archivo de texto 'DATOS_COMPAR_BLS.TXT', y organizados de tal modo que pueden estudiarse mediante la rutina de MATLAB 'EJCR.M', de la manera descrita en Capítulo 2. Esta rutina proporciona los valores ajustados de pendiente y ordenada al origen, produce una figura con la correspondiente elipse y genera un archivo de texto que contiene los valores numéricos necesarios para graficar la región elíptica mediante programas gráficos: la primera columna de este archivo contiene los valores de pendiente y la segunda y tercera los valores de ordenada al origen que corresponden a las dos mitades de la elipse. Usuarios de QB: los datos están en D_C_BLS.TXT. Ejercicios complementarios 1) Los valores siguientes conciernen a la comparación entre las predicciones efectuadas para la determinación de teofilina en sangre mediante un método espectrofotométrico, comparado con un método de inmunofluorescencia polarizada (FPIA). No se determinaron las muestras por triplicado debido a la cantidad insuficiente de muestra (sueros de pacientes pediátricos). Sin embargo, se estima que los desvíos estándar promedio para cada método son: 0.4 μg ml−1 para el método FPIA y 0.9 μg ml−1 para el espectrofotométrico. Llevar a cabo el análisis de comparación de métodos mediante la construcción de la elipse apropiada, suponiendo que los desvíos estándar anteriores son constantes para todos los datos.
40
Tabla 3.6. Concentraciones de teofilina obtenida por ambos métodos. Muestra Teofilina hallada / μg ml−1 FPIA Espectrofotométrico 1 0.0 1.4 2 6.5 5.3 3 33.2 30.6 4 9.7 12.7 5 12.2 14.9 6 14.8 17.7 7 20.1 19.9 8 15.6 18.5 9 19.3 20.4 10 16.8 22.6 11 24.2 27.1 12 28.6 29.8 13 0.0 0.0 14 3.9 1.6 15 8.0 5.7 16 11.2 14.2 17 11.4 15.3 18 14.7 17.5 19 16.5 17.6 20 16.6 19.4 21 19.8 18.7 22 19.5 18.9 23 23.0 21.2 2) En la determinación del antibiótico ciprofloxacina en orina se emplean tres métodos multivariados diferentes. La Tabla 3.7 proporciona datos para estudiar la exactitud de cada método, frente a un grupo de muestras de referencia, cuya concentración de analito es conocida. Grafique las correspondientes EJCR y comente los resultados. Note que no hay datos disponibles acerca de los desvíos estándar, por lo que deberá realizarse un análisis OLS.
41
Tabla 3.7. Concentraciones de ciprofloxacina obtenida por tres métodos. Muestra Nominal Método 1 Método 2 214 190 173 1 80 86 2 87 29 23 26 3 6 14 13 4 28 38 19 5 142 145 150 6 33 16 26 7 60 67 8 58 146 126 125 9 67 63 10 65 89 92 90 11 172 158 12 160 52 48 41 13 68 64 14 75 11 0 10 15 8 5 16 0 7 0 3 17 7 11 18 0
Método 3 208 107 46 28 50 160 47 80 146 75 120 174 61 92 26 21 30 27
Se recomienda emplear la rutina de MATLAB 'EJCR.M' (o su equivalente en QB) organizando los datos del ejercicio propuesto de la manera que se presenta en los archivos de texto correspondientes al ejercicio resuelto.
Respuesta a los ejercicios complementarios 1) La tabla de datos debe complementarse con la de los desvíos estándar. En este caso, dado que el desvío estándar para FPIA es menor que para el método espectrofotométrico, podría emplearse un análisis de tipo WLS, con los valores de desvío estándar igual a 0,9 para todos los datos de la tabla anterior. Esto último, sin embargo, es idéntico al uso de un método OLS (ver la teoría). Por lo tanto, podemos en este caso particular realizar una regresión lineal ordinaria empleando como variable y los valores provistos por el método espectrofotométrico y como variable x los provistos por el método FPIA. Los resultados del análisis OLS son: Pendiente: 0,983 Ordenada al origen: 1,35 s y/x: 2,35 La elipse correspondiente contiene, aunque marginalmente, al punto ideal (1,0):
42
Figura 3.6. Elipse construida por el
método OLS
Vale la pena destacar el resultado que se obtendría mediante un análisis BLS, esto es, considerando que tanto la variable x como la y están sujetas a incertidumbre: Pendiente: 0.996 Ordenada al origen: 1.16 s y/x: 2.39 Como puede apreciarse en la figura siguiente, el resultado final en cuanto al estudio de la comparación de los métodos es similar al hallado mediante el análisis OLS sencillo.
Figura 3.7. Elipse construida por el método BLS
43
La rutina de MATLAB 'EJCR.M', proporciona los valores ajustados de pendiente y ordenada al origen, produce una figura con la correspondiente elipse, y genera un archivo de texto que contiene los valores numéricos necesarios para graficar la región elíptica mediante programas gráficos: la primera columna de este archivo contiene los valores de pendiente y la segunda y tercera los valores de ordenada al origen que corresponden a las dos mitades de la elipse. 2) Se requiere graficar tres elipses, calculadas por OLS, que proporcionan visualmente una buena impresión de la exactitud y precisión relativas de los tres métodos probados: 30
3
n e 20 g i r o l a a 10 d a n e d r 0 O
1 2
-10 0.8
0.9
1.0
1.1
Pendiente Figura 3.8. Elipse construida para llos
3 métodos estudiados.
La conclusión es que el método más preciso es el 3 (menor tamaño de elipse), pero es muy poco exacto (alejado del punto ideal). El método 2 es el más exacto, y además es más preciso que el método 1.
Referencias 1. A. G. González, M. A. Herrador y A. G. Asuero, Intra-laboratory testing of method accuracy from recovery assays, Talanta 1999, 48 , 729-736. 2. H. C. Goicoechea, A. C. Olivieri, A. Muñoz de la Peña, Determination of theophylline in blood serum by UV spectrophotometry and partial least-squares (PLS-1) calibration, Anal. Chim. Acta 1999, 384, 95-103. 3. V. Franco, V. E. Mantovani, H. C. Goicoechea, A. C. Olivieri, Teaching chemometrics with a bioprocess: analytical methods comparison using bivariate linear regression, Chem. Educator 2002, 7, 265-269. 4. D. L. Massart, B. G. M. Vandeginste, L. M. C. Buydens, S. De Jong, P. J. Lewi y J. Smeyers-Verbeke, Handbook of Chemometrics and Qualimetrics, Elsevier, Amsterdam, 1997, Capítulo 8. 5. Los términos homoscedástico/a y homoscedasticidad existen en el contexto del "Diccionario Estadístico" que puede consultarse en http://www.estadistico.com/dic.html. También se usan, en forma equivalente, homocedástico/a y homocedasticidad. 6. J. Riu y F. X. Rius, Assessing the accuracy of analyical methods using linear regression with errors in both axes, Anal. Chem. 1996, 68 , 1851-1857. 44
CAPÍTULO 4 CALIBRACIÓN BIVARIADA Determinació n de dos analitos u sando dos senso res En la calibración multivariada, se emplean datos instrumentales medidos utilizando más de un sensor para la determinación simultánea de dos o más analitos, o de un analito en presencia de interferentes. El ejemplo típico de datos multisensoriales es un espectro de absorción electrónica, donde la señal instrumental es la absorbancia, y los sensores son las longitudes de onda. Sin embargo, las técnicas de calibración multivariada no están restringidas al uso de datos espectrales de un tipo determinado, sino que pueden extenderse a otros datos tales como fluorescencia, absorción en el infrarrojo o infrarrojo cercano, y aún datos no espectroscópicos, como voltamperogramas. El caso más simple del uso de más de un sensor para la determinación de más de un analito es el estudio de mezclas binarias de compuestos absorbentes a dos longitudes de onda, o calibración bivariada. Si bien la técnica ha caído un tanto en desuso para aplicaciones prácticas, conserva no obstante un importante valor pedagógico, ya que permite una introducción sencilla y gradual al tema más complejo de determinaciones de multianalitos utilizando un alto número de sensores. 1,2 La teoría se expone a continuación, pero se recomienda consultar paralelamente el ejemplo concreto que se analiza al final. A medida que se discutan los conceptos teóricos asociados con el uso de dos sensores, se establecerán las analogías correspondientes con la calibración multivariada utilizando múltiples sensores. La etapa de calibración Análogamente al caso univariado, la metodología bivariada consta de dos etapas: calibración y predicción. En la etapa de calibración, se requiere establecer la relación existente entre concentración y señal para cada analito calibrado, del mismo modo que cuando se estima la pendiente de una recta univariada. En el presente caso, la diferencia estriba en que las señales multivariadas son intrínsecamente menos selectivas, y es necesario un procedimiento que distinga, de algún modo, las señales que le corresponden a cada analito. Una etapa de calibración típica en análisis bivariado consiste en preparar soluciones de concentración conocida de ambos analitos, y estimar, a partir de las señales medidas a dos longitudes de onda (o, en general, a dos sensores diferentes), las respectivas relaciones señal-concentración. Las soluciones de calibrado pueden ser mezclas de ambos analitos, o más simplemente soluciones conteniendo los analitos en forma pura (si esto es experimentalmente posible). En la sección siguiente se describirá en detalle el proceso de calibración empleando la notación matricial, lo que preparará en cierta forma el camino para el uso de este tipo de herramientas matemáticas en el análisis multivariado. La calibración en notación matricial Es sumamente útil emplear la notación matricial para indicar los resultados de mediciones de mezclas binarias a dos longitudes de onda. Quienes deseen revisar conceptos básicos sobre matrices y sus operaciones, necesarios para seguir en detalle la presente capítulo, pueden consultar el Anexo I. Supongamos que en la calibración de un método bivariado se preparan dos soluciones patrón conteniendo los analitos 1 y 2, y se leen las absorbancias de estas dos 45
soluciones a las longitudes de onda 1 y 2. Las correspondientes respuestas instrumentales Y ij (absorbancias de la solución patrón i a la longitud de onda j) se reúnen en la matriz (2 ×2) de calibración Y:
⎡Y 11 Y=⎢ ⎣Y 21
Y 12 ⎤
⎥
Y 22 ⎦
(4.1)
Las concentraciones de ambos analitos en las soluciones de calibrado deben conocerse a los efectos de llevar a cabo la calibración del modelo. Estas concentraciones se agrupan en la matriz de concentraciones de calibración (2 ×2) X, cuyo elemento genérico X in es la concentración en la mezcla i del analito n: X=
⎡ X 11 ⎢ X ⎣ 21
X 12 ⎤
⎥
X 22 ⎦
(4.2)
La etapa de calibración, o sea, la determinación de las llamadas sensibilidades individuales a cada longitud de onda, se lleva a cabo suponiendo que se cumple la ley de Beer que relaciona absorbancia con concentración. La señal de la mezcla número 1 a la longitud de onda 1, por ejemplo, se obtiene a partir de la suma de las contribuciones de ambos analitos: Y 11 = X 11 S 11 + X 12 S 12
(4.3)
donde S 11 y S 21 son las sensibilidades del analito 1 y 2 respectivamente a la longitud de onda 1. Los restantes elementos de Y se obtienen mediante ecuaciones similares a (4.3). En general Y podrá escribirse entonces mediante el siguiente producto matricial: T
Y = X S
(4.4)
donde S es una matriz (2 ×2) cuyo elemento genérico S jn es la sensibilidad a la longitud de onda j del analito n. La ecuación (4.4) puede representarse en forma gráfica mediante la Figura 4.1, útil para analizar los requerimientos de tamaño de las distintas matrices (se ilustra un caso general, en que el número de analitos es N , el número de muestras de calibrado es I y el número de longitudes de onda es J ).
46
Y
I × J
S
X
I × N
=
Estos números deben coincidir ( I )
×
Estos números deben coincidir ( N )
T
N × J
Estos números deben coincidir ( J )
Figura 4.1. Esquema que muestra las relaciones de tamaño en la aplicación de la
ley
de Beer a mezclas de componentes. Si X se expresa en términos de concentraciones molares, entonces S jn es la absortividad molar a la longitud de onda j del componente n (multiplicada por el paso óptico). Sin embargo, se prefiere llamar a los elementos S jn "sensibilidades", dado que el modelo matemático no está restringido a datos de absorción. Nótese que se requiere la trasposición de la matriz S en la ecuación (4.4) para mantener la consistencia del producto matricial (Figura 4.1). La matriz S puede obtenerse a partir de la ecuación (4.4), aunque se requieren varias etapas. En primer lugar, es necesario pre-multiplicar ambos miembros por X –1 (observe que en el terreno matricial, pre-multiplicar, o multiplicar por la izquierda, no es lo mismo, en general, que pos-multiplicar o multiplicar por la derecha). A continuación es necesario trasponer ambos miembros de la igualdad obtenida, con el objeto de "despejar" la matriz S: S = (X –1 Y)T = YT (X –1)T
(4.5)
En la ecuación (4.5), al trasponer el producto matricial, se invierte el orden de aparición de las matrices. Esta última ecuación completa la calibración, lo que provee una matriz de calibración S que debe almacenarse para predicciones en muestras futuras. La obtención de S es análoga al cálculo de la pendiente de la recta de regresión en calibración univariada, en forma previa a la medición de la señal analítica de muestras incógnita. En el procedimiento bivariado más simple que se pueda concebir, las soluciones de calibración no son mezclas binarias, sino que contienen sólo analitos puros, de manera que la matriz X tiene en este caso el siguiente aspecto: X=
⎡ X 11 ⎢ X ⎣ 21
⎡ x1,cal = ⎥ ⎢ X 22 ⎦ ⎣ 0 X 12 ⎤
0 ⎤ ⎥ x2,cal ⎦
(4.6)
donde se han empleado las cantidades escalares x1,cal y x2,cal para denotar las concentraciones de los analitos 1 y 2 respectivamente en las soluciones de calibrado. Sin embargo, en un caso más general, podrían utilizarse mezclas binarias para calibración. En tal caso, es importante recalcar que las concentraciones de los patrones 47
empleados en estas mezclas binarias deben ser tales que la matriz X pueda invertirse. Este requisito es fundamental, ya que de lo contrario será imposible calcular la matriz S a través de la ecuación (4.5). Matemáticamente hablando, X podrá invertirse si su determinante es distinto de cero, y esto sucederá si sus líneas (filas o columnas) no son combinaciones lineales. Desde el punto de vista químico, esto se traduce en que las concentraciones de un analito no deben ser proporcionales a las del otro analito en las mezclas empleadas para calibrar. Este concepto cobrará mayor importancia en el
campo del análisis de múltiples analitos empleando espectros completos. Nótese que, en el caso más simple, si se cumple la ecuación (4.6), X es diagonal y por lo tanto la matriz inversa X –1 adopta una forma sencilla, de manera que S está dada por:
S = YT (X –1)T =
⎡Y 11 ⎢Y ⎣ 12
⎡ 1 Y 21 ⎤ ⎢ x1,cal ⎥⎢ Y 22 ⎦ ⎢ 0 ⎢⎣
⎤ 0 ⎥ ⎥ = 1 ⎥ x2,cal ⎥⎦
⎡ Y 11 ⎢ x ⎢ 1,cal ⎢ Y 12 ⎢ x1,cal ⎣
⎤ x2,cal ⎥ ⎥ Y 22 ⎥ x2,cal ⎥⎦ Y 21
(4.7)
donde puede reconocerse cada elemento de S como la absortividad molar de cada analito a cada longitud de onda (multiplicada por el paso óptico), obtenida dividiendo la correspondiente absorbancia de la solución de calibración por su concentración. Etapa de predicción En la etapa de predicción, se miden las señales instrumentales para una muestra incógnita, por ejemplo, dos absorbancias a las longitudes de onda a las que se realizó la calibración. Dichas señales, denominadas y1 e y2, se agrupan en el vector columna (2 ×1) y:
⎡ y1 ⎤ y=⎢ ⎥ ⎣ y2 ⎦
(4.8)
La predicción se logra recurriendo a la ley de Beer aplicada a la muestra incógnita, en forma análoga a la ecuación (4.4): y = S x
(4.9)
donde x es un vector columna que contiene los elementos buscados en el análisis: las concentraciones (desconocidas) de ambos analitos en la incógnita. Despejando x de la ecuación (4.9): –1
x = S
y
(4.10)
El vector x (2 ×1) contiene dos elementos: las concentraciones de ambos analitos en la muestra incógnita estimadas por el modelo bivariado. Estos cálculos completan, por lo tanto, la etapa de predicción. Coefic ientes de regresión La ecuación (4.10) puede interpretarse en forma gráfica mediante el siguiente esquema:
48
x1
1ra. fila de S –1 y
x2
2da. fila de S –1 =
x
2×1
–1
S
2 ×2
×
y
2 ×1
Figura 4.2.
Esquema que muestra las relaciones de tamaño en en el cálculo de concentración de un analito. Este esquema nos ayuda a establecer que la concentración de cada analito en la muestra incógnita se predice mediante el siguiente producto escalar: –1
xn = (nava fila de S
) × y
(4.11)
La nava fila de S –1, una vez traspuesta (o sea, convertida en un vector columna) cumple un papel importante en el análisis multivariado, donde corrientemente se denomina n: n
= (nava fila de S –1)T
(4.12)
La trasposicion de la nava fila de S –1 para obtener el vector columna n es una cuestión puramente formal. Con esta última definición, la ecuación (4.11) se transforma en: xn =
n
T
y = β 1n y1 + β 2n y2
(4.13)
lo cual significa que la concentración predicha es el producto escalar del vector n por el vector de respuestas instrumentales. En el terreno multivariado n se llama vector de los coeficientes de regresión. Este concepto es sumamente importante, ya que los coeficientes de regresión actúan de manera análoga a la inversa de la sensibilidad (o pendiente de la recta de regresión univariada), permitiendo definir cifras de mérito para modelos que emplean más de un sensor. Cuando se emplean múltiples sensores para la determinación de multianalitos, las ecuaciones serán similares a las arriba descritas. Los respectivos tamaños de las matrices se comparan en la Tabla 4.1, en la que I es el número de muestras de calibrado, J el número de longitudes de onda analizadas y N el número de componentes presentes en las mezclas.
49
Tabla 4.1. Tamaños en las matrices en determinaciones utilizando múltiples sensores. Matriz / Concepto Modelo bivariado Modelo multivariado vector Señales de calibrado Y 2×2 I × J Concentraciones de X 2×2 I × N calibrado S Sensibilidades 2×2 J × N Coeficientes de J ×1 2×1 n regresión Señales de la incógnita y 2×1 J ×1 Concentraciones x N ×1 2×1 estimadas en la incógnita
Colinealidad Un examen cuidadoso de la ecuación (4.10) revela que el paso crítico en la estimación de la concentración de los analitos en la muestra incógnita es la inversión de la matriz S. Una matriz es invertible si su determinante es distinto de cero, de lo contrario se dice que la matriz en cuestión es singular o no invertible, y su inversa no existe. De esto se desprende que si por alguna razón el determinante de la matriz S, aunque no sea exactamente cero, es pequeño (en comparación con el nivel de ruido instrumental), S será difícilmente invertible, en el sentido que los elementos de S –1 estarán pobremente definidos. El proceso puede ilustrarse con la inversión de un número cercano a cero; si bien el cociente existe, su valor se vuelve más y más impreciso a medida que el número decrece. La singularidad de la matriz S, y su consiguiente dificultad de inversión, están directamente relacionadas con el concepto de paralelismo o colinealidad espectral. En términos matemáticos, el determinante de S será cercano a cero si sus filas son combinaciones lineales. Podemos traducir este concepto al campo de la espectroscopía bivariada si graficamos la sensibilidad para cada analito a cada una de las dos longitudes de onda de trabajo (Figura 4.2). Uniendo los puntos correspondientes a cada analito se obtienen dos líneas rectas: cuanto más paralelas sean estas líneas rectas, más difícil será la inversión de S, y más cercano a cero su determinante. Véase la Figura 4.2, en la que se muestra una situación deseable (izquierda) y una indeseable (derecha) para el análisis a dos longitudes de onda.
50
d a d i l i b i s n e S
d a d i l i b i s n e S
λ1
λ2
λ1
λ2
Figura 4.3.
Sensibilidad en función de la longitud de onda. Izquierda: situación deseable, derecha: situación indeseable.
Cifras de mérit o Análogamente al caso univariado, pueden definirse cifras de mérito correspondientes a determinaciones usando múltiples sensores. Respecto de la sensibilidad, nótese que la ecuación (4.13) puede interpretarse como una forma particular de la ley de Beer, en la que la concentración es proporcional a la señal. Dado que la constante de proporcionalidad en este caso es la inversa de la sensibilidad (en el caso univariado c = (εb) –1 A], es natural pensar en el vector de coeficientes de regresión como midiendo la "sensibilidad inversa" para una determinación a múltiples longitudes de onda. De hecho, la definición de sensibilidad SENn para el analito n en una determinación de dos analitos en mezclas binarias a dos longitudes de onda es: SENn =
1 β12n + β 22 n
(4.14)
donde 1n y 2n son los elementos del vector n. Esta definición es análoga a la sensibilidad de calibración, que fuera definida en el contexto de la calibración univariada. Dado que la cantidad β12n + β 22 n es la "longitud" del vector n, también conocida como su norma, la ecuación para la sensibilidad adopta la siguiente forma: SENn = 1 / ||
n
||
(4.15)
donde || · || simboliza el cálculo de la norma de un vector. En la calibración para el análisis de varios componentes debe considerarse como cifra de mérito adicional la selectividad. En el caso univariado este parámetro no se toma en cuenta debido a que aquél no contempla la existencia de señales interferentes. Se puede definir la selectividad para el analito n, en presencia de otros componentes, como el cociente entre la sensibilidad dada por la ecuación (4.15), y el valor que tendría dicha sensibilidad si el analito en cuestión estuviese presente en su forma pura: SELn = SENn / || nava. columna de S ||
(4.16) 51
Puede demostrarse que SEL n es un número adimensional que varía entre 0 y 1; el cero corresponde a un sistema totalmente no selectivo para el analito n, mientras que 1 corresponde al caso totalmente selectivo, es decir específico, para el que se puede aplicar la calibración univariada. También puede definirse la sensibilidad analítica, como el cociente entre el valor de SENn y el ruido instrumental s R, obtenido a partir de replicados de una muestra blanco:
γn = SENn / s y
(4.17)
Existen también ecuaciones para la estimación de los errores estándar en la concentración predicha de cada analito, que son una extensión de la estudiada en el caso univariado. Una aproximación sencilla a la estimación de s( xn) se puede obtener ignorando el efecto de la leva, y tomando sólo en consideración el efecto de la incertidumbre en la respuesta analítica. En ese caso: s( xn) = s y / SENn
(4.18)
En relación con el límite de detección, el cálculo se complica por el hecho de que este parámetro no puede definirse para un analito sin conocer la concentración de otros analitos en una muestra dada. El lector interesado en una lectura avanzada respecto de la estimación de errores estándar y límite de detección en este caso puede consultar el documento adjunto "LECTURA ADICIONAL CAP 4.PDF", que dará una idea de la complejidad matemática del problema, aún en el caso simple de calibración bivariada. De todas maneras, si los efectos de la leva no son relevantes, en otras palabras, si la muestra incógnita no está lejos del centro de la calibración, una ecuación aproximada para el límite de detección puede obtenerse por analogía con la calibración univariada: LOD = 3,3 s y / SENn
(4.19)
Material suministrado Para este capítulo se provee la siguiente información complemetaria: LECTURA ADICIONAL CAP 4.PDF, documento de Adobe con un trabajo educativo para lectura adicional.
Ejercicio 1) Se desean analizar mezclas acuosas de permanganto y dicromato mediante mediciones a dos longitudes de onda. En la etapa de calibración, se preparan dos soluciones patrón conteniendo, respectivamente, permanganato de potasio 4,075 ×10 –4 M y dicromato de potasio 3,030 ×10 –3 M. Determinar la matriz de sensibilidades S para esta calibración. Las absorbancias medidas a 440 y 545 nm de estas soluciones se indican a continuación: Composición de la muestra KMnO4 4,075×10 –4 M K 2Cr 2O7 3,030×10 –3 M
Absorbancia a 440 nm 0,028 1,073
Absorbancia a 545 nm 0,915 0,026
2) Se miden las absorbancias de cuatro muestras incógnita a las mismas longitudes de onda de calibración, resultando los siguientes valores: 52
Tabla 4.2. Absorbancias de cuatro muestras incógnita. Muestra Absorbancia a 440 nm 1 0.364 2 0.722 3 0.153 4 0.607
Absorbancia a 545 nm 0.220 0.258 0.915 0.937
Determinar la concentración de cada analito en estas muestras. 3) Estimar las cifras de mérito del método para cada analito. Suponga que la incertidumbre típica en la señal es de 0,003 unidades de absorbancia.
Respuesta 1) Para calcular la matriz S se requiere conocer las matrices Y y X. Con los datos del ejercicio es sencillo escribir estas dos últimas matrices: ⎡0,028 0,915⎤ Y = ⎢ ⎥ ⎣1,073 0,026⎦ 0 ⎤ –4 ⎡4,075 ⎢ 0 ⎥ ×10 30 , 30 ⎣ ⎦ Luego debemos aplicar la ecuación para el cálculo de S: T –1 T S = Y (X ) Para ello, es necesario en primer lugar invertir la matriz X. Esto es sencillo, puesto que X es diagonal, de manera que: X =
−1
0 ⎤ ⎡4,075 ⎡2.454 0 ⎤ 4 X = ⎢ × 10 = ⎢ 0 30,30⎥⎦ 330⎥⎦ ⎣ 0 ⎣ Luego se multiplican las traspuestas de Y y X –1: T T ⎡0,028 0,915⎤ ⎡2.454 0 ⎤ ⎡0,028 1,073 ⎤ ⎡2.454 0 ⎤ T –1 T S = Y (X ) = ⎢ = = ⎥ ⎢ ⎢0,915 0,026⎥ ⎢ 0 330⎥⎦ 330⎥⎦ ⎣1,073 0,026⎦ ⎣ 0 ⎣ ⎦⎣ –1
⎡ 69 354⎤ =⎢ ⎥ ⎣2.245 9 ⎦ 2) Para estimar la concentración en una muestra incógnita, necesitamos la matriz inversa de S. Esta se puede calcular fácilmente recurriendo al cálculo matricial estándar, resultando en: *
*
Recuerde que la inversa de una matriz se obtiene mediante la ecuación S –1 = det(S) –1 Cof(S)T, donde det(S) es el determinante de S, y Cof(S) es la matriz cofactor, cuyo elemento i, j se obtiene multiplicando (–1)i+ j por el determinante de la matriz menor que resulta de eliminar, de la matriz original S, la fila i y la T 9 − 2.245⎤ ⎡ 5 T columna j. En el caso estudiado, det(S) = –7,95×10 , Cof(S) = ⎢ = ⎣− 354 69 ⎥⎦
53
⎡− 0,11 4,46 ⎤ = ⎢ × 10 –4 ⎥ ⎣28,27 − 0,87⎦ Las concentraciones de ambos analitos en las cuatro muestras incógnita son, por lo tanto, las siguientes: –1
S
Muestra 1: –1
x = S
y=
⎡− 0,11 4,46 ⎤ –4 ⎢28,27 − 0,87⎥ × 10 × ⎣ ⎦
⎡0,364⎤ ⎢0,220⎥ = ⎣ ⎦
⎡0,94⎤ 4 ⎢ 10,1 ⎥ × 10 ⎣ ⎦
Muestra 2: –1
x = S
y=
⎡− 0,11 4,46 ⎤ –4 ⎢28,27 − 0,87⎥ ×10 × ⎣ ⎦
⎡0.722⎤ ⎢0.258⎥ = ⎣ ⎦
⎡1,07 ⎤ 4 ⎢20,2⎥ × 10 ⎣ ⎦
⎡− 0,11 4,46 ⎤ –4 ⎢28,27 − 0,87⎥ ×10 × ⎣ ⎦
⎡0.153⎤ ⎢0.915⎥ = ⎣ ⎦
⎡4,06⎤ 4 ⎢ 3,5 ⎥ × 10 ⎣ ⎦
⎡− 0,11 4,46 ⎤ –4 ⎢28,27 − 0,87⎥ ×10 × ⎣ ⎦
⎡0.607⎤ ⎢0.937⎥ = ⎣ ⎦
⎡4,11⎤ 4 ⎢16,3⎥ × 10 ⎣ ⎦
Muestra 3: –1
x = S
y=
Muestra 4: –1
x = S
y=
La Tabla 4.3 resume los resultados. Tabla 4.3. Resultados obtenidos para cada muestra. Muestra 1
2
3
4
Señal
Concentración predicha (error estándar) M × 104
⎡0,364⎤ ⎢0,220⎥ ⎣ ⎦
x1 = 0,94(1)
⎡0,722⎤ y = ⎢ ⎥ ⎣0,258⎦
x1 = 1,07(1)
y =
⎡ 0,153⎤ ⎢0,915⎥ ⎣ ⎦
x1 = 4,06(1)
y =
⎡0,607⎤ ⎢0,937⎥ ⎣ ⎦
x1 = 4,11(1)
y =
− 354⎤ ⎡ 9 ⎢⎣− 2.245 69 ⎥⎦ , y por lo tanto, ⎡− 0,11 4,46 ⎤ ×10 –4. ⎢⎣ 28,27 − 0,87 ⎥⎦
x2 = 10,1(1)
x2 = 20,2(1)
x2 = 3,5(1)
x2 = 16,3(1)
S –1
= (–7,95×105) –1 ×
− 354⎤ ⎡ 9 ⎢⎣− 2.245 69 ⎥⎦ =
54
Los errores estándar en las concentraciones se calcularon con el modelo aproximado citado en la teoría descripta en este capítulo, esto es s( xn) = s R / SENn, con s R = 0,003 (para el cálculo de SEN n ver el punto 3). Nótese que, debido a la diferencia en sensibilidades, los errores estándar para el analito 1 son de un orden de magnitud menores que para el analito 2. 3) Cifras de mérito. La sensibilidad para cada analito se obtiene calculando la norma de la fila correspondiente de la matriz S –1: Para el analito 1: ⎡ − 1,1 × 10 −5 ⎤ 1 = ⎢ ⎥ ⎣4,46 × 10 −4 ⎦ SEN1 = || 1 || –1 = 2.243 A M –1 (A = unidades de absorbancia) ⎡2,8 ×10 −3 ⎤ 2 = ⎢ ⎥ ⎣ − 1×10 −4 ⎦ SEN2 = || 2 || –1 = 354 A M –1 Las selectividades, por otro lado, se pueden estimar dividiendo las respectivas sensibilidades por la longitud de la columna correspondiente a cada analito: SEL1 = SEN1 / || columna 1 de S || = 0,998 SEL2 = SEN2 / || columna 2 de S || = 0,998 Tener en cuenta que para llegar a estos últimos valores es necesario incluir varias cifras significativas en los cálculos intermedios. Ejercicio complementario Se desea realizar un análisis cuantitativo de una mezcla binaria realizando medidas de absorbancia a dos longitudes de onda. Se dispone de datos de absorbancia para soluciones patrón de cada analito a varias longitudes de onda, según se muestra en la siguiente tabla: Tabla 4.4. Datos de absorbancia a diferentes longitudes de onda Solución λ1 λ2 λ3 λ4 patrón Analito 1 0,550 0,610 0,720 0,850 –4 1,00×10 M Analito 2 0,510 0,505 0,710 0,800 –4 1,00×10 M
λ5 0,910 0,800
Se requiere seleccionar dos longitudes de onda para realizar el análisis, y un criterio para ello es utilizar aquellas que provean la máxima selectividad para cada analito. ¿A qué dos longitudes de onda se obtiene la mayor selectividad?
Respuesta Deben calcularse las selectividades para cada analito para todas las combinaciones posibles de longitudes de onda. Haremos el cálculo detallado para el caso de elegir λ1 y λ2 : Matrices conteniendo las señales y las concentraciones de los patrones: ⎡0,550 0,610⎤ ⎡1 0⎤ –4 R = ⎢ C = ⎥ ⎢0 1⎥ ×10 0 , 510 0 , 505 ⎣ ⎦ ⎣ ⎦ Cálculo de S y su inversa: 55
0 ⎤ ⎡5.500 5.100⎤ ⎡0,550 0,510⎤ ⎡1×10 4 ) = ⎢ = ⎢ ⎥⎢ 0 ⎥ 4⎥ 0 , 610 0 , 505 1 × 10 ⎣ ⎦⎣ ⎦ ⎣6.100 5.050⎦ ⎡− 15,14 15,29 ⎤ –1 S = ⎢ × 10 –4 ⎥ ⎣ 18,29 − 16,49⎦ Sensibilidades y selectividades: ⎡− 1,514 × 10 −3 ⎤ = 1 ⎢ 1,529 × 10 −3 ⎥ ⎣ ⎦ –1 SEN1 = || 1 || = 464; SEL1 = 464 / (5.500 2 + 6.100 2) = 0,056 ⎡ 1,829 ×10 −3 ⎤ 2 = ⎢ ⎥ ⎣− 1,649 ×10 −3 ⎦ SEN2 = || 2 || –1 = 406; SEL2 = 406 / (5.100 2 + 5.050 2) = 0,056 T
–1 T
S = R (C
Realizando este mismo análisis a todas las posibles combinaciones de dos longitudes de onda se obtienen los siguientes resultados: Tabla 4.5. Uso de distintas combinaciones de longitudes de onda Combinación SEL1 1y2 0,056 1y3 0,029 1y4 0,007 1y5 0,028 2y3 0,084 2y4 0,059 2y5 0,027 3y4 0,022 3y5 0,056 4y5 0,034
SEL2 0,056 0,029 0,007 0,028 0,084 0,059 0,027 0,022 0,056 0,034
Como puede verse, la mejor combinación de longitudes de onda, es la 2 y 3, que conduce a la máxima selectividad.
Referencias 1. G. D. Christian, Analytical Chemistry, 6a. Edición, Wiley, New York, 2003, Capítulo 16. 2. D. A. Skoog, D. M. West y F. J. Holler, Fundamentals of Analytical Chemistry, 7a. Edición, Saunders College Publishing, New York, 1996, Capítulo 20. 3. M. A. Cabezón, A. C. Olivieri, Precision in Two-Wavelength Spectroscopic Analysis of Binary Mixtures, Chem. Educator 2004, 9, 288-296.
56
CAPÍTULO 5 CALIBRACIÓN DIRECTA
MULTIVARIADA
PARTE
I:
CALIBRACIÓN
Determinación de multianalitos u sando múl tiples sensores En esta sección extenderemos los resultados presentados en el capítulo anterior al análisis de varios analitos mediante múltiples sensores. La analogía más directa del método bivariado es el llamado análisis por cuadrados mínimos clásicos o CLS (por classical least-squares ). Recomendamos especialmente la lectura del notable trabajo de Haaland sobre el tema. 1 Otras lecturas valiosas son los capítulos correspondientes de libros de quimiometría,2,3 así como los ya famosos " tutorials" de Brereton en internet. 4 La teoría se expone a continuación, pero se recomienda consultar paralelamente el ejemplo concreto que se analiza al final. El modelo CLS en notación matrici al: etapa de calibració n Continuaremos empleando la notación matricial para indicar los resultados de mediciones a varias longitudes de onda. En el caso que se desee llevar a cabo la determinación simultánea de varios analitos, es preciso preparar mezclas de patrones de dichos analitos, como mínimo en un número igual al de analitos. En general, sin embargo, se prefiere utilizar un conjunto de mezclas de calibrado compuesto por un número de mezclas mayor que el de analitos, debido a que de este modo se obtienen resultados más precisos, así como en calibración univariada se emplean varios patrones para determinar un único analito. Esto plantea inmediatamente el problema de cuáles deben ser las concentraciones de los analitos en las mezclas de calibrado, problema que se designa, en términos generales, como del diseño experimental de las mezclas. La teoría del diseño experimental queda más allá del alcance del presente curso; sólo podemos adelantar, recurriendo al sentido común, que las mezclas de calibrado deben ser representativas, en todo lo posible, de las combinaciones de concentraciones de los analitos que se espera encontrar en las mezclas incógnita. Cuántas mezclas y qué concentraciones es parte de los detalles del diseño experimental. En esta sección estudiaremos un caso simple, en el que se determinan simultáneamente dos analitos, tratando de extender las ecuaciones, donde sea posible, a la existencia de N analitos. Supongamos que se preparan varias ( I ) soluciones patrón de los analitos 1 y 2 puros, y se leen las absorbancias de estas I soluciones a J diferentes longitudes de onda. Las correspondientes respuestas instrumentales Y ij (absorbancias de la solución patrón i a la longitud de onda j) se reúnen en la matriz ( I × J ) de calibración Y:
⎡Y 11 ⎢Y 21 Y=⎢ ⎢ ... ⎢ ⎣Y I 1
Y 12 Y 22
... Y I 2
... Y 1 J ⎤ ... Y 2 J ⎥⎥ ... ... ⎥ ⎥ ... Y IJ ⎦
(5.1)
Las concentraciones de los analitos en las I soluciones de calibrado deben conocerse, tal como en el análisis a dos sensores. Aquellas se agrupan en la matriz de concentraciones de calibración ( I × N ) X, cuyo elemento genérico X in es la concentración en la mezcla i del analito n: 57
⎡ X 11 ⎢ X 21 X=⎢ ⎢ ... ⎢ ⎣ X I 1
X 12 X 22
... X I 2
... X 1 N ⎤ ... X 2 N ⎥⎥ ... ... ⎥ ⎥ ... X IN ⎦
(5.2)
Para dos analitos, la ecuación anterior se transforma en:
⎡ X 11 ⎢ X 21 X=⎢ ⎢ ... ⎢ ⎣ X I 1
X 12 ⎤
⎥ ⎥ ... ⎥ ⎥ X I 2 ⎦ X 22
(5.3)
o sea, una matriz de I ×2. La etapa de calibración, o sea, la determinación de las llamadas sensibilidades individuales a cada longitud de onda, se realiza suponiendo que se cumple la ley de Beer que relaciona absorbancia con concentración, análogamente al caso de dos longitudes de onda. Sin embargo, debe tenerse en cuenta en este caso el problema está sobredimensionado. Esto significa que el problema puede plantearse como un conjunto de ecuaciones simultáneas en el que el número de ecuaciones disponible es superior al de incógnitas. En nuestro caso, se desea relacionar la concentración con la señal a través de la sensibilidad S jn a la longitud de onda j del analito n. Si se trata de dos analitos, hay J ×2 parámetros a determinar (los valores de todos los coeficientes S jn), y un total de I × J ecuaciones; dado que en general I > 2, el problema está sobredimensionado. En estos casos el criterio que se aplica es el de obtener la solución de cuadrados mínimos, esto es, aquella que minimice el error E del siguiente modelo: Y = X ST + E
(5.4)
donde S es una matriz ( J × N ) cuyo elemento genérico S jn es la sensibilidad a la longitud de onda j del analito n. Nótese que se requiere la trasposición de la matriz S en la ecuación (5.4) para mantener la consistencia del producto matricial. Las relaciones de tamaño entre las matrices de la ecuación (5.4) se muestran en la Figura 5.1. La solución de cuadrados mínimos de la ecuación (5.4) corresponde a la obtención de la matriz S a partir de esta última, fijando E = 0 (una matriz de ceros del mismo tamaño que Y). La obtención de S a partir de la ecuación (5.4) no puede hacerse simplemente premultiplicando por X –1, dado que X no es, en general, una matriz cuadrada, y matrices no cuadradas no pueden invertirse. Para despejar S se recurre, en primer lugar, a premultiplicar ambos miembros de la ecuación (5.4) por la matriz traspuesta de X: T
T
T
X Y = X X S
(5.5)
58
T
Y
I × J
X
I × N
=
E
S
×
N × J
+
I × J
Figura 5.1. Esquema que muestra las relaciones de tamaño en la aplicación de la
ley
de Beer a mezclas de multicomponentes. Nótese que hemos fijado E = 0 en la ecuación (5.4) antes de realizar esta operación. El producto (XT X) es una matriz cuadrada (tamaño N × N ), y pre-multiplicando por su inversa ambos miembros de la ecuación (5.5): ST = (XT X) –1 XT Y
(5.6)
Trasponiendo la ecuación anterior para obtener S: T
–1
S = [(X X)
T
T
X Y]
= YT X (XT X) –1
(5.7)
La ecuación (5.7) merece varios comentarios. En primer lugar, es necesario recalcar que para que esta ecuación tenga sentido, debe poder invertirse la matriz cuadrada ( XT X). La inversión de una matriz requiere que sus líneas (filas o columnas) no sean linealmente dependientes, esto es, combinaciones lineales unas de otras. En el ejemplo que estamos analizando, esto implica, desde el punto de vista químico, que las concentraciones del analito 1 y el analito 2 en las mezclas no estén correlacionadas (por ejemplo, que no aumenten linealmente de una mezcla a otra). Diseñar un conjunto de mezclas con mínima correlación es también parte de la teoría del diseño experimental. El segundo comentario proviene de comparar la ecuación (5.7) con su análogo del Capítulo 4, en que X era cuadrada y podía invertirse directamente. Esta comparación sugiere que en la ecuación (5.7), la matriz [ X ( XT X ) –1] funciona como "una especie de inversa" de X (traspuesta, para ser más exactos). En la literatura se la ha llamado "inversa generalizada de X" o simplemente "seudoinversa de X", simbolizándola por +† X . Con esta nomenclatura, la ecuación (5.7) puede escribirse en forma más compacta: T
S = Y
(X+)T
(5.8)
†
El nombre "seudoinversa" tiene mayores implicancias en quimiometría que las discutidas aquí. En el caso de que (XT X) sea imposible o difícil de invertir, por ejemplo, porque su determinante es cero o cercano a cero, la seudoinversa aún existe, aunque la inversa generalizada no.
59
Esta última ecuación completa la calibración, lo que provee una matriz de calibración S para predicciones en muestras futuras. La obtención de S es análoga al cálculo de la absortividad molar en calibración univariada, en forma previa a la medición de la señal analítica de muestras incógnita. Como resumen de la etapa de calibrado podemos consignar los siguientes requerimientos: • El modelo CLS necesita un diseño de calibrado apropiado. • La calibración del modelo requiere conocer las concentraciones de los componentes de las mezclas de calibración. Quienes deseen revisar algunos conceptos sobre matrices y sus operaciones, que son necesarios para seguir en detalle el presente capítulo, pueden consultar el Anexo I. Etapa de predicc ión y c oeficientes de regresión En la etapa de predicción, una muestra incógnita produce J valores de la señal instrumental, por ejemplo, J absorbancias a las longitudes de onda a las que se realizó la calibración. Estas respuestas instrumentales se agrupan en el vector columna ( J ×1) y:
⎡ y1 ⎤ ⎢ y ⎥ 2 y =⎢ ⎥ ⎢ ... ⎥ ⎢ ⎥ ⎣ y J ⎦
(5.9)
La predicción se logra recurriendo a la ley de Beer aplicada a la muestra incógnita, en forma análoga a la ecuación (5.7): y = S x + e
(5.10)
donde x es un vector columna que contiene dos elementos: las concentraciones de ambos analitos en la incógnita, y e es un vector que recoge los errores del modelo lineal. Nuevamente se emplea el criterio de mínimos cuadrados para despejar x de la ecuación (5.11) (fijando e = 0). En primer lugar se debe pre-multiplicar la ecuación (5.10) por ST, de manera que se obtenga una matriz cuadrada en el segundo miembro: T
T
S y = (S S) x
(5.11)
Luego puede despejarse x pre-multiplicando por la inversa de ( ST S): T
–1
x = (S S)
T
S y
(5.12)
Nuevamente, podemos definir la seudoinversa de S de tal modo que permita obtener x directamente, pre-multiplicando a y: +
x = S y
(5.13)
60
nava fila de S
xn
+
=
x N ×1
S N × J
+
y
×
y J ×1
Figura 5.2.
Esquema que muestra las relaciones de tamaño en el cálculo de la concentración. El esquema superior muestra que la ecuación (5.13) puede interpretarse diciendo que la concentración de cada analito se predice mediante el siguiente producto escalar: xn = (nava fila de S
+
) × y
(5.14)
La nava fila de S+, una vez traspuesta (convertida en un vector columna) se conoce como el vector de los coeficientes de regresión para el componente n, n: n
= (nava fila de S+)T
(5.15)
Con esta última definición, la ecuación (5.14) se transforma en: xn =
n
T
y = β 1n y1 + β 2n y2 + ... + β Jn y J
(5.16)
lo cual significa que la concentración es el producto escalar del vector de coeficientes de regresión por el vector de respuestas instrumentales. Cifras de mérit o Análogamente al caso univariado, pueden definirse cifras de mérito correspondientes a determinaciones usando múltiples sensores. Respecto de la sensibilidad, nótese que la ecuación (5.16) puede interpretarse como una forma particular de la ley de Beer, en la que la concentración es proporcional a la señal. Dado que la constante de proporcionalidad en este caso es la inversa de la sensibilidad (en el caso univariado c = ( εb) –1 A], es natural pensar en el vector de coeficientes de regresión como midiendo la "sensibilidad inversa" para una determinación a dos longitudes de onda. De hecho, la definición de sensibilidad para cada analito en una determinación de dos analitos en mezclas binarias a dos longitudes de onda es: SENn = donde
jn son
1 β12n
+ β 22 n
+ ... + β 2
Jn
los elementos del vector
SENn = 1 / ||
n
(5.17)
n
.
||
(5.18)
donde || · || simboliza el cálculo de la norma de un vector. 61
Se puede definir la selectividad para el analito n, en presencia de otros componentes, como el cociente entre la sensibilidad dada por la ecuación (5.20), y el valor que tendría dicha sensibilidad si el analito en cuestión estuviese presente en su forma pura: SELn = SENn / || nava. columna de S || (5.19) Puede demostrarse que SEL n es un número adimensional que varía entre 0 y 1; el cero corresponde a un sistema totalmente no selectivo para el analito n, mientras que 1 corresponde al caso totalmente específico, para el que se puede aplicar la calibración univariada. También existe la sensibilidad analítica, que puede definirse como el cociente entre el valor de SENn y el ruido instrumental s y, obtenido a partir de replicados de una muestra blanco:
γn = SENn / s y
(5.20)
Existen también ecuaciones para la estimación de los errores estándar en la concentración predicha de cada analito, que son una extensión de la estudiada en el caso univariado. Una aproximación sencilla a la estimación de s(cn) se puede obtener ignorando el efecto de la leva, y tomando sólo en consideración el efecto de la incertidumbre en la respuesta analítica. En ese caso: s( xn) = s y / SENn
(5.21)
En relación con el límite de detección, el cálculo se complica por el hecho de que este parámetro no puede definirse para un analito sin conocer la concentración de otros analitos en una muestra dada. De todas maneras, si los efectos de la leva no son relevantes, en otras palabras, si la muestra incógnita no está lejos del centro de la calibración, una ecuación aproximada para el límite de detección puede obtenerse por analogía con la calibración univariada: LOD = 3,3 s y / SENn
(5.22)
El lector interesado en una lectura avanzada respecto de la estimación de errores estándar y límite de detección en este caso, puede consultar el trabajo correspondiente a la referencia número 5. Debe mencionarse también que en el marco de los modelos del tipo CLS puede obtenerse un parámetro típico de los ajustes por cuadrados mínimos: los residuos de la regresión. En el presente caso se trata del vector e de la ecuación (5.10), que contiene la incertidumbre asociada con el modelado de la señal de la muestra. Es importante calcular, para cada muestra incógnita, el desvío estándar de los residuos sres: J
∑=1 (e ) 2 j
sres =
j
J − N
(5.23)
donde e j representa cada uno de los elementos del vector e. Nótese el empleo de J – N grados de libertad en la ecuación (5.23), en atención a que la señal de la muestra proporciona J datos (las señales medidas a las J longitudes de onda), y se estiman N parámetros (las concentraciones de los N analitos en la muestra). Finalmente, es importante llevar a cabo una validación del modelo de calibrado, preparando un juego de muestras independientes, en el que los analitos estén presentes 62
en concentraciones distintas de las empleadas para calibrar el modelo, pero dentro de sus respectivos rangos lineales. La comparación de las concentraciones estimadas para este juego de validación con las nominales se lleva a cabo convenientemente mediante la prueba de la elipse discutida en el Capítulo 3. El ejercicio resuelto que acompaña este documento ilustrará el uso de los parámetros comentados en esta sección. Colinealidad espectral Análogamente al análisis de dos analitos a dos longitudes de onda, la presencia de colinealidad espectral en el modelo CLS se manifiesta a través de la dificultad en encontrar la seudoinversa S+. Específicamente, si los espectros de los analitos son colineales en un grado significativo, será difícil encontrar la inversa ( ST S) –1, y las concentraciones de los analitos estarán pobremente definidas. El resultado será una disminución en la sensibilidad, y a través de la ecuación (5.21), un aumento considerable del error de predicción. Interferentes no modelados Hemos supuesto, hasta el momento, que una muestra incógnita no debe poseer componentes que no estén presentes en la calibración, y que produzcan señal a las longitudes de onda de trabajo. En efecto, la suposición básica del método univariado es su especificidad completa. Análogamente, en el análisis multisensorial se requiere que la muestra incógnita esté compuesta por los mismos componentes que se utilizaron para calibrar. En el modelo CLS podemos sin embargo plantearnos, por primera vez, qué sucedería si una muestra incógnita estuviese compuesta por sustancias no presentes en la calibración. La respuesta es que se produciría un error significativo en la predicción, básicamente porque la ecuación (5.10) no sería correcta. En esta última ecuación, se supone que sólo existen los analitos calibrados en la muestra incógnita. Si bien es cierto que no es posible pretender que CLS estime las concentraciones correctamente en un caso como este, no es menos cierto que el modelo es capaz de "avisar" al analista que esto está ocurriendo. En un caso de interferencias no modeladas, los elementos del vector e de la ecuación (5.10) serán anormalmente grandes en relación con el nivel de ruido instrumental. De esta manera, los modelos que operan con múltiples sensores y ajustes por cuadrados mínimos son capaces de proveer información acerca de la presencia de interferentes no modelados, y a pesar de que son incapaces de corregirlos, al menos pueden informar al operador de estas anomalías. Ventajas y desvent ajas de CLS Podemos resumir las principales ventajas del modelo CLS del siguiente modo. Por un lado, se trata de un modelo matemáticamente sencillo, que puede seguirse convenientemente con el auxilio del cálculo matricial estándar, y aún mediante planillas de cálculo o programas fácilmente accesibles que realicen ajustes por cuadrados mínimos. Por otro lado, si el tipo de muestra a analizar no presenta interferencias serias de componentes desconocidos, o no se encuentran colinealidades espectrales significativas entre los analitos, el análisis CLS provee una manera rápida, simple y confiable de estimar las concentraciones en muestras de multicomponentes en forma simultánea. Las desventajas del modelo son fácilmente imaginables: es sensible a la presencia de colinealidad espectral, de manera que analitos con espectros severamente solapados no pueden estudiarse mediante esta técnica. Además, es necesario conocer los componentes químicos presentes en las mezclas incógnitas, de lo contrario, la presencia de interferentes no modelados producirá un error serio en la determinación. 63
Comparación d e métodos Vale la pena en este punto detenerse a reflexionar sobre las diferentes técnicas de calibración que hemos estudiado, y efectuar un análisis comparativo. Las propiedades analíticas que nos interesa comparar son: • Habilidad para analizar más de un analito en forma simultánea. • Conocimiento de las concentraciones de los componentes de la calibración. • Efectos de la colinealidad espectral. • Presencia de interferencias no modeladas en la calibración. Tabla 5.1. Comparación de las propiedades analíticas de los distintos métodos de calibración. Propiedad Método Univariado Bivariado Multivariado CLS Número de 1 2 Varios analitos Concentración de Conocida Conocidas Conocidas componente(s) de calibrado Efecto de la – Disminuye la Disminuye la colinealidad sensibilidad, sensibilidad, selectividad y selectividad y precisión precisión Presencia de Análisis inexacto Análisis Análisis interferentes inexacto inexacto pero con detección del problema Cifra de mérito Sensibilidad SEN = A SENn = 1 / || n || Sensibilidad γ = SEN / s y γ = SENn / s y analítica s( x) = s( xn) ≈ s y / SENn Incertidumbre en la predicción s y / x 1 1 ( xinc − x ) 2 + + A
Límite detección
de
Límite de cuantificación
n
m
Q xx
LOD = 3,3 s y / x 1 1 x 2 + + A 3 m Q xx LOQ = 10 s y / x 1 1 x 2 + + 3 m Q xx A
LOD ≈ 3,3 s y / SENn
LOQ ≈ 10 s y / SENn
En la Tabla 5.1 hemos resumido estas propiedades para los tres métodos analizados hasta el momento: univariado, bivariado y multivariado CLS. Hemos incluido, además, las definiciones de cifras de mérito más usadas en cada caso. Como puede verse, el pasaje del análisis univariado al multivariado CLS representa el logro de beneficios progresivos, respecto del número de analitos que pueden estudiarse simultáneamente, y de la detección de la presencia de interferentes. En relación con el conocimiento de las concentraciones de los componentes de las mezclas de calibración el comportamiento de los tres métodos es similar, al igual que la 64
respuesta al efecto de la colinealidad (no aplicable al caso univariado). En el próximo capítulo analizaremos un método capaz de superar estas dificultades y describiremos sus propiedades en perspectiva con las de la Tabla 5.1. Material suministrado Para este capítulo se provee la siguiente información complementaria: Archivos de texto (*.TXT) conteniendo datos típicos. Archivos (*.M) con rutinas para el entorno de programación MATLAB. Archivos (*.EXE) con programas ejecutables en QB.
Ejercicio 1) Se miden las señales instrumentales de cuatro soluciones de calibrado para dos analitos, a seis longitudes de onda distintas. La matriz de calibrado tiene la siguiente forma: Tabla 5.2. Datos de absorbancia de la matriz de calibrado. Muestra de λ1 λ2 λ3 λ4 calibrado 1 1,52 2,78 3,32 3,26 2 2,94 5,42 6,33 5,48 3 1,47 2,94 3,81 4,21 4 3,01 5,63 6,80 6,35
λ5
λ6
2,48 3,35 4,08 5,10
1,94 2,78 3,15 4,03
Las concentraciones de los dos analitos en las muestras se muestran en la Tabla 5.3:
Tabla 5.3. Concentraciones de los analitos en las muestras de calibración Ejercicio 1 x1cal x2cal Muestra de calibrado 1 2 3 4
1,00 2,00 1,00 2,00
1,00 1,00 2,00 2,00
Construir las matrices X e Y para calibrado, y calcular la matriz S de sensibilidades y los coeficientes de regresión. Informar las correspondientes cifras de mérito para el modelo. Suponga que el nivel de ruido instrumental es igual a 0,03 unidades de señal. 2) Se estudia un conjunto de cuatro muestras de validación, para las que se conocen las concentraciones nominales de ambos analitos. Las señales obtenidas a las mismas longitudes de onda que el calibrado, y las respectivas concentraciones se muestran en las tablas siguientes:
65
Tabla 5.4. Valores de las señales medidas en el Ejercicio 2. Muestra de λ1 λ2 λ3 λ4 validación 1 1,32 2,59 3,36 3,62 2 2,73 5,10 5,95 5,23 3 1,38 2,54 2,92 2,47 4 1,25 2,42 3,06 3,18
λ5
λ6
3,55 3,35 1,48 2,90
2,72 2,72 1,24 2,29
Tabla 5.5. Concentraciones de los analitos en las muestras de calibración del Ejercicio 2. Muestra de validación x1val x2val 1 2 3 4
0,89 1,86 0,93 0,83
1,69 1,05 0,40 1,34
Estimar las concentraciones de los analitos en este juego de muestras y estudiar la exactitud del método mediante la prueba de la elipse. 3) Analizar mediante el modelo CLS anterior tres muestras de prueba, para las cuales se han medido las siguientes señales a las mismas seis longitudes de onda que la calibración. Prestar atención a los residuos espectrales, ya que se sospecha que en una de estas tres muestras está presente una especie no modelada en la matriz de calibración. Tabla 5.6. Valores de las señales medidas en el Ejercicio 3. Muestra de λ1 λ2 λ3 λ4 prueba 1 1,11 2,20 2,77 2,81 2 5,54 6,71 7,02 5,83 3 2,56 4,76 5,81 5,56
λ5
λ6
2,56 3,66 4,39
2,02 2,77 3,50
Respuesta 1) Para calcular la matriz S se requiere conocer las matrices Y y X. Con los datos del ejercicio es sencillo escribir estas dos últimas matrices: ⎡1,52 2,78 3,32 3,26 2,48 1,94 ⎤ ⎢2,94 5,42 6,33 5,48 3,35 2,78⎥ ⎥ Y = ⎢ ⎢1,47 2,94 3,81 4,21 4,08 3,15 ⎥ ⎢ ⎥ ⎣ 3,01 5,63 6,80 6,35 5,10 4,03⎦
66
⎡1 ⎢2 X = ⎢ ⎢1 ⎢ ⎣2
1⎤ 1 ⎥⎥ 2⎥ ⎥ 2⎦
Luego debemos aplicar la ecuación para el cálculo de S: T T –1 S = Y X (X X) Estas operaciones son sumamente tediosas, aún para unas pocas longitudes de onda, y es preferible realizarlas con la ayuda de un programa. Para ello, los datos de calibración están organizados en los archivos de texto XCAL_E_R.TXT (concentraciones) e YCAL_E_R.TXT (señales), y pueden ser analizados convenientemente por los programas CLS_CAL.M (MATLAB) o CLS_CAL.EXE (QB). Ambos graban un archivo de texto conteniendo la matriz S. Los coeficientes de regresión pueden obtenerse a partir de las filas de la matriz S+ = (ST S ) –1 ST. Tanto el programa en MATLAB como en QB generan un archivo de texto conteniendo estos vectores de coeficientes de regresión para cada analito. Las figuras siguientes muestran en forma gráfica los espectros de calibrado y la matriz S, así como los coeficientes de regresión que serán luego útiles para la etapa de predicción.
Figura 5.3: Espectros de calibrado correspondientes al Ejercicio 1.
67
Figura 5.4: Sensibilidades y coeficientes de regresión correspondientes al Ejercicio
1. Las cifras de mérito calculadas mediante los programas para este modelo son las siguientes: Tabla 5.7. Cifras de mérito Cifra de mérito Analito 1 Analito 2 –1 Sensibilidad 4,1 Señal × concentración 1,9 Señal × concentración –1 Sensibilidad analíticaa 137 concentración –1 63 concentración –1 Selectividad 0,83 0,83 a Obtenida dividiendo la sensibilidad por el nivel de ruido instrumental (0,03 unidades). Nótese que la selectividad es idéntica para ambos analitos. En el caso de mezclas de más componentes esto no es necesariamente así. Puede apreciarse que el modelo es más sensible al analito 1 que al 2, hecho que también se ilustra en forma gráfica en la Figura 5.2. 2) Las concentraciones de ambos analitos en las cuatro muestras de validación están dadas por: + x = S y Estos cálculos pueden realizarse con ayuda de los programas CLS_PRED.M (MATLAB) o CLS_PRED.EXE (QB), organizando los datos de manera apropiada. El archivos de texto YVAL_E_R.TXT contiene la matriz de las señales de estas muestras de validación en la forma apropiada. Los resultados de la validación son los siguientes: 68
Tabla 5.8. Predicciones durante la validación Muestra Analito 1 1 2 3 4
Nominal 0,89 1,86 0,93 0,83
Predicho 0,88(1) 1,87(1) 0,93(1) 0,84(1)
Analito 2 Nominal 1,69 1,05 0,40 1,34
Predicho 1,70(1) 1,05(1) 0,40(1) 1,34(1)
Residuo espectral 0,01 0,03 0,01 0,02
Los errores estándar en las concentraciones se calcularon con el modelo aproximado citado en la teoría, esto es s( xn) = s y / SENn, con s y = 0,03. Se informan también, en la última columna de esta tabla, los residuos espectrales para cada muestra incógnita, que, como puede apreciarse, se encuentran dentro del nivel del ruido instrumental. Esto confirma que el ajuste por cuadrados mínimos para estas muestras es adecuado. Para establecer la exactitud del método, lo recomendado es analizar los datos de la tabla precedente mediante la prueba de la elipse, tal como se discutiera en el Capítulo 3. En los casos multivariados se recomienda también producir una única elipse, que recoja la comparación de las concentraciones nominales y predichas para todos los analitos. De este modo, la tabla de datos a suministrar a los programas de cálculo de la elipse será como sigue: Tabla 5.9. Datos a suministrar a los programas de cálculo de la elipse. 0,89 0,88 1,86 1,87 0,93 0,93 0,83 0,84 1,69 1,70 1,05 1,05 0,40 0,40 1,34 1,34 Dado que no se tienen resultados de réplicas de cada muestra, lo que proveería una estimación del desvío estándar de cada valor predicho, realizaremos un análisis mediante el método OLS. Se recomienda organizar los datos apropiadamente y someterlos a los programas EJCR.M (MATLAB) o EJCR.EXE (QB). El resultado se muestra en la figura siguiente, donde puede apreciarse la exactitud del método.
69
Figura 5.5: Elipse corerespondiente al ajuste de los datos de la Tabla
5.9.
El programa para el cálculo de la elipse también provee el error medio de la predicción: RMSE = 0,003 Este valor puede considerarse como sumamente satisfactorio en vista de las cifras significativas asignadas a los valores nominales de concentración, tanto de calibrado como de validación. 3) Los datos para las muestras de prueba están contenidos en el archivo de texto YPRU_E_R.TXT. Los resultados para las muestras de prueba son los siguientes: Tabla 5.10. Resultados obtenidos para el Ejercicio 3. Muestra de Analito 1 Analito 2 prueba Predicho Predicho 1 0,76(1) 1,17(1) 2 2,52(1) 0,62(1) 3 1,71(1) 1,78(1)
Residuo espectral 0,02 1,00 0,02
Evidentemente, la muestra número 2 posee una interferencia no modelada, causante de un mal ajuste. Las concentraciones de los analitos predichas para esta muestra no son confiables. Lamentablemente, el modelo CLS no puede resolver este problema, pero al menos informa al analista de su presencia. Ejercicio complementario 1) Se han recogido espectros de absorción electrónica de mezclas de dos colorantes a 281 longitudes de onda diferentes, para un conjunto de calibración compuesto por 9 muestras de calibración. Estos datos se proveen en el archivo de texto RESP_CAL.TXT, en forma de una matriz de 281 ×9. Las respectivas concentraciones (en ppm) están contenidas, en forma de matriz de 9 ×2, en el archivo de texto CONC_CAL.TXT. Los detalles experimentales de este trabajo están informados en la referencia 5. Lleve a cabo la calibración mediante el modelo CLS con el programa adecuado e informe las cifras de mérito. Suponga un nivel de ruido instrumental de 0,005 unidades de absorbancia. 70
2) También se midieron las señales de tres muestras de prueba, cuyos espectros están contenidos, en forma de matriz de 281 ×3, en el archivo de texto RESP_TST.TXT. Estimar las concentraciones de los dos analitos en estas muestras, y sus respectivos desvíos estándar.
Respuesta 1) Los resultados provistos por el programa son los siguientes: Figura con espectros de calibración:
Figura 5.6: Espectros de calibrado correspondientes al Ejercicio complementario
71
Figura con sensibilidades y coeficientes de regresión:
Figura 5.7: Sensibilidades y coeficientes de regresión correspondientes al
Ejercicio complementario. Estas figuras son provistas automáticamente por la rutina de MATLAB; los usuarios de QB pueden producirlas con cualquier programa gráfico, leyendo los datos de los correspondientes archivos de texto 'RESP_CAL.TXT' (espectros de calibración), 'S_.TXT' (sensibilidades) y 'B_.TXT' (coeficientes de regresión). Las cifras de mérito son las siguientes: Tabla 5.11. Cifras de mérito correspondientes al Ejercicio complementario. Cifra de mérito Analito 1 Analito 2 –1 Sensibilidad 0,50 A × ppm 0,76 A × ppm –1 Sensibilidad analíticaa 100 ppm –1 152 ppm –1 Selectividad 0,62 0,62 a Obtenida dividiendo la sensibilidad por el nivel de ruido instrumental (0,005 unidades de A). 2) Los resultados para las muestras de prueba son los siguientes: Tabla 5.12. Predicciones correspondientes al Ejercicio complementario. Muestra de Analito 1 Analito 2 prueba Predicho Predicho 1 2,00(1) 2,01(1) 2 1,02(1) 1,06(1) 3 4,05(1) 3,91(1)
Residuo espectral 0,009 0,01 0,007 72
La rutina de MATLAB provee la gráfica de los espectros de prueba:
Figura 5.7: Espectros de las
muestras correspondientes al Ejercicio complementario.
Referencias 1. E. V. Thomas y D. M. Haaland, Partial least-squares methods for spectral analyses. 1. Relation to other quantitative calibration methods and the extraction of qualitative information, Anal. Chem. 1988, 60, 1193-1202 2. D. L. Massart, B. G. M. Vandeginste, L. M. C. Buydens, S. De Jong, P. J. Lewi y J. Smeyers-Verbeke, Handbook of Chemometrics and Qualimetrics, Elsevier, Amsterdam, 1997, Capítulo 10. 3. R. G. Brereton, Chemometrics. Data Analysis for the Laboratory and Chemical Plant, Wiley, Chichester, 2003, Capítulo 5. http://www.chm.bris.ac.uk/org/chemometrics/pubs/chemweb.html 4. 5. D. González Gómez, A. Muñoz de la Peña, A. Espinosa Mansilla, A. C. Olivieri, Spectrophotometric Analysis of Mixtures by Classical Least-Squares Calibration: An Advanced Experiment Introducing MATLAB, Chem. Educator 2003, 8 , 187-191.
73
CAPÍTULO 6 CALIBRACIÓN MULTIVARIADA PARTE II: INVERSA (ILS Y PCR)
CALIBRACIÓN
Regresión por cuadrados mínimos inversos En este capítulo sobre calibración inversa exploraremos dos métodos para el análisis de mezclas de multianalitos: la regresión por cuadrados mínimos inversos (ILS, del inglés inverse least-squares ) y la regresión por componentes principales (PCR, del inglés principal component regression ). Las teorías de ambos métodos se exponen en este documento, pero se recomienda consultar paralelamente el ejemplo concreto que se analiza al final. Los métodos de calibración inversa reciben este nombre porque se basan en el uso de la ley de linealidad respuesta-concentración escrita en forma inversa a los métodos clásicos tales como CLS. Como se verá a continuación, los métodos inversos permiten estudiar mezclas de componentes en las que uno o más analitos son de interés, pero de los restantes componentes pueden desconocerse concentraciones, espectros e identidades químicas. De este modo, permiten superar una de las grandes desventajas de CLS: la necesidad del conocimiento de las concentraciones de todos los componentes presentes en las mezclas de calibrado. Tal como se discutió para CLS, la calibración directa implica la medida de espectros de muestras de calibración, conteniendo analitos con concentraciones conocidas, y obtención de la matriz de sensibilidades a partir de la ley "directa" por ajuste mediante cuadrados mínimos: Señal = Concentración × Sensibilidad
(6.1)
En cambio, en la calibración inversa se utiliza la ley de linealidad escrita en forma "inversa": Concentración = Señal × Coeficiente de regresión
(6.2)
donde se supone la existencia de una proporcionalidad entre la concentración de componentes calibrados y la correspondiente respuesta, a través de coeficientes de regresión que deberán calcularse mediante un modelo apropiado. Si bien el modelo CLS puede en principio interpretarse mediante una ecuación similar a la (6.2), en los métodos inversos la ecuación (6.2) se aplica cuando sólo se conoce la concentración de algunos analitos en las muestras de calibrado, pero se desconocen los restantes componentes. Este importantísimo concepto será detallado en el presente capítulo, y constituye la base sobre la que se afirman los métodos quimiométricos más provechosos para calibración multivariada. La bibliografía sobre el tema, particularmente en lo que concierne a PCR, es muy abundante. Recomendamos especialmente el texto clásico de Massart y colaboradores, 1 y el artículo de Haaland y Thomas. 2 Calibración Debemos notar que, en el campo de la calibración inversa, la literatura utiliza una notación para señales y concentraciones que es la inversa a la empleada en la discusión 74
del modelo CLS. Dado que la concentración se considera ahora la variable dependiente y la señal la variable independiente, X identificará la señal e Y la concentración. El método ILS es el más simple de los métodos inversos, y está basado en la ley de Beer inversa: Y = X BT
+E
(6.3)
donde la matriz (de tamaño I × J ) X reúne las señales instrumentales para I mezclas de calibrado, recogidas a J longitudes de onda. La matriz Y, por su parte, contiene las concentraciones de calibración en cada una de las I mezclas, de cada uno de los N analitos calibrados, y su tamaño es de I × N . En la ecuación (6.3), B es una matriz de J × N que relaciona las concentraciones con las respuestas de manera inversa a la ley de Beer, también llamada matriz de los coeficientes de regresión. Finalmente, E es una matriz de errores no modelados por la ecuación (6.3), siendo su tamaño idéntico al de Y. Para obtener la matriz B, se debe despejar ésta de la ecuación (6.3), empleando el criterio de cuadrados mínimos en el que E se considera nula. Para despejar B, se deben pre-multiplicar ambos miembros de (6.3) por XT: XT Y = (XT X) BT
(6.4)
Aquí se presenta un importante inconveniente del método ILS. Para continuar el proceso a partir de la ecuación (6.4), es preciso invertir la matriz ( XT X). Esto implica que si se han realizado mediciones a un número de sensores J que es mayor que el de mezclas I , (XT X) no puede invertirse, ya que el determinante de ( XT X) será en este caso nulo. Un ejemplo numérico aclarará el problema: supongamos que X es una matriz de 2×4 como la que se muestra en la Figura 6.1.
X=
Figura 6.1.
⎡0 1 2 1⎤ ⎢ 4 3 1 0⎥ ⎣ ⎦
⎡16 ⎢12 T X X = ⎢ ⎢4 ⎢ ⎣0
12 10 5 1
4 5 5 2
0⎤ 1⎥⎥ 2⎥ ⎥ 1⎦
Producto de una matriz por su traspuesta, generando una matriz
singular. El determinante de ( XT X) es nulo, por lo que ( XT X) es singular y no puede invertirse. La singularidad es inevitable, ya que el modo en que se produce ( XT X) hace que sus líneas sean combinaciones lineales. En el ejemplo de la Figura 6.1, una de las combinaciones lineales presentes hace que la tercera fila sea igual a (primera fila / 4 + cuarta fila × 2). En términos del modelo descrito por la ecuacion (6.3), implica que debe resolverse un sistema de ecuaciones sub-determinado, en el que el número de incógnitas es J × N (los J × N elementos de B) disponiendo solamente de I × N ecuaciones. El único modo de evitar esta singularidad es emplear menos sensores que mezclas, lo que puede considerarse como una seria limitación del método: la necesidad de contar con más mezclas de calibración que sensores. Sin embargo, este modelo ILS posee una gran ventaja en relación con CLS: logra desacoplar los componentes químicos entre sí, importantísimo concepto que ilustraremos a continuación, y que implica que sólo es necesaria la información de la 75
concentración del (o los) componente(s) de interés para calibrar el modelo. En otras palabras, se podrá cuantificar un analito en presencia de una interferencia, siempre que ésta haya sido incluida en la calibración, aunque no se conozca su concentración. En caso de que ( XT X) pueda invertirse (bajo la condición de que J < I ) es posible despejar B de la ecuación (6.4): T
–1
T
B = (X X) X Y
(6.5)
Esta última ecuación puede interpretarse diciendo que cada columna de B se obtiene por el producto de [( XT X) –1 XT] por una columna específica de Y (que contiene los datos de concentración de un componente dado en las mezclas de calibración). La Figura 6.2 muestra cómo se obtiene este vector de regresión, contando sólo con la matriz de los datos instrumentales y el vector que contiene la concentración del analito de interés (y n).
B T
n
=
–1
Y
T
(X X) X
×
yn
J × I J × 1
I × 1
Figura 6.2.
Esquema que muestra cómo es posible calcular el vector de regresión conociendo sólo la información de la concentración del analito de interés. El vector n es el producto de la matriz gris oscura por el vector yn. Por lo tanto, es posible plantear un modelo simplificado en el que no es necesario conocer la concentración de los restantes componentes de calibrado sino sólo la del analito n: n =
(XT X) –1 XT yn
(6.6)
La ecuación (6.6) ilustra lo que se conoce como "desacople" de componentes, situación que no puede lograrse en CLS, donde es preciso conocer las concentraciones de todas las especies presentes en las muestras empleadas para calibrar. En la ecuación (6.6), n representa el vector de coeficientes de regresión asociado al componente particular n, mientras que yn es un vector que contiene las concentraciones del analito n en las mezclas de calibrado. Nótese que la necesidad de invertir la matriz ( XT X) para obtener los coeficientes de regresión implica que ILS será sensible a colinealidades espectrales, tal como fuera discutido para CLS. Es preciso destacar también que la ecuación (6.6) no implica que ILS permita analizar un único analito. Si existen varios analitos de interés, además de un número no identificado de componentes adicionales, se puede plantear un modelo desacoplado como el de la ecuación (6.6) para cada analito de interés. La etapa de calibrado es análoga a la descrita en el caso del modelo CLS, excepto que en ILS los elementos del vector de regresión asociado a un componente particular se obtienen a partir de las mezclas de calibrado, ignorando las concentraciones de los restantes componentes. Esto no era posible en CLS. 76
No obstante, el precio que se paga por esta ventaja es alto: deben prepararse más mezclas de calibrado que sensores de lectura de la señal (lo cual puede ser difícil en términos de costo o tiempo experimental), o bien deben estudiarse unos pocos sensores, desperdiciando información útil que es típica de las mediciones multisensoriales. Predicción Durante la etapa de predicción se tendrá una ecuación similar a la de calibrado, la ley inversa de Beer aplicada a la muestra incógnita: yn = (
n
)T x
(6.7)
La ecuación (6.7) permite observar que n se comporta como el vector de coeficientes de regresión para el componente n, tal como fuera discutido en el caso de CLS. Si existiera más de un analito de interés, la ecuación (6.7) se aplicaría tantas veces como fuese necesario, utilizando cada vez el vector n asociado al analito n. Ventajas y desvent ajas de ILS La posibilidad de desacoplar componentes origina la principal ventaja de este método, pudiéndose estudiar mezclas complejas mediante un proceso de calibración en el que se conoce sólo la concentración del componente de interés. Su desventaja radica en que ILS sigue siendo sensible a las colinealidades espectrales discutidas en los capítulos anteriores, y que se debe usar un número reducido de sensores, con la consecuente pérdida de información y por ende de sensibilidad. Como muestra de la potencialidad de ILS, téngase en cuenta que la técnica fue originalmente desarrollada para el análisis de propiedades de polímeros o determinaciones del contenido de proteínas en semillas a partir de espectros de absorción de infrarrojo cercano (NIR). Los espectros NIR presentan bandas debidas a un enorme conjunto de especies presentes en estos materiales. En la calibración de una propiedad específica o del contenido de proteínas, sin embargo, la información acerca de las concentraciones de los componentes de estos sistemas complejos es extremadamente limitada. Aún así, ILS es capaz de proveer una respuesta inteligente a este tipo de problemas analíticos. Debido a sus desventajas, sin embargo, la práctica moderna lo ha desplazado por métodos más poderosos.
Regresión po r c omponentes princ ipales La pregunta que surge automáticamente al considerar los modelos CLS e ILS es: ¿porqué no pueden aprovecharse las ventajas de ambos a la vez?. El método de regresión en componentes principales o PCR representa uno de los intentos más simples de reunir sus principales ventajas. Emplea una calibración inversa, pero no correlaciona las concentraciones directamente con las respuestas instrumentales, sino con una matriz más pequeña, llamada de puntuaciones (en inglés scores ). Estos scores o variables latentes deben condensar de un modo eficiente la información espectral completa (las variables manifiestas) en una matriz de tamaño adecuado. Esto puede realizarse matemáticamente con ayuda de los autovectores de la matriz cuadrada ( XT X) (de tamaño J × J ). La etapa de condensación o compresión de la información contenida en X es esencial para comprender el funcionamiento del modelo PCR. Quienes deseen revisar algunos conceptos sobre matrices y sus operaciones, que son necesarios para seguir en detalle el presente capítulo, pueden consultar el Anexo I. 77
Compresión de la información La compresión de la información contenida en la matriz de señales de calibración es el paso crítico para el modelo PCR. Una técnica muy empleada en quimiometría para la eficiente compresión de datos es su "proyección" sobre los autovectores de la matriz (XTX). Existen varios algoritmos capaces de obtener dichos autovectores, entre los cuales uno muy eficiente se basa en un tipo de descomposición matricial conocido como T descomposición en valores singulares , que consiste en descomponer a la matriz X (tamaño J × I ) en el producto de otras tres matrices: XT = U W VT
(6.8)
Podemos apreciar los requerimientos de tamaño matricial en esta última ecuación a través del esquema presentado en la Figura 6.3.
T
X
U W
J × I
=
Estos números deben coincidir ( J )
J × I
×
I × I
Estos números deben coincidir ( I )
V
×
T
I × I
Estos números deben coincidir ( I )
Estos números deben coincidir ( I ) Figura 6.3.
Esquema que muestra las relaciones de tamaños de matrices en la descomposición en valores singulares. Las matrices U ( J × I ), W ( I × I ) y V ( I × I ) cumplen con las siguientes condiciones: • Las columnas de U son ortogonales, de modo que UT×U = I, así como las de V, de modo que VT×V = I (I representa una matriz identidad de tamaño apropiado). • La matriz W es diagonal y sus elementos diagonales son no negativos (los no diagonales son iguales a cero). Los elementos diagonales de W se llaman valores singulares de la matriz XT. Matemáticamente, son las raíces cuadradas de los autovalores no negativos de ( XTX); desde el punto de vista químico miden la contribución a la variación espectral que puede ser explicada por cada uno de los componentes principales de X. 78
Las columnas de U son los autovectores de ( XTX), mientras que las columnas de V son los autovectores de ( XXT). En la literatura inglesa las columnas de U se llaman corrientemente loadings; en castellano se suelen llamar factores, o también variables latentes, en oposición a las variables manifiestas, que son las experimentalmente accesibles (las latentes deben ser halladas mediante operaciones matemáticas). Reuniendo el producto W VT en una única matriz T, la ecuación de descomposición singular (6.8) se puede también escribir como: T
X
=UT
(6.9)
que es la base para la regresión en componentes principales, donde T es la matriz de scores antes mencionada. Para obtener T a partir de los datos instrumentales, esto es, despejar T de la ecuación (6.9), se requiere pre-multiplicar por UT, y luego por la inversa de ( UT U), pero esta última matriz es igual a la matriz identidad I, por lo que se obtiene, directamente: T
T = U X
T
(6.10)
Los tamaños de las tres matrices involucradas en la ecuación (6.10) son, respectivamente, I × I , I × J y J × I . Dicha ecuación puede interpretarse diciendo que los scores constituyen la proyección de la matriz original de datos en el espacio de los factores. Esta proyección es la etapa fundamental de compresión de datos, ya que logra reducir la dimensionalidad de la matriz original (de I × J ) a una matriz de scores más pequeña (de I × I ). El análisis criterioso de los scores , no obstante, permitirá discernir que estos están ordenados de un modo coherente, en orden decreciente de su contribución a la variación espectral en X. Por lo tanto, la selección de los scores significativos (estadísticamente hablando) permitirá reducir aún más el tamaño de T. Un comentario final acerca de las propiedades de la matriz de scores : ésta presenta la ventaja de estar construida con columnas que son ortogonales entre sí. La propiedad de ortogonalidad implica que el producto escalar de cualquier columna de T por cualquier otra columna es nulo: (ti )T × ti' = 0
(si i ≠ i' )
(6.11)
La consecuencia más importante de la ecuación (6.11) es que el modelo PCR está libre de los efectos de las colinealidades espectrales. Esto es así porque en PCR se correlacionan las concentraciones con los scores , que pueden considerarse como un tipo especial de "espectros". Estos espectros no muestran ningún paralelismo entre sí, debido a la propiedad estipulada en la ecuación (6.11). En la literatura existe cierta confusión respecto de a qué se llama componente principal, o simplemente componente: a veces se emplea el término refiriéndose a las columnas de U, otras veces a las columnas de T. Para evitar la confusión adicional con los componentes químicos de cada sistema, llamaremos factores a las columnas de U y scores a las de T, dejando la expresión "componente principal" para identificar a la unidad factor/score. Finalmente, es preciso mencionar que la descomposición aquí presentada no es el único método para calcular componentes principales. El más célebre, quizás, es el NIPALS (por non-linearr iterative partial least-squares), desarrollado por H. Wold. 3
79
Componentes prin cip ales y fuentes de variación espectral Se acostumbra a emplear el método de descomposición singular para identificar fuentes de variación espectral. Por fuente de variación se entiende todo fenómeno capaz de producir una variación en los espectros de una muestra a otra. Obviamente los componentes activos son fuentes de variación, pero también lo son el ruido instrumental, las derivas de la línea de base, las pérdidas de la linealidad, etc. Cuando estos últimos fenómenos son de menor importancia que la presencia de los componentes químicos espectralmente activos, se supone que hay tantas fuentes de variación como componentes químicos. Sin embargo, esto en general no se cumple, y además el número de componentes químicos de una mezcla compleja puede ser desconocido, de modo que la información acerca de las fuentes de variación que proporciona el análisis de los datos espectrales es sumamente valiosa. Hay varias maneras de estimar las fuentes de variación; una muy popular es la que resulta de considerar la contribución relativa de cada componente principal a la varianza espectral total, calculada del modo que sigue: A
∑=1 (W ) 2 i
% Varianza explicada por los primeros A factores = 100
i I
(6.12)
∑=1 (W ) 2 i
i
donde W i es cada uno de los elementos diagonales de la matriz W, y (W i)2 es el correspondiente autovalor asociado al componente principal i. Dado que cada componente principal contribuye con una porción cada vez menor de la varianza total, lo usual es tomar el número de los primeros componentes principales que, colectivamente, aportan un determinado porcentaje de la varianza total. La Figura 6.4 ilustra el comportamiento típico de un conjunto de componentes principales: mientras los valores singulares disminuyen, su contribución a la varianza total también disminuye. En el caso de la Figura 6.4, por ejemplo, los tres primeros componentes principales explican más del 99% de la varianza espectral, lo que llevaría a la conclusión de que hay tres fuentes de variación espectral en los espectros contenidos en la matriz X.
80
Figura 6.4.
Varianza explicada en función del número de componente
principal. Este método de estimación de fuentes de variación adolece de dos problemas. En primer lugar, distintos autores emplean diferentes criterios para el porcentaje óptimo de varianza explicada, y parece difícil establecer un criterio común. Por otro lado, se usan únicamente los datos de las señales instrumentales de calibración para estimar el número de factores necesarios para la reducción de la información. En el ámbito analítico, es preferible incorporar en este análisis la información disponible acerca de la concentraciones de calibrado del componente de interés. Para ello se ha diseñado el método más popular de estimación del número óptimo de factores, llamado validación cruzada. Lo discutiremos más adelante, después de explicar cómo calibrar y predecir con el modelo PCR. De todas maneras, aún cuando existen varias herramientas para estimar el número apropiado de fuentes de variación espectral en la matriz de datos X, la inspección visual del aspecto de los factores puede ser importante, en homenaje a la frase "el ojo del amo engorda el ganado". La Figura 6.5 ilustra la diferencia entre un autovector capaz de representar fenómenos físicos que llevan a la variación espectral de la matriz X, y otro que representa, básicamente, ruido instrumental. El primero tiene forma "espectral"; el segundo, de ruido al azar.
81
Figura 6.5.
Izquierda, un típico autovector que representa variaciones de señal instrumental debida a fenómenos químicos. Derecha, un autovector que representa ruido instrumental. El análisis del número de fuentes de variación es sumamente importante. Supongamos que este estudio ha indicado que el número de factores A que explican un porcentaje muy significativo de la varianza es menor que el número de mezclas de calibrado I . Dado que los primeros A factores son suficientes para explicar prácticamente todo el comportamiento espectral de la matriz X, no es necesario emplear la matriz U completa en la proyección de la ecuación (6.11), sino que pueden quitarse las columnas desde la A+1 hasta la I , quedando una matriz conformada sólo por los primeros A autovectores, que llamaremos U A (tamaño J × A) Los restantes autovectores pueden descartarse puesto que se considera que modelan únicamente el ruido espectral. De este modo, la matriz de scores puede reducir aún más su tamaño, de I × I a I × A: T
T A = U A X
T
(6.13)
En la ecuación precedente, hemos llamado T A a la matriz de scores estimada con A factores. Esta nueva matriz T A, a pesar de tener un tamaño considerablemente menor que la matriz original de espectros, cumple no obstante un papel similar, ya que la información relevante presente en X ha sido comprimida de un modo eficiente. El proceso de compresión puede ilustrarse con la serie de imágenes de la Figura 6.6, que muestran la fotografía de una flor, considerada como una matriz de puntos, que puede comprimirse utilizando componentes principales, y luego "reconstruirse" recurriendo a la ecuación (6.9) escrita en términos de los A componentes selectos, esto es XT = U A T A. La imagen que corresponde a A = 1 está reconstruida utilizando sólo el primer componente principal, que es el que más aporta a la varianza matricial. A medida que se emplean más y más componentes principales, la imagen se hace más nítida. No obstante, se aprecia que empleando unos pocos factores, la información relevante es retenida por la matriz de datos comprimida. 4
82
Imagen total
A = 1
A = 2
A = 4
A = 8
A = 16
A = 32
Figura 6.6.
Una imagen (arriba al centro), reconstruida utilizando distinto número de componentes principales (abajo) 4
Calibración En este punto reuniremos las ventajas de ILS y CLS, que era el objetivo primordial al comenzar este capítulo de calibración multivariada. Plantearemos un modelo de calibración inversa, en el que las concentraciones del analito calibrado en las muestras de calibración ( yn, tamaño I ×1) se correlacionan linealmente con los scores contenidos en T A: yn = T A vn + e
(6.14)
donde vn (tamaño A×1) es el vector de coeficientes de regresión correspondiente, y e un vector que recoge los errores de modelado. Se puede obtener el vector vn despejando de la ecuación anterior. Pre-multipicando ambos miembros por T AT se obtiene: T AT yn = T AT T A vn + e
(6.15)
Luego será necesario multiplicar por la inversa de ( T AT T A), o sea, por ( T AT T A) –1. Esta última operación de inversión es trivial, ya que ( T AT T A) es una matriz diagonal (en atención a la ortogonalidad de las l as columnas de T), y la inversión de una matriz diagonal se remite a la inversión de cada uno de sus elementos diagonales. Finalmente Finalmente se obtiene, entonces: vn
= (T AT T A) –1 T A yn
(6.16)
La inversión de ( T AT T A) en la ecuación (6.16) no presenta problemas asociados a la colinealidad, por los motivos anteriormente expuestos: las columnas de T A son ortogonales. Por analogía con el criterio empleado en el modelo CLS, podemos llamar a la matriz [(T AT T A) –1 T A] la seudoinversa de T A y denominarla T A+, con lo cual la ecuación (6.16) adopta su forma final: vn
= T A+ yn
(6.17)
Este último paso completa la calibración. La obtención de los coeficientes de regresión vn es completamente análoga al proceso realizado en CLS, y su empleo en la predicción de la concentración concentración del analito analito en muestras incógnitas es también también similar. 83
Predicción En la etapa de predicción, se registra el espectro de una muestra incógnita, y se almacenan las señales instrumentales en el vector columna x (tamaño J ×1). Antes de aplicar el modelo de predicción es necesario proyectar dicho vector en el espacio de los A factores de la matriz U A, dado que no podemos emplear los datos originales para estimar concentraciones "mezclando" el vector espectral real con los coeficientes de regresión "comprimidos" contenidos en vn. A×1) Análogamente a la ecuación (6.13), entonces, se obtiene el vector t A ( A correspondiente a la muestra incógnita: T
t A = U A x
(6.18)
Este vector t A contiene los scores de la muestra, que actuarán en calidad de "espectros" en la etapa predictiva del modelo. Esta última no es sino la repetición del modelo inverso de la ley lineal expresado anteriormente en ILS, esto es: T
yn = (vn) t A
(6.19)
en el que vn reemplaza a n y t A reemplaza a x. A partir de esta última ecuación se estima la concentración del analito en la incógnita. Validació Validació n cru zada zada La posibilidad de calibrar y predecir mediante un modelo inverso del tipo PCR ofrece la alternativa de seleccionar el número apropiado de factores ( A) mediante una combinación de información espectral y de concentraciones, que se conoce como validación cruzada. Consiste en calibrar el modelo con todas las muestras de calibración excepto una, predecir la concentración de la muestra dejada de lado, y calcular el error cometido (diferencia entre el valor nominal y el predicho). Este cálculo se realiza utilizando un número creciente de factores, desde uno hasta un cierto máximo. El máximo puede establecerse a voluntad (debe ser menor al número de mezclas de calibrado). Luego se repite el procedimiento hasta que todas las muestras hayan sido dejadas de lado una vez. En cada caso, se predicen las concentraciones del analito en cada una de las muestras dejadas de lado. Para cada número de factores, se calcula la suma de los cuadrados de los errores de predicción, que se acostumbra a llamar PRESS (por predicted error sum of squares). Luego se procede a estudiar cómo varía el PRESS así obtenido en función del número de factores mediante un procedimiento estadístico. A modo de ejemplo, supóngase que durante un procedimiento típico de validación cruzada se ha obtenido la siguiente tabla de valores de PRESS en función del número de factores (véase la Figura 6.7): Tabla 6.1. Variación del PRESS con el número de factores incluídos en la calibración. Factores PRESS 1 0,92 2 0,0217 3 4,3×10 –3 4 4,1×10 –3 5 3,7×10 –3 6 5,1×10 –3 84
Figura 6.7. Variación del PRESS en función del número de factores para un
modelo PCR típico. Se observa que, a medida que se agregan factores, el PRESS disminuye: esto se debe a que la compresión de los datos se va haciendo progresivamente más eficientes, puesto que los primeros factores contienen información relevante respecto de la variación espectral en la calibración. Si se emplean menos factores que los necesarios se obtiene una situación poco deseable llamada subajuste de los datos. Al seguir aumentando el número de factores, el PRESS parece estabilizarse y finalmente aumenta ligeramente. Esto es una fuerte indicación de que los últimos factores no están aportando información relevante sino esencialmente ruido. Emplear más factores de lo debido puede llevar a una situación también indeseable llamada sobreajuste. Intuitivamente, podría plantearse que el número óptimo de factores es aquel que lleve al mínimo PRESS. Sin embargo, estudios estadísticos cuidadosos indican que este no es el caso. Una técnica conveniente para estimar A es la descripta por Haaland. 2 Consiste en ampliar la tabla anterior, calculando los cocientes entre los distintos PRESS y el mínimo (sólo para un número de factores inferior a aquel que produce el mínimo PRESS). Estos cocientes de PRESS cumplen el papel de un cociente de varianzas, de manera que tienen asociada una probabilidad, que se estima estadísticamente con un número de grados de libertad igual al número de mezclas de calibrado I tanto tanto para el numerador como para el denominador.
85
La tabla completa sería como sigue: Tabla 6.2. Variación del PRESS con el número de factores incluídos en la calibración y cálculo de la probabilidad de la distribución F . Factores PRESS PRESS/min(PRESS) p 1 0,92 248 0,999 2 0,0217 5,88 0,997 –3 3 1,17 0,605 4,3×10 –3 4 1,12 0,576 4,1×10 –3 5 1 0,5 3,7×10 6 – – 5,1×10 –3 Haaland propone, basándose en resultados empíricos, seleccionar como número óptimo de factores el primer valor para el que la probabilidad asociada disminuye por debajo de 0,75. En la tabla anterior, este criterio llevaría a elegir A = 3. El valor de óptimo de PRESS puede emplearse para tener una idea de la bondad del ajuste de concentraciones, ya que permite acceder al llamado error medio de validación cruzada RMSECV (por root mean square error in cross-validation ), obtenido como RMSECV = [PRESS/( I )]1/2. Este parámetro debe ser del orden de la incertidumbre asociada a las concentraciones de calibrado. Residu os espectr ales Como todo método que emplea espectros completos, PCR es capaz de proveer residuos espectrales para la muestra incógnita, como la diferencia entre el espectro experimental de la muestra y el espectro estimado por el modelo. El espectro calculado por el modelo, xˆ , también llamado espectro "reconstruido", se obtiene simplemente a partir de una ecuación análoga a la ecuación (6.9), pero empleando el vector de scores de la muestra y la matriz de factores reducida U A: xˆ
= U A t A
(6.20)
Luego puede definirse el residuo espectral en forma análoga a CLS: J
∑=1 ( x
j
sres =
− xˆ j ) 2
j
J − A
(6.21)
Cifras de mérit o Las cifras de mérito se pueden calcular con ecuaciones similares a las empleadas en el modelo CLS. Para ello, se requiere el análogo de los coeficientes de regresión espectrales n, que puede obtenerse mediante una ecuación análoga a la (6.20), esto es, "reconstruyendo" el vector espectral n a partir del vector reducido vn: (6.22) n = U A vn Lo que no existe en el ámbito de PCR es la estimación del espectro del analito puro, que en CLS eran las columnas de la matriz S. Esto impide el cálculo de la selectividad en PCR mediante la aproximación discutida en el Capítulo 4. No obstante, la selectividad puede calcularse en PCR recurriendo a conceptos que están más allá del 86
alcance de este libro. Los programas suministrados con este capítulo permiten estimar todas las cifras de mérito. Ventajas y desventajas de PCR Como resumen de este capítulo, podemos enumerar las ventajas de PCR respecto de las otras técnicas multivariadas que hemos estudiado hasta el presente. PCR combina las ventajas ya analizadas de CLS con dos adicionales: 1) calibración directa, que permite ignorar las concentraciones de compuestos químicos desconocidos durante el calibrado, y 2) uso de "espectros" abstractos llamados scores , que eliminan los problemas asociados con la colinealidad espectral. En referencia a la tabla de propiedades analíticas presentada en el Capítulo 4, se mantiene, sin embargo el problema de las interferencias no modeladas. Este problema es común a la mayoría de los métodos multivariados basados en información espectral: si aparece en una muestra incógnita un compuesto no contenido en la calibración, el análisis no será exacto. Aún así, la falta de exactitud tiene su estilo en el mundo multivariado, ya que los modelos son capaces de detectarla, aunque no de corregirla. Más allá d e PCR Si PCR reúne las ventajas de CLS e ILS, y ninguna de sus desventajas, y si además su punto débil es común a todas las técnicas basadas en espectros, la pregunta lógica es: ¿qué puede ser mejor que PCR?. La respuesta es que el espacio para la mejora de los métodos multivariados es inmenso. Un defecto que puede adjudicarse a PCR es que utiliza factores calculados en base a información espectral del calibrado únicamente, sin referencia a las concentraciones de calibrado. Esta última información es valiosa, y métodos multivariados basados en la combinación de espectros y concentración para el cálculo de factores son capaces de superar a PCR en valor predictivo. El más popular es la regresión en cuadrados mínimos parciales o PLS (por partial least-squares). Material suministrado • Para este capítulo se provee la siguiente información complementaria: Archivos de texto (*.TXT) conteniendo datos típicos. Archivos (*.M) con rutinas para el entorno de programación MATLAB. Archivos (*.EXE) con programas ejecutables en QB.
Ejercicio 1) Los datos del presente ejercicio están tomados del trabajo citado como referencia número 5. Se desea determinar el contenido de un fármaco, la bromhexina, presente en muestras de jarabe para la tos. Los componentes del jarabe se conocen en forma incompleta, de manera que se preparan muestras para construir un modelo PCR. Para ello, se agregan cantidades conocidas de bromhexina a doce diferentes muestras de jarabe “blanco” (esto es, el fondo de la matriz del jarabe, sin bromhexina), y se utilizan para calibrar el modelo. Las concentraciones del analito en las muestras de calibrado son:
87
Tabla 6.3. Concentración de las muestras de calibración correspondientes al Ejercicio 1. Muestra de calibrado Concentración ×104 M 1 1.55 2 2.06 3 2.58 4 1.55 5 2.06 6 2.58 7 1.55 8 2.06 9 2.58 10 1.68 11 2.10 12 2.66 Estas concentraciones se recogen en forma de un vector de 12 ×1 en el archivo de texto BR_CON_C.TXT Los espectros de absorción de estas 12 muestras se registran a 64 diferentes longitudes de onda. Estos espectros están contenidos, en forma de matriz de 64 ×12, en el archivo de texto BR_RES_C.TXT. Informar las correspondientes cifras de mérito para el modelo. Suponga que el nivel de ruido instrumental es igual a 0,003 unidades de señal. 2) Para la validación del modelo, se prepararon 11 muestras adicionales de jarabe con contenido conocido de bromhexina, diferente al empleado para calibrar. Los espectros de estas muestras están contenidos, en forma de matriz de 64 ×11, en el archivo BR_RES_T.TXT, y las concentraciones nominales, en forma de vector de 11 ×1, en el archivo BR_CON_T.TXT. Estimar las concentraciones de los analitos en este juego de muestras y sus incertidumbres asociadas, y estudiar la exactitud del método mediante la prueba de la elipse. 3) Una muestra adicional de prueba, cuyo espectro está contenido en el archivo de texto BR_RES_P.TXT se analiza mediante el mismo modelo. Sin embargo, se sospecha que se trata de una muestra que contiene una interferencia no modelada en la calibración. ¿Qué conclusiones puede extraer al respecto del análisis mediante PCR?
Respuesta 1) El primer paso en el análisis PCR debe ser el estudio del número óptimo de factores presentes en la matriz de calibrado, que luego se emplearán para la predicción. El método más recomendado para esto es la validación cruzada, que puede implementarse mediante la rutina PCR_CV.M de Matlab o el programa PCR_CV.EXE de QB. Para ejecutar estos algoritmos, se requiere introducir un número máximo de factores de prueba. Este puede ser, como máximo, igual al número de mezclas de calibrado menos una (ya que el procedimiento consiste en calibrar con las muestras menos una), en el presente caso 11 = 12 – 1. No obstante, se supone que se han preparado más mezclas de calibración que fuentes de variación espectral, por lo que se recomienda introducir, como número máximo, un valor menor. Los resultados obtenidos para un número máximo de factores igual a ocho son los siguientes: 88
Tabla 6.4. Selección del número de factores correspondiente al Ejercicio 1. Factores PRESS PRESS/min(PRESS) p 1 0,92 248 0,999 2 5,88 0,0217 0,997 –3 3 1,17 0,605 4,3×10 –3 4 1,12 0,576 4,1×10 –3 5 1 0,5 3,7×10 6 – – 5,1×10 –3 7 – – 8,9×10 –3 8 – – 1,1×10 –2 Puede apreciarse que el PRESS disminuye al ir aumentando el número de factores, llega a un mínimo para 5 factores, y luego aumenta. El número óptimo de factores, obtenido para el primer valor de p que disminuye por debajo de 0,75 es 3. El RMSECV para 3 factores es satisfactorio (0,02) en vista de las concentraciones nominales de calibrado y sus incertidumbres asociadas (en la segunda cifra decimal). Estos primeros tres componentes principales explican más del 99,99% de la varianza de la matriz espectral. Tanto los resultados correspondientes al PRESS como la varianza explicada se observan gráficamente en la figura generada por MATLAB, figura que también puede construirse mediante los valores provistos por el programa QB correspondiente (PCR_CV.EXE).
Figura 6.8. Variación del PRESS en función del número de factores para el
Ejercicio 1. 89
Figura 6.9. Variancia explicada por los factores para el Ejercicio 1.
Una vez establecido el número óptimo de factores para la compresión de la información, se procede a calibrar el modelo, empleando los programas PCR_CAL.M (Matlab) o PCR_CAL.EXE (QB). Las cifras de mérito calculadas mediante los programas para este modelo son las siguientes: Tabla 6.5. Cifras de mérito calculadas para el Ejercicio 1. Cifra de mérito Valor Sensibilidad 1,21×104 A × M –1 Sensibilidad analíticaa 4×106 M –1 1/γ 2,5×10 –7 M Selectividad 0,46 a Obtenida dividiendo la sensibilidad por el nivel de ruido instrumental (0,003 unidades).
2) Para predecir las concentraciones de las muestras incógnitas, empleamos los programas PCR_PRED.M (Matlab) o PCR_PRED.EXE (QB), con los siguientes resultados:
90
Tabla 6.6. Concentraciones predichas correspondientes al Ejercicio 2. Muestra Residuo Concentración × 104 espectral a Nominal Predicha 1 1,96 1,97(1) 0,004 2 2,16 2,19(1) 0,002 3 0,00 –0,01(1) 0,014 4 0,82 0,84(1) 0,009 5 1,02 1,04(1) 0,006 6 1,33 1,37(1) 0,005 7 1,84 1,93(1) 0,003 8 2,35 2,43(1) 0,004 9 1,94 2,00(1) 0,004 10 2,14 2,19(1) 0,003 11 2,24 2,25(1) 0,006 a Los errores estándar en las concentraciones, calculados con el modelo aproximado citado en la teoría, esto es s( xn) = s y / SENn, con s y = 0,003, son todos iguales a 0,002. Este valor es demasiado optimista, en vista de que las concentraciones de calibrado están dadas con una incertidumbre de 0,01, por lo que se ha optado por este último valor, más conservador, en la presente tabla. Se informan también, en la última columna de esta tabla, los residuos espectrales para cada muestra incógnita, que, como puede apreciarse, se encuentran dentro del nivel del ruido instrumental. Esto confirma que el ajuste por cuadrados mínimos para estas muestras es adecuado. Dos excepciones a esta situación son las muestras número 3 y 4, cuyo residuo espectral es superior al resto. Una explicación posible para esto es que estas muestras fueron preparadas con una concentración nominal inferior a las de calibrado. En este sentido, no se trata de verdaderos outliers, que contengan interferencias no modeladas, pero se trata de muestras para las que le estamos exigiendo al modelo que realice una extrapolación hacia un ambiente para el que no está entrenado. De todas maneras, nótese que las concentraciones predichas para estas muestras son muy cercanas al valor nominal. Para establecer la exactitud del método, lo recomendado es analizar los datos de la tabla precedente mediante la prueba de la elipse, tal como se discutiera en el Capítulo 3. De este modo, la tabla de datos a suministrar a los programas de cálculo de la elipse será como sigue: Tabla 6.7. Datos a usar para construir la elipse del Ejercicio 2. 1,96 1,97 2,16 2,19 0,00 –0,01 0,82 0,84 1,02 1,04 1,33 1,37 1,84 1,93 2,35 2,43 1,94 2,00 2,14 2,19 2,24 2,25 91
Dado que no se tienen resultados de réplicas de cada muestra, lo que proveería una estimación del desvío estándar de cada valor predicho, realizaremos un análisis mediante el método OLS. Se recomienda organizar los datos apropiadamente y someterlos a los programas EJCR.M (MATLAB) o EJCR.EXE (QB). El resultado se muestra en las figuras siguientes. La primera de ellas muestra los valores predichos en función de los nominales, y la segunda la elipse.
Figura 6.10. Valores predichos en función de los nominales (Ejercicio 2).
92
Figura 6.11. Elipse obtenida a partir de la regresión de la
Figura 6.10 (Ejercicio 2).
El análisis de esta figura revela que el punto ideal (1,0) no está contenido dentro de la elipse, por lo que la validación del modelo no pasa la prueba de exactitud. Nótese que los resultados del análisis EJCR indican que los valores de la pendiente y ordenada al origen, individualmente consideradas, pasan la prueba de exactitud, ya que sus valores e intervalos de confianza son: Pendiente = 1,02 ± 0,025 Ordenada = 0,002 ± 0,004 Se observa que ambos intervalos de confianza contienen a los respectivos valores ideales (1 y 0). Sin embargo, el modelo no aprueba el test más estricto del intervalo conjunto de confianza. ¿Qué puede hacerse en un caso como el presente? Una alternativa es estudiar un número mayor de muestras de validación, incluyendo réplicas, para realizar un análisis WLS que es más cercano al real. Repetir los análisis de muestras con mayor residuo espectral, o mayor desviación del valor nominal es otro recurso. Finalmente podemos utilizar otros modelos multivariados alternativos a PCR, que no están contemplados en este curso. Si luego de estos intentos, la gráfica EJCR es similar, quizás debamos conformarnos con la falta de exactitud del modelo, y aceptar que para un problema de la complejidad del presente esta es la mejor respuesta que se puede dar, valorando que el RMSE de predicción obtenido es satisfactorio. En nuestro caso, RMSE = 0,015, que puede considerarse satisfactorio en vista de que las concentraciones nominales de calibrado y validación llevan una incertidumbre de alrededor de 0,01 unidades. 3) El análisis de la muestra contenida en el archivo BR_RES_P.TXT arroja los siguientes resultados: Concentración estimada: 2,10 Residuo espectral: 0,08 93
Aquí el residuo es significativamente mayor que el ruido espectral, lo que haría sospechar la presencia de un interferente no modelado. Ejercicio complementario Se desea modelar, mediante PCR, la determinación del antibiótico tetraciclina en suero humano. La matriz de espectros de calibración es de 101 ×50 y consiste de 50 espectros de fluorescencia registrados a 101 longitudes de onda. Esta matriz está contenida en el archivo TE_RES_C.TXT. Las concentraciones del analito en los 50 sueros empleados para calibrar están, en forma de vector de 50 ×1, en el archivo TE_CON_C.TXT. Calibrar el modelo con el número óptimo de factores, y validarlo frente a las 57 muestras de validación contenidas en el archivo TE_RES_T.TXT (espectros, matriz de 101×57) y TE_CON_T.TXT (concentraciones, vector de 57 ×1). Analizar la exactitud mediante el método EJCR. Considerar que el nivel de ruido instrumental es igual a 3 unidades de fluorescencia.
Respuesta al ejercicio complementario 1) Los resultados provistos por el programa son los siguientes: El número óptimo de factores es 4. Calibrando el modelo con cuatro factores, y prediciendo las 57 muestras incógnita se obtienen las siguientes concentraciones predichas (se informan junto con las nominales):
94
Tabla 6.8. Predicciones correspondientes al Ejercicio complementario. Muestra Concentración Concentración Muestra Concentración Concentración nominal predicha nominal predicha 1 1.25 1.08 30 3.50 3.65 2 1.25 1.29 31 3.75 3.79 3 1.50 1.43 32 3.75 3.52 4 1.50 1.45 33 4.00 3.83 5 1.50 1.35 34 1.00 1.01 6 1.50 1.37 35 1.00 1.02 7 1.50 1.38 36 1.00 1.00 8 1.50 1.51 37 1.00 1.00 9 1.50 1.39 38 0.00 0.99 10 1.75 1.77 39 0.00 0.00 11 1.75 1.73 40 0.00 –0.01 12 1.75 1.83 41 0.00 –0.03 13 2.00 2.03 42 0.60 0.62 14 2.00 2.08 43 0.60 0.67 15 2.00 1.87 44 0.60 0.67 16 2.25 2.14 45 0.60 0.67 17 2.25 2.29 46 2.00 0.21 18 2.50 2.63 47 2.00 0.19 19 2.50 2.49 48 2.00 0.17 20 2.50 2.37 49 2.00 0.22 21 2.75 2.75 50 0.40 0.38 22 2.75 2.76 51 0.40 0.39 23 2.75 2.75 52 0.40 0.41 24 3.00 2.96 53 0.40 0.36 25 3.00 3.00 54 0.80 0.81 26 3.00 2.86 55 0.80 0.82 27 3.50 3.60 56 0.80 0.78 28 3.50 3.58 57 0.80 0.79 29 3.50 3.36 Con estos datos se puede utilizar el programa EJCR.M para evaluar la exactitud, usando el método OLS:
Figura 6.12. Elipse correspondiente al Ejercicio complementario obtenida con Matlab.
95
Referencias 1. D. L. Massart, B. G. M. Vandeginste, L. M. C. Buydens, S. De Jong, P. J. Lewi y J. Smeyers-Verbeke, Handbook of Chemometrics and Qualimetrics, Elsevier, Amsterdam, 1997, Capítulos 17 y 36. 2. D. M. Haaland y E. V. Thomas, Partial Least-Squares Methods for Spectral Analysis. 1. Relation to Other Quantitative Calibration Methods and the Extraction of Qualitative Information, Anal. Chem. 1988, 60, 1193-1202. 3. H. Wold, Estimation of principal components and related models by iterative least squares, en Multivariate Analysis (Ed., P.R. Krishnaiah), Academic Press, NY, 1966, pp. 391-420. 4. Las imágenes están tomadas de la página web www.cc.gatech.edu/ people/home/adjacent/. 5. H. C. Goicoechea, A. C. Olivieri, Determination of bromhexine in cough–cold syrups by absorption spectrophotometry and multivariate calibration using partial least-squares and hybrid linear analyses. Application of a novel method of wavelength selection, Talanta 1999, 49, 793-800. 6. H. C. Goicoechea, A. C. Olivieri, Enhanced Synchronous Spectrofluorometric Determination of Tetracycline in Blood Serum by Chemometric Analysis. Comparison of Partial Least-Squares (PLS-1) and Hybrid Linear Analysis (HLA) Calibrations, Anal. Chem. 1999, 71, 4361-4368.
96
CAPÍTULO 7 CALIBRACIÓN MULTIVARIADA MÍNIMOS PARCIALES (PLS)
PARTE
III:
CUADRADOS
El método de cuadrados mínimos parciales (PLS, por partial least-squares) pretende mejorar la técnica antes descrita (PCR) introduciendo los valores de las concentraciones de calibración en el cálculo de los factores. De esta manera, en PLS se emplean factores dependientes de la concentración, mientras que en PCR los factores eran independientes de la concentración. Debemos mencionar que existen dos tipos de métodos PLS: uno denominado PLS-1, que concentra su atención en un único analito a la vez, y otro llamado PLS-2, que permite calibrar y predecir las concentraciones de varios analitos simultáneamente. Esto último puede parecer a primera vista una ventaja, ya que PLS-1 debe repetirse para cada analito diferente de interés, pero representa por otro lado una gran desventaja, ya que PLS-1 permite optimizar las condiciones de trabajo para cada analito independientemente. En general, hoy en día se prefiere utilizar PLS-1 para la mayoría de las aplicaciones, y de aquí en más nos referiremos a PLS-1 simplemente como PLS. Un algorit mo i terativo p ara PCR Como se vio en el capítulo anterior, la técnica de descomposición en valores singulares permite obtener los factores espectrales de la matriz de datos instrumentales X. Sin embargo, desde el punto de vista computacional, calcular todos los factores constituye una pérdida de tiempo, ya que usualmente sólo se requieren los primeros factores, esto es, los que más contribuyen a la varianza espectral. En general, no es aconsejable utilizar un número de factores superior a la mitad del número de mezclas de calibración, por lo que resultaría sumamente útil disponer de una herramienta que permita calcular un número determinado de factores hasta un cierto límite máximo. Existen varias técnicas computacionales iterativas que permiten realizar esta operación, entre las que se destaca el algoritmo NIPALS (por non linear iterative partial least1 squares). La posibilidad de obtener los factores uno a uno permite plantear un algoritmo iterativo para PCR, que siga estos pasos: • Calcular el factor que explica la mayor parte de la varianza de X. • Descontar de X la parte explicada por el factor anterior, obteniendo el residuo E. • Volver al primer paso y reemplazar X por E, continuando hasta obtener el número deseado A de factores. Matemáticamente, este algoritmo se expresa del modo que sigue: 1) Se calcula el primer factor espectral de X, o primer loading u1. 2) Se proyecta la matriz de datos en este factor espectral, obteniéndose el primer T score t1, a través de t1 = X u1. 3) Se substrae de X la contribución del primer factor, calculada como u1 t1T, es decir se calcula la diferencia o residuo E = X – u1 t1T. 4) Se substituye E por X en el primer paso y se continúa hasta llegar al número deseado A de factores.
97
Los vectores ua y ta encontrados a cada paso de este algoritmo se reúnen en las matrices U y T discutidas anteriormente para PCR. Un algor itm o iterativo para PLS PLS opera de manera similar al algoritmo iterativo delineado para PCR. En PLS, sin embargo, existen dos clases de factores espectrales: unos llamados weigth loading factors, contenidos en una matriz usualmente designada W, y otros llamados simplemente loadings, contenidos en una matriz designada P. Las columnas de W son ortogonales, mientras que las de P no necesariamente lo son, a diferencia de PCR. Es importante recalcar que las columnas de W no son autovectores propiamente dichos, sino factores obtenidos mediante una técnica diferente a la de PCR, cuyos elementos dependen de las concentraciones de calibración del analito de interés. La obtención de estos factores se lleva a cabo mediante un algoritmo iterativo cíclico, muy similar al descrito anteriormente para PCR. La diferencia fundamental radica en que en PCR los factores describen la máxima varianza posible en la matriz de datos instrumentales únicamente, mientras que en PLS los factores describen la máxima correlación posible entre la matriz de datos y el vector de concentraciones del analito de interés. Matemáticamente, el algoritmo PLS se resume en los siguientes pasos: 1) Se proyecta la matriz de datos X en el vector de concentraciones yn, obteniéndose el primer weigth loading factor , que luego se normaliza a longitud unitaria: w1 = X yn / (ynT yn), seguido de normalización. 2) Se obtiene el primer score: t1 = XT w1 3) Se obtiene el primer coeficiente de regresión v1: v1 = t1 yn / (t1T t1) 4) Se obtiene el primer loading p1: T T p1 = X t1 / (t1 t1) 5) Se calculan los residuos espectrales y de concentración: T T T EX = X – t1 p1 ey = yn – v1 t1 6) Se sustituyen eX y ey por X e yn respectivamente en el paso 1) y se continúa hasta llegar al número de factores deseado A. A continuación describimos los pasos anteriores de manera cualitativa, en relación con los del algoritmo correspondiente a PCR: Paso 1). En este paso del algoritmo, se supone que sólo se conocen las concentraciones de un único componente, en este caso el analito 1, en las mezclas de calibración. Se efectúa un análisis similar al de CLS, pero suponiendo que sólo está presente el analito 1. En otras palabras, w1 es una aproximación por cuadrados mínimos al espectro puro del analito 1, similar a la que se hubiese realizado en CLS suponiendo la presencia de un único componente en la calibración. A diferencia de PCR, en este paso se aprecia la introducción de información concerniente a las concentraciones contenidas en yn en el cálculo del primer factor. Recuérdese que en PCR el primer factor se calcula a partir únicamente de la matriz X, con prescindencia de las concentraciones del analito. Paso 2). Se continúa con la suposición de que únicamente está presente el analito 1, y se calcula qué contribución del primer factor w1 está presente en las mezclas de calibración. Estas "concentraciones" forman el vector t1. Nótese que si efectivamente hubiese un único componente en la calibración, los pasos 1-2 serían idénticos a los realizados mediante un método CLS. En presencia de más de un componente, PLS se desvía del método clásico de análisis. 98
Paso 3). Este paso es similar al realizado en PCR. Se calcula el coeficiente de regresión que relaciona el score t1 calculado en el paso 2) con las concentraciones de calibración. Pasos 4 y 5). Estos pasos aseguran que los vectores ta y wa subsiguientes serán ortogonales entre sí. Para ello se calculan los vectores pa, llamados loadings. Estos vectores, a diferencia de PCR, no explican la varianza espectral en la matriz X, sino que representan un intento de explicar dicha varianza, mientras simultáneamente se correlacionan los scores ta con las concentraciones yn. Calibración La etapa de calibración requiere estimar en primer lugar el número óptimo de factores A, lo que usualmente se lleva a cabo mediante la técnica de validación cruzada, tal como se describió para PCR. El resultado de la calibración es la obtención del vector de coeficientes de regresión vn, cuyos elementos ( v1,..., v A) se obtienen en cada uno de los A pasos del algoritmo cíclico anterior. Predicción En la etapa de predicción se emplean los coeficientes de regresión para estimar la concentración del analito en la muestra. El paso previo, tal como en PCR, es la obtención de los scores de la muestra, lo que se realiza con ayuda de las matrices W y P: T
–1
T
t A = (W P) W x
(7.1)
T
(7.2)
yn = (vn) t A
Residu os espectrales y cifr as de mérito En PLS también se estima el espectro de la muestra incógnita, de manera que pueden calcularse residuos espectrales, en forma análoga a PCR. La estimación del espectro de la muestra se realiza mediante la siguiente ecuación: xˆ
= P t A
(7.3)
Y luego puede definirse el residuo espectral: J
∑=1 ( x
j
sres =
− xˆ j ) 2
j
J − A
(7.4)
Ventajas y desvent ajas de PLS PLS es el método de calibración multivariada más empleado cuando la información instrumental proveniente de cada muestra es de tipo vectorial (un espectro de absorbancia es el ejemplo típico). En este sentido, su desarrollo ha superado de algún modo a PCR, incorporando información útil referida a concentraciones de calibrado durante la etapa de cálculo de las variables latentes. Sin embargo, n referencia a la tabla de propiedades analíticas presentada en el Capítulo 4, se mantiene el problema de las interferencias no modeladas. Este problema es común a la mayoría de los métodos multivariados basados en información espectral: 99
si aparece en una muestra incógnita un compuesto no contenido en la calibración, el análisis no será exacto. Más allá de PLS PLS es probablemente el más usado de los métodos quimiométricos para calibración multivariada utilizando datos vectoriales. Sin embargo, en los últimos años se han desarrollado varios competidores, desde variantes cosméticas de PLS hasta metodologías completamente disímiles. El lector interesado en algunos de estos métodos alternativos puede consultar la bibliografía específica. 2-7 Por otro lado, si desea emplear un programa completo de MATLAB, capaz de implementar varios métodos quimiométricos con una serie interesante de recursos gráficos, de preprocesamiento de los datos, etc. puede consultar la referencia reciente de nuestro grupo de trabajo, y obtener el programa de internet, junto con juegos modelo de datos.8 Material suministrado Para este capítulo se provee la siguiente información complementaria: Archivos (*.M) con rutinas para el entorno de programación MATLAB. Archivos (*.EXE) con programas ejecutables en QB.
Ejercicio 1) Los datos del presente ejercicio están tomados del trabajo número 10 de las referencias. Se desea determinar el contenido de un fármaco, la bromhexina, presente en muestras de jarabe para la tos. Los componentes del jarabe se conocen en forma incompleta, de manera que se preparan muestras para construir un modelo PLS. Para ello, se agregan cantidades conocidas de bromhexina a doce muestras de jarabe “blanco” de diferentes cpncentraciones (esto es, el fondo de la matriz del jarabe, sin bromhexina), y se utilizan para calibrar el modelo. Las concentraciones del analito en las muestras de calibrado son:
Tabla 7.1. Concentración de las muestras de calibración de jarabe de bromhexina Muestra de calibrado Concentración ×104 M 1 1.55 2 2.06 3 2.58 4 1.55 5 2.06 6 2.58 7 1.55 8 2.06 9 2.58 10 1.68 11 2.10 12 2.66 100
Estas concentraciones se recogen en forma de un vector de 12 ×1 en el archivo de texto BR_CON_C.TXT Los espectros de absorción de estas 12 muestras se registran a 64 diferentes longitudes de onda. Estos espectros están contenidos, en forma de matriz de 64 ×12, en el archivo de texto BR_RES_C.TXT. Informar las correspondientes cifras de mérito para el modelo. Suponga que el nivel de ruido instrumental es igual a 0,003 unidades de señal. 2) Para la validación del modelo, se prepararon 11 muestras adicionales de jarabe con contenido conocido de bromhexina, diferente al empleado para calibrar. Los espectros de estas muestras están contenidos, en forma de matriz de 64 ×11, en el archivo BR_RES_T.TXT, y las concentraciones nominales, en forma de vector de 11 ×1, en el archivo BR_CON_T.TXT. Estimar las concentraciones de los analitos en este juego de muestras y sus incertidumbres asociadas, y estudiar la exactitud del método mediante la prueba de la elipse. 3) Una muestra adicional de prueba, cuyo espectro está contenido en el archivo de texto BR_RES_P.TXT se analiza mediante el mismo modelo. Sin embargo, se sospecha que se trata de una muestra que contiene una interferencia no modelada en la calibración. ¿Qué conclusiones puede extraer al respecto del análisis mediante PLS?
Respuesta 1) El primer paso en el análisis PLS debe ser el estudio del número óptimo de factores presentes en la matriz de calibrado, que luego se emplearán para la predicción. El método más recomendado para esto es la validación cruzada, que puede implementarse mediante la rutina PLS_CV.M de Matlab o el programa PLS_CV.EXE de QB. Para ejecutar estos algoritmos, se requiere introducir un número máximo de factores de prueba. Este puede ser, como máximo, igual al número de mezclas de calibrado menos una (ya que el procedimiento consiste en calibrar con las muestras menos una), en el presente caso 11 = 12 – 1. No obstante, se supone que se han preparado más mezclas de calibración que fuentes de variación espectral, por lo que se recomienda introducir, como número máximo, un valor menor. Los resultados obtenidos para un número máximo de factores igual a ocho son los siguientes: Tabla 7.2. Selección del némero de factores para la calibración de bromhexina Factores PRESS PRESS/min(PRESS) p 1 0,907 249002 0,999 2 0,021 5,76 0,997 –3 3 1,139 0,587 4,15×10 –3 4 1 0,5 3,65×10 –3 5 – – 5,81×10 6 – – 9,62×10 –3 7 – – 1,64×10 –2 8 – – 2,06×10 –2 Puede apreciarse que el PRESS disminuye al ir aumentando el número de factores, llega a un mínimo para 4 factores, y luego aumenta. El número óptimo de factores, obtenido para el primer valor de p que disminuye por debajo de 0,75 es 3. 101
El RMSECV para 3 factores es satisfactorio (0,02) en vista de las concentraciones nominales de calibrado y sus incertidumbres asociadas (en la segunda cifra decimal). Estos primeros tres componentes principales explican más del 99,99% de la varianza de la matriz espectral. Tanto los resultados correspondientes al PRESS como la varianza explicada se observan gráficamente en la figura generada por MATLAB, figura que también puede construirse mediante los valores provistos por el programa QB correspondiente (PLS_CV.EXE).
Figura 7.1. Variación del PRESS en función del número de factores para la
bromhexina. Una vez establecido el número óptimo de factores para la compresión de la información, se procede a calibrar el modelo, empleando los programas PLS_CAL.M (Matlab) o PLS_CAL.EXE (QB). Las cifras de mérito calculadas mediante los programas para este modelo son las siguientes: Tabla 7.3. Cifras de mérito para el modelo PLS de la bromhexina. Cifra de mérito Valor Sensibilidad 1,21×104 A × M –1 Sensibilidad analíticaa 4×106 M –1 1/γ 2,5×10 –7 M Selectividad 0,46 a Obtenida dividiendo la sensibilidad por el nivel de ruido instrumental (0,003 unidades).
102
2) Para predecir las concentraciones de las muestras incógnitas, empleamos los programas PLS_PRED.M (Matlab) o PLS_PRED.EXE (QB), con los siguientes resultados: Tabla 7.4. Concentraciones predichas en las muestras incógnitas y residuo espectral para el modelo PLS para la bromhexina. Muestra
Concentración × 104
Residuo espectral
Nominal Predichaa 1 1,96 1,97(1) 0,004 2 2,16 2,19(1) 0,002 3 0,00 –0,01(1) 0,014 4 0,82 0,84(1) 0,009 5 1,02 1,04(1) 0,007 6 1,33 1,37(1) 0,005 7 1,84 1,93(1) 0,003 8 2,35 2,43(1) 0,004 9 1,94 1,99(1) 0,004 10 2,14 2,19(1) 0,002 11 2,24 2,25(1) 0,005 a Los errores estándar en las concentraciones, calculados con el modelo aproximado citado en la teoría, esto es s( xn) = s y / SENn, con s y = 0,003, son todos iguales a 0,002. Este valor es demasiado optimista, en vista de que las concentraciones de calibrado están dadas con una incertidumbre de 0,01, por lo que se ha optado por este último valor, más conservador, en la presente tabla. Se informan también, en la última columna de esta tabla, los residuos espectrales para cada muestra incógnita, que, como puede apreciarse, se encuentran dentro del nivel del ruido instrumental. Esto confirma que el ajuste por cuadrados mínimos para estas muestras es adecuado. Dos excepciones a esta situación son las muestras número 3 y 4, cuyo residuo espectral es superior al resto. Una explicación posible para esto es que estas muestras fueron preparadas con una concentración nominal inferior a las de calibrado. En este sentido, no se trata de verdaderos outliers que contengan interferencias no modeladas, sino de muestras para las que le estamos exigiendo al modelo que realice una extrapolación hacia un ambiente para el que no está entrenado. De todas maneras, nótese que las concentraciones predichas para estas muestras son muy cercanas al valor nominal. Para establecer la exactitud del método, lo recomendado es analizar los datos de la tabla precedente mediante la prueba de la elipse, tal como se discutiera en el capítulo 2. De este modo, la tabla de datos a suministrar a los programas de cálculo de la elipse será como sigue:
103
Tabla 7.5. Datos a suministrar a los programas de cálculo para construir la elipse 1,96 1,97 2,16 2,19 0,00 –0,01 0,82 0,84 1,02 1,04 1,33 1,37 1,84 1,93 2,35 2,43 1,94 1,99 2,14 2,19 2,24 2,25 Dado que no se tienen resultados de réplicas de cada muestra, lo que proveería una estimación del desvío estándar de cada valor predicho, realizaremos un análisis mediante el método OLS. Se recomienda organizar los datos apropiadamente y someterlos a los programas EJCR.M (MATLAB) o EJCR.EXE (QB). El resultado se muestra en las figuras siguientes. La primera de ellas muestra los valores predichos en función de los nominales, y la segunda la elipse. 3) El análisis de la muestra contenida en el archivo BR_RES_P.TXT arroja los siguientes resultados: Concentración estimada: 2,10 Residuo espectral: 0,08 Aquí el residuo es significativamente mayor que el ruido espectral, lo que haría sospechar la presencia de un interferente no modelado.
Referencias 1. H. Wold, Estimation of principal components and related models by iterative least squares, en Multivariate Analysis (Ed., P.R. Krishnaiah), Academic Press, NY, 1966, pp. 391-420. 2. H. C. Goicoechea y A. C. Olivieri, A comparison of orthogonal signal correction and net analyte preprocessing methods. Theoretical and experimental study, Chemom. Intell. Lab. Syst. 2001, 56, 73. 3. O. Svensson, T. Kourti y J. F. MacGregor, An investigation of orthogonal signal correction algorithms and their characteristics, J. Chemometrics, 2002, 16 , 176. 4. S. Wold, H. Antti, F. Lindgren y J. Öhman, Orthogonal signal correction of nearinfrared spectra, Chemom. Intell. Lab. Syst . 1998, 44, 175. 5. T. Fearn, On orthogonal signal correction, Chemom. Intell. Lab. Syst . 2000, 50, 47. 6. L. Xu e I. Schechter, A calibration method free of optimum factor number selection for automated multivariate analysis. Experimental and theoretical study, Anal. Chem. 1997, 69, 3722. 7. A.C. Olivieri, Calibración multivariada. Introducción a la programación en MATLAB, Ediciones Científicas Argentinas, Buenos Aires, 2001. 8. El programa MVC1 ( Multivariate Calibration 1) puede obtenerse libremente en www.chemometry.com
104
9.
R. Bro, Multiway Analysis in the Food Industry. Models, Algorithms, and Applications, Royal Veterinary and Agricultural University Denmark, 1998, disponible en internet en www.models.kvl.dk . 10. H. C. Goicoechea, A. C. Olivieri, Determination of bromhexine in cough–cold syrups by absorption spectrophotometry and multivariate calibration using partial least-squares and hybrid linear analyses. Application of a novel method of wavelength selection, Talanta 1999, 49, 793-800
105
CAPÍTULO 8 CALIBRACIÓN MULTIVARIADA PARA RESPUESTAS NO LINEALES RESPECTO DE LA CONCENTRACIÓN DEL ANALITO Redes neuron ales arti fic iales El desarrollo de las redes neuronales artificiales se basa en un modelo matemático que pretende imitar el funcionamiento del cerebro. Estos algoritmos contienen elementos de cálculo aritmético que equivalen a las neuronas cerebrales, las células que procesan la información en el cerebro. La red neuronal equivale, en general, a un conjunto de neuronas conectadas entre sí, por lo que se las conoce como redes neuronales artificiales (ANNs, del inglés Artificial Neural Networks).1 Las ANNs están constituidas por nodos o unidades unidas mediante conexiones. A cada conexión se le asigna un peso numérico. Los pesos constituyen el principal recurso de memoria en las redes y el aprendizaje normalmente se realiza actualizando dichos pesos. 2,3 Las unidades constan de un conjunto de conexiones de entrada provenientes de otras unidades, un conjunto de vínculos de salida que van hacia otras unidades, un nivel de activación del momento, y recursos para calcular cuál será el nivel de activación del siguiente paso, sobre la base de sus entradas y pesos respectivos. Lo importante es que en cada una de las unidades se efectúa un cálculo local basado en las entradas que le proporcionan sus vecinas, pero sin que sea necesario un control global en todo el conjunto de las unidades.1–3 Elementos si mples de cálcu lo En la Figura 8.1 se muestra una unidad típica, donde cada una de las unidades realiza un cálculo sencillo. Esta unidad recibe señales de sus vínculos de entrada y calcula el correspondiente nuevo nivel de activación, que envía a través de sus vínculos de salida. El cálculo del nivel de activación está basado en los valores de cada una de las señales de entrada que envía un nodo vecino, así como los pesos de cada uno de los vínculos de entrada. El cálculo está dividido en dos componentes. El primero es un componente lineal, denominado función de entrada ( ent i), que calcula la suma ponderada de los valores de entrada de la unidad. El segundo es un componente no lineal, conocido como función de activación ( g), que transforma la suma ponderada en el valor final, que sirve como valor de activación de la unidad ( ai). Por lo general, todas las unidades de la red utilizan la misma función de activación. 1–4 La entrada ponderada total es la suma de las activaciones de entrada, multiplicadas por sus pesos respectivos: ent i =
∑w
j ,i
T
a j = w i a i
(8.1)
j
Los pesos correspondientes al nodo i están designados por el vector wi y al conjunto de valores de entrada se lo denomina ai.
106
Figura 8.1. Representación esquemática de una unidad de activación o neurona.
El nuevo valor de activación de la unidad se calcula mediante un paso de cómputo elemental de cada una de las unidades. Para ello, se aplica la función de activación g al resultado de la función de entrada: (8.2) ai ← g (ent i ) = g ( ∑ w j ,i a j ) Al usar diferentes funciones matemáticas para g se obtienen distintos modelos. Tres de los modelos más comunes son las funciones escalón, signo y sigmoidea, cuyas representaciones esquemáticas se ilustran en la Figura 8.2. La función escalón tiene un límite t , de manera que produce un resultado igual a 1 cuando la entrada es mayor que su límite y, en caso contrario produce 0. El equivalente de la motivación biológica es que 1 representa el envío del impulso al axón, y el 0 representa que éste no es enviado. El límite, o umbral, representa la mínima entrada total ponderada necesaria para provocar la activación de la neurona. También pueden definirse versiones límite de las funciones signo y sigmoidea. Matemáticamente: Función escalón x = 1, si ent ≥ t 0, si ent < t Función signo x = +1, si ent ≥ 0 –1, si ent < 0 -ent Función sigmoidea x = 1 / ( 1 + e )
107
Figura 8.2. Representación gráfica de tres tipos
de funciones para activación de
unidades. Estructura de red óptima La elección de la estructura de la red es un hecho de gran importancia. Si la red elegida es demasiado pequeña, el modelo no podrá representar a la función deseada. Si, por otra parte, escogemos una red demasiado grande, aunque será capaz de memorizar todos los ejemplos, no será capaz de generalizar bien cuando se presenten entradas que no haya visto anteriormente. Es decir, al igual que todos los modelos estadísticos, las redes neuronales están sujetas a un sobre-ajuste cuando los parámetros (esto es, los pesos del modelo) son demasiados. Más adelante describiremos cómo se llega a obtener la arquitectura óptima en una red usada en calibración multivariada, y cómo se trata de evitar el sobre-ajuste cuando se entrena una red. Redes neuronales arti fici ales en calibración mu ltiv ariada El uso de ANNs en calibración multivariada se restringe casi exclusivamente a las llamadas redes multicapas con retroalimentación o perceptrón multicapa (MLP, del inglés multi-layer perceptron), por lo que de aquí en adelante se describirá el funcionamiento de este tipo de redes. En la Figura 3 se muestra un ejemplo de red para cuatro variables descriptoras ( x1, x2, x3 y x4) y una sola salida ( y). Las variables se presentan a las cuatro neuronas de entrada, y son luego pesadas por las conexiones (wij’). Posteriormente, las neuronas escondidas reciben las señales ponderadas y realizan dos tareas: una suma de las entradas y una proyección de esta suma sobre la función de transferencia ( f h), produciendo la activación. Luego, estas activaciones son pesadas por las conexiones ( w j’’) entre las neuronas escondidas y las de salida, que realizan un tarea similar a las neuronas escondidas, proyectando la suma sobre la función de transferencia ( f o). La salida es la respuesta estimada ( y pred), que puede expresarse según la siguiente ecuación: nh ⎡ ⎛ nd ⎞⎤ (8.3) ypred = f o ⎢θ ' '+ ∑ w' ' j f h ⎜ ∑ w'ij xi + θ ' ⎟⎥ 1 1 j i = = ⎝ ⎠ ⎣ ⎦ donde nd y nh son los números de neuronas de entrada y escondidas respectivamente. Aunque la red se considera una herramienta no paramétrica, los modelos están definidos por un conjunto de parámetros ajustables, determinados por un algoritmo. Estos parámetros ajustables son los pesos ( wij’ y w j’’) y los umbrales ( biases) (θ ´ y θ ´´ ) que actúan desplazando horizontalmente las funciones de transferencia. El procedimiento 108
iterativo que permite determinar estos parámetros se llama aprendizaje o entrenamiento (learning o training en inglés). En primer lugar, estos parámetros se eligen al azar, y se calcula un valor de la respuesta. Luego, en un paso de retroalimentación, el error entre la respuesta real y la predicha se minimiza por ajuste de los parámetros. Cada par de estos pasos se considera una época. El procedimiento se repite hasta la convergencia, es decir, hasta que se llega a un error pre-especificado o aceptablemente bajo. Resumiendo, el problema consiste en una optimización en la que se busca un mínimo de error en una superficie multidimensional caracterizada por la presencia de mínimos locales. Probablemente no se encontrará el mínimo absoluto, pero sí uno local que satisface la solución del problema. 5 El algoritmo de minimización más utilizado es el del gradiente descendente sucesivo, que se basa en la estimación de la primera derivada del error con respecto a cada peso. 6 En el esquema de la figura 8.3 puede apreciarse como se realiza el procedimiento.
Figura 8.3. Representación gráfica de una red
multicapa con retroalimentación.
Las ANNs pueden utilizarse para construir modelos de calibración multivariada empíricos de la forma y = f (X) + e. Sólo se considerarán los modelos de calibración inversa, donde X designa una matriz de medidas o señales analíticas formada por n muestras. Para una muestra dada, las mediciones son descriptas por un juego de variables descriptoras, xi, un ejemplo lo constituyen los valores de absorbancia a una longitud de onda dada (espectro de absorción). Por su parte, y es un vector conteniendo las respuestas de las muestras, por ejemplo, las concentraciones de un analito de interés en un juego de mezclas de calibración. 5 Es importante remarcar que las ANNs se deben utilizar cuando un juego de datos se sospecha o se sabe que es no lineal. Desde un punto de vista matemático, un verdadero modelo es no lineal respecto de sus parámetros. En calibración multivariada, también suele considerarse no lineal a un modelo en el que la no linealidad se presenta en la relación entre la respuesta y el descriptor. En espectroscopía, la desviación de la linealidad puede ocurrir si una muestra tiene una elevada absorción o no es homogénea, si el tamaño de la partícula no es constante en todas las muestras (para especies cristalinas) o si algunas señales están solapadas. 6 Una respuesta no lineal, o la presencia de luz parásita, introducen una curvatura en la función de respuesta de la concentración. 109
Estos efectos de no linealidad pueden ser corregidos algunas veces con el uso apropiado de técnicas de preprocesamiento de los datos, aunque estas técnicas presentan limitaciones como la pérdida de sensibilidad. 7,8 Se pueden presentar otros efectos nolineales como equilibrios químicos asimétricos, reacciones intermoleculares, reacciones intermetálicas en electroquímica, uniones puente de hidrógeno y cambios en temperatura o composición de los solventes. Desarrollo de modelos de c alibración 1. Detección de no linealidad
Como regla general, no se debería construir un modelo ANN si no se tiene seguridad de que el sistema es no lineal. Existen diversas herramientas que permiten establecer si los datos presentan un comportamiento no lineal. El método más simple, y en algunos casos suficiente, consiste en realizar una inspección de los residuos del modelo lineal versus la concentración. También ayuda mucho graficar la concentración nominal vs. la predicha por el modelo lineal.5 Algunos autores presentan métodos complejos que han sido aplicados con éxito en algunos casos. 9,10 2. Muestras de entrenamiento, monitoreo y validación.
El número de muestras disponibles suele resultar un factor limitante para la aplicación de ANNs. Si se dispone de pocas muestras, el número de parámetros a ser ajustados es tal que conduce a sobre-ajustes. El número de muestras usadas debería ser superior a 30, pero se pueden obtener buenos resultados sólo con 15. 5 Debe considerarse un aspecto importante de las ANNs, y es la necesidad de contar con un juego de muestras de monitoreo a fin de evitar el sobre-ajuste. Ocurre que las redes pueden ser entrenadas hasta un error muy bajo. Ese hecho, si bien dice que la calibración realizada es buena, implica la pérdida de capacidad de generalizar, y esto conduce a poca habilidad predictiva del modelo (sobre-ajuste). La evaluación del error obtenido en la predicción del juego de monitoreo al mismo tiempo que el de aprendizaje, cuando se desarrolla cada época, permite encontrar un punto de corte. Esto quiere decir que pueden encontrarse los pesos que originan errores de calibración y de monitoreo aceptables al mismo tiempo. La Figura 8.4 muestra cómo evolucionan los errores y cuál sería la época seleccionada de acuerdo con el criterio comentado.
Figura 8.4. Representación gráfica de la evolución típica de los errores de calibración y
monitoreo en función del número de iteraciones. 110
Resumiendo, a diferencia de los modelos de calibración como PLS, por ejemplo, las ANNs necesitan contar con muestras adicionales agrupadas en el juego de monitoreo. Algunos autores utilizan una parte del juego de muestras de calibración para monitorear. La selección se realiza al azar o aplicando herramientas como el algoritmo de Kennard-Stone, que permiten dividir los juegos de una manera inteligente, contemplando la variabilidad del sistema bajo estudio. 11 3. Arquitectura
Generalmente la arquitectura se representa con números que indican la cantidad de neuronas en cada capa. La red de la Figura 8.3 se representa entonces como 4–3–1. El número de capas escondidas es generalmente 1. No se han encontrado mejoras en los resultados de calibración para más de una capa escondida. También es posible usar conexiones directas entre las neuronas de entrada y la de salida, pero esta conexiones generalmente tienen en cuenta las relaciones lineales entre las variable y permiten un aprendizaje más rápido. Es preferible que el número de neuronas en la capa escondida sea lo más pequeño posible para evitar sobre-ajuste. El número de neuronas de entrada puede obtenerse empezando con un número deliberadamente grande, y quitando gradualmente algunas ( pruning) hasta que la calidad del ajuste deja de mejorar. Algo importante a tener en cuenta es que en calibración multivariada se necesitarían tantas neuronas como sensores hay en el espectro (por ejemplo longitudes de onda). La solución se encuentra aplicando análisis de componentes principales (PCA) y utilizando las puntuaciones ( scores ) correspondientes a los autovectores que explique la máxima varianza del sistema. El número de neuronas de salida debería ser uno, ya que se modela una respuesta a la vez. Podría existir la excepción de querer modelar respuestas correlacionadas, como las concentraciones de diferentes constituyentes de una mezcla en un sistema cerrado. Entrenamiento de una red
El entrenamiento de una red es un problema de optimización. El algoritmo más utilizado es el del gradiente descendente sucesivo ( steepest-descent ), que consiste en una minimización sucesiva en la superficie del error, ajustando los parámetros en un hiperespacio.12 La tendencia a caer en mínimos locales obliga a la utilización de un término de “momento” en la actualización de los pesos. Esto permite suavizar la superficie del error y atenuar las oscilaciones en el fondo de los valles. La velocidad del algoritmo se puede aumentar usando dos parámetros adaptativos: velocidad de aprendizaje y momento. Modelos basados en métodos l ineales Existen algunos métodos para modelar respuestas no lineales. Estos modelos se basan en la adaptación de los conocidos PCR y PLS, aunque también hay otros modelos que han sido creados especialmente para resolver el problema de la no linealidad. Entre ellos se pueden citar: regresión localmente ponderada (LWR, del inglés, locally 13 weighted regression), regresión por búsqueda de la proyección (PPR, del inglés, 14 projection pursuit regression), expectativa condicional alternante (ACE, del inglés, 14 alternating conditional expectation), regresión multivariada adaptativa basada en funciones splines (MARS, del inglés, multivariate adaptative regression splines)14 y las ya comentadas redes neuronales artificiales. Además debe considerarse que PCR y PLS pueden contemplar no linealidades por la inclusión de mayor número de variables latentes en el modelo que las requeridas para un sistema lineal o, como ya dijimos usando las versiones no lineales NL-PLS (del inglés, non linear PLS ) y SPL-PLS (del 111
inglés, spline PLS ).14 En este capítulo se describirán los dos últimos métodos, que son los que más se han usado en la literatura. NL-PLS y SPL-PLS
Como ya vimos con anterioridad, los métodos de calibración multivariada tienen como objetivo encontrar un vector de regresión ( bk ) a partir de cierta información (calibración), que permita a su vez encontrar la concentración del analito de interés ( ck ) a partir del espectro (r) de una muestra incógnita (predicción), según la siguiente ecuación: T (8.4) ck = bk r En calibración inversa, bk se obtiene a partir de la pseudoinversa de la matriz R que contiene los espectros de calibración y del vector que contiene las concentraciones del analito de interés ( k ) en esas mezclas de calibración ( ck ). La Ecuación (8.5) muestra esta relación: bk = R + ck (8.5) De cómo se obtiene la pseudoinversa, surgen los diferentes modelos, entre ellos PCR y PLS. En dicho procedimiento se calculan variables latentes y puntuaciones o scores que permiten el cálculo de bk . La principal diferencia entre PLS y PCR es que el primero incorpora información de la concentración en el cálculo de dichas variables latentes. Para el cálculo de la pseudoinversa, las respuestas de calibración son descompuestas según la siguiente ecuación: R = T P+ (8.6) donde T representa la matriz de scores y P la matriz de variables latentes. En el esquema de la Figura 8.5 puede apreciarse un diagrama de flujo en el que un vector r es afectado a las variables latentes ( p) y a los scores (t), en este caso se calcularon dos factores, para finalmente encontrar la concentración del analito:
Figura 8.5. Diagrama de flujo que muestra la transferencia de la información durante la
etapa de predicción en PCR y PLS. Para implementar las versiones no lineales, se procede a utilizar una matriz de scores aumentada con la contribución de las transformaciones cuadráticas de dichos scores. La Ecuación (8.7) muestra esta matriz aumentada ( Taum) para dos componentes, mientras que la Figura 8.6 muestra el correspondiente diagrama de flujo: 2 2 Taum = [t1 t2 t1 t2 ] (8.7)
112
Figura 8.6. Diagrama de flujo de los algoritmos PCR y
PLS no lineal durante la etapa
de predicción. Se puede incorporar cualquier función en la transformación de los vectores de scores para construir el modelo que más convenga. 14 Por otra parte, se puede evitar la especificación exacta de una transformación aplicando splines, originando el método SPL-NPLS.15 Referencias
1.
S. Russel, P. Norvig, Inteligencia artificial. Un enfoque moderno, 1ra ed., Ed. Prentice may Hispanoamericana S. A., páginas 595-627, México, 1996. 2. B. Kröse, P. van der Smagt, An introduction to neural networks, University of Ámsterdam, Ámsterdam1, 1996. 3. N.K. Kasabov, Foundation of neural netwoks, fuzzy systems and knoledge engineering, The MIT PRESS, Cadbridge, 1998. 4. J. Zupan, J. Gasteiger, Neural networks in chemistry and drug design, Second Edition, Ed. Wiley-VCH, Wenheim, 1999. 5. F. Despagne, D.L. Massart, Neural networks in multivariate calibration, Analyst 1998, 123, 157R. 6. R. Fletcher, Practical methods for optimization, Vol 1: Uncontrained optimization, Wiley, New Cor, 1980. 7. P.J. Gemperline, J.R. Long, V.G. Gregorious, Nonlinear multivariate calibration using principal componet regresión and artificial neural networks, Anal. Chem. 1991, 63, 2313. 8. L. Hadjiiski, P. Geladi, P. Hopke, A comparison of modeling nonlinear systems with artificial neural networks and partial least squares, Chemom. Intell. Lab. Syst. 1999, 49, 91. 9. V. Centner, O.E. de Noord, D.L. Massart, Detection of nonlinerity in multivariate calibration, Anal. Chim. Acta 1998, 376 , 153. 10. M.S. Collado, M.L. Satuf, H.C. Goicoechea, A.C. Olivieri, Complementary use of partial least-squares and artificial neural networks for the non-linear spectrophotometric analysis of pharmaceutical samples, Anal. Bioanal. Chem. 2002, 374, 460. 11. R.W. Kennard, L.A. Stone, Computer aided design of experiments, Technometrics 1969, 11, 137. 12. D.E. Rumelhart, J.L. McClelland, Parallel distributed processing, MIT Press, Cambridge, 1986 113
13. V. Centner, D.L. Massart, Optimization in locally weighted regression, Anal. Chem. 1998, 70, 4206. 14. S. Sekulic, M.B. Seasholtz, Z. Wang, B.R. Kowalsky, S.E. Lee, B.R. Holt, Nonlinear multivariate calibration methods in analytical chemistry, Anal. Chem. 1993, 65, 835A. 15. S. Wold, Nonlinear partial least squares modelling. II. Spline inner relation, Chemom. Intell. Lab. Syst. 1992, 14, 71. 16. M.S. Camara, F.M. Ferroni, M. De Zan, H. Goicoechea, Sustained modelling ability of artificial neural networks in the analysis of two pharmaceutical (dextropropoxyphene and dipyrone) present in unequal concentrations, Anal. Bioanal. Chem. 2003, 376 , 838-843.
114
CAPÍTULO 9 RECIENTES DESARROLLOS MULTIVARIADA
EN
CALIBRACIÓN
Nuevas tendencias en calibración multivariada. Calibración mult ivariada de orden superior a uno Debe mencionarse que los métodos para calibración multivariada descritos hasta el momento se basan en el procesamiento de datos del tipo vectorial, es decir, espectros, u otro tipo similar de datos instrumentales (voltamperogramas, por ejemplo). Una calibración basada en vectores para cada muestra se llama calibración de primer orden, debido a que un vector se considera, en lenguaje tensorial, como un tensor de primer orden.1 En este sentido, la calibración univariada se clasificaría como de orden cero.1 Existe la posibilidad de realizar una calibración empleando datos matriciales para cada muestra, por ejemplo, matrices de excitación-emisión (obtenidas fácilmente en un espectrofluorómetro convencional), matrices de absorbancia-tiempo (obtenidas a través de una reacción química en un espectrofotómetro de arreglo de diodos), etc. En este caso, la calibración se denomina de segundo orden, dado que una matriz es un tensor de segundo orden. 1 No existe límite teórico para el orden, y recientemente se han descrito en la literatura calibraciones utilizando datos de tercer orden (matrices de excitaciónemisión de fluorescencia combinadas con la cinética de una reacción química). La Figura 9.1 muestra los diferentes tipos de datos que se pueden obtener junto con las calibraciones a las que dan origen.
Figura 9.1. Esquema mostrando los tipos de
datos que se pueden obtener junto con las calibraciones a las que éstos dan origen.
La calibración de orden superior (segundo, tercero, etc.) presenta ventajas adicionales a las descritas en este curso, en particular, la llamada ventaja de segundo orden, que permite cuantificar analitos calibrados en presencia de interferencias no calibradas. Esta propiedad está ausente en los datos de primer orden, y presenta 115
inmensas posibilidades en el análisis de mezclas complejas, en particular las de origen biológico. Una descripción detallada acerca de los métodos de orden superior puede encontrarse en la tesis de R. Bro 2 y en un libro de reciente aparición que lo cuenta entre sus autores. 3 Desde un punto de vista instrumental, los datos de segundo orden pueden ser obtenidos de dos maneras diferentes: a) acoplando dos instrumentos de primer orden, lo que se conoce como técnicas acopladas o en tandem (cromatografía gaseosacromatografía gaseosa, cromatografía gaseosa-detector de masas, cromatografía líquidadetector de arreglo de diodos, etc., o b) usando un instrumento de segundo orden como un espectrofluorómetro que registra matrices de excitación emisión (EEMs) o un espectrofotómetro con arreglo de diodos que registra la evolución de una reacción cinética. La Figura 9.2 muestra una matriz de excitación emisión típica, que representa uno de los datos más utilizados por la rapidez con que pueden ser obtenidos. 4
Figura 9.2. Superficie obtenida por un
espectrofluorómetro que registra espectros de emisión y excitación en dos dimensiones (EEM).
Diferentes mod elos ut ilizados. Ventaja de segundo o rden Entre los diferentes métodos presentados en la literatura para modelar datos de segundo orden, no todos son capaces de aprovechar la seductora propiedad conocida como ventaja de segundo orden. Una revisión a la literatura científica de los últimos años muestra entre los más usados y estudiados a los siguientes métodos: PARAFAC (de análisis paralelo de factores), GRAM (de método de aniquilación del rango), MCRALS (de resolución muiltivariada de curvas con cuadrados mínimos alternados), N-PLS (de PLS multidimensional), U-PLS (de PLS desdoblado o unfolded ), BLLS/RBL (de cuadrados mínimos bilineales seguidos de bilinearización residual) y PLS/RBL (de PLS seguido de bilinearización residual). Es importante resaltar que tanto N-PLS como UPLS no son capaces de explotar la ventaja de segundo orden, mientras que de los restantes, los que usan RBL, operan de manera diferente a los que no la usan. 2 Consideraremos como ejemplo al modelo PARAFAC para intentar explicar la manera en que este método permite predecir analitos en muestras, algunos de cuyos componentes no fueron modelados en las mezclas de calibración. El modelo PARAFAC agrega la matriz con los datos de la muestra incógnita al juego de datos de calibración. De esta manera, el modo I (correspondiente al número de muestras de calibración) se ve incrementado en una unidad ( I +1 en la Figura 9.3). Los modos restantes ( J y K ) permanecen iguales. Estos últimos modos pueden corresponder 116
a las longitudes de onda de emisión y de excitación, si se tratara de EEMs. Así es como se agrega información sobre las posibles interferencias de la muestra en la matriz de calibración, aunque no es necesario conocer de que interferencias se trata, ni su concentración.
Figura 9.3. Esquema que representa como se realiza la
calibración en PARAFAC.
La Figura 9.4 trata de explicar esquemáticamente como se descompone el arreglo de tres vías según un modelo PARAFAC, generando tres matrices loadings (A, B y C), cuyos elementos son aif , b jf , y ckf respectivamente, donde f representa el número de factores usados en la descomposición, en este caso 2. La matriz de datos X puede redefinirse según cada uno de sus elementos, como: ∧
X ijk =
F
∑a b c
if jf kf
(9.1)
f =1
=
Figura 9.4. Esquema que muestra como PARAFAC realiza la descomposición de los
datos. Los perfiles obtenidos en la descomposición de X, corresponden a los espectros en ambas dimensiones de los N analitos puros presentes en la mezcla. Estos perfiles son utilizados para construir un modelo de regresión lineal simple para cada analito en particular; lo que transforma finalmente la calibración multivariada en una calibración pseudo-univariada. 117
A partir de la ecuación 9.1, una vez resuelta la descomposición matricial y los perfiles de cada uno de los componentes ha sido calculado, los scores para las muestras de calibración asociados con un determinado analito son linealmente relacionados con las concentraciones nominales del analito, según: (an1 | ani |…….| anI ) = k PARAFAC y
(9.2)
donde n corresponde a un determinado componente, y es un vector que contiene las concentraciones nominales del analito n en las I mezclas de calibración, y k PARAFAC la constante de proporcionalidad calculada. Como el modelo usado para la predicción es un modelo de regresión lineal simple, se predice la concentración del analito de interés en la muestra ( yu), a partir del score de la muestra (an( I +1)): yu = an( I+1) / k PARAFAC
(9.3)
Como puede verse, se obtiene un modelo muy simple, y calibrando con un número elevado de muestras, no es necesario realizar una prueba de validación cruzada para evaluar la capacidad de predicción, tal como debe hacerse cuando se modelan datos de primer orden. Finalmente debemos remarcar que para la aplicación de todos los métodos que permiten modelar datos de segundo orden, existe una serie de programas los que pueden ser obtenidos gratuitamente de Internet. 5-7
Referencias 1. K.S. Booksh, B.R. Kowalski, The Theory of Analytical Chemistry, Anal. Chem. 1994, 66 , 782A-791A. 2. R. Bro, Multiway Analysis in the Food Industry. Models, Algorithms, and Applications, Royal Veterinary and Agricultural University Denmark, 1998, disponible en internet en www.models.kvl.dk 3. A. Smilde, R. Bro, P. Geladi, Multiway Analysis. Aplications in the Chemical Sciences, Ed. Wiley, Chichester, 2004. 4. A. Muñoz de la Peña, A. Espinosa Mansilla, D. Gonzalez Gomez, A. Olivieri, H. C. Goicoechea, Interference-free analysis using three-way fluorescence data and the parallel factor model. Determination of fluoroquinolone antibiotics in human serum, Anal Chem. 2003, 75, 2640-2646. 5. www.models.kvl.dk/source/ 6. http://www.ub.es/gesq/mcr/mcr.htm 7. El programa MVC2 ( Multivariate Calibration 2) puede obtenerse libremente en www.chemometry.com
118
ANEXO A NEXO I
REPASO REPASO DE ÁLGEBRA ÁL GEBRA MATRICIAL MATRICIAL Matriz Una matriz es un arreglo rectangular de números, generalmente expresados entre corchetes o paréntesis. Una matriz está formada por arreglos en columnas, llamados vectores columna y arreglos en filas, llamados ll amados vectores fila. La matriz se representa con una letra mayúscula en negrita. Por ejemplo, la matriz X con dimensiones ( I × × J ) tendrá la siguiente estructura: Matriz X Fila
X
I : J :
número de filas número de columnas
Columna En este caso, se dice que la matriz X es rectangular, ya que el número de filas es distinto al de columnas. En el caso en que I = = J , se tiene una matriz cuadrada . Los vectores fila y columna se representan como con una letra minúscula en negrita. Por ejemplo, la j-ésima columna de la matriz X se puede representar como el vector columna x j con dimensión ( I × × 1). Por otra parte, podemos identificar a la i-ésima fila de la matriz X como un vector fila yiT con dimensión (1 × J ). ). Observar que para representar un vector fila, se utiliza el signo T como superíndice, indicando que corresponde a un vector columna traspuesto. La trasposición es una operación común en álgebra matricial, y consiste en transformar las filas en columnas y viceversa. También se pueden trasponer las matrices. Los elementos de la matriz se simbolizan como escalares con subíndices indicando la fila y columna a las que pertenecen: xij. Un escalar puede ser considerado un caso especial de matriz de dimensión (1 × 1). Los escalares se representan con una letra minúscula cursiva. Por ejemplo el escalar c = 5. Como ejemplo, dada la siguiente matriz :
⎡1 ⎢− 1 X=⎢ ⎢4 ⎢ ⎣0
2⎤ 3⎥⎥ 2⎥ ⎥ 1⎦
podemos decir que X es una matriz rectangular de 4 ×2 (4 filas y 2 columnas), que x1 es el primer vector columna:
119
⎡1⎤ ⎢− 1⎥ x1 = ⎢ ⎥ ⎢4⎥ ⎢ ⎥ ⎣0⎦ y que y3 es la tercera fila de la matriz X: y3 =
⎡ 4⎤ ⎢ 2⎥ ⎣ ⎦
ó
T
y3
= [4 2]
Por otro lado, la trasposición de X lleva a: X
⎡1 − 1 4 0 ⎤ = ⎢ ⎥ ⎣2 3 2 1⎦
T
Operaciones con matrices: suma, resta y multiplicación Para sumar o restar matrices, estas deben poseer igual dimensión. La matriz obtenida de la suma o resta de dos matrices es una nueva matriz con igual dimensión a las anteriores y con sus elementos resultando de la suma o resta de los elementos de las matrices originales. Ejemplo:
⎡ 1 5⎤ Dadas: A = ⎢ ⎥ y B = 2 3 ⎣ ⎦
⎡2 9 ⎤ ⎢0 2 ⎥ , ⎣ ⎦
A + B = C =
⎡ (1 + 2) (5 + 9) ⎤ ⎢(2 + 0) (3 + 2)⎥ ⎣ ⎦
El producto de dos matrices sólo existe si el número de columnas de la primera matriz es igual que el de filas de la segunda. El producto de matrices es distributivo y asociativo, pero no es, en general, conmutativo: (A + B) × C = A × C + B × C (A × B) × C = A × (B × C) A × B ≠ B × A (en general)
⎡1 2 5⎤ Ejemplo: dadas: A = ⎢ ⎥ y 3 4 1 ⎣ ⎦
A × B = C =
⎡3 4 ⎤ ⎢ ⎥ B = 5 6 ⎢ ⎥ ⎢⎣1 0⎥⎦
⎡(1 × 3 + 2 × 5 + 5 × 1) (1 × 4 + 2 × 6 + 5 × 0)⎤ ⎢(3 × 3 + 4 × 5 + 1 × 1) (3 × 4 + 4 × 6 + 1 × 0)⎥ ⎣ ⎦
Es decir, se obtiene una matriz con igual número de filas que la primera e igual número de columnas que la segunda. Cada elemento C ijij de la matriz producto C es igual al producto escalar del vector fila i de la matriz A por el vector columna j de la matriz B. Resulta interesante ver las distintas posibilidades de productos de vectores y matrices, y los resultados obtenidos. El siguiente esquema muestra el producto de un vector fila por un vector columna, originando un escalar: 120
T
x y = z
= z x
T
z =
∑ xi yi
y
Nótese que xT debe tener la misma cantidad de columnas que filas y. Si el producto es de un vector columna por un vector fila, el resultado es una matriz.: = Z
xy
T
= Z
T
y
x
Matri Matri ces especiales especiales Una matriz nula es una matriz cuyos elementos son ceros. La matriz identidad representa la unidad en en el entorno matricial. Es una matriz cuadrada en cuya diagonal hay unos, y ceros en el resto ( xij = 1 si i = j, y xij = 0 si i ≠ j:):
⎡1 0 0⎤ ⎢ ⎥ I = 0 1 0 ⎢ ⎥ ⎢⎣0 0 1⎥⎦ Se dice que una matriz es singular cuando cuando por lo menos hay una fila o columna que resulte de combinaciones lineales entre filas o columnas de la matriz. Una matriz es simétrica cuando xij = x ji para todos los i y j:
⎡1 5 2 ⎤ ⎢ ⎥ A = 5 8 0 ⎢ ⎥ ⎢⎣2 0 5⎥⎦ Matriz triangular es es cuando xij = 0 para 0 para todo i > j o para todo i < j:
⎡1 0 0 ⎤ ⎢ ⎥ A = 5 8 0 ⎢ ⎥ ⎢⎣2 0 5⎥⎦ La traza de una matriz cuadrada es la suma de los elementos diagonales de una matriz:
121
Si
⎡1 0 0 ⎤ ⎢ ⎥ A = 5 8 0 ⎢ ⎥ ⎢⎣2 0 5⎥⎦
tr(A) = 1 + 8 + 5 = 14
Norma o longit ud de un vector
⎡ x ⎤ Dado el vector x = ⎢ 11 ⎥ , este se puede representar en un eje cartesiano de la ⎣ x12 ⎦ siguiente manera: x2 x x12
x11
x1
El cuadrado de la distancia desde el extremo del vector al origen se calcula como: ║x║2 = x112 + x122, que equivale al producto xTx. Como puede apreciarse, el resultado es un escalar. La raíz cuadrada, ║x║, se llama longitud o norma del vector x, y se puede generalizar a vectores de n dimensiones:
║x║2 = x112 + x122 + ... + x1n2 Normalización d e un vector Un vector normalizado es igual al vector dividido por su norma. Ya que la norma es un escalar, esto corresponde a dividir cada elemento por la norma. Dado el vector x, el vector u normalizado es: u = x /║x║
Una propiedad interesante de los vectores normalizados es que tienen longitud igual a uno. Además, todos los elementos de dicho vector tienen valor absoluto menor que 1. Determinante Una matriz cuadrada X puede ser caracterizada a través de un escalar llamado determinante, que se simboliza como det( X). Para el caso de una matriz X de (2 × 2) se tiene: det(X) = x11 x22 – x12 x21 Para matrices de órdenes superiores existen complejos procedimientos que se facilitan usando rutinas de computación fácilmente disponibles. 122
Inversa La inversa de una matriz cuadrada X se simboliza por X –1, y es una matriz tal que X X –1 = X –1 X = I. Si X es singular, X –1 no existe, ya que el cálculo de la inversa implica un paso en el que se debe dividir por el determinante y este resulta igual a cero para matrices singulares. Invers a Generalizada La inversa generalizada es el análogo de la matriz inversa en el mundo de las matrices rectangulares. Dada la matriz rectangular X de dimensión ( m × n), se puede obtener una matriz cuadrada pre- o post multiplicándola por su traspuesta. En el primer caso se obtiene la matriz XTX (n × n), en el segundo XXT (m × m). En el caso en que m > n, el producto por la trsnspuesta conduce a la compresión y expansión respectivamente de la matriz X. En ambos casos, se tendrán las siguientes expresiones, útiles para despejar una matriz rectangular X de una ecuación matricial: (XTX) –1 XT X = I T T –1 X X (XX ) = I Como puede verse, dependiendo de la necesidad del caso, se debe invertir la matriz cuadrada obtenida por compresión o la obtenida por expansión. En el primer caso, la operación es posible siempre y cuando la matriz obtenida no presente colinealidades. En el segundo caso la inversión no es posible porque el determinante de (XXT) es igual a cero. Esto último resulta del hecho que la expansión genera filas y columnas que son combinaciones lineales de las filas y columnas originales de la matriz rectangular. Seudoinversa Si bien el concepto de seudoinversa se ha desarrollado para matrices singulares, hoy en día está generalizado a todo tipo de matrices, y es importante comprender cómo funciona en diferentes casos. En general, dada una matriz X, se entiende por X+ a una matriz que cumple con la condición: X X+ X = X
Existen varias posibilidades. 1) Si la matriz X es cuadrada y no singular, entonces la seudoinversa es la propia inversa, esto es, X+ = X –1. Dado que X X –1 = X –1 X = I, es evidente que se cumple la ecuación de la seudoinversa: X X –1 X = X
123
2) Si la matriz X es rectangular, entonces una de las matrices cuadradas (ya sea T X X ó XX ) es necesariamente singular (aquella cuyo tamaño es mayor). Sólo si la matriz de menor tamaño es no singular, la seudoinversa es igual a la inversa generalizada. Supongamos que la matriz cuadrada ( XT X) es la de menor tamaño y es no singular, entonces X+ = (XT X) –1 XT, ya que X+ X = I, y por lo tanto: T
X [(XT X) –1 XT] X = X
Si por el contrario, la matriz cuadrada ( X XT) es la de menor tamaño y es no singular X+ = XT (X XT) –1, ya que X X+ = I, y por lo tanto: X [XT (X XT) –1] X = X
3) Si la matriz X es cuadrada y singular, o es rectangular y las matrices cuadradas XTX y XXT son ambas singulares, entonces X+ es, en general, un tipo de matriz seudoinversa que debe obtenerse mediante un procedimiento matemático especial. Entre ellos, se destaca la descomposición en valores singulares (que conduce a la llamada seudoinversa de Moore-Penrose). En este caso no se calcula directamente la inversa de ninguna matriz, sino que se busca una matriz X+ que cumpla con la condición: +
X X X = X
Existen muchas seudoinversas capaces de cumplir la ecuación anterior, de las cuales la de Moore-Penrose es quizás la más popular. A modo de ejemplo, sea la matriz:
⎡1 ⎢1 X = ⎢ ⎢1 ⎢ ⎣1
1⎤ 1⎥⎥ 1⎥ ⎥ 1⎦
La seudoinversa de Moore-Penrose es: X
⎡0,125 0,125 0,125 0,125⎤ = ⎢ ⎥ ⎣0,125 0,125 0,125 0,125⎦
+
Verifique que se cumple la ecuación X X+ X = X. Sin embargo, no es la única posibilidad, por ejemplo, puede verificarse que la siguiente matriz: X
⎡1 0 0 0 ⎤ = ⎢ ⎥ ⎣0 0 0 0 ⎦
+
es también una posible seudoinversa. De hecho, cualquier matriz cuyos elementos genéricos cumplan con la siguiente condición: X
⎡a b c d ⎤ = ⎢ ⎥ , tal que a + b + c + d + e + f +g + h = 1 e f g h ⎦ ⎣
+
124
será también seudoinversa de X. La seudoinversa de Moore-Penrose se aplica fundamentalmente a sistemas de ecuaciones inconsistentes o que tienen varias soluciones posibles. En este sentido, el uso de la seudoinversa para resolver este tipo de ecuaciones lleva a soluciones llamadas mínimas. Por ejemplo, en la resolución de la ecuación Ax = b, donde A es una matriz y tanto x como b son vectores, si el sistema tiene un conjunto de soluciones posibles, la solución mínima es aquella para la cual la norma || Ax – b || es mínima. Para más detalles, véase la obra Análisis Numérico, de D. Kincaid y W. Cheney.1 Autovectores y autovalor es Dada una matriz cuadrada M (n × n), un problema de autovectores y autovalores implica resolver la siguiente ecuación: M u = λ u
donde λ es un autovalor y u es su autovector correspondiente. La ecuación anterior puede resolverse si introducimos la matriz identidad I del mismo tamaño que M: I M u = λ I u I M u – λ I u = 0 (M – λ I) u = 0
Esta última ecuación se cumple si det( M – λI) = 0.La resolución de la última ecuación genera un polinomio de grado n, cuyas n raíces corresponden a los n autovalores, que a su vez tienen n autovectores asociados. Es importante tener en cuenta las siguientes propiedades: 1) La suma de los n autovalores es igual a la suma de los elementos de la diagonal principal de la matriz o traza. 2) Cada autovalor tiene un autovector asociado. El valor absoluto del autovalor indica el grado de importancia del autovector, en términos de su contribución a la variación espectral contenida en la matriz M. 4) Los autovectores son ortogonales entre si, es decir que su producto escalar es cero: T ui u j = 0 (si i ≠ j) 5) Los autovectores generalmente se normalizan, dividiendo cada uno por su norma: T ui ui = 1 Supongamos que la matriz M está dada por:
⎡1 1⎤ ⎢1 2⎥ ⎣ ⎦ La ecuación a resolver para encontrar los autovalores y autovectores es: M =
det(M – λI) = 0 que puede ser escrita en forma detallada como sigue:
125
⎡1 − λ 1 ⎤ det( ⎢ ⎥ ) = (1 – λ)(2 – λ) – 1 = 0 1 2 − λ ⎣ ⎦ Esta ecuación cuadrática tiene dos soluciones, los autovalores: λ1 = 2,618
y
λ2 = 0,382
Para encontrar los autovectores, recurrimos a la definición del problema: M u = λ u Para el primer autovalor, por ejemplo, encontramos:
⎡1 1⎤ ⎢1 2⎥ u = 2,618 u ⎣ ⎦ Si el autovector u contiene elementos u1 y u2:
⎡1 1⎤ ⎡ u1 ⎤ ⎡ u1 + u 2 ⎤ ⎡ u1 ⎤ = = 2,618 ⎢1 2⎥ ⎢u ⎥ ⎢u + 2 u ⎥ ⎢u ⎥ ⎣ ⎦ ⎣ 2⎦ ⎣ 1 2⎦ ⎣ 2⎦ Por lo tanto, deben cumplirse simultáneamente las ecuaciones siguientes: u1 + u2 = 2,618 u1 u1 + 2 u2 = 2,618 u2
Este sistema de dos ecuaciones con dos incógnitas está indeterminado. Cualquier par de valores u1 y u2 que cumpla con la condición u2 / u1 = 1,618 es solución. En definitiva, el autovector puede escribirse como:
⎡ 1 ⎤ ⎥ ⎣1,618⎦
u = t ⎢
donde t es una constante. Si normalizamos este vector, la indeterminación desaparece, ya que: 1 ⎡ ⎤ ⎢ (t 2 + 1,618 2 t 2 )1 / 2 ⎥ u = t ⎢ ⎥ = 1 , 618 ⎢ ⎥ ⎢⎣ (t 2 + 1,618 2 t 2 )1 / 2 ⎥⎦
1 ⎡ ⎤ ⎢ (1 + 1,618 2 )1 / 2 ⎥ ⎢ ⎥ 1 , 618 ⎢ ⎥ ⎢⎣ (1 + 1,618 2 )1 / 2 ⎥⎦
Y finalmente: u =
⎡0,5257⎤ ⎢0,8507 ⎥ ⎣ ⎦
que es el primer autovector normalizado de M, asociado al autovalor λ1. De idéntica manera se puede encontrar el segundo autovector:
126
⎡− 0,8507 ⎤ ⎢ 0,5257 ⎥ ⎣ ⎦ Puede comprobarse que, además de estar normalizados, los autovectores son ortogonales: u' =
T
u
T
⎡0,5257⎤ ⎡0,5257⎤ ⎡0,5257 ⎤ × u = ⎢ [ ] = 0 , 5257 0 , 8507 ⎥ ⎢0,8507 ⎥ ⎢0,8507⎥ = 1 0 , 8507 ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ T
⎡− 0,8507 ⎤ ⎡− 0,8507⎤ ⎡− 0,850⎤ u'T × u' = ⎢ − [ ] = 0 , 8507 0 , 5257 ⎥ ⎢ 0,5257 ⎥ ⎢ 0,5257 ⎥ = 1 0 , 5257 ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ T
u
T
⎡− 0,8507 ⎤ ⎡0,5257 ⎤ ⎡− 0,8507 ⎤ × u' = ⎢ [ ] = 0 , 5257 0 , 8507 ⎥ ⎢ 0,5257 ⎥ ⎢ 0,5257 ⎥ = 0 0 , 8507 ⎣ ⎦ ⎣ ⎦ ⎣ ⎦
Descomposi ción en valores s ingulares Cualquier matriz Z (de n×m) puede descomponerse de acuerdo a la siguiente ecuación: Z n ×m
=
U n×m
×
W m×m
×
T
V m×m
En la descomposición en valores singulares, las matrices U, W y V cumplen con las siguientes condiciones. W es una matriz cuadrada diagonal. Los elementos de la diagonal principal son distintos de cero y se llaman valores singulares; usualmente se ordenan de mayor a menor. Los demás elementos son todos iguales a cero. Los valores singulares elevados al cuadrado son los autovalores de la matriz Z, de manera que W puede representarse por:
⎡(λ1 )1 / 2 0 0 ⎢ (λ 2 )1 / 2 0 ⎢ 0 W = ⎢ 0 (λ 3 )1 / 2 0 ⎢ 0 0 ⎢ 0 ⎢ 0 0 0 ⎣
0 0 ⎤ ⎥ 0 0 ⎥ 0 0 ⎥ ⎥ ... 0 ⎥ 0 (λ m )1 / 2 ⎥⎦
donde (λ1)1/2 ≥ (λ2)1/2 ≥ ... ≥ (λm)1/2. De este modo, la matriz diagonal de autovalores matriz W:
⎡λ 1 ⎢ ⎢0 =⎢ 0 ⎢ ⎢0 ⎢0 ⎣
0 0 0 0 λ2 0 0 0 0 λ3 0 0 0 0 ... 0 0 0 0 λm
es igual al cuadrado de la
⎤ ⎥ ⎥ ⎥ = W2 ⎥ ⎥ ⎥ ⎦ 127
Desde el punto de vista quimiométrico, los valores singulares (y los autovalores) explican de algún modo las variaciones contenidas en la matriz Z. Dado que los valores singulares están ordenados en forma decreciente, cada uno explica sucesivamente menos información. La suma de los autovalores es igual a la suma de la varianza total, y por lo tanto, una medida de la contribución de cada autovalor a la varianza está dada por: Varianza explicada por el autovalor λi (%) = 100
λi m
∑=1 λ
i
i
Las matrices U y V contienen vectores columna llamados vectores singulares izquierdos y derechos respectivamente. Cada columna está asociada al autovalor respectivo. Tanto U como V son matrices ortonormales, en el sentido que están constituidas por columnas que son ortogonales entre sí, y además están normalizadas: T
U U = I T V V = I
Las columnas de U son los autovectores de la matriz cuadrada ( Z ZT), ya que: (Z ZT) = (U W VT) (U W VT)T = U W VT V W UT = U W W UT = U W2 UT = = U UT de donde, pos-multiplicando ambos miembros de la ecuación anterior por U, y usando la propiedad de que UT U = I. (Z ZT) U = U que tiene la misma forma que un problema de autovalores y autovectores para la matriz (Z ZT). Del mismo modo, se puede demostrar que las columnas de V son los autovectores de la matriz cuadrada ( ZT Z). Supongamos como ejemplo que la matriz Z es la siguiente:
⎡1 2 ⎤ ⎢ ⎥ Z = 4 5 ⎢ ⎥ ⎢⎣5 4⎥⎦ La descomposición en valores singulares provee las siguientes matrices:
⎡− 0,2307 − 0,5513 ⎤ T 0 ⎤ ⎡− 0,6944 0,7196 ⎤ ⎡9,2481 ⎢ ⎥ Z = − 0,6894 − 0,4889 × ⎢ × ⎢ ⎥ ⎣ 0 1,2138⎥⎦ ⎢⎣− 0,7196 − 0,6944 ⎥⎦ ⎢⎣− 0,6867 0,6761 ⎥⎦ Puede comprobarse que estas matrices cumplen con las condiciones anteriormente comentadas. Desde el punto de vista quimiométrico, interesa comentar las siguientes propiedades. 128
Los autovalores están en relación (9,2481) 2/(1,2138)2 = 58,05. Esto significa que el primer autovector contribuye a la varianza de Z alrededor de 58 veces más que el segundo. De la varianza total, entonces, el primer autovector es capaz de explicar el siguiente porcentaje:
(9,2481)2 = 98,3 % Explicado por el primer autovector = 100 × (9,2481)2 + (1,2138)2 La segunda conclusión importante es que, dado que este autovector es capaz de explicar un porcentaje muy significativo de la variación en Z, esta matriz puede aproximarse utilizando sólo este componente principal en la descomposición matricial:
⎡− 0,2307 ˆ = ⎢ − 0,6894 Z ⎢ ⎢⎣− 0,6867
⎤ ⎥ × 9,2481 × ⎥ ⎥⎦
⎡− 0,6944 ⎢− 0,7196 ⎣
⎡1,48 1,54 ⎤ ⎤ ⎢ ⎥ ⎥ = ⎢4,43 4,59⎥ ⎦ ⎢⎣ 4,41 4,57⎥⎦ T
donde Zˆ representa la matriz reconstruida con un número limitado de componentes, en este caso uno. Puede notarse que la reconstrucción no es completa, pero se aproxima a la matriz original. En el caso de matrices de mayor tamaño, la aproximación es mejor a medida que se agregan componentes principales.
Referencias 1. D. Kincaid. W. Cheney, Análisis Numérico, Addison-Wesley, 1994, capítulo 5. 2. D. L. Massart, B. G. M. Vandeginste, L. M. C. Buydens, S. De Jong, P. J. Lewi y J. Smeyers-Verbeke, Handbook of Chemometrics and Qualimetrics, Elsevier, Amsterdam, 1997, Capítulo 29.
129
ANEXO II
PROGRAMAS DE COMPUTACIÓN En el CD adjunto existe una serie de programas de computación que permiten la aplicación de todas las técnicas desarrolladas en el presente libro. Programas ejecut ables escr ito s en Quick Basic (QB) Los programas ejecutables escritos en QB que se proveen son muy sencillos. Para operar con ellos deben seguirse los siguientes pasos: 1) Copiar el programa ejecutable en una carpeta del disco rígido. 2) Copiar en esa carpeta los datos organizados apropiadamente. En el caso del programa LR_CAL.EXE, los datos deben estar contenidos en un archivo de texto con las siguientes columnas: Columna 1: concentración de cada patrón Columnas 2 ... : respuesta de las réplicas de cada patrón En el caso del programa LR_PRED.EXE, los datos deben estar contenidos en un archivo de texto con las siguientes columnas: Columnas 1 ... : respuesta de las réplicas de cada muestra incógnita. 3) Ejecutar el programa haciendo doble click en su nombre desde el explorador de Windows. 4) Suministrar al programa los datos que pide. En el caso del programa LR_CAL.EXE, pide el número de niveles de concentración, número de réplicas y nombre del archivo de texto con los datos. En el caso del programa LR_PRED.EXE, pide el número de muestras incógnita, el número de réplicas y nombre del archivo de texto con los datos. Notas importantes:
1) Los programas ejecutables escritos en QB que se proveen son muy sencillos y no tienen soporte gráfico. 2) Estos programas no admiten nombres largos de archivo. Sólo se permiten nombres con no más de ocho caracteres. 3) Al organizar los datos en un archivo de texto, no separar los números usando el tabulador, sino mediante espacios en blanco. Instrucciones específicas para el Capítulo 2:
1) Ubicar el archivo comprimido "QB CLASE 1.ZIP" en un directorio del disco rígido. 2) Descomprimir el archivo anterior. 3) Editar los archivos de texto (doble click sobre el ícono del archivo) y examinar su contenido. Los archivos "D_E_R_C.TXT" y "D_E_R_L.TXT" contienen datos de calibración organizados del modo que sigue: Columna 1: concentración de cada patrón Columnas 2 ... : respuesta de las réplicas de cada patrón 130
El archivo "D_E_R_T.TXT" contiene datos para la predicción, organizados así: Columnas 1 ... : respuesta de las réplicas de cada muestra incógnita. 4) Tomar nota del número de filas y columnas de cada archivo. 5) Para calibrar un sistema, ejecutar el programa LR_CAL.EXE, haciendo doble click en su nombre desde el explorador de Windows, e ingresar la información solicitada. Este programa genera un archivo llamado "CALIBR_.TXT" en el que se guardan datos de la calibración para futuras predicciones. También genera un archivo "RESIDUOS.TXT", donde se guardan los residuos de la regresión para estudios de linealidad. Editarlo para examinar su contenido, y graficar los residuos en función de la concentración de cada muestra. 6) Luego de efectuada la calibración, predecir muestras incógnitas ejecutando el programa LR_PRED.EXE, e ingresando la información solicitada. Este programa genera un archivo de texto "RESULT_.TXT" donde se guardan las predicciones de cada muestra, incluyendo intervalo de confianza, desvío estándar y desvío estándar relativo. Instrucciones específicas para el Capítulo 3:
1) Ubicar el archivo comprimido "QB CLASE 2.ZIP" en un directorio del disco rígido. 2) Descomprimir el archivo anterior. 3) Editar los archivos de texto (doble click sobre el ícono del archivo) y examinar su contenido. El archivo "D_E_WLS.TXT" contiene datos para estudios de exactitud, organizados del modo que sigue: Columna 1: concentración de cada patrón. Columna 2: concentración predicha por el método. Columna 3: desvío estándar en la concentración predicha. El archivo "D_C_BLS_T.TXT" contiene datos para la comparación de métodos, organizados así: Columna 1: concentración predicha por el método 1. Columna 2: desvío estándar en la concentración del método 1. Columna 3: concentración predicha por el método 2. Columna 4: desvío estándar en la concentración del método 2. NOTA: si no existen datos de desvíos estándar, el archivo de texto contendrá solamente dos columnas, la primera con las concentración de cada patrón, y la segunda con las concentraciones predichas por el método. 4) Tomar nota del número de filas y columnas de cada archivo. 5) Para estudios de exactitud o comparación, ejecutar el programa EJCR.EXE, haciendo doble click en su nombre desde el explorador de Windows, e ingresar la información solicitada. Este programa genera un archivo llamado "EJCR_.TXT" en el que se guardan datos para graficar la elipse de confianza conjunta. 6) Editar el archivo "EJCR_.TXT" para confirmar su estructura. 7) Para graficar la elipse, se debe utilizar un programa gráfico. Deben superponerse dos gráficas, una que tenga como eje x a la primera columna del archivo "EJCR_.TXT" y como eje y la segunda columna, y otra que tenga como eje x la primera columna y como eje y la tercera columna. 8) En EXCEL, por ejemplo, se pueden seguir estos pasos: 131
Ir a "Archivo", "Abrir". Seleccionar el archivo "EJCR_.TXT" y abrirlo. Elegir la opción "tipo de archivo" como "delimitado". Clickear "Siguiente". Seleccionar únicamente la opción "Coma", deseleccionando otras opciones. Clickear "Siguiente". Clickear "Finalizar". Las tres columnas deben quedar en la pantalla. Seleccionar las tres columnas. Ir a "Insertar", "Gráfico". Seleccionar "XY Dispersión" Seleccionar la opción "Dispersión con líneas suavizadas y sin marcadores de datos". Clickear en "Siguiente". Clickear en "Siguiente". Agregar título y etiquetas en los ejes. Clickear en "Finalizar". Instrucciones específicas para el Capítulo 5:
1) Ubicar el archivo comprimido "QB CLASE 4.ZIP" en un directorio del disco rígido. 2) Descomprimir el archivo anterior. 3) Editar los archivos de texto (doble click sobre el ícono del archivo) y examinar su contenido. El archivo "XCAL_E_R.TXT" contiene las concentraciones de dos analitos para calibrar un modelo CLS para 4 muestras, en una matriz de 4 ×2. El archivo "YCAL_E_R.TXT" contiene las respuestas instrumentales de las 4 muestras de calibración, organizados en una matriz de 6 filas y 4 columnas. El archivo "XVAL_E_R.TXT" contiene las concentraciones para 4 muestras, en una matriz de 6 ×4. El archivo "YVAL_E_R.TXT" contiene las respuestas instrumentales de las 4 muestras de validación, organizados en una matriz de 6 filas y 4 columnas. El archivo "YPRU_E_R.TXT" contiene las respuestas instrumentales de 3 muestras de prueba, organizado en una matriz de 6 ×3. También se proveen archivos similares para otro sistema diferente, con los nombres "CONC_CAL.TXT" (matriz de 9 ×2), "RESP_CAL.TXT" (matriz de 281 ×9) y "RESP_TST_T.TXT" (matriz de 281×3). 4) Para calibrar el modelo CLS, ejecutar el programa CLS_CAL.EXE, e ingresar la información solicitada. Este programa genera los archivos "S_.TXT" en el que se guardan las sensibilidades, en forma de matriz de J × N (longitudes de onda × número de analitos), y "B_.TXT" en el que se guardan los coeficientes de regresión, en forma de matriz de J× N (longitudes de onda × número de analitos). 5) Las sensibilidades y coeficientes de regresión se pueden graficar utilizando un programa gráfico, por ejemplo EXCEL, siguiendo pasos análogos a los descriptos en el capítulo 2. 6) Para predecir muestras incógnita, luego de haber calibrado el modelo, ejecutar el programa CLS_PRED.EXE, e ingresar la información solicitada. Instrucciones específicas para el Capítulo 6:
1) Ubicar el archivo comprimido "QB CLASE 5.ZIP" en un directorio del disco rígido. 2) Descomprimir el archivo anterior. 3) Editar los archivos de texto (doble click sobre el ícono del archivo) y examinar su contenido. 132
El archivo "BR_CON_C.TXT" contiene las concentraciones para calibrar el modelo PCR para 12 muestras conteniendo bromhexina, en una única columna. El archivo "BR_RES_C.TXT" contiene las respuestas instrumentales de las 12 muestras de calibración, organizados en una matriz de 64 filas y 12 columnas. El archivo "BR_CON_T.TXT" contiene las concentraciones para validar el modelo PCR para 11 muestras, en una única columna. El archivo "BR_RES_T.TXT" contiene las respuestas instrumentales de las 11 muestras de validación, organizados en una matriz de 64 filas y 11 columnas. El archivo "BR_RES_P.TXT" contiene las respuestas instrumentales de una muestra de prueba, organizado en un vector de 64 filas. También se proveen archivos similares para otro sistema diferente, conteniendo tetraciclina, con los nombres "TE_CON_C.TXT" (50 valores), "TE_RES_C.TXT" (matriz de 101×50), "TE_CON_T.TXT" (57 valores) y "TE_RES_T.TXT" (matriz de 101×57). 4) Tomar nota del número de filas y columnas de cada archivo. 5) Para realizar la validación cruzada, ejecutar el programa PCR_CV.EXE, haciendo doble click en su nombre desde el explorador de Windows, e ingresar la información solicitada. Este programa genera un archivo llamado "PRESS.TXT" en el que se guarda el informe estadístico de la validación cruzada. 6) Editar el archivo "PRESS.TXT" para confirmar su estructura. 7) Para calibrar el modelo PCR, ejecutar el programa PCR_CAL.EXE, e ingresar la información solicitada. Este programa genera archivos llamados "U_1.TXT", "U_2.TXT", etc., en el que se guardan los autovectores o factores. 8) Para graficar los factores se debe utilizar un programa gráfico, por ejemplo EXCEL, siguiendo pasos análogos a los descriptos en el capítulo 2. 9) Para predecir muestras incógnita, luego de haber calibrado el modelo, ejecutar el programa PCR_PRED.EXE, e ingresar la información solicitada. Se genera el archivo "RESULT_.TXT" conteniendo los resultados de la predicción. Instrucciones específicas para el Capítulo 7:
1) Ubicar el archivo comprimido "QB CLASE 6.ZIP" en un directorio del disco rígido. 2) Descomprimir el archivo anterior. 3) Editar los archivos de texto (doble click sobre el ícono del archivo) y examinar su contenido. El archivo "BR_CON_C.TXT" contiene las concentraciones para calibrar el modelo PCR para 12 muestras conteniendo bromhexina, en una única columna. El archivo "BR_RES_C.TXT" contiene las respuestas instrumentales de las 12 muestras de calibración, organizados en una matriz de 64 filas y 12 columnas. El archivo "BR_CON_T.TXT" contiene las concentraciones para validar el modelo PCR para 11 muestras, en una única columna. El archivo "BR_RES_T.TXT" contiene las respuestas instrumentales de las 11 muestras de validación, organizados en una matriz de 64 filas y 11 columnas. El archivo "BR_RES_P.TXT" contiene las respuestas instrumentales de una muestra de prueba, organizado en un vector de 64 filas. También se proveen archivos similares para otro sistema diferente, conteniendo tetraciclina, con los nombres "TE_CON_C.TXT" (50 valores), "TE_RES_C.TXT" (matriz de 101×50), "TE_CON_T.TXT" (57 valores) y "TE_RES_T.TXT" (matriz de 101×57). 4) Tomar nota del número de filas y columnas de cada archivo. 5) Para realizar la validación cruzada, ejecutar el programa PLS_CV.EXE, haciendo doble click en su nombre desde el explorador de Windows, e ingresar la información 133
solicitada. Este programa genera un archivo llamado "PRESS.TXT" en el que se guarda el informe estadístico de la validación cruzada. 6) Editar el archivo "PRESS.TXT" para confirmar su estructura. 7) Para calibrar el modelo PLS, ejecutar el programa PLS_CAL.EXE, e ingresar la información solicitada. Este programa genera los archivos "W_1.TXT", ..., conteniendo los factores espectrales de calibración. También se generan archivos auxiliares (TEMP.TXT, R1.TXT, ...) para uso posterior en predicciones de incógnitas. 8) Para predecir muestras incógnita, luego de haber calibrado el modelo, ejecutar el programa PLS_PRED.EXE, e ingresar la información solicitada. Se genera el archivo "RESULT_.TXT" conteniendo los resultados de la predicción. Instrucciones específicas para el Capítulo 8:
1) Ubicar el archivo comprimido "QB CLASE EXTRA.ZIP" en un directorio del disco rígido. 2) Descomprimir el archivo anterior. 3) Editar los archivos de texto (doble click sobre el ícono del archivo) y examinar su contenido. El archivo "BR_CON_C.TXT" contiene las concentraciones para calibrar el modelo PCR para 12 muestras conteniendo bromhexina, en una única columna. El archivo "BR_RES_C.TXT" contiene las respuestas instrumentales de las 12 muestras de calibración, organizados en una matriz de 64 filas y 12 columnas. El archivo "BR_CON_T.TXT" contiene las concentraciones para validar el modelo PCR para 11 muestras, en una única columna. El archivo "BR_RES_T.TXT" contiene las respuestas instrumentales de las 11 muestras de validación, organizados en una matriz de 64 filas y 11 columnas. El archivo "BR_RES_P.TXT" contiene las respuestas instrumentales de una muestra de prueba, organizado en un vector de 64 filas. También se proveen archivos similares para otro sistema diferente, conteniendo tetraciclina, con los nombres "TE_CON_C.TXT" (50 valores), "TE_RES_C.TXT" (matriz de 101×50), "TE_CON_T.TXT" (57 valores) y "TE_RES_T.TXT" (matriz de 101×57). 4) Tomar nota del número de filas y columnas de cada archivo. 5) Para realizar la validación cruzada, ejecutar el programa PLS_CV.EXE, haciendo doble click en su nombre desde el explorador de Windows, e ingresar la información solicitada. Este programa genera un archivo llamado "PRESS.TXT" en el que se guarda el informe estadístico de la validación cruzada. 6) Editar el archivo "PRESS.TXT" para confirmar su estructura. 7) Para calibrar el modelo PLS, ejecutar el programa PLS_CAL.EXE, e ingresar la información solicitada. Este programa genera los archivos "W_1.TXT", ..., conteniendo los factores espectrales de calibración. También se generan archivos auxiliares (TEMP.TXT, R1.TXT, ...) para uso posterior en predicciones de incógnitas. 8) Para predecir muestras incógnita, luego de haber calibrado el modelo, ejecutar el programa PLS_PRED.EXE, e ingresar la información solicitada. Se genera el archivo "RESULT_.TXT" conteniendo los resultados de la predicción. Programas ejecut ables escr ito s en MATLAB Cualquiera sea la versión de MATLAB que se utilice, guarde los archivos con los datos a procesar y las rutinas de trabajo en un mismo directorio (los usuarios avanzados también pueden guardar las rutinas en un directorio asignado como "permanente" para MATLAB y los datos en directorios temporarios). 134
Organice los datos en archivos de texto, que pueden generarse mediante el "Block de notas", o exportando datos desde planillas gráficas (tipo EXCEL) en formato texto. Utilice como modelo los archivos *.TXT que acompañan el capítulo. Pueden usarse nombres largos de archivo, pero hay ciertas restricciones para los nombres de archivo válidos en MATLAB: No empiece el nombre con números, por ejemplo "2_DAT.TXT" no es válido. No utilice signos matemáticos en el nombre, por ejemplo "DAT-1.TXT" o "DAT*1.TXT" no son válidos. El guión bajo puede usarse sin problemas, por ejemplo "DAT_1.TXT" es válido. Sólo incluya números en los archivos; no agregue texto de ninguna naturaleza. Instrucciones específicas para el Capítulo 2:
1) Ubicar el archivo comprimido "MATLAB CLASE 1.ZIP" en un directorio del disco rígido. 2) Descomprimir el archivo anterior. 3) Editar los archivos de texto (doble click sobre el ícono del archivo) y examinar su contenido. Los archivos "DATOS_EJ_RES_COMPLETOS.TXT" y "DATOS_EJ_RES_LINEAL.TXT" contienen datos de calibración organizados del modo que sigue: Columna 1: concentración de cada patrón Columnas 2 ... : respuesta de las réplicas de cada patrón El archivo "DATOS_EJ_RES_TEST.TXT" contiene datos para la predicción, organizados así: Columnas 1 ... : respuesta de las réplicas de cada muestra incógnita. 4) Para calibrar un sistema, ejecutar rutina LR_CAL.M, de la forma que se explica más adelante en este documento. Ingresar el nombre del archivo de datos entre comillas simples. Este programa genera un archivo llamado "CALIBRACION", de uso interno de MATLAB, en el que se guardan datos de la calibración para futuras predicciones. 5) Luego de efectuada la calibración, predecir muestras incógnitas ejecutando la rutina LR_PRED.M, e ingresando el nombre del archivo de datos de predicción entre comillas simples. 6) Junto con las rutinas LR_CAL.M y LR_PRED.M también se encuentran las rutinas auxiliares PROBABI.M y STUD.M, que realizan cálculos de probabilidad y coeficiente de student para las rutinas de calibración y predicción. Instrucciones específicas para el Capítulo 3:
1) Ubicar el archivo comprimido "MATLAB CLASE 2.ZIP" en un directorio del disco rígido. 2) Descomprimir el archivo anterior. 3) Editar los archivos de texto (doble click sobre el ícono del archivo) y examinar su contenido. El archivo "D_EXACT_WLS.TXT" contiene datos para estudios de exactitud, organizados del modo que sigue: Columna 1: concentración de cada patrón. Columna 2: concentración predicha por el método. Columna 3: desvío estándar en la concentración predicha.
El archivo "D_COMPAR_BLS_T.TXT" contiene datos para la comparación de métodos, organizados así: Columna 1: concentración predicha por el método 1. Columna 2: desvío estándar en la concentración del método 1. Columna 3: concentración predicha por el método 2. Columna 4: desvío estándar en la concentración del método 2. 4) Para estudios de exactitud o comparación, ejecutar la rutina EJCR.M desde MATLAB, tal como se explicó en el capítulo 1 para otras rutinas. Ingresar el nombre del archivo entre comillas simples. Este programa genera un archivo llamado ("EJCR_"+Nombre del archivo de datos +".TXT") en el que se guardan datos para graficar la elipse de confianza conjunta, y que puede utilizarse para realizar el gráfico con otros programas. Si no existen datos de desvíos estándar, el archivo de texto contendrá solamente dos columnas, la primera con las concentración de cada patrón, y la segunda con las concentraciones predichas por el método. Instrucciones específicas para el Capítulo 5:
1) Ubicar el archivo comprimido "MATLAB CLASE 4.ZIP" en un directorio del disco rígido. 2) Descomprimir el archivo anterior. 3) Editar los archivos de texto (doble click sobre el ícono del archivo) y examinar su contenido. El archivo "XCAL_E_R.TXT" contiene las concentraciones de dos analitos para calibrar un modelo CLS para 4 muestras, en una matriz de 4 ×2. El archivo "YCAL_E_R.TXT" contiene las respuestas instrumentales de las 4 muestras de calibración, organizados en una matriz de 6 filas y 4 columnas. El archivo "XVAL_E_R.TXT" contiene las concentraciones para 4 muestras, en una matriz de 6 ×4. El archivo "YVAL_E_R.TXT" contiene las respuestas instrumentales de las 4 muestras de validación, organizados en una matriz de 6 filas y 4 columnas. El archivo "YPRU_E_R.TXT" contiene las respuestas instrumentales de 3 muestras de prueba, organizado en una matriz de 6 ×3. También se proveen archivos similares para otro sistema diferente, con los nombres "CONC_CAL.TXT" (matriz de 9 ×2), "RESP_CAL.TXT" (matriz de 281 ×9) y "RESP_TST_T.TXT" (matriz de 281×3). 4) Para calibrar el modelo CLS, ejecutar la rutina CLS_CAL.M, e ingresar la información solicitada. Este programa genera los archivos "S_.TXT" en el que se guardan las sensibilidades, en forma de matriz de J × N (longitudes de onda × número de analitos), y "B_.TXT" en el que se guardan los coeficientes de regresión, en forma de matriz de J× N (longitudes de onda × número de analitos). 5) Las sensibilidades y coeficientes de regresión también se pueden graficar utilizando un programa gráfico, por ejemplo EXCEL, siguiendo pasos análogos a los descriptos en el capítulo 2. 6) Para predecir muestras incógnita, luego de haber calibrado el modelo, ejecutar la rutina CLS_PRED.M, e ingresar la información solicitada. Instrucciones específicas para el Capítulo 6:
1) Ubicar el archivo comprimido "MATLAB CLASE 5.ZIP" en un directorio del disco rígido. 2) Descomprimir el archivo anterior. 136
3) Editar los archivos de texto (doble click sobre el ícono del archivo) y examinar su contenido. El archivo "BR_CON_C.TXT" contiene las concentraciones para calibrar el modelo PCR para 12 muestras conteniendo bromhexina, en una única columna. El archivo "BR_RES_C.TXT" contiene las respuestas instrumentales de las 12 muestras de calibración, organizados en una matriz de 64 filas y 12 columnas. El archivo "BR_CON_T.TXT" contiene las concentraciones para validar el modelo PCR para 11 muestras, en una única columna. El archivo "BR_RES_T.TXT" contiene las respuestas instrumentales de las 11 muestras de validación, organizados en una matriz de 64 filas y 11 columnas. El archivo "BR_RES_P.TXT" contiene las respuestas instrumentales de una muestra de prueba, organizado en un vector de 64 filas. También se proveen archivos similares para otro sistema diferente, conteniendo tetraciclina, con los nombres "TE_CON_C.TXT" (50 valores), "TE_RES_C.TXT" (matriz de 101×50), "TE_CON_T.TXT" (57 valores) y "TE_RES_T.TXT" (matriz de 101×57). 4) Para realizar la validación cruzada, ejecutar la rutina PCR_CV.M. 5) Para calibrar el modelo PCR, ejecutar la rutina PCR_CAL.M, e ingresar la información solicitada. Este programa genera el archivo "U_.TXT" en el que se guardan los autovectores o factores, en forma de matriz de J ×A (longitudes de onda × número de factores). 6) Los factores pueden graficarse en cualquier programa gráfico, por ejemplo EXCEL, siguiendo pasos análogos a los descriptos en el capítulo 2. 7) Para predecir muestras incógnita, luego de haber calibrado el modelo, ejecutar la rutina PCR_PRED.M, e ingresar la información solicitada. Instrucciones específicas para el Capítulo 7:
1) Ubicar el archivo comprimido "MATLAB CLASE 6.ZIP" en un directorio del disco rígido. 2) Descomprimir el archivo anterior. 3) Editar los archivos de texto (doble click sobre el ícono del archivo) y examinar su contenido. El archivo "BR_CON_C.TXT" contiene las concentraciones para calibrar el modelo PLS para 12 muestras conteniendo bromhexina, en una única columna. El archivo "BR_RES_C.TXT" contiene las respuestas instrumentales de las 12 muestras de calibración, organizados en una matriz de 64 filas y 12 columnas. El archivo "BR_CON_T.TXT" contiene las concentraciones para validar el modelo PLS para 11 muestras, en una única columna. El archivo "BR_RES_T.TXT" contiene las respuestas instrumentales de las 11 muestras de validación, organizados en una matriz de 64 filas y 11 columnas. El archivo "BR_RES_P.TXT" contiene las respuestas instrumentales de una muestra de prueba, organizado en un vector de 64 filas. También se proveen archivos similares para otro sistema diferente, conteniendo tetraciclina, con los nombres "TE_CON_C.TXT" (50 valores), "TE_RES_C.TXT" (matriz de 101×50), "TE_CON_T.TXT" (57 valores) y "TE_RES_T.TXT" (matriz de 101×57). 4) Para realizar la validación cruzada, ejecutar la rutina PLS_CV.M. 5) Para calibrar el modelo PLS, ejecutar la rutina PLS_CAL.M, e ingresar la información solicitada. Este programa genera el archivo "W_.TXT" en el que se guardan los factores, en forma de matriz de J ×A (longitudes de onda × número de factores). 137