How do I plot points and the model's line with glm gamma link=log?

Question

I used glm with Gamma distribution and link=log to create my regressions. I am trying to create the model's line with visreg, but it creates just for the response and it doesn't show the points. I want to create a graph with all my points and the model's line. Is it possible using visreg? If it isn't, Is there another package that I could use to do that?

I have already tried to create the model's line using 'abline' and it didn't work.

MODEL=glm(y~x, family=Gamma(link='log'), data)
visreg(MODEL, scale='response')

Answer 1

I am unaware of the visreg package implementation, but creating a line from scratch is not that difficult. Sadly abline wont work as the line wont be linear on the form y = ax + b. Instead we can use 4 small steps:

1. Create your fit
2. Sort your old data, or better yet generate some new data
3. Predict the model with your sorted or new data
4. Plot the prediction against the data

I assume the visreg package takes care of step 2-4.

Using the gamma model example from help(glm), I'll illustrate how one can achieve this with a minimal reproducible example:

Lets start by making a model.

#minimum reproducible example 
#From help(glm)
clotting <- data.frame(
    u = c(5,10,15,20,30,40,60,80,100),
    lot1 = c(118,58,42,35,27,25,21,19,18))
summary(fit1 <- glm(lot1 ~ log(u), data = clotting, family = Gamma(link = "log")))

Now that we've fit our model, we need to actually visualize the probability. For this i want some realistic data, so i just use 100 points, evenly spaced between the minimum and maximmum value of my variable u

rn <- range(clotting$u)
newdata <- data.frame(u = seq(rn[1], rn[2], length.out = 100)) #note the data.frame

To get some plot-able data, I use the predict function.

plotData <- xy.coords(x = newdata$u, y = predict(fit1, newdata = newdata))

All that is left is to plot your data.

plot(plotData, ylab = "log( lot1 | u )", main = "Model plot on link scale")

Now we've got a plot as the one below:

In order to add the points, one can use points. Note that you may have to change the xlim and ylim argument as i have done below (these set the limits on the x and y axis respectively).

plot(plotData, ylab = "log( lot1 | u )", main = "Model plot on link scale",
     ylim = range(clotting$lot1),
     xlim = rng)
points(clotting$u, clotting$lot1)

which yields the plot below:

Bonus: Plot on original scale

If you instead wanted to plot it in the value of lot1 (or y in your case) simply change

plotData <- xy.coords(x = newdata$u, y = predict(fit1, newdata = newdata))

to

plotData <- xy.coords(x = newdata$u, y = predict(fit1, newdata = newdata, type = "response"))