Regresión lineal en una muestra que abarca muchos órdenes de magnitud.

La ley de Beer de la química dice que la absorbancia de un líquido es proporcional a la concentración , entonces: Lo que hay que hacer entonces es preparar un conjunto de soluciones con concentraciones conocidas, medir la absorbancia para formar una 'curva estándar '(una curva de calibración básicamente), y haga una regresión lineal simple en esos datos para obtener la proporcionalidad (que luego puede usarse para predecir las concentraciones de soluciones desconocidas).$A$$C$

$$A = kC$$

Una manera fácil de hacer esto es comenzar con una concentración conocida y realizar una dilución en serie, lo que le daría una dilución 2x, 4x, 8x, 16x .... etc. Es decir, si comienza con una solución de obtendría soluciones con , , etc.$100\mu \mathrm g/\mathrm {mL}$$50\mu \mathrm g/\mathrm {mL}$$25\mu \mathrm g/\mathrm {mL}$$12.5\mu \mathrm g/\mathrm {mL}$

Ahora, cuando realiza la regresión lineal, tiene un conjunto de datos con muchos puntos de datos a bajas concentraciones y muy pocos a concentraciones más altas. Parece mucho más natural representar este problema en una escala logarítmica. Mi pregunta es, entonces, ¿debería estar haciendo una regresión lineal de vs. o vs. ? Cuando comparo los modelos, parecen dar respuestas que están en el mismo orden de magnitud pero en el orden del 30% diferente.$A$$C$$\log A$$\log C$

Deja que la física (del experimento y el aparato de medición) te guíe.

En última instancia, la absorción se determina midiendo las cantidades de radiación que pasan a través del medio y esas mediciones se reducen a contar fotones. Cuando el medio es macroscópico, las fluctuaciones termodinámicas en la concentración son insignificantes, por lo que la principal fuente de error radica en el recuento. Este error (o "ruido de disparo" ) tiene una distribución de Poisson . Esto implica que el error es relativamente grande a altas concentraciones cuando pasa poca radiación.

Con suficiente cuidado en el laboratorio, las concentraciones generalmente se miden con extrema precisión, por lo que no me preocuparé por errores en las concentraciones.

La absorbancia en sí está directamente relacionada con el logaritmo de la radiación medida . Tomar el logaritmo iguala la cantidad de error en todo el rango posible de concentraciones. Solo por esta razón, es mejor analizar la absorbancia en términos de sus valores habituales en lugar de volver a expresarlos. En particular, debemos evitar tomar registros de absorbancia, aunque eso simplificaría la expresión de la ley de Beer-Lambert.

También debemos estar atentos a posibles no linealidades. La derivación de la Ley de Beer-Lambert sugiere que la curva de absorbancia vs concentración se volverá no lineal a altas concentraciones. Se necesita alguna forma de detectar o probar esto.

Estas consideraciones sugieren un procedimiento simple para analizar una serie de de concentraciones y absorbancias medidas:$(C_i, A_i)$

• Estime el coeficiente como la media aritmética de , .$\kappa$$A/C$$\hat{\kappa} = \sum_i \frac{A_i}{C_i}$

• Predecir la absorbancia en cada concentración en términos del coeficiente estimado:$\hat{A}(C) = \hat{\kappa}C.$

• Verifique los residuos aditivos para tendencias no lineales en .$A_i - \hat{A_i}$$C_i$

Por supuesto, todo esto es teórico y algo especulativo, no tenemos datos reales para analizar, pero es un lugar razonable para comenzar. Si la experiencia de laboratorio repetida sugiere que los datos se apartan de los comportamientos estadísticos descritos aquí, entonces se requerirían algunas modificaciones de estos procedimientos.

---

Para ilustrar estas ideas, he creado una simulación que implementa los aspectos clave de la medición, incluido el ruido de Poisson y posiblemente las respuestas no lineales. Al ejecutarlo muchas veces, podemos observar el tipo de variación que es probable que se encuentre en el laboratorio. Aquí están los resultados de una ejecución de simulación. (Se pueden realizar otras simulaciones simplemente cambiando la semilla inicial en el siguiente código y modificando varios parámetros según se desee).

Este experimento simulado midió la absorbancia a concentraciones de hasta . Los márgenes verticales en valores aparentes en el diagrama de dispersión muestran los efectos de (a) ruido de disparo en las mediciones de transmisión y (b) ruido de disparo en la medición de transmisión inicial a concentración cero. (Observe cómo esto realmente crea algunos valores de absorbancia negativos ). Aunque los errores resultantes no van a tener exactamente las mismas distribuciones en cada concentración, los spreads más o menos iguales son evidencia empírica de que las distribuciones son lo suficientemente cercanas como para ser lo mismo que necesitamos. No te preocupes por eso. En otras palabras, no hay necesidad de ponderar las absorbancias de acuerdo con las concentraciones.$1$$1/32$

La línea diagonal roja se ha estimado a partir de las 50 simulaciones. Tiene una pendiente de , que difiere ligeramente de la pendiente físicamente correcta de que se usó en las simulaciones. Esta desviación es muy grande porque supuse que había muy poca radiación para medir; el recuento máximo de fotones fue de solo . En la práctica, los recuentos máximos podrían ser muchos órdenes de magnitud mayores que esto, lo que llevaría a estimaciones de pendiente muy precisas, ¡pero entonces no aprenderíamos mucho de esta cifra!$\hat{\kappa}=2.13$$2$$1000$

El histograma de los residuos no se ve bien: está sesgado a la derecha. Esto indica algún tipo de problema. Ese problema no proviene de la asimetría en los residuos en cada concentración; más bien, proviene de una falta de ajuste. Eso es evidente en los diagramas de caja a la derecha: aunque los primeros cinco se alinean casi horizontalmente, el último, en la concentración más alta, difiere claramente en su ubicación (es demasiado alta) y escala (es demasiado larga) . Esto resulta de una respuesta no lineal que construí en la simulación. Aunque la no linealidad está presente en todo el rango completo de concentraciones, tiene un efecto apreciable solo en las concentraciones más altas. Esto es más o menos lo que sucedería también en el laboratorio. Sin embargo, con solo una ejecución de calibración disponible, no pudimos dibujar tales diagramas de caja. Considere analizar múltiples ejecuciones independientes si la no linealidad puede ser un problema.

---

La simulación se realizó en R. Sin embargo, los cálculos con datos reales son simples de realizar a mano o con una hoja de cálculo: solo asegúrese de verificar que los residuos no sean lineales.

#
# Simulate instrument responses:
#   `concentration` is an array of concentrations to use.
#   `kappa` is the Beer-Lambert law coefficient.
#   `n.0`   is the largest  expected photon count (at 0 concentration).
#   `start` is a tiny positive value used to avoid logs of zero.
#   `beta`  is the amount of nonlinearity (it is a quadratic perturbation
#           of the Beer-Lambert law).
# The return value is a parallel array of measured absorbances; it is subject
# to random fluctuations.
#
observe <- function(concentration, kappa=1, n.0=10^3, start=1/6, beta=0.2) {
  transmission <- exp(-kappa * concentration - beta * concentration^2)
  transmission.observed <- start + rpois(length(transmission), transmission * n.0)
  absorbance <- -log(transmission.observed / rpois(1, n.0))
  return(absorbance)
}
#
# Perform a set of simulations.
#
concentration <- 2^(-(0:5)) # Concentrations to use
n.iter <- 50                # Number of iterations
set.seed(17)                # Make the results reproducible
absorbance <- replicate(n.iter, observe(concentration, kappa=2))
#
# Put the results into a data frame for further analysis.
#
a.df <- data.frame(absorbance = as.vector(absorbance))
a.df$concentration <- concentration # ($ interferes with TeX processing on this site)
#
# Create the figures.
#
par(mfrow=c(1,3))
#
# Set up a region for the scatterplot.
#
plot(c(min(concentration), max(concentration)), 
     c(min(absorbance), max(absorbance)), type="n",
     xlab="Concentration", ylab="Absorbance",
     main=paste("Scatterplot of", n.iter, "iterations"))
#
# Make the scatterplot.
#
invisible(apply(absorbance, 2, 
                function(a) points(concentration, a, col="#40404080")))
slope <- mean(a.df$absorbance / a.df$concentration)
abline(c(0, slope), col="Red")
#
# Show the residuals.
#
a.df$residuals <- a.df$absorbance - slope * a.df$concentration # $
hist(a.df$residuals, main="Histogram of residuals", xlab="Absorbance Difference") # $
#
# Study the residual distribution vs. concentration.
#
boxplot(a.df$residuals ~ a.df$concentration, main="Residual distributions",
        xlab="Concentration")

Ninguno de sus modelos ajustados puede ser correcto para su ecuación original, y su ecuación original no puede ser un modelo correcto para las variables aleatorias que observe.

Abordemos algunos problemas uno a la vez.

1) Su ecuación original es , pero, por supuesto, los datos no se observan en los valores de la población (de lo contrario, solo necesitaría un valor e , ya que . Claramente, este no es el modelo correcto para los datos o no estaría tratando de ajustar las regresiones. Veremos cómo reescribirlo en breve.$A=kC$$x$$y$$k=y/x$

2) Si ambas variables se observan con error, necesita técnicas más especializadas que los modelos de regresión lineal simples.

3) Si se observa sin error (o muy bajo error relativo a ), presumiblemente se refiere a un modelo algo así como decir .$C$$A$$E(A|C) = kC$

La pregunta es "cómo modelar la variación sobre la media": necesitamos una distribución para en cada$A$$C$, o al menos alguna información relacionada con la distribución.

a) Como sugirió, podría ajustar un modelo lineal directamente a los datos en la escala original:

$$\text{E}(A|C)=kC$$

Esto podría ser adecuado si la variabilidad es casi constante en esta escala ($\text{Var}(A|C)$ casi constante).

Tenga en cuenta que este modelo, basado en su modelo de población, no tiene término de intercepción; pasa por el origen.

Alternativamente, si tiene algún otro modelo para la variación sobre la línea, como $\text{Var}(A|C)\propto C^p$, entonces podría ajustarse a una regresión ponderada.

b) Como usted menciona, otro posible modelo puede ajustarse en la escala logarítmica:

$$\text{E}(\log_{10}A|\log_{10}C)=\log_{10}k+\log_{10}C$$

Esto podría ser adecuado si la variación es casi constante en la escala logarítmica ($\text{Var}(\log_{10}A|\log_{10}C)$ casi constante) ... que tendería a ver si la varianza-potencia en (a) fuera 2 o cercana a 2.

Tenga en cuenta que este modelo tiene una intersección, pero tiene un coeficiente de pendiente de 1. Una forma de ajustar este modelo es ajustar:

$$\text{E}(\log_{10}A|\log_{10}C)-\log_{10}C=\log_{10}k$$

(Dicho esto, es posible que desee verificar un modelo más general que el original, como por ejemplo considerar una pendiente no unitaria en (b).

Los ajustes en (a) y (b) ponderan los datos de manera diferente (aunque con $p=2$en (a) los resultados serán bastante cercanos), por lo que no darán respuestas idénticas. En sus datos difieren en aproximadamente un 26%, lo que ilustra que para tal muestra la elección es bastante importante.

Si no tiene ningún conocimiento previo para guiar la elección, el análisis residual podría ser una forma de llegar a un modelo (especialmente si tiene otro conjunto de datos en el que puede basar esa elección); alternativamente, podrías hacer$p$ Un parámetro de un modelo más grande.

(Como sucede, un pequeño análisis residual me sugiere que puede haber problemas potenciales tanto con (a) con varianza constante como con el modelo en (b); ni el modelo, ni siquiera, quizás, los modelos más generales (a) - with0intercept y (b) -with-non-unit pendiente son necesariamente adecuados (hay una sugerencia de curvatura sobre la línea). Una de las preocupaciones es el apalancamiento en (a), que se deriva de los 'muchos órdenes de magnitud')

4) Tenga en cuenta que hay muchos otros modelos que pueden instalarse.

Por ejemplo, considere $E(A^q|C) = k^qC^q$ para alguna constante especificada $q$ (el modelo de registro corresponde a $q=0$)