如何通过 nls() 计算拟合值的置信区间
How to compute confidence interval of the fitted value via nls()
我的数据由两列组成——时间和累计数,如下所示:
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)
我的 objective 是使用非线性回归对我的数据建模,使用 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
的预测值以及下限和上限吗?
非常感谢您的帮助!!!!
在你的例子中,模型似乎不太适合数据,样本量也很小。通常,这意味着出现问题,您应该在进行任何进一步分析之前修改您的模型。但是我还是提供了一些通过bootstrap方法计算"confidence interval"的方法,虽然在这种情况下可能无效。
这些是我们需要的数据:
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结果的矩阵。
每个 bootstrap 个样本与原始样本具有相同的大小(本例中为 14 个),并且 bootstrap 个样本是从原始样本中通过带放回的简单随机抽样生成的。因此 sample.int(14, replace = TRUE))
用于为 bootstrap 个样本生成索引。 estimate
函数用于为每个随机生成的索引生成拟合值向量和预测值。
由于 predict1
是 bootstrapped 拟合值和预测,我从 bootstrapped 估计中计算出 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 的建议进行一些编辑以提高可读性并说明如何将代码用于一般用途。
我的数据由两列组成——时间和累计数,如下所示:
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)
我的 objective 是使用非线性回归对我的数据建模,使用 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
的预测值以及下限和上限吗?
非常感谢您的帮助!!!!
在你的例子中,模型似乎不太适合数据,样本量也很小。通常,这意味着出现问题,您应该在进行任何进一步分析之前修改您的模型。但是我还是提供了一些通过bootstrap方法计算"confidence interval"的方法,虽然在这种情况下可能无效。
这些是我们需要的数据:
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结果的矩阵。
每个 bootstrap 个样本与原始样本具有相同的大小(本例中为 14 个),并且 bootstrap 个样本是从原始样本中通过带放回的简单随机抽样生成的。因此 sample.int(14, replace = TRUE))
用于为 bootstrap 个样本生成索引。 estimate
函数用于为每个随机生成的索引生成拟合值向量和预测值。
由于 predict1
是 bootstrapped 拟合值和预测,我从 bootstrapped 估计中计算出 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 的建议进行一些编辑以提高可读性并说明如何将代码用于一般用途。