r - R中的箱线图比较不同样本量的95％置信区间的结果

Question

我试图通过使用自举来显示 gamma 分布的样本大小对 95% 置信区间边界的影响。现在，我需要将 4 个不同样本大小的结果放在一个箱线图中。R代码如下：

y <- rgamma(30,1,1) + rnorm(30,0,0.01)
y60 <- rgamma(60,1,1) + rnorm(60,0,0.01)
y100 <- rgamma(100,1,1) + rnorm(100,0,0.01)
y200 <- rgamma(200,1,1) + rnorm(200,0,0.01)
minusL <- function(params, data) {
-sum(log(dgamma(data, params[1], params[2])))
}
fit <- nlm(minusL, c(1,1), data=y)
fit
gammamedian<-function(data) {
fit <- nlm(minusL, c(1,1), data=data)
qgamma(.5, fit$estimate[1], fit$estimate[2])
}
gammamedian(y)
gammamedian(y60)
gammamedian(y100)
gammamedian(y200)
gengamma<- function(data, params){
rgamma(length(data), params[1], params[2])}
library(boot)
pbootresults<-boot(y, gammamedian, R=1000, sim="parametric",            
ran.gen=gengamma, mle=fit$estimate)
pbootresults
boot.ci(pbootresults, type=c("basic", "perc", "norm"))
pbootresults<-boot(y60, gammamedian, R=1000, sim="parametric",  
ran.gen=gengamma, mle=fit$estimate)
pbootresults
boot.ci(pbootresults, type=c("basic", "perc", "norm"))
pbootresults<-boot(y100, gammamedian, R=1000, sim="parametric",     
ran.gen=gengamma, mle=fit$estimate)
pbootresults
boot.ci(pbootresults, type=c("basic", "perc", "norm"))
pbootresults<-boot(y200, gammamedian, R=1000, sim="parametric",  
ran.gen=gengamma, mle=fit$estimate)
pbootresults
boot.ci(pbootresults, type=c("basic", "perc", "norm"))

 [An Excel image example ][1]


 [1]: https://i.stack.imgur.com/JXR6P.jpg

Answer 1

您必须将min和max值存储在数据表中，然后使用geom_rect()

我为其中两个案例提供了一个简单的示例，您可以为其余案例做类似的事情，并且还可以使用美学来获得更好的视野。

library(boot)
pbootresults<-boot(y, gammamedian, R=1000, sim="parametric",            
                   ran.gen=gengamma, mle=fit$estimate)
pbootresults
temp2 <-boot.ci(pbootresults, type=c("basic", "perc", "norm")) # store teh results
pbootresults<-boot(y60, gammamedian, R=1000, sim="parametric",  
                   ran.gen=gengamma, mle=fit$estimate)
pbootresults
boot.ci(pbootresults, type=c("basic", "perc", "norm"))
pbootresults<-boot(y100, gammamedian, R=1000, sim="parametric",     
                   ran.gen=gengamma, mle=fit$estimate)
pbootresults
 boot.ci(pbootresults, type=c("basic", "perc", "norm"))
pbootresults<-boot(y200, gammamedian, R=1000, sim="parametric",  
                   ran.gen=gengamma, mle=fit$estimate)
pbootresults
temp <- boot.ci(pbootresults, type=c("basic", "perc", "norm")) # store the results

datatable <- temp$percent %>% as.data.table()
datatable <- rbind(datatable, temp2$percent %>% as.data.table())
datatable[, group :=.I]

 ggplot(data = datatable) + geom_rect(mapping =  aes(xmin = group - 0.25, xmax = group+ 0.25, ymin = V4, ymax = V5))