Curso de quimiometría

Calibraciones univariada y multivariada de primer orden

Alejandro C. Olivieri

Departamento de Química Analítica, Facultad de Ciencias Bioquímicas y Farmacéuticas, Universidad Nacional de Rosario,

Suipacha 531, Rosario (S2002LRK), Argentina. E-mail: aolivier@fbioyf.unr.edu.ar

2

Tabla de Contenidos

CLASE 1 3 REGRESIÓN LINEAL 3 Material suministrado con la clase 1 3 Parte 1: calibración univariada 3 Determinación del extremo superior del rango lineal 4 Preparación de patrones 4 Medición de la respuesta de los patrones 5 Estimación de los parámetros de la regresión 5 Predicción en muestras incógnita 6 Cifras de mérito del método 7 Sensibilidad de calibración 7 Sensibilidad analítica 7 Límite de detección 8 Límite de cuantificación 9 Rango dinámico 9 Rango lineal 9 Programas de computación 12 EJERCICIO RESUELTO 12 RESPUESTA DETALLADA 14 EJERCICIOS PROPUESTOS 18 CLASE 2 20 REGRESIÓN LINEAL 20 Material suministrado con la clase 2 20 Parte 2: exactitud y comparación de métodos analíticos 20 Exactitud de un método analítico 21 Región de confianza en el caso homoscedástico 22 Regresión ponderada 24 Región de confianza en el caso heteroscedástico 25 Comparación de métodos analíticos 26 Programas de computación 27 EJERCICIO RESUELTO 27 RESPUESTA DETALLADA 27 EJERCICIOS PROPUESTOS 32 CLASE 3 34 CALIBRACIÓN BIVARIADA 34 Material suministrado con la clase 3 34 Determinación de dos analitos usando dos sensores 34 La etapa de calibración 35 La calibración en notación matricial 35 Etapa de predicción 37 Coeficientes de regresión 38 Colinealidad 38 Cifras de mérito 39 EJERCICIO RESUELTO 40 RESPUESTA DETALLADA 41 EJERCICIO PROPUESTO 43 CLASE 4 44 CALIBRACIÓN MULTIVARIADA 44 Material suministrado con la clase 4 45 Determinación de multianalitos usando múltiples sensores 45 El modelo CLS en notación matricial: etapa de calibración 45

Etapa de predicción y coeficientes de regresión48 Cifras de mérito 49 Colinealidad espectral 50 Interferentes no modelados 50 Ventajas y desventajas de CLS 50 Comparación de métodos 51 EJERCICIO RESUELTO 52 RESPUESTA DETALLADA 53 EJERCICIO PROPUESTO 56 CLASE 5 58 CALIBRACIÓN MULTIVARIADA 58 Material suministrado con la clase 5 58 Regresión por cuadrados mínimos inversos 58

Calibración 59 Predicción 61 Ventajas y desventajas de ILS 62

Regresión por componentes principales 62 Compresión de la información 62 Componentes principales y fuentes de variación espectral 64 Calibración 67 Predicción 67 Validación cruzada 68 Residuos espectrales 70 Cifras de mérito 70 Ventajas y desventajas de PCR 70 Más allá de PCR 71

EJERCICIO RESUELTO 71 RESPUESTA DETALLADA 72 EJERCICIO PROPUESTO 77 CLASE 6 78 CALIBRACIÓN MULTIVARIADA 78 Material suministrado con la clase 6 78 Regresión por cuadrados mínimos parciales 78

Un algoritmo iterativo para PCR 79 Un algoritmo iterativo para PLS 79 Calibración 80 Predicción 80 Residuos espectrales y cifras de mérito 81 Ventajas y desventajas de PLS 81

Más allá de PLS 81 EJERCICIO RESUELTO 82 RESPUESTA DETALLADA 83 EJERCICIO PROPUESTO 86 RESOLUCIONES A LOS EJERCICIOS PROPUESTOS 87 RESPUESTAS A LOS EJERCICIOS PROPUESTOS EN LA CLASE 1 87 RESPUESTA A LOS EJERCICIOS PROPUESTOS EN LA CLASE 2 88 RESPUESTA AL EJERCICIO PROPUESTO EN LA CLASE 3 90 RESPUESTA AL EJERCICIO PROPUESTO EN LA CLASE 4 91 RESPUESTA AL EJERCICIO PROPUESTO EN LA CLASE 5 93 REFERENCIAS 95

3

La estadística es un método sistemático para llegar a la conclusión incorrecta con un 95% de confianza.

Clase 1

Regresión lineal

"Camino recto", fotografía, tomada de www34.brinkster.com.

Material suministrado con la clase 1 Para esta clase se proveen los siguientes archivos: • Archivos de texto (*.TXT) conteniendo datos típicos. • Rutinas (*.M) para el entorno de programación MATLAB. • COMO OPERAR CON MATLAB.PDF, documento de Adobe que explica el empleo

del entorno MATLAB. • Programas ejecutables en QB (*.EXE). • COMO OPERAR CON QB.PDF, documento de Adobe que explica el uso de los

programas en QB.

Parte 1: 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 en este

4

documento, pero se recomienda consultar paralelamente el ejemplo concreto que se analiza en la sección Ejercicio Resuelto.

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 • 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 Gardiner2 y Miller y Miller.3

Determinación del extremo 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 de patrones 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 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 el hecho de que la recta de 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 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 100 mg de reactivo, preparar soluciones en matraces calibrados de al menos 100 mL, tomar alícuotas con pipetas aforadas calibradas, etc.

5

Medición de la respuesta de los patrones 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 de los parámetros de la 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:

A =

∑

∑

=

=

−

−−= m

ii

m

iii

xx

xy

xx

yyxx

QQ

1

2

1

)(

))(( (1)

B = xAy − (2) donde xi es la concentración de cada uno de los m patrones de calibrado, x es el promedio de las concentraciones de calibrado, yi es la respuesta en cada punto e y es el promedio de las respuestas de los patrones de 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 siguientes ecuaciones:

sA = xx

xy

Qs / (3)

sB = xx

xy Qx

ms

2

/1

+ (4)

En las ecuaciones precedentes, el parámetro sy/x es el desvío estándar de los residuos de la regresión y está dado por:

sy/x = 2

)ˆ(1

2

−

−∑=

m

yym

iii

(5)

donde yi es la respuesta experimental de cada patrón de calibrado e iy representa la respuesta estimada en cada punto, esto es, iy = A xi + B. En la ecuación (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 sy/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 – iy ) 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.

6

Predicción en muestras incó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 (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) = xx

xy

QAyy

mnAs

2

2inc/ )(11 −

++ = xx

xy

Qxx

mnAs 2

inc/ )(11 −++ (7)

donde sy/x es el desvío estándar de los residuos de la regresión dado por la ecuación (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 Qxx fue definido en la ecuación (1).

La ecuación (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 (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 (7), y la segunda con los términos ⎟⎟⎠

⎞⎜⎜⎝

⎛ −+

xxQxx

m

2inc )(1 , que

colectivamente reciben el nombre 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 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 (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 (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) = sy/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 (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).

7

Cifras de mérito del método 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. Esta última 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 (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.

Sensibilidad 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 / sy (9) donde sy 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 (sy) a través de los desvíos de las réplicas de las mediciones de calibrado respecto de sus promedios:

sy = pm

yyp

i

r

jiij

−

−∑∑= =1

2

1)(

(10)

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

iy es el promedio de las respuestas de las réplicas para cada nivel de concentración. En la ecuación (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 iy . Este cálculo se ilustra en forma detallada en el ejercicio resuelto que acompaña al presente documento.

En el segundo método de estimación del nivel de ruido, se lo estima como el desvío estándar de los residuos de la regresión lineal, el parámetro ya definido sy/x [véase la ecuación (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.

8

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 (7), escrita del modo siguiente:

s(xinc) = xx

xy

Qxx

mnAs 2

inc/ )(11 −++ (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 (11) se reduce a:

s0 = xx

xy

Qx

mAs 2/ 1

31

++ (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 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 (LC en la Figura 1), a partir de la cual se toman decisiones respecto de la detección del analito. Para concentraciones superiores a LC, 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 1, la probabilidad de cometer este error de tipo I está dada por la zona sombreada de azul (á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 LC 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 de rojo en la Figura 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 LC a la concentración correspondiente a dicho valor de β es el producto del coeficiente tβ,ν por 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 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 × t0,05,m–2 × s0 (13) definición que ha sido adoptada también por IUPAC6 e ISO.7 En la práctica, dado que m es un número relativamente grande, el valor de (2×t0,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 = 3sbl / A, donde sbl 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 t0,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.

9

Figura 1. Prueba de significación empleada para estimar el límite de detección. LC 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 (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, adaptada de la 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 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). Una manera conveniente de medir el cumplimiento de la linealidad es a través de la relación que existe entre la variancia de la regresión, medida por (sy/x)2 [ecuación (5)], y la del ruido instrumental, medida por (sy)2 [ecuación (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 (sy/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:

0

αβ

Hipótesis alternativa: analito presente a este nivel

Predicción LC

Hipótesis nula: analito ausente

LOD

(tα,ν + tβ,ν) s0

10

Concentración

Res

pues

ta

Rango lineal

Rango dinámico

LODLOQ

Extremo superiordel rango lineal

Pérdida de la relaciónrespuesta-concentración

Fexp = ( )( )2

2/

y

xy

s

s (15)

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 Fexp < F, se acepta que los datos se comportan linealmente. Alternativamente, se calcula la probabilidad pF asociada a este valor de Fexp, 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

Figura 2. Rangos dinámico y lineal de un método analítico.

11

Concentración

Res

iduo

s0

Res

iduo

s

0

Res

iduo

s

0

A

B

C

Figura 3. 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.

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 Fexp 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 3 ilustra casos representativos al respecto. En el caso A), el comportamiento es lineal: se espera que la

12

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 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 computación Los métodos descritos en esta clase 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 Gardiner2 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.

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.PDF', 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.PDF'.

Ejercicio resuelto 1) La Tabla 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.

13

Tabla 1. Concentraciones y respuestas para un rango en el que se sospecha que existen desvíos de la linealidad.

Concentración del patrón

Respuesta 1 Respuesta 2 Respuesta 3

0,00 1,00 2,00 3,00 4,00 5,00 6,00 7,00 8,00 9,00 10,00

0,06 1,44 2,82 4,15 5,29 6,61 7,79 8,89 10,03 10,84 11,87

0,08 1,56 2,76 4,20 5,46 6,54 7,70 8,97 9,88 10,91 11,81

–0,06 1,41 2,90 4,08 5,52 6,69 7,69 8,83 9,77 10,65 11,90

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 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 muestra los mismos datos que la Tabla 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).

Tabla 2. Concentraciones y respuestas para un rango en el que existe linealidad.

Concentración del patrón

Respuesta 1 Respuesta 2 Respuesta 3

0,00 1,00 2,00 3,00 4,00 5,00

0,06 1,44 2,82 4,15 5,29 6,61

0,08 1,56 2,76 4,20 5,46 6,54

–0,06 1,41 2,90 4,08 5,52 6,69

Usuarios de MATLAB: los datos de la Tabla 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'.

14

Calcule los valores de la pendiente y ordenada al origen para la recta ajustada con los datos de la Tabla 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 3 muestra los valores de la respuesta para cuatro muestras incógnita, todos por

triplicado.

Tabla 3. Respuestas para cuatro muestras incógnita. Muestra Respuesta 1 Respuesta 2 Respuesta 3

1 2 3 4

0,69 2,20 3,55 4,82

0,65 2,13 3,41 4,71

0,75 2,05 3,52 4,70

Los datos de la Tabla 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 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 detallada 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 Fexp de 8,88, con una probabilidad asociada pF de 0,001. La gráfica de los residuos es informativa al respecto:

15

2) Los valores estimados, dados por las ecuaciones (1) y (2) son, para el ejemplo de la Tabla 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 sy/x = 0,1, sA = 0,01 y sB = 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 3 encontrará un resumen de todos los cálculos intermedios necesarios para

estimar A, B y sus errores estándar.

Tabla 3. Parámetros necesarios para el cálculo de A, B, sA y sB. i xi xi – x yi yi – y (xi – x )2 (xi – x ) (yi – y ) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

0,00 1,00 2,00 3,00 4,00 5,00 0,00 1,00 2,00 3,00 4,00 5,00 0,00 1,00 2,00 3,00 4,00 5,00

–2,50 –1,50 –0,50 0,50 1,50 2,50 –2,50 –1,50 –0,50 0,50 1,50 2,50 –2,50 –1,50 –0,50 0,50 1,50 2,50

0,06 1,44 2,82 4,15 5,29 6,61 0,08 1,56 2,76 4,20 5,46 6,54 –0,06 1,41 2,90 4,08 5,52 6,69

–3,36 –1,98 –0,60 0,73 1,87 3,19 –3,34 –1,86 –0,66 0,78 2,04 3,12 –3,48 –2,01 –0,52 0,66 2,10 3,27

6,25 2,25 0,25 0,25 2,25 6,25 6,25 2,25 0,25 0,25 2,25 6,25 6,25 2,25 0,25 0,25 2,25 6,25

8,39 2,97 0,30 0,37 2,81 7,98 8,34 2,79 0,33 0,39 3,06 7,81 8,69 3,01 0,26 0,33 3,15 8,18

Total Qxx = 52,5 Qxy = 69,17 Promedio x = 2,50 y = 3,42

4) Los valores de predicción se muestran en la Tabla 4. Tabla 4. Predicciones en muestras incógnita. Muestra Respuesta

promedio (yinc) Concentración predicha (xinc)

Desvío estándara s(xinc)

DSR = 100 s(xinc) / xinc (%)

1 2 3 4

0,70 2,13 3,49 4,74

0,44 1,52 2,56 3,51

0,05 0,05 0,05 0,05

12 3,3 1,9 1,4

a A partir de la ecuación (6), insertando sy/x = 0,1; A = 1,32; n = 3; m = 18; yinc de la columna 2 de la Tabla 4, y = 3,42 y Qxx = 52,5. Note que los valores pueden aproximarse por s(xinc) = sy/x / (A n1/2), tal como se dijo en la parte teórica.

16

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 4 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 4 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 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 No. 4 en la Tabla 4:

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. 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: Fexp = 1,58, pF = 0,21. La Tabla 5 ilustra el cálculo detallado de sy 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, p = 6, r = 3, sy = 0,08 (véase la Tabla 5 para el detalle del cálculo).

17

Tabla 5. Parámetros requeridos para el cálculo de sy. i j yij iy (yij – iy )2 1

1 2 3

0,06 0,08 –0,06

0,03

0,0009 0,0025 0,0081

2

1 2 3

1,44 1,56 1,41

1,47

0,0009 0,0081 0,0036

3

1 2 3

2,82 2,76 2,90

2,83

0,0001 0,0049 0,0049

4

1 2 3

4,15 4,20 4,08

4,14

0,0001 0,0036 0,0036

5

1 2 3

5,29 5,46 5,52

5,42

0,0169 0,0016 0,0100

6

1 2 3

6,61 6,54 6,69

6,61

0,0000 0,0049 0,0064

Total

=−∑∑

= =

p

i

r

jiij yy

1

2

1)( 0,081

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, sy/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 (10), como γ = SEN / sy/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×t0,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 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 No. 1 de la Tabla 4 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 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 1, y se comprueba que para este caso, si se incluyen todos los datos, Fexp = 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 Fexp y sus pF asociadas, se obtienen los resultados informados en la Tabla 6.

18

Tabla 6. Rangos de concentración y estudio de la linealidad mediante la prueba F. Rango de concentración Fexp pF

0-10 0-9 0-8 0-7 0-6 0-5

8,88 6,69 4,62 3,50 2,73 1,58

0,001 0,001 0,001 0,007 0,031 0,214

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 7 resume las cifras de mérito calculadas.

Tabla 6. Cifras de mérito. Cifra de mérito Valor (unidades)

Sensibilidad Sensibilidad analítica Límite de detección Límite de cuantificación Rango dinámico Rango lineal

SEN = 1,32 (Unidades de respuesta)×(Unidades de concentración)–1 γ = SEN / sy/x = 13 (Unidades de concentración)–1

LOD = 0,2 (Unidades de concentración) LOQ = 0,6 (Unidades de concentración) 0,2-10,0 (Unidades de concentración) 0,6-6,0 (Unidades de concentración)

Ejercicios propuestos 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:

Concentraciones de patrones y respuestas obtenidas mediante dos métodos analíticos. Concentración

del patrón Método A Método B

Respuesta 1

Respuesta 2

Respuesta 3

Respuesta 1

Respuesta 2

Respuesta 3

0,000 0,100 0,200 0,300 0,400 0,500

0,01 0,17 0,32 0,48 0,64 0,79

0,02 0,17 0,33 0,48 0,64 0,79

0,02 0,17 0,32 0,48 0,64 0,79

2,0 17,4 32,5 47,8 63,2 78,4

1,9 17,4 32,6 47,8 63,3 78,5

1,9 17,3 32,6 48,0 63,3 78,4

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:

19

Método A Método B

Respuesta 1

Respuesta 2

Respuesta 3

Respuesta 1

Respuesta 2

Respuesta 3

0,25 0,26 0,25 25,2 25,1 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 Ramana 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?

Caso A

Muestra Concentración Respuesta 1 Respuesta 2 Respuesta 3 1 2 3 4 5 6

0,000 0,198 0,392 0,583 0,769 0,950

0,78 3,38 5,75 8,53 10,97 13,40

0,80 3,44 6,16 8,51 11,04 13,08

0,82 3,51 6,01 8,68 10,89 13,37

Caso B

Muestra Concentración Respuesta 1 Respuesta 2 Respuesta 3 1 2 3 4 5 6 7 8

0,000 0,198 0,392 0,583 0,769 0,950 1,130 1,310

0,01 1,96 3,75 5,59 7,30 9,07 10,83 12,08

0,03 1,88 3,75 5,52 7,35 8,95 10,71 12,11

0,04 1,90 3,80 5,56 7,27 9,03 10,46 12,21

20

El 42,57 % de toda la estadística está equivocado.

Clase 2 Regresión lineal

"Elliptical viewpoint", escultura, tomada de www.sculpture-design.com.

Material suministrado con la clase 2 Para esta clase se proveen los siguientes archivos: • LECTURA ADICIONAL CLASE 2.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.

Parte 2: exactitud y comparación de métodos analíticos En este segundo 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 en este documento, pero se recomienda consultar paralelamente el ejemplo concreto que se analiza en la sección Ejercicio Resuelto.

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,9 así como en la obra clásica de Massart y colaboradores.10

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.

21

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 la presente clase.

Exactitud 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 de la Clase 1.

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:

s(yi) = 1

)(1

2

−

−∑=

n

yyn

jiij

(1)

En la ecuación (1), yij indica la concentración para el patrón i en la réplica j, e iy 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 variancia de la variable y sea aproximadamente constante, u homoscedástica.11 La Figura 1 muestra las diferencias entre una variancia homoscedástica y otra 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ó en la Clase 1, 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 en esta clase no es homoscedástica. En estos casos, se recomienda realizar una regresión lineal mediante cuadrados mínimos ponderados (WLS, por weighted least-squares) y no una regresión ordinaria (OLS, por ordinary least-squares) como la empleada en la Clase 1.

Dado que el método WLS es más complicado que el OLS, lo recomendable es previamente verificar si efectivamente la variancia no es constante, para utilizar el primero en los casos en los que es estrictamente necesario. Una prueba de constancia de la variancia (o prueba de la homoscedasticidad) puede realizarse mediante el uso del parámetro estadístico F, calculando el valor "experimental" Fexp definido por el cociente entre el máximo y el mínimo valor de las variancias en las réplicas de los patrones [se toma como medida de cada variancia el valor de s(yi)2]:

22

Fexp = [ ][ ]2

2

)(min)(max

i

i

ysys (2)

Este valor se compara luego con el valor crítico de tablas para n – 1 y n – 1 grados de libertad (usualmente con el 95% de confianza). Si Fexp > Fcrit 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.

Figura 1. Arriba, variancia homoscedástica; abajo, variancia heteroscedástica.

Región de confianza en el caso homoscedástico Si se ha podido aplicar el método OLS descrito en la Clase 1, debido a que las variancias

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

23

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:9

2,22

/1

22

1

2 2)())((2)( −==

=−α+−β−α+−β ∑∑ qxy

q

ii

q

ii FsxAxBABq (3)

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 F2,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 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 sy/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 sy/x 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 para aclarar este punto.

Nótese que el valor de sy/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):

RMSE = q

yy∑ − 2nominalpredicho )(

(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.

24

Pendiente

1.0 1.1

Ord

enad

a al

orig

en

-0.2

0.0

0.2

Pendiente

1.0 1.1

Ord

enad

a al

orig

en

-0.2

0.0

0.2

Figura 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).

Figura 3. Distintos tipos de elipses, de acuerdo con la exactitud y precisión: verde, exacta y precisa; celeste, exacta e imprecisa; amarilla, inexacta e imprecisa; naranja, inexacta y precisa. El cuadrado negro marca el punto ideal (1,0).

Regresión ponderada 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):

SC = ∑=

−q

iiii yyw

1

2)ˆ( (5)

Pendiente

Ord

enad

a al

orig

en ■

25

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 iy 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 variancia de la variable en el punto i:

wi = 2)(

1

iys (6)

El efecto concreto del pesado de los datos en forma inversamente proporcional a su variancia 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:

A = 2

1

1

)(

))((

w

q

iii

q

iwiwii

xxw

yyxxw

−

−−

∑

∑

=

= (7)

B = wy – A wx (8) donde xi es la concentración de cada uno de los q patrones de validación, y los parámetros wx e wy son las coordenadas del centro de gravedad pesado por donde pasa la recta ajustada, que están dadas por:

wx =

∑

∑

=

=q

ii

q

iii

w

xw

1

1 (9)

wy =

∑

∑

=

=q

ii

q

iii

w

yw

1

1 (10)

En el método WLS el parámetro sy/x (el desvío estándar de los residuos de la regresión) está dado por:

sy/x = 2

)ˆ(1

2

−

−∑=

q

yywq

iiii

(11)

donde yi es la respuesta experimental, e iy representa la respuesta estimada en cada punto, esto es, iy = 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 la Clase 1.

Región de confianza 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:

26

2,22

/1

22

11

2 2)())((2)( −===

=−α+−β−α+−β ∑∑∑ qxy

q

iii

q

iii

q

ii FsxwAxwBAwB (12)

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 variancias. 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 en este caso difiere 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 least-squares).12

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 (5), excepto que los pesos son una función de las variancias en ambas variables:

[ ] 1222 )()( −+= iii xsAysw (13) En otras palabras, los pesos de la regresión "doblemente ponderada" BLS se eligen como

inversamente proporcionales a una combinación de las variancias 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 (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 (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 variancia 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 (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 el trabajo que se adjunta (LECTURA ADICIONAL CLASE 2.PDF), en el que se ilustran los peligros de no emplear el método correcto de regresión para la comparación de métodos 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.

27

Programas de computació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.

Ejercicio resuelto 1) La Tabla 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 1. Concentraciones nominales de patrones, y valores hallados por un método analítico (con sus desvíos estándar).

Muestra Nominal Hallada (promedio de cinco réplicas)

Desvío estándar

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,06 0,05 0,04 0,02 0,03 0,04 0,06

2) La Tabla 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.

Tabla 2. Concentraciones halladas por dos métodos analíticos con sus desvíos estándar.

Muestra Método 1 Desvío estándar

Método 2 Desvío estándar

1 2 3 4 5 6 7

0,05 5,16 9,91 14,90 19,80 24,90 30,00

0,03 0,02 0,02 0,01 0,02 0,01 0,03

0,06 5,02 10,00 15,20 19,90 25,00 30,00

0,06 0,05 0,04 0,02 0,03 0,04 0,06

Respuesta detallada 1) En primer lugar debemos determinar si los datos de la Tabla 1 son homoscedásticos.

Para ello calculamos el cociente:

Fexp = [ ][ ] 9

)02,0()06,0(

)(min)(max

2

2

2

2==

i

i

ysys

28

Dado que este último valor es mayor que el de tabla [Fcrit (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:

wi =

∑=

−

−

q

ii

i

ys

yqs

1

2

2

)(

)(

De esta manera, se consigue que la suma de los pesos sea igual a q, lo que facilita los cálculos.

Tabla 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:

w1 = 33,0

)06,0(1

)04,0(1

)03,0(1

)02,0(1

)04,0(1

)05,0(1

)06,0(1

)06,0(7

2222222

2=

++++++

Luego debemos calcular los valores de los diferentes productos de variables y pesos, que se muestran en la Tabla 4.

Tabla 4. Cálculos parciales para el método WLS. i ii xw 2

ii xw ii yw iii yxw 1 2 3 4 5 6 7

0,0167 2,4839 7,4538 44,8281 26,4756 18,7285 10,0287

0,0008 12,8169 73,8671 667,9384 524,2178 466,3399 300,8596

0,0201 2,4165 7,5215 45,7307 26,6094 18,8037 10,0287

0,0010 12,4692 74,5380 681,3868 526,8653 468,2128 300,8596

Total 110,0153 2.046,0405 111,1304 2.064,3327 Con los resultados anteriores, calculamos:

wx = 110,0153 / 7 = 15,72

wy = 111,1304 / 7 = 15,88

29

A = 2

1

1

)(

))((

w

q

iii

q

iwiwii

xxw

yyxxw

−

−−

∑

∑

=

= =

= 0022,1)72,15(70405,046.2

88,1572,1573327,064.22

2

1

2

1 =×−

××−=

−

−

∑

∑

=

=

w

q

iii

ww

q

iiii

xvxw

yxvyxw

B = wy – A wx = 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:

sA = xx

xy

Qs /

sB = xx

wxy Q

xm

s2

/1

+

donde sy/x se determina mediante la ecuación apropiada para datos pesados (WLS), tal como se describió en la parte teórica:

sy/x = 2

)ˆ(1

2

−

−∑=

q

yywq

iiii

= 0,16

Por su parte, Qxx está dado por:

Qxx = 2

1

2w

q

iii xqxw −∑

=

= 316,2

A partir de estos parámetros, se obtiene (redondeando a una cifra significativa): sA = 0,01 sB = 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

iii xw

1110,0153

=∑=

q

iii xw

1

2 2.046,0405

=2/ xys 0,026

=−2,2 qF 8,6 Por lo tanto, la ecuación de la elipse estará dada por:

44,0)1(0405,046.2)1,0)(1(0306,220)1,0(7 22 =−α+−β−α+−β

30

La ecuación anterior tiene la siguiente forma: 4

232

21 )())(()( aBaBAaAa =−β+−β−α+−α

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

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):

[ ] 0)())(()( 42

122

3 =−−α+−β−α+−β aAaBAaBa Luego calculamos los valores de (β – B) a partir de la resolvente de segundo grado:

(β – B) = [ ]

3

42

1322

22

2)(4)()(

aaAaaAaAa −−α−−α±−α−

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:

[ ] 0)(4)( 42

1322

2 =−−α−−α aAaaAa de donde se pueden calcular los límites superior e inferior de (α – A) como:

LIM(α – A) = ± 13

22

43

44

aaaaa

+− = ± 0,0373

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 + [ ]

3

42

1322

22

2))(4)()(

aaAaaAaAa −−α−−α±−α−

Ejemplos de pares de valores de α y β calculados con la ecuación anterior son:

α – A α β –0,0373 –0,0273 –0,0173 –0,0073 0,0027 0,0127 0,0227 0,0327

0,9627 0,9727 0,9827 0,9927 1,0027 1,0127 1,0227 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

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):

31

Pendiente (α)

0.96 0.98 1.00 1.02 1.04 1.06

Ord

enad

a al

orig

en (β

)

-1

0

1 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 contenidos 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 la Clase 1. 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 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) Nótese 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 la Clase 1. Esta rutina proporciona

32

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 propuestos 1) Los valores siguientes corresponden 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.

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 que sigue 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.

33

Muestra Nominal Método 1 Método 2 Método 3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

190 87 23 13 38 150 26 58 125 65 90 160 48 75 0 0 0 0

173 80 26 6 19 142 33 67 146 63 89 158 41 64 10 5 3 11

214 86 29 14 28 145 16 60 126 67 92 172 52 68 11 8 7 7

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.

34

Un estadístico tenía sus pies sobre hielo

y su cabeza en un horno encendido. Al preguntársele cómo se sentía, respondió:

"en promedio, me siento bien".

Clase 3 Calibración bivariada

Material suministrado con la clase 3 Para esta clase se proveen los siguientes archivos: • MATRICES PARTE 1. PDF, documento de Adobe conteniendo conceptos básicos

sobre álgebra matricial. • LECTURA ADICIONAL CLASE 3.PDF, documento de Adobe con un trabajo

educativo para lectura adicional.

Determinación de dos analitos usando dos sensores 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

35

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.13,14

La teoría se expone en este documento, pero se recomienda consultar paralelamente el ejemplo concreto que se analiza en la sección Ejercicio Resuelto. 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 estos 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 clase, pueden consultar el documento adjunto 'MATRICES PARTE 1.PDF'.

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 soluciones a las longitudes de onda 1 y 2. Las correspondientes respuestas instrumentales Yij (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:

⎥⎦

⎤⎢⎣

⎡=

2221

1211

YYYY

Y (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 Xin es la concentración en la mezcla i del analito n:

36

⎥⎦

⎤⎢⎣

⎡=

2221

1211

XXXX

X (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:

Y11 = X11 S11 + X12 S12 (3) donde S11 y S21 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 (3). En general Y podrá escribirse entonces mediante el siguiente producto matricial:

Y = X ST (4) donde S es una matriz (2×2) cuyo elemento genérico Sjn es la sensibilidad a la longitud de onda j del analito n.

La ecuación (4) puede representarse en forma gráfica mediante la Figura 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).

Figura 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 Sjn 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 Sjn "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) para mantener la consistencia del producto matricial (Figura 1).

La matriz S puede obtenerse a partir de la ecuación (4), aunque se requieren varias etapas. En primer lugar, es necesario pre-multiplicar ambos miembros por X–1 (nótese 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 (5) En la ecuación (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

Estos números deben coincidir (N)

Estos números deben coincidir (I)

Estos números deben coincidir (J)

I×J = I×N × N×J

Y

X

ST

37

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:

⎥⎦

⎤⎢⎣

⎡=⎥

⎦

⎤⎢⎣

⎡=

,cal2

,cal1

2221

1211

00

xx

XXXX

X (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 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 (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 (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 = ⎥⎦

⎤⎢⎣

⎡

2212

2111

YYYY

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

cal

cal

x

x

,2

,110

01

=

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

calcal

calcal

xY

xY

xY

xY

,2

22

,1

12

,2

21

,1

11

(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:

⎥⎦

⎤⎢⎣

⎡=

2

1

yy

y (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):

y = S x (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 (9):

x = S–1 y (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.

38

Coeficientes de regresión La ecuación (10) puede interpretarse en forma gráfica mediante el siguiente esquema:

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: xn = (nava fila de S–1) × y (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 (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 (11) se transforma en: xn = βn

T y = β1n y1 + β2n y2 (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 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.

Tabla 1. Tamaños en las matrices en determinaciones utilizando múltiples sensores.

Matriz / vector

Concepto Modelo bivariado Modelo multivariado

Y Señales de calibrado 2×2 I×J X Concentraciones de

calibrado 2×2 I×N

S Sensibilidades 2×2 J×N βn Coeficientes de regresión 2×1 J×1 y Señales de la incógnita 2×1 J×1 x Concentraciones

estimadas en la incógnita 2×1 N×1

Colinealidad Un examen cuidadoso de la ecuación (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

x = S–1 × y 2×1 2×2 2×1

1ra. fila de S–1x1

y

2da. fila de S–1x2

39

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 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 2, en la que se muestra una situación deseable (izquierda) y una indeseable (derecha) para el análisis a dos longitudes de onda.

Figura 2. Sensibilidad en función de la longitud de onda. Izquierda: situación deseable,

derecha: situación indeseable. Los datos en rojo representan el analito 1, los datos en azul el analito 2.

Cifras de mérito 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 (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 = 22

21

1

nn ββ + (14)

λ1 λ2 λ1 λ2

Sen

sibi

lidad

Sen

sibi

lidad

40

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 2

221 nn ββ + 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 || (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 (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 || (16) Puede demostrarse que SENn 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 puede definirse la sensibilidad analítica, como el cociente entre el valor de SENn y el ruido instrumental sR, obtenido a partir de replicados de una muestra blanco:

γn = SENn / sy (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) = sy / SENn (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 CLASE 3.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 sy / SENn (19)

Ejercicio resuelto 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 muestran en la siguiente tabla:

Composición de la muestra Absorbancia a 440 nm Absorbancia a 545 nm KMnO4 4,075×10–4 M 0,028 0,915 K2Cr2O7 3,030×10–3 M 1,073 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:

41

Muestra Absorbancia a 440 nm Absorbancia a 545 nm 1 0,364 0,220 2 0,722 0,258 3 0,153 0,915 4 0,607 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 detallada 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:

Y = ⎥⎦

⎤⎢⎣

⎡026,0073,1915,0028,0

X = ⎥⎦

⎤⎢⎣

⎡30,300

0075,4×10–4

Luego debemos aplicar la ecuación para el cálculo de S: S = YT (X–1)T Para ello, es necesario en primer lugar invertir la matriz X. Esto es sencillo, puesto que X

es diagonal, de manera que:

X–1 = 1

30,3000075,4 −

⎥⎦

⎤⎢⎣

⎡×104 = ⎥

⎦

⎤⎢⎣

⎡3300

0454.2

Luego se multiplican las traspuestas de Y y X–1:

S = YT (X–1)T = T

026,0073,1915,0028,0

⎥⎦

⎤⎢⎣

⎡T

33000454.2

⎥⎦

⎤⎢⎣

⎡ = ⎥

⎦

⎤⎢⎣

⎡026,0915,0073,1028,0

⎥⎦

⎤⎢⎣

⎡3300

0454.2 =

= ⎥⎦

⎤⎢⎣

⎡9245.2

35469

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:1 1 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 columna j. En el

caso estudiado, det(S) = –7,95×105, Cof(S)T = T

69354245.29

⎥⎦⎤

⎢⎣⎡−

− = ⎥⎦

⎤⎢⎣⎡−

−69245.23549

, y por lo tanto,

S–1 = (–7,95×105)–1 × ⎥⎦⎤

⎢⎣⎡−

−69245.23549

= ⎥⎦⎤

⎢⎣⎡

−−

87,027,2846,411,0

×10–4.

42

S–1 = ⎥⎦

⎤⎢⎣

⎡−

−87,027,28

46,411,0 × 10–4

Las concentraciones de ambos analitos en las cuatro muestras incógnita son, por lo tanto, las siguientes:

Muestra 1:

x = S–1 y = ⎥⎦

⎤⎢⎣

⎡−

−87,027,28

46,411,0 × 10–4 × ⎥

⎦

⎤⎢⎣

⎡220,0364,0

= ⎥⎦

⎤⎢⎣

⎡1,10

94,0 × 104

Muestra 2:

x = S–1 y = ⎥⎦

⎤⎢⎣

⎡−

−87,027,28

46,411,0×10–4 × ⎥

⎦

⎤⎢⎣

⎡258.0722.0

= ⎥⎦

⎤⎢⎣

⎡2,20

07,1 × 104

Muestra 3:

x = S–1 y = ⎥⎦

⎤⎢⎣

⎡−

−87,027,28

46,411,0×10–4 × ⎥

⎦

⎤⎢⎣

⎡915.0153.0

= ⎥⎦

⎤⎢⎣

⎡5,306,4

× 104

Muestra 4:

x = S–1 y = ⎥⎦

⎤⎢⎣

⎡−

−87,027,28

46,411,0×10–4 × ⎥

⎦

⎤⎢⎣

⎡937.0607.0

= ⎥⎦

⎤⎢⎣

⎡3,16

11,4 × 104

La tabla que sigue resume los resultados.

Muestra Señal Concentración predicha (error

estándar)

M × 104

1 y = ⎥

⎦

⎤⎢⎣

⎡220,0364,0

x1 = 0,94(1)

x2 = 10,1(1)

2 y = ⎥

⎦

⎤⎢⎣

⎡258,0722,0

x1 = 1,07(1)

x2 = 20,2(1)

3 y = ⎥

⎦

⎤⎢⎣

⎡915,0153,0

x1 = 4,06(1)

x2 = 3,5(1)

4 y = ⎥

⎦

⎤⎢⎣

⎡937,0607,0

x1 = 4,11(1)

x2 = 16,3(1)

43

Los errores estándar en las concentraciones se calcularon con el modelo aproximado citado en la teoría, esto es s(xn) = sR / SENn, con sR = 0,003 (para el cálculo de SENn ver el punto 3). Nótese que, debido a la diferencia en sensibilidades, los errores estándar para el analito 1 son 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 = ⎥⎦

⎤⎢⎣

⎡××−

−

−

4

5

1046,4101,1

SEN1 = || β1 ||–1 = 2.243 A M–1 (A = unidades de absorbancia)

β2 = ⎥⎦

⎤⎢⎣

⎡×−×

−

−

4

3

101108,2

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 Téngase en cuenta que para llegar a estos últimos valores es necesario incluir varias cifras

significativas en los cálculos intermedios.

Ejercicio propuesto 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:

Solución patrón

λ1 λ2 λ3 λ4 λ5

Analito 1 1,00×10–4 M

0,550

0,610 0,720 0,850 0,910

Analito 2 1,00×10–4 M

0,510 0,505 0,710 0,800 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?

44

Celebrar cumpleaños es saludable. La estadística demuestra que la gente

que más cumpleaños celebra vive más.

Clase 4

Calibración multivariada

45

Material suministrado con la clase 4 Para esta clase se proveen los siguientes archivos: • MATRICES PARTE 2. PDF, documento de Adobe conteniendo conceptos básicos

sobre álgebra matricial. • LECTURA ADICIONAL CLASE 4.PDF, documento de Adobe con un trabajo

educativo para lectura adicional. • 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.

Determinación de multianalitos usando múltiples sensores En esta sección extenderemos los resultados presentados en la clase 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 trabajo clásico de Haaland sobre el tema.15 Otras lecturas valiosas son los capítulos correspondientes de libros de quimiometría,16,17 así como los ya famosos "tutorials" de Brereton en internet.18

La teoría se expone en este documento, pero se recomienda consultar paralelamente el ejemplo concreto que se analiza en la sección Ejercicio Resuelto.

El modelo CLS en notación matricial: 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 Yij (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:

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

IJII

J

J

YYY

YYYYYY

...............

...

...

21

22221

11211

Y (20)

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

46

calibración (I×N) X, cuyo elemento genérico Xin es la concentración en la mezcla i del analito n:

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

INII

N

N

XXX

XXXXXX

...............

...

...

21

22221

11211

X (21)

Para dos analitos, la ecuación anterior se transforma en:

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

21

2221

1211

......

II XX

XXXX

X (22)

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 Sjn 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 Sjn), 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 (23) donde S es una matriz (J×N) cuyo elemento genérico Sjn 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 (4) para mantener la consistencia del producto matricial. Las relaciones de tamaño entre las matrices de la ecuación (4) se muestran en la Figura 1.

La solución de cuadrados mínimos de la ecuación (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 (4) no puede hacerse simplemente pre-multiplicando 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 pre-multiplicar ambos miembros de la ecuación (4) por la matriz traspuesta de X:

XT Y = XT X ST (24)

47

Figura 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 (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):

ST = (XT X)–1 XT Y (25) Trasponiendo la ecuación anterior para obtener S: S = [(XT X)–1 XT Y]T = YT X (XT X)–1 (26) La ecuación (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 (7) con su análogo de la Clase 3, en que X era cuadrada y podía invertirse directamente. Esta comparación sugiere que en la ecuación (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 (7) puede escribirse en forma más compacta:

S = YT (X+)T (27) 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.

† 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.

I×J = I×N × N×J + I×J

ST

E

Y

X

48

Quienes deseen revisar algunos conceptos sobre matrices y sus operaciones, que son necesarios para seguir en detalle la presente clase, pueden consultar el documento adjunto 'MATRICES PARTE 2.PDF'.

Etapa de predicción y coeficientes 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:

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

=

Jy

yy

...2

1

y (28)

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 (7):

y = S x + e (29) 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 (11) (fijando e = 0). En primer lugar se debe pre-multiplicar la ecuación (10) por ST, de manera que se obtenga una matriz cuadrada en el segundo miembro:

ST y = (ST S) x (30) Luego puede despejarse x pre-multiplicando por la inversa de (ST S): x = (ST S)–1 ST y (31) Nuevamente, podemos definir la seudoinversa de S de tal modo que permita obtener x

directamente, pre-multiplicando a y: x = S+ y (32)

El esquema superior muestra que la ecuación (13) puede interpretarse diciendo que la

concentración de cada analito se predice mediante el siguiente producto escalar: xn = (nava fila de S+) × y (33) 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 (34) Con esta última definición, la ecuación (14) se transforma en: xn = βn

T y = β1n y1 + β2n y2 + ... + βJn yJ (35) lo cual significa que la concentración es el producto escalar del vector de coeficientes de regresión por el vector de respuestas instrumentales.

x = S+ × y N×1 N×J J×1

nava fila de S+xn y

49

Cifras de mérito 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 (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 = 22

221 ...

1

Jnnn βββ +++ (36)

donde βjn son los elementos del vector βn. SENn = 1 / || βn || (37)

donde || · || simboliza el cálculo de la norma de un vector. 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 (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 || (38) Puede demostrarse que SENn 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 sy, obtenido a partir de replicados de una muestra blanco:

γn = SENn / sy (39) 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) = sy / SENn (40) 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 sy / SENn (41) 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 CLASE 4.PDF".

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 (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:

50

sres = NJ

eJ

jj

−

∑=1

2)( (42)

donde ej representa cada uno de los elementos del vector e. Nótese el empleo de J – N grados de libertad en la ecuación (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 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 la Clase 2. 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 (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 (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 (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 desventajas 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

51

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.

Comparación de 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 1. Comparación de las propiedades analíticas de los distintos métodos de calibración. Propiedad Método

Univariado Bivariado Multivariado CLSNúmero de analitos

1 2 Varios

Concentración de componente(s) de calibrado

Conocida Conocidas Conocidas

Efecto de la colinealidad

– Disminuye la sensibilidad, selectividad y

precisión

Disminuye la sensibilidad, selectividad y

precisión Presencia de interferentes

Análisis inexacto Análisis inexacto Análisis inexacto pero con

detección del problema

Cifra de mérito Sensibilidad SEN = A SENn = 1 / || βn || Sensibilidad analítica

γ = SEN / sy γ = SENn / sy

Incertidumbre en la predicción

s(x) =

xx

xy

Qxx

mnAs 2

inc/ )(11 −++

s(xn) ≈ sy / SENn

Límite de detección LOD = 3,3

xx

xy

Qx

mAs 2/ 1

31

++LOD ≈ 3,3 sy / SENn

Límite de cuantificación LOQ = 10

xx

xy

Qx

mAs 2/ 1

31

++ LOQ ≈ 10 sy / SENn

52

En la Tabla 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 respuesta al efecto de la colinealidad (no aplicable al caso univariado). En la Clase 5 analizaremos un método capaz de superar estas dificultades y describiremos sus propiedades en perspectiva con las de la Tabla 1.

Ejercicio resuelto 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:

Muestra de calibrado

λ1 λ2 λ3 λ4 λ5 λ6

1 2 3 4

1,52 2,94 1,47 3,01

2,78 5,42 2,94 5,63

3,32 6,33 3,81 6,80

3,26 5,48 4,21 6,35

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 son las siguientes:

Muestra de calibrado

x1cal x2cal

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:

Muestra de validación

λ1 λ2 λ3 λ4 λ5 λ6

1 2 3 4

1,32 2,73 1,38 1,25

2,59 5,10 2,54 2,42

3,36 5,95 2,92 3,06

3,62 5,23 2,47 3,18

3,55 3,35 1,48 2,90

2,72 2,72 1,24 2,29

53

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.

Muestra de

prueba λ1 λ2 λ3 λ4 λ5 λ6

1 2 3

1,11 5,54 2,56

2,20 6,71 4,76

2,77 7,02 5,81

2,81 5,83 5,56

2,56 3,66 4,39

2,02 2,77 3,50

Respuesta detallada 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:

Y =

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

03,410,535,680,663,501,315,308,421,481,394,247,178,235,348,533,642,594,294,148,226,332,378,252,1

X =

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

22211211

Luego debemos aplicar la ecuación para el cálculo de S: S = YT X (XT X)–1 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.

54

Figura 1: Espectros de calibrado.

Figura 2: Sensibilidades y coeficientes de regresión.

55

Las cifras de mérito calculadas mediante los programas para este modelo son las siguientes:

Cifra de mérito Analito 1 Analito 2

Sensibilidad Sensibilidad analíticaa Selectividad

4,1 Señal × concentración–1

137 concentración–1 0,83

1,9 Señal × concentración–1 63 concentración–1

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 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:

Muestra Analito 1 Analito 2 Residuo espectral

Nominal Predicho Nominal Predicho 1 2 3 4

0,89 1,86 0,93 0,83

0,88(1) 1,87(1) 0,93(1) 0,84(1)

1,69 1,05 0,40 1,34

1,70(1) 1,05(1) 0,40(1) 1,34(1)

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) = sy / SENn, con sy = 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 la Clase 2. 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:

0,89 1,86 0,93 0,83 1,69 1,05 0,40 1,34

0,88 1,87 0,93 0,84 1,70 1,05 0,40 1,34

56

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.

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:

Muestra de prueba

Analito 1 Analito 2 Residuo espectral

Predicho Predicho 1 2 3

0,76(1) 2,52(1) 1,71(1)

1,17(1) 0,62(1) 1,78(1)

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 propuesto 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 el documento 'LECTURA ADICIONAL CLASE 4.PDF'.

57

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.

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.

58

¿Cuántos estadísticos hacen falta para cambiar una lamparita?

Uno ± dos

Clase 5 Calibración multivariada

El monte Nipals, ubicado en el norte de Suecia, homónimo del algoritmo desarrollado por H. Wold para el cálculo de componentes principales.

Material suministrado con la clase 5 Para esta clase se proveen los siguientes archivos: • MATRICES PARTE 3. PDF, documento de Adobe conteniendo conceptos básicos

sobre álgebra matricial. • LECTURA ADICIONAL CLASE 5.PDF, documento de Adobe con un trabajo

educativo para lectura adicional. • 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.

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 en la sección Ejercicio Resuelto.

59

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 (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 (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 mediante un modelo apropiado.

Si bien el modelo CLS puede en principio interpretarse mediante una ecuación similar a la (2), en los métodos inversos la ecuación (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 la presente clase, 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,19 y el artículo de Haaland y Thomas.15

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 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 B + E (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 (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 (3), siendo su tamaño idéntico al de Y.

Para obtener la matriz B, se debe despejar ésta de la ecuación (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 (3) por XT:

XT Y = (XT X) B (4) Aquí se presenta un importante inconveniente del método ILS. Para continuar el proceso a

partir de la ecuación (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

60

numérico aclarará el problema: supongamos que X es una matriz de 2×4 como la que se muestra en la Figura 1.

⎥⎦

⎤⎢⎣

⎡=

01341210

X XT X =

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

12102554151012041216

Figura 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 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 (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 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 (4):

B = (XT X)–1 XT Y (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 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 (yn).

61

Figura 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) La ecuación (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), β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) 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) para cada analito de interés.

La etapa de calibrado es análoga a la descrita en el caso del modelo CLS de la clase anterior, 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.

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 (7) La ecuación (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 (7) se aplicaría tantas veces como

fuese necesario, utilizando cada vez el vector βn asociado al analito n.

= × (XT X)–1 XT yn βn

J × 1 J × I

I × 1

Y B

62

Ventajas y desventajas 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 las clases 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 por componentes principales 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 la presente clase, pueden consultar el documento adjunto 'MATRICES PARTE 3.PDF'.

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 descomposición en valores singulares, que consiste en descomponer a la matriz XT (tamaño J×I) en el producto de otras tres matrices:

XT = U W VT (8) Podemos apreciar los requerimientos de tamaño matricial en esta última ecuación a través

del esquema presentado en la Figura 3.

63

Figura 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.

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)T en una única matriz T, la ecuación de descomposición singular (8) se puede también escribir como:

XT = U TT (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 (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 = (UT XT)T = X U (10)

Estos números deben coincidir (I)

Estos números deben coincidir (J)

Estos números deben coincidir (I)

J×I = J×I × I×I × I×I

XT

U W

VT

Estos números deben coincidir (I)

64

Los tamaños de las tres matrices involucradas en la ecuación (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') (11)

La consecuencia más importante de la ecuación (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 (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-linear linear iterative partial least-squares), desarrollado por H. Wold.20

Componentes principales 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 variancia espectral total, calculada del modo que sigue:

% Variancia explicada por los primeros A factores =

∑

∑

=

=I

ii

A

ii

W

W

1

2

1

2

)(

)(100 (12)

donde Wi es cada uno de los elementos diagonales de la matriz W, y (Wi)2 es el correspondiente autovalor asociado al componente principal i.

65

Dado que cada componente principal contribuye con una porción cada vez menor de la variancia total, lo usual es tomar el número de los primeros componentes principales que, colectivamente, aportan un determinado porcentaje de la variancia total. La Figura 4 ilustra el comportamiento típico de un conjunto de componentes principales: mientras los valores singulares disminuyen, su contribución a la variancia total también disminuye. En el caso de la Figura 4, por ejemplo, los tres primeros componentes principales explican más del 99% de la variancia 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.

Figura 4. Variancia 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 variancia 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 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,

66

básicamente, ruido instrumental. El primero tiene forma "espectral"; el segundo, de ruido al azar.

Figura 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 variancia 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 (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 UA (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: TA = X UA (13) En la ecuación precedente, hemos llamado TA a la matriz de scores estimada con A

factores. Esta nueva matriz TA, 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, 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 (9) escrita en términos de los A componentes selectos, esto es XT = UA TA

T. 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 variancia 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.21

67

Imagen total

A = 1 A = 2 A = 4 A = 8 A = 16 A = 32

Figura 6. Una imagen (arriba al centro), reconstruida utilizando distinto número de componentes principales (abajo).

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 TA:

yn = TA vn + e (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 TA

T se obtiene: TA

T yn = TAT TA vn + e (15)

Luego será necesario multiplicar por la inversa de (TAT TA), o sea, por (TA

T TA)–1. Esta última operación de inversión es trivial, ya que (TA

T TA) es una matriz diagonal (en atención a la ortogonalidad de las 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 se obtiene, entonces:

vn = (TAT TA)–1 TA yn (16)

La inversión de (TAT TA) en la ecuación (16) no presenta problemas asociados a la

colinealidad, por los motivos anteriormente expuestos: las columnas de TA son ortogonales. Por analogía con el criterio empleado en el modelo CLS, podemos llamar a la matriz [(TA

T TA)–1 TA] la seudoinversa de TA y denominarla TA

+, con lo cual la ecuación (16) adopta su forma final:

vn = TA+ yn (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 del analito en muestras incógnitas es también similar.

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 UA, 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.

68

Análogamente a la ecuación (13), entonces, se obtiene el vector tA (A×1) correspondiente a la muestra incógnita:

tA = UAT x (18)

Este vector tA 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:

yn = (vn)T tA (19) en el que vn reemplaza a βn y tA reemplaza a x. A partir de esta última ecuación se estima la concentración del analito en la incógnita.

Validación cruzada 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 se 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 7):

Factores PRESS

1 2 3 4 5 6

0,92 0,0217

4,3×10–3 4,1×10–3 3,7×10–3 5,1×10–3

69

Figura 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.¡Error! Marcador no definido. 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 variancias, 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 para el numerador como para el denominador.

70

La tabla completa sería como sigue:

Factores PRESS PRESS/min(PRESS) p 1 2 3 4 5 6

0,92 0,0217

4,3×10–3

4,1×10–3 3,7×10–3 5,1×10–3

248 5,88 1,17 1,12

1 –

0,999 0,997 0,605 0,576 0,5 –

Haaland propone, basándose en resultados empíricos, seleccionar el 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)]1/2. Este parámetro debe ser del orden de la incertidumbre asociada a las concentraciones de calibrado.

Residuos espectrales 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 (9), pero empleando el vector de scores de la muestra y la matriz de factores reducida UA:

x = UA tA (20) Luego puede definirse el residuo espectral en forma análoga a CLS:

sres = AJ

xxJ

jjj

−

−∑=1

2)ˆ( (21)

Cifras de mérito 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 (20), esto es, "reconstruyendo" el vector espectral βn a partir del vector reducido vn:

βn = UA vn (22) 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 la Clase 4. No obstante, la selectividad puede calcularse en PCR recurriendo a conceptos que están más allá del alcance de este curso. Los programas suministrados con esta clase permiten estimar todas las cifras de mérito.

Ventajas y desventajas de PCR Como resumen de esta clase, podemos enumerar las ventajas de PCR respecto de las otras

técnicas multivariadas que hemos estudiado hasta el presente. PCR combina las ventajas ya

71

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 la Clase 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á de 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 achacarse 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).

Ejercicio resuelto 1) Los datos del presente ejercicio están tomados del trabajo que acompaña la presente

clase (LECTURA ADICIONAL CLASE 5.PDF). 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:

Muestra de calibrado

Concentración×104 M

1 2 3 4 5 6 7 8 9 10 11 12

1.55 2.06 2.58 1.55 2.06 2.58 1.55 2.06 2.58 1.68 2.10 2.66

72

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 detallada 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:

Factores PRESS PRESS/min(PRESS) p

1 2 3 4 5 6 7 8

0,92 0,0217

4,3×10–3

4,1×10–3 3,7×10–3 5,1×10–3

8,9×10–3

1,1×10–2

248 5,88 1,17 1,12

1 – – –

0,999 0,997 0,605 0,576 0,5 – – –

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 variancia de la matriz espectral.

73

Tanto los resultados correspondientes al PRESS como la variancia 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).

74

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:

Cifra de mérito Valor

Sensibilidad Sensibilidad analíticaa

1 / γ Selectividad

1,21×104 A × M–1

4×106 M–1 2,5×10–7 M

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:

Muestra Concentración × 104 Residuo espectral

Nominal Predichaa 1 2 3 4 5 6 7 8 9 10 11

1,96 2,16 0,00 0,82 1,02 1,33 1,84 2,35 1,94 2,14 2,24

1,97(1) 2,19(1) –0,01(1) 0,84(1) 1,04(1) 1,37(1) 1,93(1) 2,43(1) 2,00(1) 2,19(1) 2,25(1)

0,004 0,002 0,014 0,009 0,006 0,005 0,003 0,004 0,004 0,003 0,006

a Los errores estándar en las concentraciones, calculados con el modelo aproximado citado en la teoría, esto es s(xn) = sy / SENn, con sy = 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 la Clase 2. De este modo, la tabla de datos a suministrar a los programas de cálculo de la elipse será como sigue:

75

1,96 2,16 0,00 0,82 1,02 1,33 1,84 2,35 1,94 2,14 2,24

1,97 2,19 –0,01 0,84 1,04 1,37 1,93 2,43 2,00 2,19 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.

76

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

77

Aquí el residuo es significativamente mayor que el ruido espectral, lo que haría sospechar la presencia de un interferente no modelado.

Ejercicio propuesto 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.

78

Un estadístico cometió un delito y fue encarcelado.

Ahora tiene cero grados de libertad.

Clase 6 Calibración multivariada

"PLS", un grafitti tomado de www.home.aone.net.au/byzantium/aerosol/images1.html

Material suministrado con la clase 6 Para esta clase se proveen los siguientes archivos: • Archivos (*.M) con rutinas para el entorno de programación MATLAB. • Archivos (*.EXE) con programas ejecutables en QB.

Regresión por cuadrados mínimos parciales 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

79

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 algoritmo iterativo para PCR Como se vio en la Clase 5, 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 variancia 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 least-squares).20 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 variancia 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 score t1, a través de t1 = XT u1. 3) Se substrae de X la contribución del primer factor, calculada como u1 t1

T, es decir se calcula la diferencia o residuo E = X – u1 t1

T. 4) Se substituye E por X en el primer paso y se continúa hasta llegar al número deseado A de factores. 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 algoritmo 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 llamada W, y otros llamados simplemente loadings, contenidos en una matriz llamada 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 estriba en que en PCR los factores describen la máxima variancia 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:

80

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 / (t1

T t1) 4) Se obtiene el primer loading p1: p1 = XT t1 / (t1

T t1) 5) Se calculan los residuos espectrales y de concentración: eX

T = XT – t1 p1T

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 en este caso 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. 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,..., vA) 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:

tA = (WTP)–1 WT x (1)

81

yn = (vn)T tA (2)

Residuos espectrales y cifras 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 tA (3) Y luego puede definirse el residuo espectral:

sres = AJ

xxJ

jjj

−

−∑=1

2)ˆ( (4)

Ventajas y desventajas 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.

En referencia a la tabla de propiedades analíticas presentada en la Clase 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.

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.22-26

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.27

Debe mencionarse que los métodos para calibración multivariada descritos en este curso 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. En este sentido, la calibración univariada se clasificaría como de orden cero. Existe la posibilidad de realizar una calibración empleando datos matriciales para cada muestras, 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. 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ón-emisión de fluorescencia combinadas con la cinética de una reacción química).

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

82

cuantificar analitos calibrados en presencia de interferencias no calibradas. Esta propiedad está ausente en los datos de primer orden, y presenta 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.28

Ejercicio resuelto 1) Los datos del presente ejercicio están tomados del trabajo que acompaña la clase 5

(LECTURA ADICIONAL CLASE 5.PDF). 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 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:

Muestra de calibrado

Concentración×104 M

1 2 3 4 5 6 7 8 9 10 11 12

1.55 2.06 2.58 1.55 2.06 2.58 1.55 2.06 2.58 1.68 2.10 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 PLS?

83

Respuesta detallada 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:

Factores PRESS PRESS/min(PRESS) p

1 2 3 4 5 6 7 8

0,907 0,021

4,15×10–3 3,65×10–3 5,81×10–3 9,62×10–3 1,64×10–2 2,06×10–2

249002 5,76 1,139

1 – – – –

0,999 0,997 0,587 0,5 – – – –

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.

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 variancia de la matriz espectral.

Tanto los resultados correspondientes al PRESS como la variancia 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).

84

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:

Cifra de mérito Valor

Sensibilidad Sensibilidad analíticaa

1 / γ Selectividad

1,21×104 A × M–1

4×106 M–1 2,5×10–7 M

0,46 a Obtenida dividiendo la sensibilidad por el nivel de ruido instrumental (0,003 unidades).

85

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:

Muestra Concentración × 104 Residuo

espectral Nominal Predichaa 1 2 3 4 5 6 7 8 9 10 11

1,96 2,16 0,00 0,82 1,02 1,33 1,84 2,35 1,94 2,14 2,24

1,97(1) 2,19(1) –0,01(1) 0,84(1) 1,04(1) 1,37(1) 1,93(1) 2,43(1) 1,99(1) 2,19(1) 2,25(1)

0,004 0,002 0,014 0,009 0,007 0,005 0,003 0,004 0,004 0,002 0,005

a Los errores estándar en las concentraciones, calculados con el modelo aproximado citado en la teoría, esto es s(xn) = sy / SENn, con sy = 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 la Clase 2. De este modo, la tabla de datos a suministrar a los programas de cálculo de la elipse será como sigue:

1,96 2,16 0,00 0,82 1,02 1,33 1,84 2,35 1,94 2,14 2,24

1,97 2,19 –0,01 0,84 1,04 1,37 1,93 2,43 1,99 2,19 2,25

86

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.

Ejercicio propuesto Se desea modelar, mediante PLS, 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.

87

Resoluciones a los ejercicios propuestos

Respuestas a los ejercicios propuestos en la clase 1

1) Empleando las ecuaciones de regresión lineal y cálculo de cifras de mérito expuestas en la teoría de la clase 1, se obtienen los siguientes resultados respecto de la sensibilidad:

Método Sensibilidad de calibración Sensibilidad analítica

A B

1,552 153.0

4,3×102 1,8×103

Nótese 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: Método Concentración (desvío

estándar) A B

0,153(1) 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: Caso Sensibilidad Sensibilidad

analítica 1/γ LOD LOQ Rango

lineal A 13,2 114,1 0,009 0,02 0,06 0,06-0,95 B 9,3 91,2 0,011 0,03 0,08 0,08-1,31

Debe notarse que el caso A posee efectivamente un blanco significativo, ya 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.

88

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.

Respuesta a los ejercicios propuestos en la clase 2 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 de la clase 2). 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 sy/x: 2,35 La elipse correspondiente contiene, aunque marginalmente, al punto ideal (1,0):

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 sy/x: 2.39

89Pendiente0.8 0.9 1.0 1.1

Ord

enad

a al

orig

en

-10

0

10

20

30

1

2

3

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.

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:

90

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.

Respuesta al ejercicio propuesto en la clase 3 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:

R = ⎥⎦

⎤⎢⎣

⎡505,0510,0610,0550,0

C = ⎥⎦

⎤⎢⎣

⎡1001

×10–4

Cálculo de S y su inversa:

S = RT (C–1)T = ⎥⎦

⎤⎢⎣

⎡505,0610,0510,0550,0

⎥⎦

⎤⎢⎣

⎡×

×4

4

10100101

= ⎥⎦

⎤⎢⎣

⎡050.5100.6100.5500.5

S–1 = ⎥⎦

⎤⎢⎣

⎡−

−49,1629,18

29,1514,15 × 10–4

Sensibilidades y selectividades:

β1 = ⎥⎦

⎤⎢⎣

⎡××−

−

−

3

3

10529,110514,1

SEN1 = || β1 ||–1 = 464 SEL1 = 464 / (5.5002 + 6.1002) = 0,056

β2 = ⎥⎦

⎤⎢⎣

⎡×−

×−

−

3

3

10649,110829,1

SEN2 = || β2 ||–1 = 406 SEL2 = 406 / (5.1002 + 5.0502) = 0,056 Realizando este mismo análisis a todas las posibles combinaciones de dos longitudes de

onda se obtienen los siguientes resultados: Combinación SEL1 SEL2

1 y 2 0,056 0,056 1 y 3 0,029 0,029 1 y 4 0,007 0,007 1 y 5 0,028 0,028 2 y 3 0,084 0,084 2 y 4 0,059 0,059 2 y 5 0,027 0,027 3 y 4 0,022 0,022 3 y 5 0,056 0,056 4 y 5 0,034 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.

91

Respuesta al ejercicio propuesto en la clase 4 1) Los resultados provistos por el programa son los siguientes: Figura con espectros de calibración:

Figura con sensibilidades y coeficientes de regresión:

92

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:

Cifra de mérito Analito 1 Analito 2 Sensibilidad Sensibilidad analíticaa Selectividad

0,50 A × ppm–1

100 ppm–1 0,62

0,76 A × ppm–1 152 ppm–1

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:

Muestra de prueba

Analito 1 Analito 2 Residuo espectral

Predicho Predicho 1 2 3

2,00(1) 1,02(1) 4,05(1)

2,01(1) 1,06(1) 3,91(1)

0,009 0,01 0,007

La rutina de MATLAB provee la gráfica de los espectros de prueba:

93

Respuesta al ejercicio propuesto en la clase 5 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):

Muestra Concentración

nominal Concentración

predicha Muestra Concentración

nominal Concentración

predicha 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

1.25 1.25 1.50 1.50 1.50 1.50 1.50 1.50 1.50 1.75 1.75 1.75 2.00 2.00 2.00 2.25 2.25 2.50 2.50 2.50 2.75 2.75 2.75 3.00 3.00 3.00 3.50 3.50 3.50

1.08 1.29 1.43 1.45 1.35 1.37 1.38 1.51 1.39 1.77 1.73 1.83 2.03 2.08 1.87 2.14 2.29 2.63 2.49 2.37 2.75 2.76 2.75 2.96 3.00 2.86 3.60 3.58 3.36

30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

3.50 3.75 3.75 4.00 1.00 1.00 1.00 1.00 0.00 0.00 0.00 0.00 0.60 0.60 0.60 0.60 2.00 2.00 2.00 2.00 0.40 0.40 0.40 0.40 0.80 0.80 0.80 0.80

3.65 3.79 3.52 3.83 1.01 1.02 1.00 1.00 0.99 0.00 –0.01 –0.03 0.62 0.67 0.67 0.67 0.21 0.19 0.17 0.22 0.38 0.39 0.41 0.36 0.81 0.82 0.78 0.79

Con estos datos se puede utilizar el programa EJCR.M para evaluar la exactitud, usando

el método OLS:

94

95

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. 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.

10. 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.

11. 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.

12. 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.

13. G. D. Christian, Analytical Chemistry, 6a. Edición, Wiley, New York, 2003, Capítulo 16.

14. 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.

15. 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

16. 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.

17. R. G. Brereton, Chemometrics. Data Analysis for the Laboratory and Chemical Plant, Wiley, Chichester, 2003, Capítulo 5.

18. http://www.chm.bris.ac.uk/org/chemometrics/pubs/chemweb.html 19. 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.

20. 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.

96

21. Las imágenes están tomadas de la página web www.cc.gatech.edu/

people/home/adjacent/. 22. 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.

23. O. Svensson, T. Kourti y J. F. MacGregor, An investigation of orthogonal signal correction algorithms and their characteristics, J. Chemometrics, 2002, 16, 176.

24. S. Wold, H. Antti, F. Lindgren y J. Öhman, Orthogonal signal correction of near-infrared spectra, Chemom. Intell. Lab. Syst. 1998, 44, 175.

25. T. Fearn, On orthogonal signal correction, Chemom. Intell. Lab. Syst. 2000, 50, 47. 26. 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.

27. El programa MVC1 (Multivariate Calibration 1) puede obtenerse libremente en www.chemometry.com

28. 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.