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

Posted 2023-02-19

【中文标题】在 R 中使用 Weibull 链接函数对数据进行建模

【英文标题】：Modelling data with a Weibull link function in R

【发布时间】：2013-01-24 12:05:12

【问题描述】：

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

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

# 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 函数建模，因为数据符合 sigmoid 曲线关系。这是我对数据建模并生成代表性图的尝试：

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()）。因此，下面的一些 cmets 可能没有意义。 .....

【问题讨论】：

我最初将您指向weibreg()，但这似乎是一个红鲱鱼。我很抱歉。 weibreg() 显然只处理 Weibull 回归生存模型（通常用 Weibull 建模） - 但心理物理学似乎是独一无二的，因为它们使用 Weibull 链接函数对非生存数据建模 其他人都会使用 logit 或 probit 的地方。但是，看起来VGAM 包中的vglm() 函数可能会起作用：rss.acs.unt.edu/Rdoc/library/VGAM/html/weibull.html 如果您可以将dput(dframe) 的输出添加到您的帖子中，我会尽力提供更多帮助。

感谢 Stephan，这对我来说是一次学习经历！我在我的问题中添加了“dput()”。任何有关如何运行该功能的建议将不胜感激。

好吧，我当然希望你有三个以上的观察结果！我猜您的 p 值来自多个观察结果，所以我建议您将它们全部放入数据框中。然后我将使用model <- vglm(p~size,family=weibull,data=dframe) 拟合模型（您需要告诉vglm() 什么是因变量，什么是自变量）并使用summary(model) 检查结果。您的警告消息意味着 ML 估计产生无效的形状参数；它可能会随着更多数据而消失。但我当然不会说我很了解vglm；也许其他人可以提供帮助？

好的，我可以从您的示例中看到，您的自变量似乎遵循累积威布尔形状。但是：观察值的统计特性是什么？它们是正态分布的吗？它们是否成比例，在这种情况下它们可能是 beta 分布的？需要知道这一点以适合统计模型...我查看了 cornea.berkeley.edu/pubs/148.pdf ，看起来您的数据可能是是/否比例？为了正确地做到这一点，我们可能需要分母（即每个点的试验次数）。

下渐近线是 0.5 而不是 1 似乎也很有趣……你能解释一下吗？

【参考方案1】：

好的，我只是晚了几个月才发现这个，但你也可以使用 带有 glm 的 psyphy 包中的 mafc.cloglog 链接。如果 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)

【参考方案2】：

这是我的解决方案，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()

【讨论】：

非常感谢本本。在我的一些试验中，我的演讲超过了 40 次。我面临以下选择：a) 忽略 40 日之后收集的数据，或 b) 修改“m1”的计算以考虑超过 40 次演示的试验。虽然它可能对结果几乎没有什么影响，但我想知道是否有办法合并这些额外的数据？我已经设法在最后一步合并了一个变量“n_presentations”，但不知道如何生成一个允许在每个数据中使用不同样本大小的 p_frame。

您当然应该能够考虑不同的样本量：只需确保上述模型中的 y 是成功次数，size 是实际试验次数（可以是向量，当然）。由于您正在尝试预测概率，我认为您可以将任何您想要的内容放入n_presentations。尝试n_presentations=1 的列，看看是否有效。否则，手动生成预测应该不会太难。

谢谢。使用mle2 中生成的模型预测“y”值时似乎出现了问题。如果我输入一个向量n_presentations 作为size= 参数，pframe$y <- predict(m1,pframe) 行不知道如何处理它。据推测，由于这条线试图从 9 个输入值中推断出 101 个点，因此它不知道每个点使用什么“大小”（即使 n_presentations 对每个数据都是“40”，这也会失败）......因为每个点的试验次数没有“趋势”，模型肯定不可能知道如何缩放 y 的每个值？

【参考方案3】：

您也可以使用 drc-package (dose-response-modelling)。

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

这里我拟合了一个四参数 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))