2017-06-01 78 views
2

我的數據由兩列的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) 

我的目標是使用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) 

我可以具有較低和沿着計算new.time的預測值上限?

非常感謝您的幫助!

+0

請正確格式化您的問題。 –

+0

在nls()函數中它應該是cum.num not cum.vul –

+0

通常'predict()'將接受interval參數,您希望將其設置爲「confidence」,但幫助頁面?predict.nls'告訴我們「這個論點被忽略了」,所以你不能做你想做的事情,至少用'stats :: predict()' – Nate

回答

3

在您的示例中,似乎該模型不適合數據很好,並且樣本量很小。通常情況下,這意味着出了問題,在進行任何進一步分析之前,應該修改模型。但是我仍然提供了一些通過引導方法計算「置信區間」的方法,雖然在這種情況下可能無效。

這些數據,我們需要:

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的建議一般用途。

+0

謝謝你...太多了.. @一致性 –

+0

我*想*你的回答是有幫助的準確......但我無法圍繞如何編寫函數來匹配我的數據。有沒有辦法以更一般的形式寫出來?或解釋條款?我不想問一個新問題,因爲我認爲你的答案接近我所需要的。 – user3386170

+1

@ user3386170我很高興答案是有幫助的。我在答案中增加了更多解釋,試圖說明bootstrap方法的基本思想。 – Consistency