2015-05-14 75 views
1

我正在用R中的mle2進行非線性迴歸,並且我想圍繞最佳擬合曲線生成95%的逐點置信區間。我試圖做的一個簡化例子如下所示,我嘗試使用mle2擬合對象的預測。有關如何做到這一點的任何建議?mle2迴歸線的置信區間

library(bbmle) 

# Fabricated data 
e.u <- function(x, k) { exp(-k * x) } 
n <- 40 
t.bio <- 1:n 
bio <- 10*e.u(t.bio,log(2)/10) + rnorm(n,0,sqrt(e.u(t.bio,log(2)/10))) 

#Use mle2 to estimate the parameters 
intake.guess <- 10 
rc.guess <- 0.07 

n.log.like <- function(intake,k) { 
    sum.y <- 0 
    for (i in 1:length(bio)) { 
    x <- intake * e.u(t.bio[i],k) 
    y <- bio[i] 
    sum.y <- sum.y + log(dnorm(y,x,0.1*sqrt(y))) } 
    return(-sum.y) 
} 

b <- mle2(n.log.like, 
     start=list(
      intake=intake.guess, 
      k=rc.guess), 
     data=list(
      t.bio=t.bio, 
      bio=bio), 
     method="Nelder-Mead", 
     skip.hessian=FALSE) 

intake <- coef(summary(b))[1,1]   
rc <- coef(summary(b))[2,1]   
summary(b) 

#Scatter plot 
bio.p <- numeric(n) 
x <- 1:n 
for (i in 1:n) { bio.p[i] <- intake * e.u(x[i],rc) } 
plot(x,bio.p,type="l",log="xy",main="", 
    xlab="Days After Intake",ylab="Excretion") 
points(t.bio,bio) 

#I want to generate a confidence interval on the regression line 
bio.hat <- predict(b) 

回答

0

你可能想試試這個資源:

https://stat.ethz.ch/R-manual/R-patched/library/stats/html/predict.lm.html

有一個se.fit參數與predict()通過。如果TRUE,它會計算SE。我沒有嘗試過。

另一種方法是使用library(ggplot2),其中一個參數可以自動覆蓋圖上的置信度。一個例子,

c <- ggplot(mtcars, aes(qsec, wt)) 
c + stat_smooth(se = TRUE) + geom_point() 

這將圍繞劇情放一個樂隊,以表明置信區間。默認情況下,se設置爲TRUE。

參考: http://svitsrv25.epfl.ch/R-doc/library/ggplot2/html/stat_smooth.html

第三,你可以計算置信區間,並通過將參數疊加圖(新= TRUE)。這可能很麻煩但仍然可行。

+0

這對於通用的MLE配合效果不佳。 Bootstrapping可能是最簡單的方法,除非問題真的很大/很難。 –

+0

我按照Ben Bolker的建議使用bootstrap。 – user1930565