visreg not working with geeglm function {in geepack package} or gee function {in gee package}
Closed this issue · 1 comments
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,
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.