Моделирование данных с помощью функции связи Вейбулла в R

Я пытаюсь моделировать некоторые данные, которые следует за отношением сигмовидной кривой. В моей области работы (психофизика) функция Вейбулла обычно используется для моделирования таких отношений, а не пробита.

Я пытаюсь создать модель, использующую 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()

Это должно быть смоделировано с помощью функции Вейбулла, поскольку данные соответствуют сигмовидной кривой. Вот моя попытка смоделировать данные и создать репрезентативный сюжет:

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()). Следовательно, некоторые из комментариев ниже могут не иметь смысла. .....

Ответы

Ответ 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))

Постройте кумулятивный Weibull, который по определению составляет от 0.5 до 1.0:

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

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

Установите вейбулл (соответствующие параметры логарифмической шкалы) с биномиальным изменением:

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()

Ответ 2

Вы также можете использовать drc-пакет (моделирование доза-реакция).

Я на самом деле нуб для таких моделей, но perhabs это помогает как-то...

Здесь я установил четыре параметра 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))

Однако я согласен с Бен Болкером в том, что другие модели могут быть лучше подходят.

Я знаю только такие модели из экотоксикологии (модели дозозависимости, где вас интересует концентрация, где у нас 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))

Ответ 3

Хорошо, я только что настал через несколько месяцев, но вы также можете использовать ссылка mafc.cloglog из пакета psyphy с 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)