机器算法验证 - 检验非线性回归模型系数之间的差异 - 吾爱随笔录

检验非线性回归模型系数之间的差异

机器算法验证回归广义线性模型非线性回归

2022-03-23 05:41:58

让我们考虑以下数据，显示两种不同化合物（蓝色和红色）的 S 型剂量依赖性：

我想知道比较蓝色和红色“曲线”的半激活点（= 半最大效应，这里分别在浓度 4 和 5 左右）的最佳方法。科学问题是“两种不同化合物在相同浓度下达到半最大活化的零假设的 p 值是多少”？

在这种特殊情况下，这些点没有以任何方式链接（即，设计不是配对的；不能只在每组中拟合五个 sigmoid，然后比较半激活系数，这在配对设计中是可能的）。

我知道使用线性回归模型，可以添加一个交互项来获得偏移或斜率差异的 p 值。但是，我不确定如何将这种方法扩展到非线性回归。

（为了完整起见，我知道在某些领域，这种情况的默认方法是 2-way ANOVA，但由于多种原因，我认为这不是一个好方法）。

编辑 1：描述数据的模型的形式为，其中是 x 轴位置，是半最大激活点，是sigmoid 的斜率。 $1 - 1/(1+(x/a)^b)$ $x$ $a$ $b$

编辑2：澄清数据点之间的关系：可以使用以下过程生成数据：

100 个独立* 细胞被随机分为 20 组（2 个化合物 x 10 个先验选择浓度），并放置在单独的室中。
测量每个细胞的特性（例如离子电流的密度）
每个细胞都用给定浓度的化合物 A（蓝色）或 B（红色）处理。这可能具有将浴液从无化合物更改为化合物 A 或化合物 B 的形式，或者通过在每个孔中移取浓缩原液，这可能会引入稍微不同的噪音模式。
重复步骤 2，并为每个点绘制对应于治疗效果的值（对于离子电流阻断，这可能是 $1-\frac{value-after-treatment}{value-before-treatment}$

这个问题是针对真正独立的细胞提出的。细胞在实践中是否真正独立可能取决于确切的实验设置 - 通常这可能在某种程度上被违反，需要对批次效应偏差进行调查。

3个回答

我可能会使用近似排列测试来解决这个问题。

一般的想法是，如果数据来自同一个函数（请注意，这将假设 $a_1 = a_2$ 和 $b_1 = b_2$ 对于两组数据，以及两个数据集的随机分量来自相同的分布），任何给定数据点分配给数据的哪个子集（在您的情况下为蓝色或红色）是无关紧要的。因此，我们可以打乱分配，计算红色和蓝色参数估计值的差异，然后一遍又一遍地重复这个过程，以估计如果原假设为真，两个参数估计值的差异有多大。在更多的统计术语中，我们可以从产生与我们使用原始数据计算的估计值一样大或更大的估计值的置换样本的分数中计算置换 p 值。

置换检验的两个优点是它们是无分布的，只要求观察值在原假设下是可交换的，并且您可以使用您想要的任何检验统计量 - 您不限于分布在原假设是已知的或可以计算的。它们也是渐近最强大的数据条件。

这是一个演示，使用似乎或多或少模仿您的构建数据。首先，构建数据：

# Create data; 50 observations for each of sample 1 and 2
x1 <- ceiling(10*runif(50))
x2 <- ceiling(10*runif(50))

a1 <- 4; b1 <- 8
a2 <- 5; b2 <- 6

func <- function(x, a, b) {
  p <- 1 - 1/(1 + (x/a)^b)
  y <- rbeta(length(p), p*15, (1-p)*15)
  y
}

df1 <- data.frame(y = func(x1, a1, b1), x = x1)
df2 <- data.frame(y = func(x2, a2, b2), x = x2)
df_H0 <- rbind(df1, df2)

构建的数据如下所示：

plot(y~x, data=df1, col="blue", pch=19, cex=1.25, main="Sample points")
points(y~x, data=df2, col="red", bg="red", pch=24, cex=1.25)

接下来，我们计算参数估计之间的差异 $a$ 在两个样本上：

fm <- as.formula("y ~ 1 - 1/(1+ (x/a)^b)")
m1 <- nls(fm, data=df1, start=list(a=5, b=5))
m2 <- nls(fm, data=df2, start=list(a=5, b=5))
stat_21 <- m2$m$getPars()[1] - m1$m$getPars()[1]
stat_21
       a 
1.153752

现在进行排列计算：

boot_21 <- rep(0, 1000)
for (i in seq_along(boot_21)) {
  indx1 <- sample(nrow(df_H0))[1:nrow(df1)]
  m1 <- nls(fm, data=df_H0[indx1,], start=list(a=5, b=5))
  m2 <- nls(fm, data=df_H0[-indx1,], start=list(a=5, b=5))
  boot_21[i] <- as.numeric(m2$m$getPars()[1] - m1$m$getPars()[1])
}

hist(boot_21)
abline(v=stat_21)

这给了我们以下情节：

最右边的垂直条等于实际统计数据。置换 p 值为：

mean(abs(boot_21) >= stat_21)
[1] 0

这意味着我们的 1,000 个排列样本中没有一个在参数估计中的差异与我们实际观察到的一样大。

在数据集中定义一个新变量red，红色点为 1，蓝色点为 0。c然后在回归方程中包含一个额外的系数：

1−1/(1+(x/(a+c*red))^b)

c表示拟合到蓝色和红色点的曲线之间的偏移量。如果c与零显着不同，则可以认为两条响应曲线不同。

另外，请参阅 R 中的drc包。

评论部分充斥着迂腐的注释和关于数据中噪音来源的问题，我感觉有点糟糕。所以我会用一个比那些小评论更体面的答案来弥补它。

这个答案将证明一个简单的（非线性）最小二乘拟合（也可以用 glm* 完成）将如何高估两种化合物之间差异的重要性。发生这种高估的原因是噪声不是均匀的，而用于计算 p 值的模型假设所有点都具有相同的噪声。实际上，在模型之间产生很大差异的主要是中点周围的几个点。因此，有效自由度远小于计算假设的。

这使得来自@jbowman 答案的排列测试成为一种更强大的方法。置换测试的问题是您需要足够数量的相关数据点才能使测试变得强大。此外，由于特定的实验程序导致错误不独立，因此置换测试不能避免数据点相关的情况。其他一些方法，例如使用对整个过程进行建模并包含所有错误的贝叶斯模型可能会更好。

底线是没有单一的解决方案可以应用于所有这些类型的曲线。您不能将一组数据连同一个函数而没有其他信息一起提供给统计学家。了解数据的来源和生成方式非常重要。这需要纳入统计模型。

在此示例中，我们生成数据点，不是通过向曲线添加噪声，而是通过更改每个数据点的参数。

a \sim Unif (4, 6) b = 3 y = 1 - \frac{1}{1 + (x / a)^{b}}

$a \sim \text{Unif}(4,6) \qquad b = 3 \\ y = 1- \frac{1}{1+(x/a)^b}$

然后我们会得到一个这样的图来生成两次 80 个数据点。

曲线的可变性不像添加到某个“真实”值的噪声，而是因为曲线本身的变化每次都不同。这可能与使用的（细菌？）细胞的变化相对应。它们中的每一个都是不同的，并且可能根据不同的系数表现 $a$ 和 $b$ . 结果是噪声不均匀。大多在中间 $x=5$ 变化最大的地方。

如果我们重复这个模拟 $10^4$ 次并使用 glm 模型或非线性最小二乘拟合计算 p 值，那么 p 值的分布如下所示

因此，您不会得到 p 值的均匀分布，并且模型会高估出现某种差异的概率。

原因是变化主要发生在几个点，而 p 值的计算假设它均匀地发生在所有点。这高估了自由度。特别是包含点 $x=0$ 完全没用，因为系数的变化对结果的价值没有影响。

示例代码：

sigmoid = function(x,a=5,b=3) {
  1-(1+(x/a)^b)^-1
}

simulate_test = function(plot = TRUE) {
  ### create data
  x = rep(0:15,10) 
  y = sigmoid(x, a = runif(length(x),4,6))
  compound = c(rep(0,16*5), rep(1,16*5))
  
  ### make a plot if desired
  if (plot == TRUE) {
    plot(x, y, pch = 21, col = 1, bg = 0 + compound * 2, cex = 0.7)
  }
  
  ### Perform fitting
  ###
  ### we can do the fitting with nls but the lines below shows that glm works as well
  ### For glm points with x=0 need to be removed because we use log(x), but these do not add information anyway
  modnls = nls(y ~ sigmoid(x,a+c*compound,b), start = c(a=5,b=3,c=0),
               control = nls.control(minFactor = 10^-4,warnOnly = TRUE))
  modglm = glm(y[-which(x==0)] ~ log(x[-which(x==0)]) + compound[-which(x==0)], family = gaussian(link = "logit"))
  
  #lines below demonstrate how you would get the (same) coefficients with the two methods glm vs nls
  #coefficients(modnls)
  #coefficients(modglm)
  #exp(-modglm$coefficients[1]/modglm$coefficients[2])
  #exp((-modglm$coefficients[1]-modglm$coefficients[3])/modglm$coefficients[2])-exp(-modglm$coefficients[1]/modglm$coefficients[2])

  return(list(nls_p = summary(modnls)$coefficients[3,4],
              glm_p = summary(modglm)$coefficients[3,4]))

}


set.seed(1)
x = replicate(10^4, as.numeric(simulate_test(plot = FALSE)))
x = list(nls_p = x[1,], 
         glm_p = x[2,]) 

hist(x$glm_p, breaks = seq(0,1,0.01))
hist(x$nls_p, breaks = seq(0,1,0.01))

simulate_test(plot = TRUE)

^{*拟合也可以使用广义线性模型来完成，因为我们可以找到一个链接函数，使得转换后的值是回归量的线性函数： $\log \left(\frac{y}{1-y}\right) = b \log(x) - b \log(a)$ . 这在代码中得到了证明。由于计算错误，系数和估计值之间存在细微差异，而 p 值则因为 glm 模型必须使用 $\log(x)$ 并排除该值 $x=0$ . 但这些价值观 $x=0$ 无论如何不添加任何信息，并使 glm 模型实际上比 nls 模型更精确。}

其它你可能感兴趣的问题

上一篇是否可以使用机器学习作为学习统计数据的方法，而不是反之亦然？下一篇卡方和假设检验的正确使用？