Espectro de Covarianzas (longitudes de onda - parámetro)

El espectro de covarianza nos indica como dos medidas varian conjuntamente. Ver esto de una manera gráfica nos ayuda mejor a encontrar las bandas de absorción de los analitos de interes.

En el siguiente gráfico están conjuntamente los espectros de correlación y de covarianza respecto al parámetro "Sustancia Seca" (100 - Humedad). El espectro de covarianza (esta autoescalado para verlo mejor comparativamente respecto al de correlación) nos marca de manera más nítida  la banda combinacional de absorción del agua que esta a 1940 nm que el de correlación, también el primer sobretono por la zona de 1450 nm.

MLR: ¿Porque seleccionar 2072 nm?

En la entrada anterior optamos por seleccionar la longitud de onda de 2072 nm:

Se ordenan las muestras por valor de "Índice de Hidroxilo", restando el espectro de la de mayor concentración, de la de menor contenido.

El espectro "diferencia", lo aumentamos multiplicándolo por 5, para ver mejor las bandas.

Se observa que la mayor diferencia esta a una longitud de onda aproximada de 2072.

Por otra parte, en las tablas de absorción NIR para los ROH, en los que se incluye el índice de hidroxilo, las bandas de absorción están entorno a esta longitud de onda.

http://es.wikipedia.org/wiki/Grupo_hidroxilo

MLR: Espectros de Correlación y Sensibilidad

En la gráfica se representan los espectros extremos en valor cuantitativo para un determinado constituyente (índice de hidroxilo) en rojo (mínimo) y azul (máximo), entre ellos estarían el resto.

Los espectros están es 2ª derivada, ´por lo que el valor en el eje Y disminuye al aumentar la concentración, por lo que tenemos que tenerlo en cuenta a la hora de interpretar la correlación y sensibilidad y seleccionar la banda de absorción. Tenemos por tanto que buscar valores de correlación  (en marrón) negativos y lo mas próximos a -1.

También que la sensibilidad (en naranja) sea lo mas próxima a cero, lo que hará que la ecuación sea más insensible a pequeños cambios.

En este caso se optó por seleccionar la longitud de onda de 2072 nm, pudiendo haber elegido otras zonas:

1416 a 1450 nm, 1548-1578 nm o la de 2056 a 2102 nm.

Se recomienda leer bibliografía adecuada para seleccionar la primera longitud de onda en una calibración MLR, seleccionando después longitudes de onda (por adición) de otros absorbentes que interfieran con el analito de interés o (por división) longitudes de onda que eliminen el efecto multiplicativo a causa de fuentes de variación como el tamaño de partícula.

Software usado (VISION - FOSS NIRSystems)

Ecuaciones con variables indicadoras

En la entrada anterior hemos visto un posible caso para aplicar una ecuación con variables indicadoras, para ello unimos las dos bases de datos (X matriz de espectros) e (Y constituyente: Humedad). Ahora generamos un nuevo constituyente a modo de variable indicadora “Ind1”. En este constituyente damos valores “0” a los espectros con dato de laboratorio del NIR A y “1” a los espectros con datos de laboratorio del NIR B. En este caso el software utilizado (Win ISI) nos permite hacer la calibración ajustada  al primer grupo (valor 0), al segundo (valor 1) o a la media.

En azul la ecuación con varianle indicadoras, desarrollada con los mismos tratamientos matemáticos.

El estadístico SECV (Cross Validation Error), se realizó con la opción “leave one out”.

Corrección de Bias

En la primera gráfica se observa una recta de regresión desarrollada con un equipo (NIR A) para humedad en harina de ave (ESP), y en la siguiente los estadísticos de otro conjunto de validación del mismo producto, pero pertenecientes a una población totalmente distinta (U.S.A), y pasados por un equipo del mismo modelo, pero distinto número de serie  (NIR B),y los datos de laboratorio, por supuesto de otro diferente.

Obviamente observamos un Bias, y podíamos actuar de diferentes formas:

a) Poner la ecuación del equipo A en el equipo ajustando el Bias.

b) Mezclar los datos de los dos y desarrollar una ecuación para ambos.

c) Estandarizar el NIR B respecto al A.

d) Usar ficheros de repetibilidad

e) Hacer una ecuación con variables indicadoras.

Por supuesto, depende de las circunstancias, hay que tener en cuenta que los equipos estan a mucha distancia y realizar la estandarización o la calibración con los ficheros de repetibilidad resulta complicado. De una manera provisional podríamos aplicar la opción "a", pero la "e" me parece interesante, pues podemos generar una ecuación adaptada al NIR A que aplique la nueva variabilidad de la base de datos del NIR B, o viceversa, o ir a un término medio.

Las ecuaciones con variables indicadoras son bastante desconocidas, pero pienso que este es un buen caso para aplicarlas.

¿Que opináis vosotros?,¿Añadir vuestros comentarios?

¿Local or Global?

Recuerdo un interesante artículo de Tom Fearn con este título.

En estas gráficas se ven los estadísticos de validación con un conjunto 230 muestras de semilla de girasol (sunflower seed) elegidas al azar, usando una calibración LOCAL (gráfico superior) y una calibración GLOBAL PLS (gráfico inferior) para el constituyente "humedad" (moisture).

Se utilizaron los mismos tratamientos matemáticos para ambas. (Software utilizado Win ISI 4).

La validación muestra como las muestras de bajo contenido en humedad, así como las de alto contenido se ajustan mejor a la recta de regresion, permitiendo al modelo extrapolar mejor, que la ecuación global. Se necesitan una cantidad importante de muestras y que estas cubran toda la variabilidad que se nos presentará en rutina. Si una muestra analizada no está representada por la base de datos LOCAL para un constituyente concreto, no obtendremos resultado de predicción para dicho constituyente.

Las ecuaciones locales se desarrollan en "Win ISI" con ficheros RED. Accede a "Learning Win ISI" para mas información sobre creación de ficheros RED.

Distancia Euclidiana_001

Hemos visto como se representan los vectores en un espacio multidimensional (un punto). Ahora queremos conocer la distancia entre esos dos puntos. Una de los cálculos de distancia mas usados es el de la "distancia euclidiana".  Se calcula usando el teorema de Pitágoras.

En el caso de dos puntos (P1 y P2) representados en un plano de ejes X e Y, con coordenadas (X1 Y1) y (X2 Y2), la distancia euclidiana la calculamos como:

d(P1 P2) = SQRT[(X1-X2)^2 - (Y1-Y2)^2].

Si aplicamos este concepto a un espacio de N dimensiones, para calcular la distancia entre dos vectores A y B:

Distancia Euclidiana = SQRT (sum (A - B)^2)
La distancia euclidiana da el mismo peso a todas las direcciones del espacio que representa y se podría representar como un círculo entorno a un punto.

En espectroscopía este cálculo se utiliza en ocasiones para saber si los espectros son redundantes, al estar muy próximos entre sí, dejando unicamente uno para representar a los próximos.

Uno de los problemas de la distancia euclidiana es que no tiene en cuenta si las variables estan correlacionadas entre si, que es uno de los problemas que tenemos  en la espectroscopía NIR.
Una alternativa a la distancia euclidiana en NIR es la distancia de Mahalanobis en el espacio de los componentes principales, en este caso la representación es por medio de una elipse, y no existe el problema de la correlación entre las variables como en la distancia euclidiana..

PCA ( 3ª Parte )

La matriz “X”,  la descomponemos en dos nuevas matrices (una de loadings “P” y otra de scores “T”), quedándonos un residual de varianza no explicada representado por la matriz “E”.

La forma de hacerlo es mediante la fórmula:

                                                     X = T.P^t + E

Para multiplicar matrices debe de cumplirse la propiedad que el nº de columnas de una sea igual al nº de filas de la otra.

Hemos hecho la transpuesta de P, para que se pueda  dar la condición para hacer la multiplicación de estas matrices.

La matriz T, describe como sus variables se relacionan con las variables originales de X, por otra  P define las combinaciones lineales.

PCA (2ª Parte)

Este espectro en el rango de 400 a 2500 nanómetros, con datos de log(1/R) cada 2 nm,

dispone de 1050 variables, pero lo acotaremos para mayor comodidad al rango de 1100 a

2500 nm, lo que representan 700 variables. En este rango es donde están las bandas de los

analitos de interés (proteína, grasa y humedad). El espectro acotado queda así: 

La matriz de datos original consta de 66 filas (equivalen al número de muestras) y de 700
columnas (nº de longitudes de onda). El espectro en negro es el de desviación estándar y nos
da una idea de donde esta la variabilidad en este conjunto de datos.     

Ahora vamos a ver el efecto del centrado quitando a cada espectro, el espectro medio:

 No están todas las muestras, pero las representadas, nos dan una idea mejor de donde esta

la variabilidad en el espectro, y vemos que coincide con el espectro de desviación

estándar anterior.
Estas muestras se pueden ubicar en un espacio multidimensional de tantas dimensiones

como variables (longitudes de onda), como un punto. La muestra promedio de todas esta 

en el centroide de dicho espacio. Como hemos comentado antes, este espacio  tiene

correlacionadas las variables entre sí, siendo necesario encontrar nuevas variables no correlacionadas y que sean combinaciones lineales de todas las demás (con mayor o

menor peso). Afortunadamente los programas quimiométricos hacen estas funciones,

quedando en  nuestro caso del siguiente modo para los tres primeros componentes

principales  de los diez que selecciona:

 Podemos elegir las combinaciones que queramos de componentes principales. En este
caso es el PC1, PC2 y PC3, pero podemos seleccionar el 1-2-7, el 1-4-8,....,etc.

Para representar el 100% de la varianza, hubiésemos necesitado las 700 dimensiones, pero
en este caso con diez componentes principales representamos el 99,5%.
Hemos dicho que los componentes principales son combinaciones lineales de las
variables originales y que en el caso del primero representa la mayor fuente de varianza
(en este caso el 82,6%).
Veamos en la siguiente figura como es este primer componente principal:

 Se aprecian bandas importantes en entorno a 1940 y 1450 nm (bandas de humedad),
también aparecen otras bandas, de las que hablaremos mas adelante. Todo parece indicar
que este primer componente principal esta representando como parte importante las
bandas de humedad así como de algún otro parámetro.
Comentamos anteriormente que podemos representar estos componentes principales,
los unos respecto a los otros formando planos.
Vamos a mirar el plano formado por el PC1 y PC2:

Se han dividido las muestras en 3 tercios acorde a su valor de humedad: Las rojas
representan al tercio con menor humedad, las blancas al intermedio y las azules a las de
mayor humedad. Al proyectar las muestras sobre el eje del PC1, se observa que hay
una gran correlación (0,9) respecto a los valores de humedad.Las proyecciones sobre el
segundo tienen ya muy poca correlación con la humedad.
El primero también tiene cierta correlación con la grasa y la proteína (cerca de 0,3), esto lo
justifican las otras bandas que aparecen en el espectro del primer componente principal,
asociadas a estos dos constituyentes

En los softwares quimiométricos, a partir de la matriz original “X”, creamos dos nuevas
matrices a las que llamaremos: ”P”  (loadings o cargas) y “T” (scores o puntuaciones).
Estas dos nuevas matrices combinadas en una determinada fórmula  se acercarán en
mayor o menor medida a dar un resultado similar a “X”. Todo dependerá de la matriz de
residuales “E”, que representa a la varianza no explicada.

En el caso de que existan tantos loadings como variables en X, se explicará toda la varianza
(100%), siendo por tanto el residual de cero.

En el caso del software Win ISI, la matriz de loadings es el fichero PCA y la matriz de
scores el fichero LIB. Ambos ficheros deben de estar en la misma carpeta y
conjuntamente representan un cubo entorno a un centroide en el que confluyen los
ejes que representan a los loadings.

Como es lógico solo se podrán representar gráficamente en “3D”, pero podemos hacerlo
con una combinación de diferentes loadings. Un fichero LIB solo lo podremos
representar gráficamente, si en su mismo proyecto está el fichero “PCA”.

El fichero PCA nos aporta los coeficientes de cada longitud de onda respecto a los
distintos componentes principales (PC1, PC2,…..). Por tanto donde se produzcan las
mayores fuentes de varianza, sus longitudes de onda tendrán altos coeficientes en los
primeros componentes principales. En muchas ocasiones estas fuentes de varianza
están relacionadas con características físicas (tamaño de partícula, color, viscosidad,….),
o con determinados analitos incluidos en la muestra y que tienen un amplio rango de
concentración y presentan bandas de absorción en la región NIR (humedad, grasa, proteína,….).

Cada componente principal no representa únicamente a un analito, puede representar a
varios en mayor o menor medida. Para ello es importante fijarse en las los coeficientes en
todas las longitudes de onda.

Puede que la información, en lo que concierne a ciertos analítos, este en los primeros

componentes principales o en los últimos, por lo que decidir sobre el número de componentes principales a seleccionar  será uno de los temas a tratar

PCA (1ª Parte)

v  Se generará una base de datos con dos matrices: "Y" datos de Lab y "X" datos espectrales. (Una columna por parámetro y por longitud de onda).

v  Que correlacionen los datos espectrales con el dato de referencia es esencial para obtener un buen modelo predictivo.

v  La matriz espectral "X",la podemos representar en un espacio multidimensional (tantas dimensiones como long. de onda), como un punto.Se trata de una matriz NxM, donde N es el número de muestras y M el número de longitudes de onda.

v  Al representar los espectros de un mismo producto en este espacio, aparecerá una nube de puntos. En espectroscopía NIR tiene forma elipsoidal.

v  Cada espectro se representa como un vector que parte del centro del espacio multidimensional hasta el punto que lo representa en la nube.

v  CENTRADO: Es como representar multidimensionalmente los espectros restándole a cada uno el espectro medio.

v  Los espectros contienen información a muchas long. de onda, lo que complica la interpretación de estas nubes espectrales. Se debe comprimir.

v  Ademas, las columnas de la matriz X están correlacionadas, debido a que contienen información de sobretonos de las bandas fundamentales (sobretonos 1, 2),y a que la información de las longitudes de onda contiguas es similar. Todo esto hace mas necesaria la compresión de la matriz X. Hay que buscar nuevas variables no correlacionadas.

v  Desarrollaremos en los próximos puntos la manera de comprimir esta información por "Análisis de Componentes Principales", a los que denominaremos PCA.

v  Imaginemos una nube de puntos en 3D. Del centro de la nube parten los 3 ejes X1,X2,X3 correspondientes a 3 longitudes de onda (centrado).

v  Trazaremos una nueva variable PC1 combinación lineal de X1,X2,X3, que pase por el centroide en la dirección de la máxima varianza.

v  Este PC1 es un vector que explicará una determinada varianza (autovalor de PC1), lo que queda por explicar es la varianza residual.

v  La varianza explicada por PC1 (autovalor) se calcula sobre la proyección de los puntos de la nube sobre  PC1 y respecto al centroide.(Como sabemos, la varianza es la desviación estándar al cuadrado). Se puede expresar en %.(Ej.: 62%). No habrá ningún otro vector que explique más varianza que él. La varianza residual será el resto hasta el 100%, es decir el 38%.

v  El valor de esta variable para cualquier muestra en particular, se calcula proyectando el punto que la representa sobre este eje (perpendicularmente), siendo esta la distancia mas corta posible entre el punto y el nuevo eje. Al valor de este punto proyectado se le denomina "puntuación" o "score".

v  Se trazará un segundo componente principal PC2, ortogonal con el primero, que pase a su vez por el centroide y en la dirección que explique la máxima variabilidad residual que dejo el primer componente .Este componente PC2 formará un plano con PC1.Hay que tener en cuenta que podrémos trazar como máximo tantos componentes principales como variables originales. Se proyectarán de igual modo los puntos sobre este nuevo eje, calculando los “scores” y el autovalor.

v  Otro nuevo componente PC3, ortogonal con los otros dos, explicara el resto de la varianza que dejaron los otros dos, llegando entre los tres a explicar el 100%, ya que disponíamos de 3 variables originales. Otro nuevo componente PC3, ortogonal con los otros dos, explicara el resto de la varianza que dejaron los otros dos, llegando entre los tres a explicar el 100%, ya que disponíamos de 3 variables originales. Este nuevo componente formará junto a PC1 y PC2 dos nuevos planos.

v  La situación queda de modo que disponemos de un nuevo espacio dimensional, cuyo centroide coincide con el del espacio original una vez centrado. Se trata por tanto de una rotación de los ejes en nuevas direcciones, que van explicando de una manera ordenada (de mayor a menor) la varianza. Conocemos los ángulos que forman con los componentes originales pudiendo expresarlos como combinaciones lineales de ellos.

v  MAPAS DE “SCORES”: La representación de estos planos, en forma gráfica es lo que se conoce como “Mapas de Scores”, y nos ayudan a interpretar de una manera muy útil los datos. Hablaremos de ellos más adelante.

v  Los scores se representan en una nueva Matriz a la que llamaremos T, y en ella estarán representadas las muestras con sus valores de scores en las filas y las nuevas variables (componentes principales) en las columnas. Como se ha comentado anteriormente estas nuevas variables son ortogonales entre sí.

v  VECTOR DE “SCORES”: Se conoce como vector de scores a cada una de las columnas de la matriz de scores, que representan las proyecciones de las muestras sobre cada uno de  los componentes principales. Estos valores son en los que se fundamentan los cálculos de varianza para calcular los autovalores.
      
     Hacemos un paréntesis para aplicar lo hasta ahora visto a la espectroscopía NIR.