如何通过 nls() 计算拟合值的置信区间答案

【问题标题】：How to compute confidence interval of the fitted value via nls()如何通过 nls() 计算拟合值的置信区间
【发布时间】：2017-06-01 16:43:14
【问题描述】：

我的数据由两列组成 - 时间和累计数，如下所示：

time <- c(1:14)
cum.num <- c(20, 45, 99, 195, 301, 407, 501, 582, 679, 753, 790, 861, 1011, 1441)

我的非线性函数是：

B/(B*C*exp(-A*B*time) + 1)

我的目标是使用nls() 使用非线性回归对我的数据进行建模，并找到拟合值的置信区间。我已经尝试了以下

m1 <- nls(cum.num ~ B/((B*C)*exp(-A*B*time) + 1),start=list(A=0.001,B=1000,C=0.5))

我尝试了以下方法来计算模型的拟合值：

predict(m1,interval="predict")

我只得到了没有上下置信区间的拟合值：

[1]  116.9912  145.7954  181.1951  224.4367  276.8663  339.8665  414.7550
[8]  502.6399  604.2369  719.6632  848.2417  988.3638 1137.4632 1292.1377

我的问题是：

a) 有什么方法可以计算拟合值的下限和上限？（通常lm()函数默认产生拟合值、下限和上限）

b) 假设我有新的时间，比如：

new.time<-c(15:20)

我可以计算cum.num 在new.time 的预测值以及上下限吗？

非常感谢您的帮助！！！！

【问题讨论】：

请正确格式化您的问题。
在 nls() 函数中应该是 cum.num 而不是 cum.vul
通常predict() 将接受间隔参数，您希望将其设置为“置信度”，但帮助页面?predict.nls 告诉我们“此参数被忽略”所以你不能做什么你想要，至少用stats::predict()

标签： r nls

【解决方案1】：

在您的示例中，模型似乎不太适合数据，并且样本量很小。通常，这意味着出现问题，您应该在进行任何进一步分析之前修改您的模型。但是我仍然提供了一些通过引导方法计算“置信区间”的方法，尽管在这种情况下它可能无效。

这些是我们需要的数据：

time <- c(1:14)
cum.num <- c(20, 45, 99, 195, 301, 407, 501, 582, 679, 753, 790, 861, 1011, 1441)
new.time <- c(15:20)
all.time <- c(time, new.time)

我们可能会给它们起其他名字，这有助于更一般的使用：

y=cum.num # the dependent variable values from data
x=time # the independent variable values from data
new.x=all.time # the independent variable values over which we want to predict

这是本例中使用的非线性最小二乘模型，将在方程中使用，但需要修改以用于一般情况：

nls(y ~ B/((B*C)*exp(-A*B*x) + 1), start=list(A=0.001,B=1000,C=0.5),
    control = nls.control(maxiter = 500, warnOnly = TRUE))

基于模型，我们可以定义一个estimate 函数，用于为每个随机生成的索引生成拟合值和预测向量。函数的参数应该是一些样本索引，并且在函数中，基于具有输入索引的样本拟合模型，并从拟合模型中生成拟合值和预测的向量（因为在问题 a需要拟合值和预测的 CI）。

estimate <- function(ind){
    x <- x[ind]
    y <- y[ind]
    m1 <- nls(y ~ B/((B*C)*exp(-A*B*x) + 1), start=list(A=0.001,B=1000,C=0.5),
          control = nls.control(maxiter = 500, warnOnly = TRUE))
    predict(m1, newdata = list(x = new.x))
}


m1 <- nls(cum.num ~ B/((B*C)*exp(-A*B*time) + 1),start=list(A=0.001,B=1000,C=0.5))
predict0 <- predict(m1, newdata = list(time = all.time))
predict1 <- replicate(1000, estimate(sample.int(14, replace = TRUE)))
intervals <- apply(predict1, 1, quantile, probs = c(0.05, 0.95))
rbind(predict0, intervals)

predict1 是一个存储引导结果的矩阵。每个 bootstrap 样本与原始样本（本例中为 14 个）具有相同的大小，并且 bootstrap 样本是从原始样本通过带放回的简单随机抽样生成的。所以sample.int(14, replace = TRUE)) 用于生成引导样本的索引。 estimate 函数用于为每个随机生成的索引生成拟合值和预测向量。

由于predict1 是自举拟合值和预测值，因此我根据自举估计计算 90% CI。在 bootstrap 过程中，nls 函数有很多警告，这意味着数值上有问题，这与小样本量和失拟模型一致。最终的结果是这样的：

> rbind(predict0, intervals)
[,1]      [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]     [,9]
predict0 116.99118 145.79538 181.1951 224.4367 276.8663 339.8665 414.7550 502.6399 604.2369
5%        39.22272  67.34464 111.2190 173.7619 231.7736 289.7346 358.8469 436.2569 524.8187
95%      162.92948 190.60295 224.2462 266.1298 314.1032 392.3228 504.1270 611.3698 704.2803
[,10]    [,11]     [,12]     [,13]     [,14]     [,15]     [,16]     [,17]    [,18]
predict0 719.6632 848.2417  988.3638 1137.4632 1292.1377 1448.4271 1602.2033 1749.5981 1887.374
5%       627.1981 739.8984  822.7940  838.2366  846.9043  851.8955  854.2859  855.8558  856.873
95%      799.1904 923.1220 1068.4667 1231.6091 1416.4405 1631.2212 1900.6581 2220.5415 2617.839
[,19]     [,20]
predict0 2013.1701 2125.5890
5%        857.4619  857.8027
95%      3072.8531 3594.9036
>

编辑：根据@user3386170 的建议进行一些编辑以提高可读性并说明如何将代码用于一般用途。

【讨论】：

谢谢你...非常感谢..@Consistency
我认为您的回答很有帮助且准确......但我无法思考如何编写函数以匹配我的数据。有没有办法可以用更一般的形式写出来？或者解释一下条款？我不想问一个新问题，因为我认为你的答案接近我所需要的。
@user3386170 我很高兴这个答案很有帮助。我在答案中添加了更多解释，以尝试说明引导方法的基本思想。
感谢您提供更多信息。我已经对其进行了重新设计并设法使其与我的数据一起使用。我已经编写了一个通用案例方程式作为单独的答案，并带有一些额外的解释，如果您认为它是正确的，您可以使用这些解释来替换您的一些代码。我不想直接编辑你的答案，因为变化太大了。
@user3386170 如果百分比适中，例如 50%。那么新的“bootstrap”样本的样本量比原始样本的样本量要小得多。这很重要，因为样本量对置信区间很重要。如果您仅对 50% 的数据进行抽样，您的置信区间通常会更宽。如果有更宽的置信区间是可以的，那么只抽取中等百分比（如 50% 或 60%）的样本而不进行替换是没有问题的。