I have used the vglm() function from the VGAM package to obtain a Weibull model for some data. I would like to plot the curve that this model is represented by. How do I do this?

Here's my input and output (I am comparing 'p_seen' to 'calculated_logmar', which exist in the dataframe 'dframe1'):

> my_model <- vglm(formula = p_seen ~ calculated_logmar, family = weibull, data = dframe1)

Pearson Residuals:
  log(shape) log(scale)
1    0.78928  -0.093136
2   -1.21092  -0.978527
3    0.12519   0.020018

                     Estimate Std. Error   z value
(Intercept):1      1.4797e+01 1.0803e+00  13.69619
(Intercept):2      4.5824e-06 2.3095e-05   0.19841
calculated_logmar -1.5324e+01 1.4448e+00 -10.60582

Number of linear predictors:  2 

Names of linear predictors: log(shape), log(scale)

Dispersion Parameter for weibull family:   1

Log-likelihood: 9.46132 on 3 degrees of freedom

Number of iterations: 30 

I understand that the coefficients must relate somehow to the model itself, but I have no idea how to turn these coefficients into meaningful plots.

PS: I would like to use ggplot2 to visualise my data if possible.

Edit: I have tried using the command:

plotvgam(my_model, se = TRUE, lcol = dframe1$calculated_logmar[1], scol = dframe1$calculated_logmar[1], which.term = 1)

This gives me a plot of 'calculated_logmar' against 'partial for calculated_logmar'. What I'm looking for is the actual curve of the model, showing calculated_logmar plotted against p_seen. If it is possible, I'd like to render this using ggplot2. Any help would be appreciated.

  • $\begingroup$ It would be much easier to help if you will describe your variables, show real data or make a script to create a fictional data with given properties. $\endgroup$ – Yuriy Petrovskiy Feb 12 '13 at 16:17
  • $\begingroup$ No problem - see the image here; this is the same thing (stress is what I'm calling calculated_logmar, and probability of failure is what I'm calling p_seen). Otherwise, it's all the same >> opticalengineering.spiedigitallibrary.org/data/Journals/OPTICE/… ... Essentially I want to model those data, and plot on the same axes. $\endgroup$ – CaptainProg Feb 12 '13 at 16:41

I'm not sure if this is what you mean, but you can plot the fitted values from the model against the predictor variable in the normal way:

## Generate example data
n <- 100
x <- runif(n)
d <- data.frame(
  x = x, 
  y = rweibull(n, shape=exp(1+x), scale=exp(-0.5)))

## Fit the model and create a new data set with the fitted values
m <- vglm(y ~ x, data=d, family=weibull)
p <- data.frame(
  x = d$x, 
  y.fit = fitted(m))

## Create the plot using ggplot (and mimicking the style in OP's comment link) 
ggplot(p, aes(x, y.fit)) + geom_smooth() + geom_point(size=4, shape=1, alpha=0.6)

enter image description here

You might get a better response for this kind of question on SO, and my disclaimer is that I'm not familiar with this kind of model or the VGAM package, but hopefully this will at least point you in the right direction. A better CV answer would explain how all these things work...


According to package documentation there is a function called plotvgam for plotting vglm results:

fit2.ms <- vglm(mstatus ~ poly(age, 2) + foo(age) + age,
+ family = multinomial(refLevel = 2), constraints = clist,
+ data = nzmarital)
par(mfrow = c(2, 2))
plotvgam(fit2.ms, se = TRUE, scale = 12, lcol = mycol[1],
+ scol = mycol[1], which.term = 1)
plotvgam(fit2.ms, se = TRUE, scale = 12, lcol = mycol[2],
+ scol = mycol[2], which.term = 2)
plotvgam(fit2.ms, se = TRUE, scale = 12, lcol = mycol[3],
+ scol = mycol[3], which.term = 3)
  • $\begingroup$ Please see question edit above. Thanks for the tip, this has at least allowed me to plot something... $\endgroup$ – CaptainProg Feb 12 '13 at 14:29

Not the answer you're looking for? Browse other questions tagged or ask your own question.