Dati di modellazione con una funzione di collegamento Weibull in R

Question

Sto provando a modellare alcuni dati che seguono una relazione curva sigmoidale. Nel mio campo di lavoro (psicofisica), una funzione di Weibull viene solitamente utilizzata per modellare tali relazioni, piuttosto che il probit.

Sto provando a creare un modello utilizzando R e sto lottando con la sintassi. So che ho bisogno di utilizzare la funzione vglm() dal pacchetto VGAM, ma non riesco a ottenere un modello ragionevole. Ecco il mio dati:

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

Ecco un grafico dei dati in dframe1:

library(ggplot2) 

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

Questo dovrebbe essere in grado di essere modellato da una funzione di Weibull, dal momento che i dati in forma di un relazione curva sigmoidea. Qui è il mio tentativo di modellare i dati e generare una trama rappresentante:

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

Come potete vedere, questo non rappresenta miei dati originali a tutti. Sto generando il mio modello in modo errato o sto generando la mia trama del modello in modo errato. Che cosa sto facendo di sbagliato?

Nota: Ho modificato questa domanda per renderlo più comprensibile; in precedenza avevo utilizzato completamente la funzione sbagliata (weibreg()). Quindi, alcuni dei commenti sottostanti potrebbero non avere senso. .....

Answer 1

Ecco la mia soluzione, con bbmle.

dati:

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

Costruire una Weibull cumulativo che va 0,5-1,0 per definizione:

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

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

montare un Weibull (parametri rilevanti di log-in scala), con variazione binomio:

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

Generare previsioni:

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

È anche possibile utilizzare il pacchetto drc (modellazione dose-risposta).

Sono in realtà un noob per questo tipo di modelli, ma perhabs aiuta in qualche modo ...

Qui ho montato un quattro parametri Weibull, con parametri fissi per gli asintoti (altrimenti l'asintoto superiore sarebbe leggermente superiore 1, non so se questo è un problema per te). Dovevo anche trasformare la variabile indipendente (+0.2) in modo che fosse> = 0, a causa dei problemi di convergenza.

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

Tuttavia concordo con Ben Bolker sul fatto che altri modelli potrebbero essere più adatti.

Conosco solo questo tipo di modelli dall'ecotossicologia (modelli dose-risposta, dove si è interessati alla concentrazione in cui si ha una mortalità del 50% [= EC50]).

Aggiornamento Un modello log-logistico a quattro parametri si adatta anche abbastanza bene (più piccolo AIC e RSE poi Weibull): Ancora una volta ho fissato qui il parametro asintoto e trasformato il 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)) 

Answer 3

OK, ho appena imbattuto in questo parecchi mesi di ritardo, ma si potrebbe anche utilizzare il link mafc.cloglog dal pacchetto psyphy con glm. Se x segue il cloglog, il log (x) seguirà una funzione psicometrica weibull. Il fermo come con le risposte precedenti è che è necessario il numero di prove per la proporzione corretta. L'ho impostato su 100 in modo da fornire un numero intero di prove ma è necessario correggerlo in modo che corrisponda ai numeri effettivamente utilizzati da . Ecco il codice per farlo.

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)