Régression linéaire sur un échantillon couvrant plusieurs ordres de grandeur

La loi de Beer de la chimie dit que l'absorbance d'un liquide est proportionnelle à la concentration , donc: La chose standard à faire est alors de préparer un ensemble de solutions avec des concentrations connues, de mesurer l'absorbance pour former une `` courbe standard '(une courbe d'étalonnage essentiellement), et effectuez une régression linéaire simple sur ces données pour obtenir la proportionnalité (qui peut ensuite être utilisée pour prédire les concentrations de solutions inconnues).$A$$C$

$$A = kC$$

Un moyen facile de le faire est de commencer avec une concentration connue et d'effectuer une dilution en série, ce qui vous donnerait une dilution 2x, 4x, 8x, 16x .... etc. Autrement dit, si vous commencez avec une solution de vous obtiendrez des solutions avec , , 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}$

Maintenant, lorsque vous effectuez la régression linéaire, vous disposez d'un ensemble de données avec beaucoup de points de données à de faibles concentrations et très peu à des concentrations plus élevées. Il semble beaucoup plus naturel de représenter ce problème sur une échelle logarithmique. Ma question est alors, dois-je faire une régression linéaire de vs ou vs ? Lorsque je compare les modèles, ils semblent donner des réponses qui sont du même ordre de grandeur mais de l'ordre de 30% différentes.$A$$C$$\log A$$\log C$

Laissez-vous guider par la physique (de l'expérience et de l'appareil de mesure).

En fin de compte, l'absorption est déterminée en mesurant les quantités de rayonnement traversant le milieu et ces mesures se résument au comptage des photons. Lorsque le milieu est macroscopique, les fluctuations thermodynamiques de concentration sont négligeables, la principale source d'erreur réside donc dans le comptage. Cette erreur (ou "bruit de grenaille" ) a une distribution de Poisson . Cela implique que l'erreur est relativement importante à des concentrations élevées lorsque peu de rayonnement passe.

Avec suffisamment de soin en laboratoire, les concentrations sont généralement mesurées avec une extrême précision, donc je ne vais pas m'inquiéter des erreurs de concentrations.

L'absorbance elle-même est directement liée au logarithme du rayonnement mesuré . Prendre le logarithme uniformise la quantité d'erreur sur toute la plage de concentrations possible. Pour cette seule raison, il est préférable d'analyser l'absorbance en fonction de ses valeurs habituelles plutôt que de les ré-exprimer. En particulier, il faut éviter de prendre des logs d'absorbance, même si cela simplifierait l'expression de la loi de Beer-Lambert.

Nous devons également être attentifs aux éventuelles non-linéarités. La dérivation de la loi de Beer-Lambert suggère que la courbe absorbance vs concentration deviendra non linéaire à des concentrations élevées. Un moyen de détecter ou de tester cela est nécessaire.

Ces considérations suggèrent une procédure simple pour analyser une série de paires de concentrations et d'absorbances mesurées:$(C_i, A_i)$

• Estimez le coefficient comme la moyenne arithmétique de , .$\kappa$$A/C$$\hat{\kappa} = \sum_i \frac{A_i}{C_i}$

• Prédisez l'absorbance à chaque concentration en fonction du coefficient estimé:$\hat{A}(C) = \hat{\kappa}C.$

• Vérifiez les résidus additifs pour les tendances non linéaires dans .$A_i - \hat{A_i}$$C_i$

Bien sûr, tout cela est théorique et quelque peu spéculatif - nous n'avons pas de données réelles à analyser - mais c'est un point de départ raisonnable. Si des expériences répétées en laboratoire suggèrent que les données s'écartent des comportements statistiques décrits ici, alors certaines modifications de ces procédures seraient nécessaires.

---

Pour illustrer ces idées, j'ai créé une simulation qui met en œuvre les aspects clés de la mesure, y compris le bruit de Poisson et éventuellement les réponses non linéaires. En l'exécutant plusieurs fois, nous pouvons observer le type de variation susceptible d'être rencontrée en laboratoire. Voici les résultats d'une simulation. (D'autres simulations peuvent être effectuées simplement en changeant la graine de départ dans le code ci-dessous et en modifiant divers paramètres comme vous le souhaitez.)

Cette expérience simulée a mesuré l'absorbance à des concentrations de à . Les écarts verticaux dans les valeurs apparentes dans le diagramme de dispersion montrent les effets (a) du bruit de grenaille dans les mesures de transmission et (b) du bruit de grenaille dans la mesure de transmission initiale à concentration nulle. (Remarquez comment cela crée en fait des valeurs d'absorbance négatives .) Bien que les erreurs résultantes ne vont pas avoir exactement les mêmes distributions à chaque concentration, les écarts à peu près égaux sont la preuve empirique que les distributions sont assez proches d'être les mêmes que celles dont nous avons besoin '' t vous en souciez. En d'autres termes, il n'est pas nécessaire de pondérer les absorbances en fonction des concentrations.$1$$1/32$

La ligne diagonale rouge a été estimée à partir des 50 simulations. Il a une pente de , qui diffère légèrement de la pente physiquement correcte de qui a été utilisée dans les simulations. Cet écart est si important parce que je supposais qu'il y avait très peu de rayonnement à mesurer; le nombre maximal de photons n'était que de . Dans la pratique, les dénombrements maximaux pourraient être de plusieurs ordres de grandeur supérieurs à cela, conduisant à des estimations de pente très précises - mais nous n'apprendrions pas grand-chose de ce chiffre!$\hat{\kappa}=2.13$$2$$1000$

L'histogramme des résidus ne semble pas bon: il est biaisé vers la droite. Cela indique une sorte de problème. Ce problème ne vient pas de l'asymétrie dans les résidus à chaque concentration; il s'agit plutôt d'un manque d'ajustement. Cela est évident dans les boîtes à moustaches à droite: bien que les cinq premiers d'entre eux s'alignent presque horizontalement, le dernier - à la concentration la plus élevée - diffère clairement par l'emplacement (il est trop élevé) et l'échelle (il est trop long) . Cela résulte d'une réponse non linéaire que j'ai intégrée à la simulation. Bien que la non-linéarité soit présente dans toute la gamme des concentrations, elle n'a un effet appréciable qu'aux concentrations les plus élevées. C'est plus ou moins ce qui se passerait aussi en laboratoire. Cependant, avec un seul cycle d'étalonnage disponible, nous n'avons pas pu dessiner de tels diagrammes. Envisagez d'analyser plusieurs exécutions indépendantes si la non-linéarité peut être un problème.

---

La simulation a été réalisée en R. Les calculs avec des données réelles, cependant, sont simples à effectuer à la main ou avec un tableur: assurez-vous simplement de vérifier la non-linéarité des résidus.

#
# 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")

Aucun de vos modèles ajustés ne peut être correct pour votre équation d'origine, et votre équation d'origine ne peut pas être un modèle correct pour les variables aléatoires que vous observez.

Abordons certains problèmes un par un.

1) Votre équation d'origine est $A=kC$, mais bien sûr, les données ne sont pas observées aux valeurs de la population (sinon vous n'en auriez besoin que d'un $x$ et $y$ valeur, puisque $k=y/x$. De toute évidence, ce n'est pas le bon modèle pour les données ou vous n'essaieriez pas d'ajuster les régressions. Nous verrons comment le réécrire sous peu.

2) Si les deux variables sont observées avec erreur, vous avez besoin de techniques plus spécialisées que de simples modèles de régression linéaire.

3) Si $C$ est observé sans erreur (ou très faible erreur par rapport à $A$), vous voulez probablement dire un modèle quelque chose comme dire $E(A|C) = kC$.

La question est alors "comment modéliser la variation de la moyenne" - nous avons besoin d'une distribution pour $A$ à chaque $C$, ou au moins certaines informations relatives à la distribution.

a) Comme vous l'avez suggéré, vous pouvez adapter un modèle linéaire directement aux données sur l'échelle d'origine:

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

Cela pourrait convenir si la variabilité est presque constante sur cette échelle ($\text{Var}(A|C)$ presque constante).

Notez que ce modèle, basé sur votre modèle de population, n'a pas de terme d'interception; il passe par l'origine.

Alternativement, si vous avez un autre modèle pour la variation de la ligne, tel que $\text{Var}(A|C)\propto C^p$, vous pouvez alors ajuster une régression pondérée.

b) Comme vous le mentionnez, un autre modèle possible peut être monté sur l'échelle logarithmique:

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

Cela pourrait convenir si la variation est presque constante sur l'échelle logarithmique ($\text{Var}(\log_{10}A|\log_{10}C)$ presque constante) ... que vous auriez tendance à voir si la variance-puissance en (a) était de 2 ou proche de 2.

Notez que ce modèle a une ordonnée à l'origine, mais a un coefficient de pente de 1. Une façon d'ajuster ce modèle consiste à ajuster:

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

(Cela dit, vous voudrez peut-être vérifier un modèle plus général que votre modèle d'origine, par exemple en considérant une pente non unitaire en (b) par exemple.

Les ajustements en (a) et (b) pondèrent les données différemment (mais avec $p=2$dans (a) les résultats seront assez proches), donc ils ne donneront pas de réponses identiques. Sur vos données, elles diffèrent d'environ 26%, ce qui montre que pour un tel échantillon, le choix est assez important.

Si vous n'avez aucune connaissance préalable pour guider le choix, l'analyse résiduelle peut être un moyen d'arriver à un modèle (surtout si vous avez un autre ensemble de données sur lequel vous pouvez baser ce choix); vous pouvez également faire$p$ un paramètre d'un modèle plus grand.

(En l'occurrence, une petite analyse résiduelle me suggère qu'il peut y avoir des problèmes potentiels à la fois (a) à variance constante et avec le modèle en (b); ni modèle, ni même, peut-être, les modèles plus généraux (a) - with0intercept et (b) -with-non-unit pente sont nécessairement appropriés (il y a une suggestion de courbure sur la ligne) .L'une des préoccupations est l'effet de levier en (a), qui découle des «nombreux ordres de grandeur»)

4) Notez qu'il existe de nombreux autres modèles qui pourraient être installés.

Par exemple, considérez $E(A^q|C) = k^qC^q$ pour une constante spécifiée $q$ (le log-modèle correspond à $q=0$).