pbreheny / visreg

Visualization of regression functions
http://pbreheny.github.io/visreg/
61 stars 18 forks source link

Partial residual plots and beta regression #122

Closed luroy closed 1 month ago

luroy commented 2 months ago

Dear all, I have just started to use visreg package to perform partial residual plots. I therefore apologize in advance if my following comment is due to my unfamiliarity with this function. I was wondering if it was possible (and correct) to obtain a partial residuals plot from a beta regression model. Initially, I ran my model using the betareg package as follows:

mod1<-betareg(y ~ scale(x1) + scale(x2) + scale (x3), data=data, link="logit", na.action = na.fail) But I couldn't get a graph of the partial residuals: visreg(mod1, "x1")

Error in plot.window(...) : valeurs finies requises pour 'ylim'
In addition: Warning messages:
1: In plot.visreg(v, ...) :
  The generic function residuals() is not set up for this type of model object.  To plot partial residuals, you will need to define your own residuals.betareg() function.
2: In min(x) : no non-missing arguments to min; returning Inf
3: In max(x) : no non-missing arguments to max; returning –Inf

I wonder if this isn't due to the fact that the residuals() function automatically calculates sweighted2 residuals for a betareg model Consequently, I've written the same model using the gam() function, and specifying family="betar" and link="logit mod2<-gam(y ~ scale(x1) + scale(x2) + scale (x3), data=data, family="betar", link="logit", na.action = na.fail). This seems to be equivalent to mod1, since I get very close estimates and p-values. In this case, I obtain a partial residuals plot from the following function: visreg(mod2, "x1") Firstly, I have a silly question: is it normal for the abline to take the form of a linear regression? Secondly, why do I get a graph of partial residuals from mod2, but not from mod1? Thank you for your attention to this matter, All the best, Léa

pbreheny commented 2 months ago

Hi Léa,

TLDR: don't use scale()

The issue is that scale() is not a genuine transformation in the sense of log/sqrt/etc: it depends on the data that is included. In particular, visreg depends on calling predict() using different sets of data (changing x1 and fixing the rest, etc.). When you run

visreg(mod1, "x1")

visreg fixes x2 and x3; since they are constant, it is impossible to scale them and we just get a bunch of NaN values when predict() is called.

I should write some code that checks for the use of scale() and communicates this to the user.

You have two options (use either one and everything should work fine):

To answer your other questions:

Q: Is it normal for the abline to take the form of a linear regression?

A: Depends on the scale, but yes -- all regression models have this property (unless you add nonlinear terms)

Q: Why do I get a graph of partial residuals from mod2, but not from mod1?

A: I have no idea why mod2 works; neither one should work since multiple uses of scale() should cause everything to blow up.

luroy commented 1 month ago

@pbreheny Thanks a lot!