我可以走大部分路，但不能走完全程。主要步骤是写一个predict()方法用于dr4pl对象：
predict.dr4pl <- function (object, newdata=NULL, se.fit=FALSE, level, interval) {
xseq <- if (is.null(newdata)) object$data$Dose else newdata$x
pred <- MeanResponse(xseq, object$parameters)
if (!se.fit) {
return(pred)
}
qq <- qnorm((1+level)/2)
se <- sapply(xseq,
function(x) car::deltaMethod(object,
"UpperLimit + (LowerLimit - UpperLimit)/(1 + (x/IC50)^Slope)")[["Estimate"]])
return(list(fit=data.frame(fit=pred,lwr=pred-qq*se,
upr=pred+qq*se), se.fit=se))
}
我提供了一种通过 delta 方法计算置信区间的稍微有点古怪的方法 - 这可能不太可靠（引导会更好......）
它对您的数据工作正常（有点）（将名称更改为dd因为有时给数据命名很冒险data (fortunes::fortune("dog"))).
dd <- data.frame(dose = c(0.078125,0.156250,0.312500,0.625000,1.25,
2.50,5.0,10.0,20.0),
POC = c(1.05637425, 0.87380081, 0.79171200,
0.83166848, 0.77361290, 0.35199288,
0.19404609, 0.09079221, 0.09850658))
library(dr4pl)
ggplot(dd, aes(dose,POC)) + geom_point() +
geom_smooth(method="dr4pl",se=TRUE) + coord_trans(x="log10")
- 置信区间很糟糕，将其关闭
se=FALSE
-
dr4pl默认情况下，将 x 轴置于 log10 刻度上，但标准scale_x_log10()搞砸了，因为它已被应用before拟合和预测，所以我用coord_trans(x="log10")反而。
- 然而，
coord_trans()如果轴处于非常宽的对数刻度上，则效果不太好 - 我尝试了上面的示例sample_data_1数据包中的数据，但它不起作用。
但恐怕我现在已经在这上面花了足够的时间了。
使用该方法会更加稳健predict上面的方法分别生成您想要的值，在您想要的范围内，然后使用geom_line() + geom_ribbon()将信息添加到绘图中......
如果您愿意先拟合模型（在外部geom_smooth）你可以这样做（这是使用sample_data_1 from dr4pl包 - 它来自第一个示例?dr4pl)
model2 <- dr4pl(dose = sample_data_1$Dose,
response = sample_data_1$Response)
ggplot(sample_data_1, aes(Dose,Response)) + geom_point() +
stat_function(fun=function(x) predict(model2,newdata=data.frame(x=x))) +
scale_x_log10()
它对 x 轴缩放/取消缩放的顺序不太敏感。
改进但缓慢的引导 CI：
predictdf.dr4pl <- function (model, xseq, se, level, nboot=200) {
pred <- MeanResponse(xseq, model$parameters)
if (!se) {
return(base::data.frame(x=xseq, y=pred))
}
## bootstrap residuals
pred0 <- MeanResponse(model$data$Dose, model$parameters)
res <- pred0-model$data$Response
bootres <- matrix(nrow=length(xseq), ncol=nboot)
pb <- txtProgressBar(max=nboot,style=3)
for (i in seq(nboot)) {
setTxtProgressBar(pb,i)
mboot <- dr4pl(model$data$Dose,
pred0 + sample(res, size=length(pred0),
replace=TRUE))
bootres[,i] <- MeanResponse(xseq, mboot$parameters)
}
fit <- data.frame(x = xseq,
y=pred,
ymin=apply(bootres,1,quantile,(1-level)/2),
ymax=apply(bootres,1,quantile,(1+level)/2))
return(fit)
}
print(ggplot(dd, aes(dose,POC))
+ geom_point()
+ geom_smooth(method="dr4pl",se=TRUE) + coord_trans(x="log10")
)
