Modellierung von Daten mit einer Weibull-Link-Funktion in R

Question

Ich versuche, einige Daten zu modellieren, die einer sigmoiden Kurvenbeziehung folgen. In meinem Arbeitsgebiet (Psychophysik) wird normalerweise eine Weibull-Funktion verwendet, um solche Beziehungen zu modellieren, anstatt Probit.

Ich versuche ein Modell mit R zu erstellen und kämpfe mit Syntax. Ich weiß, dass ich die vglm() Funktion aus dem VGAM Paket verwenden muss, aber ich bin nicht in der Lage, ein vernünftiges Modell heraus zu bekommen. Hier ist meine Daten:

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

Hier ist eine grafische Darstellung der Daten in dframe1:

library(ggplot2) 

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

Dies sollte durch eine Weibull-Funktion modelliert werden können, da die Daten ein passen sigmoide Kurvenbeziehung. Hier ist mein Versuch, die Daten zu modellieren und generieren eine repräsentative Darstellung:

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

Wie Sie sehen können, ist dies nicht meine ursprünglichen Daten überhaupt darstellen. Ich erzeuge entweder mein Modell falsch oder ich erzeuge mein Diagramm des Modells falsch. Was mache ich falsch?

Hinweis: Ich habe diese Frage bearbeitet, um sie verständlicher zu machen; zuvor hatte ich die falsche Funktion komplett benutzt (weibreg()). Daher sind einige der folgenden Kommentare möglicherweise nicht sinnvoll. .....

Answer 1

Hier ist meine Lösung, mit bbmle.

Daten:

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

eine kumulative Weibull-Konstrukt, das von 0,5 bis 1,0 per Definition lautet:

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

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

Fit ein Weibull (log-Skala relevanten Parameter), mit binomischen Variation:

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

Vorhersagen generieren:

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

Sie könnten auch das DRC-Paket (Dosis-Antwort-Modellierung) verwenden.

Ich bin eigentlich ein noob für diese Art von Modellen, aber es irgendwie hilft Eventuell kann ...

Hier I ausgerüstet ist ein Vier-Parameter-Weibull, mit festen Parametern für die Asymptoten (sonst wird die obere Asymptote etwas größer sein würde, 1, weiß nicht, ob dies ein Problem für dich ist). Ich musste auch die unabhängige Variable (+0.2) so transformieren, dass sie wegen Konvergenzproblemen> = 0 ist.

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

Allerdings stimme ich mit Ben Bolker, dass andere Modelle besser geeignet sein können.

Ich kenne nur diese Art von Modellen aus Ökotoxikologie (Dosis-Wirkungs-Modelle, wo man an der Konzentration interessiert ist, wo wir 50% Mortalität haben [= EC50]).

aktualisieren Ein Vier-Parameter log-logistischen Modell passt auch ganz gut (kleiner AIC und RSE dann weibull): Wieder festen ich hier die Asymptote Parameter und die IV umgewandelt.

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

Answer 3

OK, ich kam gerade über diese mehrere Monate zu spät, aber Sie auch den mafc.cloglog Link aus dem psyphy Paket mit GLM nutzen könnten. Wenn das x dem cloglog folgt, folgt log (x) einer weiblich-psychologischen Funktion. Der Haken wie bei den obigen Antworten ist , dass Sie die Anzahl der Versuche für den richtigen Anteil benötigen. Ich habe es nur auf 100 gesetzt, so dass es eine ganze Zahl von Versuchen geben würde, aber Sie sollten das beheben, um den Zahlen zu entsprechen, die Sie tatsächlich verwenden. Hier ist der Code, um es zu tun.

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)