search

pbreheny / visreg

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

visreg not working with geeglm function {in geepack package} or gee function {in gee package} #100

Closed adusabimana closed 2 years ago

adusabimana commented 2 years ago

Hello Pbreheny, Thanks for bring this nice package. However it is not working for me. I have been checking the issues discussed but I could not find how to fix it for geeglm. I am fitting the following gee model: gee.log.poisson <- geeglm(formula = outcome ~ var1 + var2 + var3, data = data, family = poisson(link = "log"), id = interaction(site,id), corstr = "exchangeable")

When I use visreg, I get this error :

Error in XRinv^2 %*% rep(res.var, p) : requires numeric/complex matrix/vector arguments

gee.log.poisson <- gee(formula = outcome ~ var1 + var2 + var3, data = data, family = poisson(link = "log"), id = interaction(site,id), corstr = "exchangeable") When I use visreg I get this error Error in if (object$df.residual > 0) { : argument is of length zero

Do have an idea how I can handle this? Thanks,

pbreheny commented 2 years ago

Sorry for the delay in getting back to you. The visreg package relies on printing the output from a model's predict() method. In the case of geeglm() and gee(), neither of these models provides such a method:

# Setup
library(geepack)
library(gee)
data(dietox)
dietox$Cu     <- as.factor(dietox$Cu)
mf <- formula(Weight ~ Cu * (Time + I(Time^2) + I(Time^3)))

# gee
fit <- gee(breaks ~ tension, id=wool, data=warpbreaks, corstr="exchangeable")
predict(fit, dietox[3:4,], se.fit=TRUE)
# Error in sqrt(dispersion) : non-numeric argument to mathematical function

# geepack
fit <- geeglm(mf, data=dietox, id=Pig, family=poisson("identity"), corstr="ar1")
predict(fit, dietox[3:4,], se.fit=TRUE)
# Error in XRinv^2 %*% rep(res.var, p) : requires numeric/complex matrix/vector arguments

Fundamentally, this is more of a gee and geepack problem than a visreg problem, at least in my opinion: the authors of these packages should provide a working predict() function. Without that, there's not a whole lot one can do. One (not that great) workaround would be to simply use predict.glm() for prediction:

predict.geeglm <- function(object, newdata, ...) {predict.glm(object, newdata)}
visreg(fit, 'Time')

This can plot the fitted model and partial residuals, but doesn't work for standard errors/confidence bands.