跨越多个数量级的样本的线性回归

机器算法验证 回归
2022-03-25 12:36:56

化学中的比尔定律说液体成正比,因此: 标准的做法是准备一组已知浓度的溶液,测量吸光度以形成“标准曲线” '(基本上是一条校准曲线),并对该数据进行简单的线性回归以获得比例(然后可用于预测未知溶液的浓度)。AC

A=kC

一种简单的方法是从已知浓度开始并执行连续稀释,这将使您获得 2 倍稀释、4 倍、8 倍、16 倍......等。即如果你从的解决方案开始,你会得到的解决方案, ,等...100μg/mL50μg/mL25μg/mL12.5μg/mL

现在,当您进行线性回归时,您会得到一个数据集,其中包含大量低浓度数据点,而高浓度数据点则很少。在对数尺度上表示这个问题似乎更自然。那么我的问题是,我应该对进行线性回归吗?当我比较这些模型时,它们似乎给出了相同数量级但相差 30% 的答案。AClogAlogC

以线性比例绘制的数据 在对数刻度上绘制的数据

2个回答

让(实验和测量仪器的)物理学指导您。

最终,吸收是通过测量穿过介质的辐射量来确定的,而这些测量归结为光子计数。当介质是宏观的时,浓度的热力学波动可以忽略不计,因此误差的主要来源在于计数。该误差(或“散粒噪声”)具有泊松分布这意味着当很少有辐射通过时,在高浓度时误差相对较大。

在实验室中,如果有足够的注意,通常可以非常准确地测量浓度,因此我不会担心浓度误差。

吸光度本身与测量辐射的对数直接相关取对数可以平衡整个可能浓度范围内的误差量。仅出于这个原因,最好根据其通常值来分析吸光度,而不是重新表达它们。 特别是,我们应该避免记录吸光度,即使这会简化 Beer-Lambert 定律的表达。

我们还应该警惕可能的非线性。Beer-Lambert 定律 的推导表明,吸光度浓度曲线在高浓度下将变为非线性。需要某种方法来检测或测试这一点。

这些考虑提出了一个简单的程序来分析一系列对浓度和测量的吸光度:(Ci,Ai)

  • 估计系数作为κA/Cκ^=iAiCi

  • 根据估计系数预测每个浓度的吸光度:A^(C)=κ^C.

  • 检查中非线性趋势的加性残差AiAi^Ci

当然,所有这些都是理论上的,有些推测——我们没有任何实际数据可供分析——但这是一个合理的起点。如果重复的实验室经验表明数据与此处描述的统计行为不同,则需要对这些程序进行一些修改。


为了说明这些想法,我创建了一个模拟来实现测量的关键方面,包括泊松噪声和可能的非线性响应。通过多次运行,我们可以观察到实验室中可能遇到的那种变化。这是一次模拟运行的结果。(其他模拟可以简单地通过更改下面代码中的起始种子并根据需要修改各种参数来执行。)

数据

的吸光度散点图中明显的垂直扩散值显示了 (a) 透射测量中的散粒噪声和 (b) 零浓度下初始透射测量中的散粒噪声的影响。(注意这实际上是如何产生一些吸光度值的。)虽然由此产生的误差在每个浓度下不会具有完全相同的分布,但大致相等的分布是经验证据,表明分布足够接近我们需要的相同'不用担心。换言之,无需根据浓度对吸光度进行加权。11/32

红色对角线是根据所有 50 次模拟估算的。它的斜率为 ,与模拟中使用略有不同。这个偏差如此之大,因为我认为要测量的辐射非常少;最大光子数仅为在实践中,最大计数可能比这大许多数量级,从而导致高度精确的斜率估计——但是我们不会从这个数字中学到很多东西!κ^=2.1321000

残差直方图看起来不太好:它向右倾斜。这表明某种麻烦。问题不在于每个浓度的残差不对称;相反,它来自于不合身。这在右边的箱线图中很明显:虽然前五个几乎水平排列,但最后一个 - 在最高浓度 - 位置(太高)和规模(太长)明显不同. 这是我在模拟中内置的非线性响应造成的。尽管非线性存在于整个浓度范围内,但仅在最高浓度时才具有明显的影响。这或多或少也将在实验室中发生。然而,只有一次校准运行可用,我们无法绘制这样的箱线图。 如果非线性可能是一个问题,请考虑分析多个独立运行。


模拟在R. 但是,使用实际数据的计算很容易手动或使用电子表格进行:只需确保检查残差的非线性即可。

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

对于您的原始方程,您的拟合模型都不能正确,并且您的原始方程不能是您观察到的随机变量的正确模型。

让我们一次解决一些问题。

1)您的原始方程是,但当然在总体值上没有观察到数据(否则您只需要一个值,因为。显然这不是数据的正确模型否则你不会试图拟合回归。我们很快就会看到如何重写它。A=kCxyk=y/x

2)如果观察到两个变量都有错误,您需要比简单的线性回归模型更专业的技术。

3)如果没有错误(或相对于错误非常低),您可能是指一个模型,例如CAE(A|C)=kC

那么问题是“如何对均值的变化进行建模”——我们需要在每个分配一个分布,或者至少需要一些与分布相关的信息。AC

a)正如您所建议的,您可以将线性模型直接拟合到原始比例的数据:

E(A|C)=kC

如果可变性在这个尺度上接近恒定(接近恒定),这可能是合适的。Var(A|C)

请注意,此模型基于您的人口模型,没有截距项;它通过原点。

或者,如果您有一些其他模型用于线的变化,例如,那么您可以拟合加权回归。Var(A|C)Cp

b)正如你所提到的,另一个可能的模型可能适合对数尺度:

E(log10A|log10C)=log10k+log10C

如果变化在对数尺度上接近恒定(接近恒定),这可能适合... (a) 中的方差功率为 2 或接近 2。Var(log10A|log10C)

请注意,此模型具有截距,但斜率系数为 1。拟合此模型的一种方法是拟合:

E(log10A|log10C)log10C=log10k

(也就是说,您可能想要检查比原始模型更通用的模型,例如考虑(b)中的非单位斜率。

(a) 和 (b) 中的拟合对数据的权重不同(尽管,结果将非常接近),因此它们不会给出相同的答案。根据您的数据,它们相差约 26%,这说明对于这样的样本,选择非常重要。p=2

如果您没有任何先验知识来指导选择,残差分析可能是获得模型的一种方式(特别是如果您有另一个数据集可以作为该选择的基础);或者,您可以使成为更大模型的参数。p

(碰巧,一点残差分析向我表明,(a)具有恒定方差和(b)中的模型可能存在潜在问题;既不是模型,也可能是更一般的模型(a)- with0intercept 和 (b)-with-non-unit slope 一定是合适的(有关于线的曲率的建议)。其中一个问题是(a)中的杠杆作用,它源于“许多数量级”)

4) 请注意,还有许多其他型号可以安装。

例如,考虑某些指定常数(对数模型对应于)。E(Aq|C)=kqCqqq=0