r - 在 R 中使用 Weibull 链接函数对数据进行建模

Question

我正在尝试对一些遵循 sigmoid 曲线关系的数据进行建模。在我的工作领域（心理物理学）中，通常使用 Weibull 函数来模拟这种关系，而不是概率。

我正在尝试使用 R 创建一个模型并且正在努力解决语法问题。我知道我需要使用包中的vglm()功能VGAM，但我无法得到一个合理的模型。这是我的数据：

# Data frame example data
dframe1 <- structure(list(independent_variable = c(0.3, 0.24, 0.23, 0.16, 
0.14, 0.05, 0.01, -0.1, -0.2), dependent_variable = c(1, 1, 
1, 0.95, 0.93, 0.65, 0.55, 0.5, 0.5)), .Names = c("independent_variable", 
"dependent_variable"), class = "data.frame", row.names = c(NA, 
-9L))

这是 dframe1 中的数据图：

library(ggplot2)

# Plot my original data
ggplot(dframe1, aes(independent_variable, dependent_variable)) + geom_point()

这应该可以通过 Weibull 函数建模，因为数据符合 S 型曲线关系。这是我对数据建模并生成代表性图的尝试：

library(VGAM)

# Generate model
my_model <- vglm(formula = dependent_variable ~ independent_variable, family = weibull, data = dframe1)

# Create a new dataframe based on the model, so that it can be plotted
model_dframe <- data.frame(dframe1$independent_variable, fitted(my_model))

# Plot my model fitted data
ggplot(model_dframe, aes(dframe1.independent_variable, fitted.my_model.)) + geom_point()

如您所见，这根本不代表我的原始数据。我要么错误地生成了我的模型，要么错误地生成了我的模型图。我究竟做错了什么？

注意：我已编辑此问题以使其更易于理解；以前我完全使用了错误的功能（weibreg()）。因此，下面的一些评论可能没有意义。......

Answer 1

这是我的解决方案，使用bbmle.

数据：

dframe1 <- structure(list(independent_variable = c(0.3, 0.24, 0.23, 0.16, 
0.14, 0.05, 0.01, -0.1, -0.2), dependent_variable = c(1, 1, 
1, 0.95, 0.93, 0.65, 0.55, 0.5, 0.5)), .Names = c("independent_variable", 
"dependent_variable"), class = "data.frame", row.names = c(NA, 
-9L))

根据定义，构造一个从 0.5 到 1.0 的累积 Weibull：

wfun <- function(x,shape,scale) {
    (1+pweibull(x,shape,scale))/2.0
}

dframe2 <- transform(dframe1,y=round(40*dependent_variable),x=independent_variable)

拟合 Weibull（对数尺度相关参数），具有二项式变化：

library(bbmle)
m1 <- mle2(y~dbinom(prob=wfun(exp(a+b*x),shape=exp(logshape),scale=1),size=40),
     data=dframe2,start=list(a=0,b=0,logshape=0))

生成预测：

pframe <- data.frame(x=seq(-0.2,0.3,length=101))
pframe$y <- predict(m1,pframe)

png("wplot.png")
with(dframe2,plot(y/40~x))
with(pframe,lines(y/40~x,col=2))
dev.off()

Answer 2

好的，我只是晚了几个月才遇到这个问题，但是您也可以将 psyphy 包中的 mafc.cloglog 链接与 glm 一起使用。如果 x 遵循 cloglog，则 log(x) 将遵循 weibull 心理测量函数。与上述响应一样，您需要试验次数才能使比例正确。我只是将它设置为 100，因此它会提供整数次试验，但您应该将其修复为与您实际使用的数字相对应。这是执行此操作的代码。

dframe1 <- structure(list(independent_variable = c(0.3, 0.24, 0.23, 0.16, 
0.14, 0.05, 0.01, -0.1, -0.2), dependent_variable = c(1, 1, 
1, 0.95, 0.93, 0.65, 0.55, 0.5, 0.5)), .Names = c("independent_variable", 
"dependent_variable"), class = "data.frame", row.names = c(NA, 
-9L))

library(psyphy)

plot(dependent_variable ~ independent_variable, dframe1)
fit <- glm(dependent_variable ~ exp(independent_variable), 
    binomial(mafc.cloglog(2)), 
    data = dframe1, 
    weights = rep(100, nrow(dframe1)))  # assuming 100 observations per point
xx <- seq(-0.2, 0.3, len = 100)
pred <- predict(fit, newdata = data.frame(independent_variable = xx), type = "response")
lines(xx, pred)

Answer 3

您还可以使用 drc-package (dose-response-modeling)。

我实际上是这种模型的菜鸟，但也许它以某种方式有所帮助......

这里我拟合了一个四参数 Weibull，渐近线的参数是固定的（否则上渐近线会略大于 1，不知道这对你来说是否有问题）。由于收敛问题，我还必须转换自变量 (+0.2) 使其 >= 0。

require(drc)
# four-parameter Weibull with fixed parameters for the asymptote, added +0.2 to IV to overcome convergence problems
mod <- drm(dependent_variable ~ I(independent_variable+0.2), 
           data = dframe1, 
           fct = W1.4(fixed = c(NA, 0.5, 1, NA)))

# predicts
df2 <- data.frame(pred = predict(mod, newdata = data.frame(idenpendent_variable = seq(0, 0.5, length.out=100))), 
                  x = seq(0, 0.5, length.out=100))

ggplot() + 
  geom_point(data = dframe1, aes(x = independent_variable + 0.2, y = dependent_variable)) +
  geom_line(data = df2, aes(x = x, y = pred))

然而，我同意 Ben Bolker 的观点，即其他模型可能更适合。

我只知道生态毒理学中的这类模型（剂量反应模型，人们对我们有 50% 死亡率的浓度感兴趣 [=EC50]）。

更新 四参数对数逻辑模型也非常适合（更小的 AIC 和 RSE，然后是 weibull）：我再次在这里固定了渐近线参数并转换了 IV。

# four-parameter log-logistic with fixed parameters for the asymptote, added +0.2 to IV to overcome convergence problems
mod1 <- drm(dependent_variable ~ I(independent_variable+0.2), 
           data = dframe1, 
           fct = LL2.4(fixed=c(NA, 0.5, 1, NA)))
summary(mod1)

# predicts
df2 <- data.frame(pred = predict(mod1, newdata = data.frame(idenpendent_variable = seq(0, 0.5, length.out=100))), 
                  x = seq(0, 0.5, length.out=100))

ggplot() + 
  geom_point(data = dframe1, aes(x = independent_variable + 0.2, y = dependent_variable)) +
  geom_line(data = df2, aes(x = x, y = pred))