在您的示例中,似乎該模型不適合數據很好,並且樣本量很小。通常情況下,這意味着出了問題,在進行任何進一步分析之前,應該修改模型。但是我仍然提供了一些通過引導方法計算「置信區間」的方法,雖然在這種情況下可能無效。
這些數據,我們需要:
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))
基於該模型,我們可以定義一個函數用於爲每個隨機生成的索引生成擬合值和預測的向量。該函數的參數應該是一些樣本索引,並且在該函數中,基於具有輸入索引的樣本的模型被擬合,並且從擬合模型生成擬合值和預測的向量(因爲在問題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
是存儲引導結果的矩陣。 每個引導程序樣本與原始樣本(本例中爲14)的大小相同,並且引導程序樣本是從原始樣本生成的,並使用簡單的隨機抽樣替換。因此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的建議一般用途。
請正確格式化您的問題。 –
在nls()函數中它應該是cum.num not cum.vul –
通常'predict()'將接受interval參數,您希望將其設置爲「confidence」,但幫助頁面?predict.nls'告訴我們「這個論點被忽略了」,所以你不能做你想做的事情,至少用'stats :: predict()' – Nate