r - R のワイブル リンク関数を使用したデータのモデリング

Question

シグモイド曲線の関係に従ういくつかのデータをモデル化しようとしています。私の専門分野 (心理物理学) では、ワイブル関数は通常、プロビットではなく、そのような関係をモデル化するために使用されます。

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())。したがって、以下のコメントの一部は意味をなさない場合があります。.....

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 になる累積ワイブルを作成します。

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

Answer 2

OK、数か月遅れてこの問題に遭遇しましたが、psyphy パッケージの mafc.cloglog リンクを glm で使用することもできます。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)

Answer 3

drc-package (dose-response-modelling) を使用することもできます。

私は実際にはこの種のモデルの初心者ですが、おそらくそれが何らかの形で役立ちます...

ここでは、漸近線の固定パラメーターを使用して、4 つのパラメーターのワイブルを適合させました (そうしないと、上限の漸近線はわずかに 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] になる濃度に関心がある場合)。

更新 4 パラメーターの対数ロジスティック モデルも非常によく適合します (小さい AIC と RSE、次にワイブル): ここでも漸近線パラメーターを修正し、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))