Linealidad: método del Lack-of-fit

Lo prometido es deuda y he aquí una de las técnicas prometidas en la entrada anterior para el estudio de linealidad. Aunque realmente me gusta dejar claro que esta prueba tiene como cometido principalmente comprobar cuanto error puede aportar a las predicciones la falta de ajuste del modelo. Para hablar de linealidad creo que es preferible comprobar a cada nivel de concentración, porque un modelo que puede parecer globalmente lineal puede no serlo para ciertos niveles de concentración. Pero esa es mi opinión, y por suerte la de muchos.

Debéis disculparme que para esta entrada me ponga en plan técnico más que divulgativo, pero esta entrada tiene fines docentes, y prefiero escribirlo de esta forma.

Partimos de un conjunto de n_i puntos de calibración (x_i, y_i) que presentan una relación lineal aparente, donde se considera que los cada valor de x_i está exento de error  y los valores de y_i están sujetos a errores de medida pero son homocedásticos. Los valores x_i deben aparecer replicados, con lo que estarán agrupados en j (1 a n_j) niveles con k (1 a n_k) replicados en cada nivel, con lo que n_i = n_j·n_k. Los datos pueden ajustarse a un modelo con c parámetros de ajuste (en regresión lineal c = 2) dando lugar a una ecuación del tipo (1). Siendo y la variable dependiente,  x la variable independiente, b₁ la pendiente de la recta de regresión y b₀ la ordenada en el origen. La función se obtiene fácilmente mediante el método de mínimos cuadrados minimizando la suma de cuadrados de residuales (2), que considerando la replicación se puede escribir como (3), con grados de libertad  (4).

La prueba  F de falta de ajuste o de Lack-of-fit se basa en un análisis de la varianza (ANOVA) de residuales en el que la suma de residuales, SSE, se descompone en componentes de falta de ajuste. SLF, y error puro SPE, (5) y (6) comprobando la influencia del primero en el error de los residuales. En (7) y (8) se muestra la descomposición de los grados de libertad de cada término.

Los cuadrados medios o varianzas pueden obtenerse fácilmente dividiendo las sumas de cuadrados entre sus correspondientes grados de libertad. De este modo se puede tener la siguiente tabla de ANOVA de falta de ajuste.

Tabla de ANOVA de falta de ajuste (Lack-of-fit)

La prueba de falta de ajuste consiste en calcular un valor F como el cociente entre la varianza (o cuadrado medio) de falta de ajuste y la varianza de error puro.

La hipótesis nula ha de ser entonces que no existe falta de ajuste, el modelo ajusta de forma adecuada a los datos. Si el valor de F calculado es menor que un F tabulado para α = 0.05 (95% de nivel de confianza), n_j-c grados de libertad para el numerador y n_i-n_j grados de libertad para el denominador, se acepta la hipótesis nula. En caso contrario se rechaza y se dice que existe falta de ajuste.

Aunque os parezca tedioso, el cálculo es muy simple y se puede hacer en hoja de cálculo:

1) Obtener SSE a partir de cálculos de regresión o según la fórmula tabla.
2) Obtener SPE según la fórmula de la tabla.
3) Obtener SLF despejando en (5).
4) Dividir entre los grados de libertad indicados en la tabla para obtener MSLF y MSPE.
5) Obtener el F calculado según (9)
6) Comparar con el valor tabulado de F.

Una ventaja de este método es que está implementado en muchos entornos y software como Minitab, R, Statistica, entre otros.

Espero sea de utilidad. Más adelante hablaré del método %RE-plot, propuesto por nosotros, como ya he indicado.

CONTINUARA

Sobre linealidad en el rango de calibración analítica

Desde hace unos días he mencionado varias veces un trabajo sobre linealidad en el rango de calibración analítica [1]. Creo que es necesario explicar algunas cuestiones sobre este parámetro de calidad de un método. Una comparación detallada de los distintos métodos para evaluar la linealidad se puede encontrar en un trabajo previo [2] de uno de los autores, Francisco Raposo, y, de hecho, fue el punto de partida para despertar mi interés en este tema. Tengo que agradecer al Dr. Raposo por las discusiones previas que me llevó a plantearme este trabajo y fructificó en una colaboración, por haber insistido en intentarlo en una gran revista como Talanta y por haberse encargado él de la última parte de la revisión (muy luchada) en estos momentos en los que me daba por vencido. El tiempo nos dirá si nuestra propuesta es útil o no.

Cuando somos estudiantes y comenzamos a realizar ajustes lineales en nuestras prácticas de laboratorio nos gusta mucho que el coeficiente de correlación (r) o el de determinación (r^2) sean próximos a la unidad. Si eso es así, nuestra respuesta es lineal... Falso, simplemente si nos aceramos a la unidad podemos asegurar que el modelo empleado para ajustarse a nuestros datos es bueno, pero no que sea el mejor. Este parámetro informa sobre como es de pequeña la suma de cuadrados de residuales, pero no indica como se distribuyen los mismos. Por eso no se recomienda el empleo del coeficiente de correlación o determinación como sinónimo de linealidad. Habría que observar la distribución de los residuales.

Este valor de r^2 puede darse por válido en una práctica, pero sería un error. Si se observa con detalle, los puntos se distribuyen de forma no lineal

Otro parámetro muy usado, que yo he usado mucho para no recurrir al valor de r, es el de linealidad on-line [3], que se calcula como 100 menos la desviación estándar relativa de la pendiente. La desviación estándar de la pendiente aumenta cuando aumenta el error estándar de residuales, con lo que estamos en una caso parecido a hablar de coeficientes de correlación, es decir, te da cuenta del porcentaje de ajuste de acuerdo al error de la pendiente, pero no te dice como se distribuyen los puntos alrededor de la recta de calibración. 

Fórmula de la linealidad on-line

¿Cuál es la solución? Como siempre digo a mis alumnos, lo primero que ha de hacerse cuando se tienen datos en el laboratorio es pintarlos en un papel. Hay que hacer una gráfica, siendo a veces la inspección visual lo que da idea, si los datos se distribuyen bien alrededor de la función ajustada, de si el modelo elegido es el adecuado. Supongamos el modelo lineal, a veces es difícil ver la distribución sobre la misma recta y es preferible usar un gráfico de residuales. Este tipo de gráficos es muy interesante, porque puede indicar que la función elegida no es adecuada simplemente con observar la distribución de los residuales. También puede indicar que la distribución no es homocedástica, lo que implicaría el uso  de un modelo ponderado. Pero si tenemos un gráfico de residuales ¿cómo establecemos el límite que pueda indicar  que un punto está fuera de la tendencia de los demás. Para eso es interesante trabajar con residuos estandarizados. Hay muchas formas de estandarizar, pero nosotros proponemos los residuos Studentizados. Si un residuo supera el valor de 1.96 (aproximamos a 2), se considera que el valor es sospechoso se denomina valor extremo (outlier) si supera el 3. Para aceptar un modelo no debe haber valores sospechosos a lo largo del rango de calibración.

Gráficos de residuales Studentizados correspondiente a la recta anterior. Aunque todos están dentro de los límites se observa una distribución de los mismos que indican la idoneidad de un modelo no lineal

Otra opción es el gráfico de linealidad o de factor de respuesta, que se hizo conocido gracias a Huber [4]. Desde mi punto de vista, este método da cuenta del la variación del factor de respuesta (señal debida al analito dividida entre su cantidad o concentración) a lo largo del rango de calibración, pero poco más. Huber establece unos límites del 5% por encima y por debajo de la mediana de los factores de respuesta y establece la linealidad en el rango en que los puntos se mantienen entre esos límites. Esto funciona generalmente bien en métodos cromatográficos, pero no es así en métodos espectroscópicos. Además no tiene en cuenta que las desviaciones de los factores de respuesta pueden ser muy grandes si trabajamos a concentraciones muy pequeñas. no obstante es un método muy útil para observar cambios en la respuesta del método. 

Se puede observar que el gráfico de factor de respuesta (C) falla en este caso 

El método del lack-of-fit (falta de ajuste)  no lo incluimos en el trabajo. Podéis ver una explicación del mismo en aquellos apuntes de Excel que compartí hace años. Aunque esta entrada en el blog es una declaración de intenciones y quiero trabajar en un par de entradas más sobre linealidad. En la primera describiré exclusivamente el método de lack-of-fit y en la segunda el procedimiento que hemos denominado %RE-plot, el gráfico de errores relativos recalculados.

CONTINUARA...

Referencias

[1] J. M. Jurado, A. Alcázar, R. Muñiz-Valencia, S. G. Ceballos-Magaña, F. Raposo, Some practical considerations for linearity assessment of calibration curves as function of concentration levels according to the fitness-for-purpose approach, Talanta, 2017, 172, 221-229.

[2] F. Raposo, Evaluation of analytical calibration based on least-squares linear regression for instrumental techniques: A tutorial review, TrAC Trends in Analytical Chemistry, 2016, 77, 167-185.

[3] L. Cuadros Rodríguez , A. M. García-Campaña , C. Jiménez, M Román, Estimation of performance characteristics of an analytical method using the data set of the calibration experiment, Analytical Letters, 1993, 26,1243-1258.

[4] L. Huber, Validation of analytical methods: review and strategy, LC-GC Europe,1998, 11,  96-105.

Sensibilidad de un método analítico y límite de detección

La sensibilidad de un método analítico se define de acuerdo a la IUPAC como la pendiente de la curva de calibración, tratándose de una característica del método que depende sólo del proceso de medida. Así definida, la sensibilidad no es otra cosa que el factor de respuesta, o lo que es lo mismo, el cociente entre la variación de señal asociada a un determinado analito y la variación de su concentración o cantidad (si tuviésemos un solo patrón simplemente es el cociente entre señal y concentración del mismo). La definición es muy clara pero, sin embargo, cuando se presentan las características de un método analítico en publicaciones académicas y científicas se opta normalmente por usar el límite de detección. Pocos autores hablan de sensibilidad dando el valor de la pendiente de calibración, y aquí quiero explicar el motivo de ello.  

El límite de detección (LOD, limit of detection), expresado como la cantidad o concentración, proviene de  la señal más pequeña que puede detectarse con razonable certeza en un determinado procedimiento analítico (IUPAC). Según esto, se trata de la cantidad asociada a la mínima señal que pueda atribuirse al analito, es decir, que sea distinguible de la señal del blanco de medida.  Esa señal mínima se suele definir como la señal del blanco más tres veces la desviación estándar del blanco:

Si en esta expresión se usase la señal del blanco más diez veces su desviación estándar estaríamos defiendo la señal correspondiente al límite de cuantificación (LOQ, limit of quantification), la mínima cantidad cuantificable. Pero en esta entrada nos centraremos en el límite de detección.

Imaginemos que tenemos una recta de calibrado externo del tipo Y = b·X + a, donde Y es la señal correspondiente a una concentración de analito X, b es la pendiente y a es la ordenada en el origen. Si se ha corregido la señal del blanco de los patrones de calibración en la recta anterior, la ordenada en el origen tendrá un valor de cero (a = 0). Ahora debemos sustituir la señal del límite de detección en la ecuación de la recta, pero será la señal corregida, es decir la señal correspondiente al límite de detección menos la señal del blanco:

Esta expresión para el límite de detección es muy común en técnicas de espectroscopia atómica, en las que a baja concentración puede asumirse una ordenada en el origen prácticamente nula. existen otras formas de expresar los límites de detección, pero en todo caso siempre se trata de un parámetro con unidades de señal (en nuestro caso la desviación estándar del blanco) dividido entre el factor de respuesta (la pendiente de la recta de calibrado), cuyas unidades son de señal dividido por concentración. De esta forma siempre se tiene un limite de detección con unidades similares a las de los patrones de calibrado. En otras ocasiones se usa, en lugar de la desviación estándar de la señal del blanco, la desviación estándar de residuales, la de la ordenada en el origen o la desviación estándar de la linea base. Lo importante es que el LOD quede expresado en unidades de concentración, y tras los siguientes ejemplos veremos el porqué.  

a) Comparación de la sensibilidad de una técnica para distintos analitos

En esta comparación voy a usar datos publicados en el trabajo Direct determination of copper, lead and cadmium in aniseed spirits by electrothermal atomic absorption spectrometry, publicado en Food Chemistry en el año 2007, como fruto de mi tesis doctoral: Caracterización analítica de aguardientes anisados.  Vamos a considerar dos ejemplos, un caso donde la diferencia de sensibilidad es obvia y otro en que parece obvia, pero no lo es tanto.

Ejemplo 1. Determinación de cobre y cadmio mediante espectroscopia de absorción atómica con atomización electrotérmica.

El  método propuesto para ambos elementos incluyen una dilución de la muestra en una mezcla de agua / etanol / ácido nítrico (58:40:2) y determinación directa tras un programa optimizado de secado, mineralización y atomización, previa adición del modificador de matriz adecuado. Por lo tanto el calibrado se realiza en una mezcla de los tres componentes en esa misma proporción. 

En el caso del cobre se añade una disolución de nitrato de paladio como modificador de matriz, se seca en dos etapas a 80 y 110 ºC, se mineraliza a 1300 ºC y se atomiza y mide a 2300 ºC. Para el cadmio se emplea una mezcla d nitrato de magnesio y paladio, secando de igual modo que para el cobre. En este caso la mineralización ocurre a 700 ºC y la atomización a 1500 ºC. Las rectas de calibrado tienen las siguientes ecuaciones:

Cu: Y = 0.012·X - 0.001

Cd: Y=0.11·X - 0.004

Siendo Y la señal obtenida y X la concentración en µg/L. Como puede observarse la pendiente de calibración del Cd es casi 10 veces mayor que la del Cu, luego la técnica ETAAS (de electrothermal atomic absorption spectroscopy) es casi diez veces más sensible para Cd que para Cu en estas condiciones de trabajo. El límite de detección debe ser una diez veces más bajo para Cd que para Cu. Los valores calculados como se explica al comienzo de esta entrada  son 0.04 µg/L para cadmio y 0.6 µg/L para cobre, algo más de diez veces menor.  

Ejemplo 2. Determinación de cobre y plomo mediante espectroscopia de absorción atómica con atomización electrotérmica.

La determinación de cobre ya se ha explicado, para el plomo se emplea una mezcla de nitratos de magnesio y paladio, con el mismo secado que para los otros dos elementos, mineralización a 900 ºC y atomización a 1800 ºC. Si se comparan las rectas de calibrado:

Cu: Y = 0.012·X - 0.001

Pb: Y = 0.07·X + 0.000

A primera vista la técnica parece casi el doble más sensible para cobre que para plomo. Pero, ¿que hemos olvidado? Hemos olvidado el error asociado a la pendiente, que esta relacionado con el ruido de fondo que afecta a toda medida. En el caso del cobre el error es de ± 0.004 (33%) y para plomo ± 0.001 (14%). Esto se debe a que la señal del fondo para cobre es más acuciada que para plomo.

Señales obtenidas para cobre y plomo, así como señal de fondo. La señal que se emplea en el calibrado es el área de pico, cuyas unidades son unidades de absorbancia por segundo (u.a.·s) 

Esto hace que cuando se obtienen los límites de detección se encuentre 0.6 µg/L para cobre y 0.7 µg/L para plomo. En realidad no son tan diferentes, puesto que la desviación estándar del blanco es mayor en el caso del cobre y esto hace que, aunque su pendiente sea mayor, no se aprecie un límite de detección mucho menor. Por ese motivo es necesario conocer el error de la pendiente si se quieren comparar sensibilidades o bien realizar el cálculo de los límites de detección.

  

b) Comparación de la sensibilidad de dos técnica diferentes

Para este ejemplo compararemos la determinación de cobre en aguardientes mediante espectroscopia de absorción atómica con atomización electrotérmica (ETAAS) y la determinación del mismo elemento en tequila mediante espectroscopia de emisión atómica de plasma acoplado inductivamente (ICP-AES, de inductively coupled plasma atomic emission spectrometry). Los datos del segundo caso han sido obtenidos de la tesis de Silvia Ceballos Magaña, titulada "Caracterización analítica de destilados de Agave Tequilana mediante técnicas de análisis multivariante". En este caso se digiere una cantidad de muestra en un microondas usando ácido nítrico como oxidante y pentóxido de divanadio como catalizador. Las rectas de calibrado tienen las siguientes ecuaciones:

Cu (ETAAS): Y = 0.012·X - 0.001

Cu (ICP-AES): Y = 2.52·X - 0.13

En principio puede parecer que la técnica ICP-AES es unas 200 veces más sensible que la ETAAS para la determinación de cobre, pero, ¿qué hemos olvidado? Hemos olvidado decir cuales son las unidades de dichas pendientes. En el caso de ETAAS, como la señal viene en u.a.·s y la concentración en µg/L, las unidades de la pendiente son u.a.·s·L/µg. En el caso del ICP-AES, las señales en el equipo de medida vienen dadas en kilocuentas (kcuentas) y las concentraciones de la recta de calibrado en mg/L, así la pendiente viene dada en kcuentas·L/mg. Obviamente no se puede comparar las pendientes con unidades diferentes. Por eso se hace necesario el cálculo de los límites de detección en las mismas unidades. 

Cuando se obtienen los límites de detección se observa que para ETAAS el LOD del cobre es de 0.6 µg/L y para ICP-AES es de 3 µg/L. Es decir, para cobre medido mediante ICP-AES el límite de detección es aproximadamente  5 veces mayor que cuando se mide mediante ETAAS. La espectroscopia de absorción atómica con atomización electrotérmica es más sensible que la espectroscopia de emisión atómica de plasma acoplado inductivamente para la determinación de cobre en estas condiciones experimentales.

Conclusión

Aunque la IUPAC define la sensibilidad como la pendiente de la recta de calibrado (o factor de respuesta), existen situaciones donde es mejor usar los límites de detección para comparar este parámetro de calidad del método. De hecho, es lo que recomiendo siempre. 

Espero que estos ejemplos hayan sido ilustrativos.

Regresión lineal múltiple con Excel: resolución de mezclas en espectroscopia molecular

Allá por julio de 2015 proponía una entrada sobre la resolución de sistemas de ecuaciones con Excel en la que se usaba como ejemplo la cuantificación en una mezcla de dos sustancias, previa medida de patrones de ambas sustancias por separado. Aquel era un ejemplo simplista, fácil de encontrar en manuales de Química Analítica. A veces la realidad es otra, y para resolver una mezcla de dos o más sustancias no basta con medir los patrones de cada una por separado para obtener unos coeficientes de absortividad molar y resolver así el sistema de ecuaciones. Generalmente se suele medir una serie de patrones, mezcla de los componentes a determinar, registrando la absorbancias a varias longitudes onda. Con estos datos se puede obtener un modelo de regresión lineal múltiple que permita relacionar mediante una función la concentración de cada analito con las absorbancias medidas a las distintas longitudes de onda y cuantificarlos así en una muestra.

El siguiente ejemplo es una simulación hecha en Excel para tres componentes (C1, C2 y C3) en concentraciones molares, midiendo la absorbancia (A1, A2, A3) a tres longitudes de onda. En la simulación se ha empleado un error aleatorio para las señales de un 2% de media. 

C1	C2	C3	A1	A2	A3
0.0075	0.0075	0.0075	1.076	1.08	0.646
0.0025	0.0075	0.0075	0.55	0.965	0.64
0.0075	0.0025	0.0075	0.981	0.615	0.591
0.0025	0.0025	0.0075	0.475	0.55	0.63
0.0075	0.0075	0.0025	1.031	0.96	0.346
0.0025	0.0075	0.0025	0.54	0.91	0.325
0.0075	0.0025	0.0025	0.936	0.55	0.316
0.0025	0.0025	0.0025	0.465	0.435	0.285
0.005	0.005	0.005	0.736	0.775	0.465
0.005	0.005	0.005	0.731	0.785	0.465

Los datos se introducen en la matriz A1:G11, incluyendo encabezados de columna y fila.

Datos de calibración para el ejemplo propuesto

 Tendremos que construir tres modelos de regresión lineal múltiple, uno por cada analito, para relacionar las absorbancias medidas (variables independientes en nuestro modelo) con las concentraciones (variables dependientes). Para ello empleamos la fórmula matricial =ESTIMACION.LINEAL(). 

En el caso de C1, seleccionamos la matriz K2:N6 e introducimos la fórmula =ESTIMACION.LINEAL() desde el menú Formulas/ Insertar función. Como valores de Conocido_y introducimos la matriz B2:B11, que se corresponde con los valores de C1. Como Conocido_x se introducen los valores para A1, A2 y A3, es decir, la matriz E2:G11. Se emplea Constante 1 (VERDADERO) si el modelo contempla un término independiente. En principio lo dejaremos así, si se quisiese obviar se introduce 0 (FALSO). En Estadística introducimos 1 para que calcule errores de los coeficientes, coeficiente de determinación, error de residuales y otros parámetros ya explicados en la entrada sobre regresión lineal en Excel, como el valor F de Fisher, los grados de libertad, suma de cuadrados de regresión y suma de cuadrados de residuales.  Si se prefiere se puede escribir directamente  =ESTIMACION.LINEAL(B2:B11,E2:G11,1,1) con la matriz K2:N6 seleccionada previamente. De cualquiera de las formas, pulsar al mismo tiempo "Ctrl+Shift+Enter".

Formulario de entrada para la función =ESTIMACION.LINEAL()

Si todo ha ido bien, la matriz K2:N6 queda rellena con una serie de valores. En la siguiente figura aparecen dichos valores con unos encabezados explicativos. C1 (en la celda J1) se refiere a la especie 1. C_A1, C_A2 y C_A3 se refiere a los coeficientes que relacionan cada absorbancia con C1. Constant se refiere al termino independiente. En la matriz K2:N2 están los valores de los coeficientes y en la matriz K3:N3 sus errores. El coeficiente de determinación aparece en K4 y el error de residuales en L4. El valor de F, grados de libertad, suma de cuadrados de regresión y suma de cuadrados de residuales aparecen en K5, L5, K6 y L6, respectivamente. Además, más abajo, se incluye el valor calculado de t (valor del coeficiente dividido entre su error) y la probabilidad p de que ese coeficiente valga cero. La explicación de esta prueba la podéis encontrar en una entrada anterior. En este ejemplo, t se calcula en valor absoluto.

Resultados del modelo de regresión lineal para la especie C1

De acuerdo a la prueba t, el coeficiente para A3 (para una probabilidad de 0.05) podría ser obviado, y el modelo recalculado con menos parámetros. Para no alargar la entrada, las pruebas eleminando coeficientes no han sido realizadas.

Algo que conviene recordar es el hecho que los coeficientes par A1, A2 y A3 aparecen en la matriz de resultados en orden inverso a como estén ordenandas las columnas en los datos.

Los resultados para C2 y C3 se realizarían con la misma función. En mi hoja de cálculo lo hice en K18:N22 y K32:N36, respectivamente.

Resultados para C2 y C3

De esta forma se tienen las ecuaciones:

C1 = 0.0103*A1- 0.0019*A2 - 0.00035*A3 + 0.0009

C2 = - 0.002*A1 + 0.0119*A2 - 0.0032*A3 +0.0006

C3 = 7.15*10^(-5)*A1 - 0.0013*A2 + 0.0164 *A3 + 0.002

Los coeficientes de correlación son de 0.996, 0.992 y 0.992 para los tres modelos. El ajuste parece ser adecuado.

Para comprobar la calidad del modelo se han simulado las absorbancias para tres muestras con concentraciones nominales conocidas de la misma forma que se hicieron los patrones. Con las ecuaciones anteriores se calcularon las concentraciones experimentales. Se calculan recuperaciones como (Valor calculado/valor nominal *100). Se observa que los valores de recuperación oscilan entre 93% y 112%, debido al error aleatorio que se le introdujo a las señales. Estos resultados podrían mejorarse si se incluyese un mayor número de medidas de absorbancias y un mayor número de patrones. Pero ese no es el objeto de esta entrada.

Muestras simuladas con sus concentraciones nominales, señales, concentraciones calculadas con el modelo y recuperaciones.

Hemos decidido calcular el modelo con la función matricial =ESTIMACION.LINEAL(), pero el mismo cálculo se podría haber hecho empleando la herramienta de Regresión del complemento Análisis de datos. Por ejemplo, para el compuesto C1:

Herramienta Regresión

Entrada de datos para el compuesto C1

 Hemos seleccionado las matrices de entrada de datos incluyendo el encabezado. Cuando se hace eso es necesario seleccionar Rótulos en el formulario. Ademas se ha seleccionado Residuos y Gráfico de residuales, por si alguien quiere analizar los mismos.

Resultados par C1

La ventaja es que, si se seleccionan los rótulos en el formulario, en la matriz de resultados queda claramente establecido que coeficiente corresponde a cada variable. Además de poder ver el ANOVA y los gráficos de residuales (y la representación de los valores reales y los estimados para los patrones, en caso de seleccionar la Curva de regresión ajustada. Otra ventaja es que se presentan los resultados de significación de cada coeficiente directamente (Probabilidad tras el valor del Estadístico t). En este caso, t no se obtiene en valor absoluto, pero la probabilidad si se obtiene para el valor positivo de t.

Resultados para C1, detalle de los coeficientes

En cuanto a los residuales, es otra ventaja el no tener que calcularlos a mano. En este caso se observa una distribución aleatoria de los mismos.

Detalle de los gráficos de residuales para C1

El mismo procedimiento podría llevarse a cabo para las concentraciones de C2 y C3. La herramienta a emplear es elección de quien realiza los cálculos. Aquí no continuaremos con ello, pero el lector puede comprobar los resultados por sí mismo.

Nota: En este ejemplo hemos generado tres modelos (uno para cada sustancia) que relacionan la concentración de la sustancia con absorbancias medidas a varias longitudes de onda (tres en este caso). Un ejemplo parecido a este, con esta misma forma de proceder se puede encontrar en Miller y Miller, 2002. Está forma de relacionar las variables facilita mucho el cálculo posterior en la muestra, al obtener directamente la concentración de cada analito mediante una función. Si se hubiesen relacionado las absorbancias con las concentraciones podríamos haber obtenido los coeficientes de absortividad molar para cada sustancia a cada longitud de onda. En ese caso, al medir cada muestra nos quedaría  un sistema con tres ecuaciones (tantas como absorbancias medidas) de las que habría que despejar las concentraciones. Esto complica el cálculo, pues primero habría que solucionar el ajuste lineal múltiple y luego el sistema de ecuaciones. Por eso parece más lógico relacionar directamente la concentración de cada sustancia con las absorbancias medidas. Ver: J. N. Miller, J. C. Miller, Estadística y Quimiometría para Química Analítica, Prentice Hall, Madrid, 2002, pp. 239-242

¿Son mis coeficientes de ajuste significativamente distintos de cero?

No es la primera vez ni será la última en la que me encuentre a científicos que incluyen el punto (0, 0) en una curva de calibración (me refiero a calibración lineal en toda la entrada). Yo siempre lo desaconsejo, pues para mí el calibrado es válido solo entre los puntos que se incluyen de forma experimental. Además, es muy común que a concentraciones bajas existan desviaciones de la supuesta linealidad del calibrado. En algunas técnicas, como la espectroscopia de absorción atómica con atomización electrotémica (ETAAS) es fácil asumir que el punto (0, 0), que se obtiene poniendo el equipo a cero cuando se mide el blanco, podría ser incluido, porque ciertamente hay buena linealidad a concentraciones muy bajas para esta técnica. Pero en fin, en el fondo es cuestión de escuelas de pensamiento...

Hoy no pretendo hablar de esto, aunque si de algo relacionado. Porque una cosa es incluir el punto (0, 0) en un calibrado cuando se ha medido el blanco y patrones de muy baja concentración (ng/mL, en el ejemplo de ETAAS), y otra es asumir ese valor sin haber comprobado lo que ocurre a concentraciones bajas. Y eso es lo que hace mucha gente cuando "obliga" a la recta de calibración a pasar por el origen de coordenadas. A veces, un valor muy distinto al cero puede ser significativamente igual al mismo, y un valor muy próximo a cero no serlo en absoluto. Aquí repasaremos el test estadístico más habitual para comprobar si un coeficiente es significativamente igual a cero, lo que puede ser utilizado para cualquier tipo de ajuste.

Una serie de datos de calibración, dos opciones de ajuste

Antes de empezar decir que estos datos son simulados, y que en un ajuste real, posiblemente la mayor variabilidad de los resultados hagan que no sea tan simple tomar decisiones. En mi opinión, tampoco es tan imperante eliminar la ordenada en el origen de una regresión lineal simple, pues la ecuación resultante es sencilla para realizar posteriores operaciones. No suelo emplear este procedimiento salvo que estuviésemos comprobando varias variables (cada una con su coeficiente) en un ajuste múltiple, o queramos eliminar algún orden superior de un polinomio. Otra advertencia es que este test es extremadamente sensible al nivel de errores aleatorios del sistema de medida, es decir, una mayor variabilidad puede eliminar un coeficiente sin necesidad y una poca variabilidad mantener un coeficiente innecesario. Pero al menos tenemos unas reglas que se pueden aplicar para tomar decisiones.

Imaginemos los siguientes datos de señal (Y) y de concentración (X). Calculemos la ecuación de la recta de mejor ajuste mediante la fórmula matricial =ESTIMACION.LINEAL(B2:B7,A2:A7,1,1). Este procedimiento se explica en la entrada del blog Cálculo de regresión en Excel 2007, que es perfectamente extrapolable a cualquier otra versión de Excel. Se observa una pendiente b=0.0244 ± 0.0001 y una ordenada en el origen a= 0.0027 ± 0.0006, con un coeficiente de determinación R^2=0.99984. La ordenada en el origen es muy pequeña, con lo que uno puede pensar en eliminarla. Pero, ¿sería correcto? Si obtuviésemos la ecuación de la recta haciendo cero la ordenada en el origen (=ESTIMACION.LINEAL(B2:B7,A2:A7,0,1)), el nuevo coeficiente de determinación sería R^2=0.99983. Casi el mismo valor, con lo que uno se piensa si merece la pena eliminar la ordenada en el origen del ajuste. 

Introducción de los datos del primer ejemplo y cálculo de la recta de mejor ajuste, con ordenada en el origen.

Pero no es esa la forma correcta de comprobarlo. Lo habitual en la mayoría de los paquetes estadísticos, y Excel no es una excepción, es mostrar los resultados con una prueba t de Student asociada que compara el valor del coeficiente con el cero (la hipótesis nula es que el valor del coeficiente es estadísticamente igual a cero). Es muy simple, porque el valor de t se obtiene dividiendo el coeficiente entre su error y se compara este valor con el t crítico para una probabilidad α y n-k grados de libertad (n es el número de puntos del calibrado y k el número de parámetros que se estiman en el modelo). En las versiones más recientes de Excel se emplea la fórmula =INV.T.2C(probabilidad,grados_de_libertad) para obtener el valor de t crítico (en versiones antiguas =DISTR.T.INV(), que aún funciona en las nuevas versiones). Esta es la forma que prefiero personalmente para comprobarlo, calcular los valores de t de los parámetros y el valor crítico, y compararlos directamente. Si el valor calculado es mayor que el crítico, se rechaza la hipótesis nula y se dice que el coeficiente es significativo. En caso contrario, el coeficiente es igual a cero, desde un punto de vista estadístico, para la probabilidad seleccionada (generalmente α=0.05).

 Aunque en la mayoría de los paquetes estadísticos no se suele calcular el valor crítico de t y compararlo directamente con el t calculado para el parámetro, sino que se calcula la probabilidad de que  t calculado sea menor que t crítica, o lo que es lo mismo, que el coeficiente sea igual a cero. En Excel se puede usar la función =DISTR.T.2C(x,grados_de_libertad) para obtener esta probabilidad, siendo x el valor de t calculado para el parámetro. En las siguientes figuras se ve como se introducen estas fórmulas en nuestro ejemplo y como quedarán los resultados.

Introducción de los datos del ejemplo para comprobar la significación de los coeficientes 

Resultados de la comprobación

Como puede verse, ambos valores de t son mayores que el valor crítico, o bien ambas probabilidades (p) son inferiores a 0.05. Es decir, los coeficientes no son significativamente iguales a cero y no se deben eliminar del modelo.

Esto mismo lo hace Excel empleando la función Regresión del complemento Análisis de datos del menú Datos. El complemento hay que activarlo en Archivo/Opciones/Complementos. Esta función se explica en  Cálculo de regresión en Excel 2007 y también se puede ver en el tutorial de ajuste en Excel publicado en la revista Educación Química. El formulario de esta función, que aparece en la siguiente figura, genera una hoja nueva en el libro de la que podemos sacar la misma información que  he indicado antes.

Formulario de la función Regresión del complemento Análisis de Datos

Resultados para la función Regresión del complemento análisis de datos.

En las celdas B17 y B18 aparecen los valores de ordenada en el origen (intercepción o intercepto) y pendiente, respectivamente. En las celdas C17 y C18 aparecen sus errores. Los valores de t calculado aparecen en las celdas D17 y D18 y la probabilidad de que el coeficiente sea igual a cero en las celdas E17 y E18. Además calcula unos límites de confianza para los coeficientes como (Coeficiente ± error del coeficiente* t calculado) para un nivel de confianza dado. Como se observa, los resultados son similares a los que se han obtenido mediante fórmulas.

El segundo ejemplo lo dejo a modo de ejercicio. Es curioso como ahora que tenemos una ordenada en el origen de 1.3 ± 1.0, el coeficiente es estadísticamente igual a cero. Como he dicho, todo depende de los errores del parámetro...

Segundo ejemplo, para que lo haga aquel que esté interesado

Un tutorial sobre ajuste de datos con Excel

En muchas entradas del blog he ido introduciendo problemas de ajuste de datos con Excel. Con lo cual he ido sacando algunas ideas que puestas en práctica y con ayuda de mis colegas ha quedado reflejado en una publicación en la revista Educación Química, de la Universidad Nacional Autónoma de México. 

En el tutorial se revisan las distintas opciones que presenta Excel para llevar a cabo el ajuste de funciones a datos, como son la obtención gráfica de la ecuación de ajuste, la herramienta de regresión del menú de datos, las funciones de estimación lineal y logarítmica y el empleo de Solver, incluyendo la estimación de errores.

El trabajo es de acceso libre, así que podéis verlo en el enlace Ajustando datos químicos con Excel: un tutorial práctico. Espero os sea útil.

Encabezado del trabajo

Caracterización analítica y diferenciación geográfica de pimentón mediante técnicas de reconocimiento de patrones

El pasado 11 de noviembre de 2015, una investigadora de nuestro grupo, Análisis Aplicado (FQM-347) defendió la tesis doctoral  homónima a esta entrada del blog. La tesis la realizó la Dra. Ana Palacios Morillo bajo la dirección del Profesor Fernando de Pablos y un servidor. Os desgrano un poco de la misma, aunque creo que es preferible que le echéis un vistazo en el repositorio de la Universidad de Sevilla:  Caracterización analítica y diferenciación geográfica de pimentón mediante técnicas de reconocimiento de patrones (Ana Palacios, noviembre 2015), por estar el texto completo y con una introducción muy amena.

Primera página de la tesis

La hipótesis de partida fue pensar que como en España disponemos de dos denominaciones de origen de pimentón, que se elaboran con pimientos cultivados en áreas concretas (La Vera en Extremadura y Murcia) y bajo unos estándares de producción diferentes, éstas podían dar lugar a productos diferenciables. El tema es interesante, si no, ¿a qué proteger estos productos con una marca del tipo Denominación de Origen Protegida?

El primer paso fue trabajar con parámetros químicos relacionados con el suelo. Sí señor, hablamos de los metales. Se desarrolla entonces un método de digestión y análisis para este tipo de muestras mediante espectroscopia de emisión atómica de plasma acoplado inductivamente (ICP-AES). Se determina un total de 14 elementos: Al, B, Ca, Cu, Fe, K, Mg, Mn, Na,  Ni, Sr, P, Pb y Zn. Para poder crear modelos de clasificación se estudian las muestras separadas en tipos (pimentón dulce, agridulce y picante) y en conjunto. 

Las técnicas lineales habituales, como el análisis discriminante lineal (LDA) consigue resultados  favorables en el caso de pimentones dulces y picantes, no así en agridulces y por tanto globalmente. 

Distribución de las muestras respecto a la función discriminante calculada mediante LDA para diferenciación geográfica de pimentón de la Vera y de Murcia. A) Modelo global, B) dulces, C) agridulces, D) picantes.

Se probaron otras técnicas como las máquinas de vectores soporte (SVM) y  modelado suave independiente por analogía de clases, sin mejora de los resultados (90-955 de eficacia de clasificación). Estos resultados mejoran notablemente al emplear redes neuronales artificiales de perceptrones multicapa (MLP-ANN), con un porcentaje de eficacia del 99 % para el modelo de diferenciación global. Todo esto da lugar a una publicación titulada: Geographical characterization of Spanish PDO paprika by multivariate analysis of multielemental content (Palacios-Morillo et al., Talanta, 2014, 128,15-22)

Resumen gráfico de la publicación en Talanta

Otro de los objetivos era la diferenciación en base a parámetros relacionados con el color (color ASTA y coordenadas CIELab de un extracto en acetona) o bien a partir de espectros de absorbancia en el UV-Vis (de 380-780 nm), previamente combinados mediante análisis en componentes principales (PCA). En este caso el PCA sirve para reducir el número de variables y evitar trabajar con datos correlacionados, hecho común en los espectros del visible, como se observa en la figura de abajo.

Espectros promedio y desviaciones estándar de muestras de pimentón agridulce. Las bandas anchas implican mucha correlación entre variables.

Tras probar modelos de LDA, SVM y MLP-ANN, se consigue diferenciar con estos últimos con eficacias de en torno al 95% tanto al partir de los parámetros del color como de los espectros reducidos con PCA. Esto da lugar a otra publicación titulada: Differentiation of Spanish paprika from Protected Designation of Origin based on color measurements and pattern recognition (Palacios-Morillo et al., Food Control, 2016,62, 243-249).

Resumen gráfico de la publicación en Food Control

Finalmente se desarrolla un método de determinación de capsacinoides mediante cromatografía líquida de alta eficacia (HPLC) con detector UV-Vis. Se desarrollan además modelos que permiten predecir el grado de picor de pimentón picante a partir de los espectros de absorción (280-800 nm). Para ello se aplicaron técnicas de regresión lineal múltiple, regresión en componentes principales (PCR), regresión por mínimos cuadrados parciales (PLS) y redes neuronales artificiales (ANN). El modelo que ofrece los mejores resultados está basado en la combinación de PLS y ANN, con una eficacia del 80%.

Os dejo sin más dilación con el texto completo de la tesis: Caracterización analítica y diferenciación geográfica de pimentón mediante técnicas de reconocimiento de patrones (Ana Palacios, noviembre 2015).

Nota: Enhorabuena Ana, otra vez.

Resolviendo sistemas de ecuaciones en Excel

Para ilustrar este procedimiento me van a permitir que use un problema sencillo de los que aparecen en algunos libros de texto de análisis instrumental.

  Se tiene una mezcla de dos especies, A y B, que presentan espectros de absorción parcialmente superpuestos. Calcule la concentración de A y B en la mezcla a partir de los siguientes datos. El espesor de la cubeta de muestra es de 1 cm.

	C (M)	A (254 nm)	A (550 nm)
Patrón de A	1.50 x 10 -4	0.003	0.975
Patrón de B	1.80 x 10 -4	0.589	0.017
Mezcla		0.857	0.909

Para resolver el problema, primero se calculan los coeficientes de absortividad molar de A y B a ambas longitudes de onda usando los datos de los patrones por separado. Ya que para cada longitud de onda se cumple la Ley de Beer, se tendrá:

Patrón de especie A:

Patrón de especie B:

Para la mezcla se tendrá el siguiente sistema de ecuaciones:

Si se hubiese llamado Y a la absorbancia, X1 al producto del paso de luz por la absortividad molar  para el compuesto A  y X2 al producto del paso de luz por la absortividad molar de B, el problema se reduce a una calibración lineal múltiple con dos niveles para X1 y X2 (uno por cada longitud de onda a la que se mide) con unos coeficiente a ajustar que se corresponderán con las concentraciones de A y B en la mezcla. Así el modelo lineal múltiple a resolver será:

Aquí es donde entra la función ESTIMACION.LINEAL() de Excel. Para solucionar el sistema de ecuaciones se escribe en una primera columna (columna B en la imagen) los valores de absorbancia de la mezcla a cada longitud de onda. En una segunda columna (C) se escriben los valores para X1 a cada longitud de onda y en la tercera columna (D) lo mismo para X2.

Siguiendo la disposición de las celdas de la figura anterior se selecciona el rango B8:D8 y se inserta la función =ESTIMACION:LINEAL(). En el formulario de entrada se selecciona el rango B2:B3 para la Y y el rango C2:D3 para los valores de X (que son X1 y X2). El cuadro de constante debe aparecer con el valor cero o falso, puesto que el modelo propuesto no presenta término independiente. Esto es lógico si en las medidas de los espectros se ha hecho el cero con el blanco. El apartado de estadística lo pondremos con el valor lógico verdadero, aunque para este tipo de problemas no nos sirve para mucho, puesto que la solución del sistema es única. 

Para terminar pulsamos a la vez Ctrl + shift + enter y apareceran los valores ce CB en B6 y CA en C6. 

¡Cuidado, que si ponemos el orden de las columnas de datos con  la especie A primero y B después, Excel invierte el orden y devuelve primero el coeficiente (concentración en nuestro caso) de B y luego la de A!

Es decir, si ordenamos las columnas de los valores de X poniendo primero X1 y a su derecha X2, Excel devuelve en la matriz de resultados primero el coeficiente de X2 y luego el de X1.