Modelling data with a Weibull link function in R

Question

I am trying to model some data that follows a sigmoid curve relationship. In my field of work (psychophysics), a Weibull function is usually used to model such relationships, rather than probit.

I am trying to create a model using R and am struggling with syntax. I know that I need to use the vglm() function from the VGAM package, but I am unable to get a sensible model out. Here's my data:

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

Here is a plot of the data in dframe1:

library(ggplot2)

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

This should be able to be modelled by a Weibull function, since the data fit a sigmoid curve relationship. Here is my attempt to model the data and generate a representative plot:

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

As you can see, this doesn't represent my original data at all. I'm either generating my model incorrectly, or I'm generating my plot of the model incorrectly. What am I doing wrong?

Note: I have edited this question to make it more understandable; previously I had been using the wrong function entirely (weibreg()). Hence, some of the comments below may not make sense. .....

Answer 1

OK, I just came across this several months late, but you could also use the mafc.cloglog link from the psyphy package with glm. If the x follows the cloglog then log(x) will follow a weibull psychometric function. The catch as with the above responses is that you need the number of trials for the proportion correct. I just set it to 100 so it would give an integer number of trials but you should fix this to correspond to the numbers that you actually used. Here is the code to do it.

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)