我可以走大部分路，但不是一路走来。主要步骤是为dr4pl对象写一个predict()方法：
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() 将其搞砸了，因为它是在 拟合和预测之前应用的，所以我改用 coord_trans(x="log10")。
- 但是，如果坐标轴在非常宽的对数刻度上，
coord_trans() 的效果就不那么好了 - 我使用包中的 sample_data_1 数据尝试了上面的示例，但它不起作用。
但恐怕我现在已经在这方面花费了足够的时间。
使用上面的predict 方法在你想要的范围内分别生成你想要的值会更健壮，然后使用geom_line() + geom_ribbon() 将信息添加到绘图中...... .
如果您愿意先拟合模型（在geom_smooth 之外），您可以这样做（这是使用来自dr4pl 包的sample_data_1 - 它来自?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")
)