使用R中的Weibull链接函数对数据建模

Question

我试图模拟一些遵循S形曲线关系的数据.在我的工作领域(心理物理学)中,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

好的，我几个月后才发现这个问题，但您也可以使用 glm 的 psyphy 包中的 mafc.cloglog 链接。如果 x 遵循 cloglog，则 log(x) 将遵循威布尔心理测量函数。与上述响应相同的问题是，您需要正确比例的试验次数。我只是将其设置为 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)