pbreheny / visreg

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

Error plotting glmmTMB model with weights argument #67

Closed erinsauer closed 1 year ago

erinsauer commented 5 years ago

First, thanks for making an awesome package and keeping up with user issues! I'm having some trouble plotting a glmmTMB model with a weights argument. When I remove the weights argument from model and run the same code, the visreg code works fine. There is no missing data in the in.var column. Thanks!

ENDmodel<-glmmTMB(ln.or~host.bio1*Temp+logdose+lifestage+
              (1|genus.and.species)+( 1|ID)+ (1|bd.ID),
                family="gaussian", data=END, weights=in.var, 
                dispformula = ~1, doFit=TRUE)
 visreg(ENDmodel, "Temp", by="host.bio1.10", breaks=2,   
        points=list(col="black", cex=1.5, pch=21), 
        line.par=list(lty=2, col="black", lwd=3), 
        strip=F, par.strip.text = list(col = "transparent"), 
        xlab="Laboratory temperature (C)", 
        ylab="Mortality relative to controls (ln odds ratio)", 
        scales=list(y=list(at=c(0,1,2,3,4,5),limits=c(0, 5)), 
        alternating=FALSE, tck=c(1,0)), 
        jitter=TRUE, fill.par=list(col=c("grey")))
Error in model.frame.default(weights = in.var, drop.unused.levels = TRUE,  : 
  variable lengths differ (found for '(weights)')
pbreheny commented 5 years ago

Hi Erin,

This would seem to be a problem with glmmTMB, not with visreg:

library(glmmTMB)
w <- runif(nrow(Salamanders))
fit <- glmmTMB(count~ mined + (1|site), zi=~mined, family=poisson, data=Salamanders, weights=w)
predict(fit, newdata = Salamanders[1:5,])

produces

Error in model.frame.default(weights = w, drop.unused.levels = TRUE, formula = ~mined +  : 
  variable lengths differ (found for '(weights)')

In other words, it seems that predict.glmmTMB() does not work if the model was fit with weights. You should let the maintainer of the glmmTMB know about this.

erinsauer commented 5 years ago

Hi Patrick,

I actually tested predict before posting this issue, thinking that this might be the problem. But when I use predict with my glmmTMB model it seems to work fine (no error message). You may be getting that error because w is not part of your data frame so it has a different length than Salamanders[1:5,] which is only 5 rows.

I tried running your code without limiting the length of Salamanders and predict worked fine. However, when I try plotting fit in visreg I get they same error I mentioned in my initial comment.

Thanks for the speedy response!

Erin

pbreheny commented 5 years ago

Hmm, interesting. The issue is that visreg needs to be able to make predictions not just for rows of the data that appear in the data frame, but for hypothetical rows that could appear in the data frame -- this is the line being plotted. It seems to me that there are two possibilities:

  1. predict.glmmTMB doesn't actually need weights in order to make predictions. If this is the case, it shouldn't output an error if weights aren't supplied.
  2. predict.glmmTMB does need weights in order to make predictions. If this is the case, it seems a little hopeless since visreg isn't going to know how to set up those weights.

I guess there is perhaps a third possibility:

  1. predict.glmmTMB could potentially depend on weights, but a reasonable default would just be to give everything a default weight of 1 if weights aren't supplied.

I have no idea what glmmTMB is, so I don't know which one of these is correct, but I guess I'm sort of inclined to stand by my original conclusion that predict.glmmTMB should still work if weights aren't supplied, and that's something that should be fixed on the glmmTMB end.

All they would have to do is add the line

mf$weights <- rep(1, nrow(mf$data))

with perhaps some sort of check to only do this if weights aren't supplied or something.

As a workaround, not sure how you feel about this, but you can write your own predict.glmmTMB function. Just copy the existing one, add one line:

mf$na.action <- na.action   # This line exists already
mf$weights <- rep(1, nrow(mf$data))   # This is the line you're adding

then run

environment(predict.glmmTMB) <- asNamespace('glmmTMB')

And then visreg() works fine.

mattjohnsonhsu commented 1 year ago

I get the same error with visreg using a simple glm model with weights (not glmmTMB):

forest_rsfdeermcp <- glm(case~forest, family = "binomial", data = ua_95, weights = weight) visreg(forest_rsf_deermcp) Error in eval(extras, data, env) : object 'weight' not found

pbreheny commented 1 year ago

Is there a variable called weight in the ua_95 data set? If so, the code will work. For example:

library(visreg)
df <- data.frame(y=rbinom(1000, 1, 0.5), x=runif(1000), w=1+rpois(1000, 5))
fit <- glm(y~x, data=df, family=binomial, weights=w)
visreg(fit)
mattjohnsonhsu commented 1 year ago

Yes, there is a variable called weight.

Fitting the mode works, but visreg does not.

Matt


Matt Johnson Professor, Wildlife Habitat Ecology Webpage https://wildlife.humboldt.edu/people/matthew-johnson Director, HSI-STEM grant https://hsistem.humboldt.edu/hsi-stem-grant Co-PI, Échale Ganas USDA grant https://echaleganas.humboldt.edu/ Department of Wildlife Cal Poly Humboldt (formerly Humboldt State University) Arcata, CA 95521 Phone: (707) 826-3218 FAX : (707) 826-4060 email: @.***


On Mon, Aug 1, 2022 at 1:54 PM Patrick Breheny @.***> wrote:

Is there a variable called weight in the ua_95 data set? If so, the code will work. For example:

library(visreg)df <- data.frame(y=rbinom(1000, 1, 0.5), x=runif(1000), w=1+rpois(1000, 5))fit <- glm(y~x, data=df, family=binomial, weights=w) visreg(fit)

— Reply to this email directly, view it on GitHub https://github.com/pbreheny/visreg/issues/67#issuecomment-1201709527, or unsubscribe https://github.com/notifications/unsubscribe-auth/A2APW5CODGBRTGN5GAZVURTVXA2O7ANCNFSM4G3TC3YQ . You are receiving this because you commented.Message ID: @.***>

pbreheny commented 1 year ago

Without a reproducible example, there's nothing I can do to help. The code itself would seem to work -- see my example -- so presumably the problem lies in some code/data you haven't provided.

mattjohnsonhsu commented 1 year ago

Thanks for the quick reply.

I'll distill the data & code down to something I can send that's reproducible. One moment...

Matt


Matt Johnson Professor, Wildlife Habitat Ecology Webpage https://wildlife.humboldt.edu/people/matthew-johnson Director, HSI-STEM grant https://hsistem.humboldt.edu/hsi-stem-grant Co-PI, Échale Ganas USDA grant https://echaleganas.humboldt.edu/ Department of Wildlife Cal Poly Humboldt (formerly Humboldt State University) Arcata, CA 95521 Phone: (707) 826-3218 FAX : (707) 826-4060 email: @.***


On Wed, Aug 3, 2022 at 3:28 PM Patrick Breheny @.***> wrote:

Without a reproducible example https://myweb.uiowa.edu/pbreheny/reproducible.html, there's nothing I can do to help. The code itself would seem to work -- see my example -- so presumably the problem lies in some code/data you haven't provided.

— Reply to this email directly, view it on GitHub https://github.com/pbreheny/visreg/issues/67#issuecomment-1204547994, or unsubscribe https://github.com/notifications/unsubscribe-auth/A2APW5EKCDAHP6SSQMQ7VDTVXLXCRANCNFSM4G3TC3YQ . You are receiving this because you commented.Message ID: @.***>

mattjohnsonhsu commented 1 year ago

Hi Patrick,

Attached is the data.

And here's the code distilled to the visreg error I can't seem to figure out. Any help would be appreciated.

Thanks

reproducible visreg error with weights

install.packages("glmmTMB") library(glmmTMB) library(visreg)

goats.data <- read.csv('data/goats.data_scaled_wtd.csv') gm.1_rand_slopeELEV_wt <- glmmTMB(STATUS ~ poly(ELEVATION,2)+SLOPE+ET+poly(HLI,2) + (1|ID) + (0+ELEVATION |ID), family=binomial(), data = goats.data, weights = weight) summary(gm.1_rand_slopeELEV_wt)

visreg(gm.1_rand_slopeELEV_wt, scale="response")

###########################################

Thanks,

Matt


Matt Johnson Professor, Wildlife Habitat Ecology Webpage https://wildlife.humboldt.edu/people/matthew-johnson Director, HSI-STEM grant https://hsistem.humboldt.edu/hsi-stem-grant Co-PI, Échale Ganas USDA grant https://echaleganas.humboldt.edu/ Department of Wildlife Cal Poly Humboldt (formerly Humboldt State University) Arcata, CA 95521 Phone: (707) 826-3218 FAX : (707) 826-4060 email: @.***


On Wed, Aug 3, 2022 at 3:28 PM Patrick Breheny @.***> wrote:

Without a reproducible example https://myweb.uiowa.edu/pbreheny/reproducible.html, there's nothing I can do to help. The code itself would seem to work -- see my example -- so presumably the problem lies in some code/data you haven't provided.

— Reply to this email directly, view it on GitHub https://github.com/pbreheny/visreg/issues/67#issuecomment-1204547994, or unsubscribe https://github.com/notifications/unsubscribe-auth/A2APW5EKCDAHP6SSQMQ7VDTVXLXCRANCNFSM4G3TC3YQ . You are receiving this because you commented.Message ID: @.***>

pbreheny commented 1 year ago
  1. First of all, contrary to your claim, this works fine with glm (your data, your model, just replacing glmmTMB with glm).
  2. I still think that fundamentally, this is a problem with predict.glmmTMB() in the sense that it should allow the user to make predictions for new observations without requiring weights.
  3. However, I just committed a change that allows an easy workaround to this limitation in glmmTMB. You can directly set the weights for the new predictions to be 1 (or something else, I guess):
    visreg(gm.1_rand_slopeELEV_wt, scale="response", cond=list(weight=1))

    This should now function as you want it to (keep in mind that you need to run remotes::install_github("pbreheny/visreg") first).

This functionality will appear in visreg 2.7.1.

mattjohnsonhsu commented 1 year ago

Hi Patrick,

I had a min back at my computer and ran what you suggested -- worked great! Many thanks!

Cheers,

Matt


Matt Johnson Professor, Wildlife Habitat Ecology Webpage https://wildlife.humboldt.edu/people/matthew-johnson Director, HSI-STEM grant https://hsistem.humboldt.edu/hsi-stem-grant Co-PI, Échale Ganas USDA grant https://echaleganas.humboldt.edu/ Department of Wildlife Cal Poly Humboldt (formerly Humboldt State University) Arcata, CA 95521 Phone: (707) 826-3218 FAX : (707) 826-4060 email: @.***


On Thu, Aug 4, 2022 at 10:05 AM Patrick Breheny @.***> wrote:

  1. First of all, contrary to your claim, this works fine with glm (your data, your model, just replacing glmmTMB with glm).
  2. I still think that fundamentally, this is a problem with predict.glmmTMB() in the sense that it should allow the user to make predictions for new observations without requiring weights.
  3. However, I just committed a change that allows an easy workaround to this limitation in glmmTMB. You can directly set the weights for the new predictions to be 1 (or something else, I guess):

visreg(gm.1_rand_slopeELEV_wt, scale="response", cond=list(weight=1))

This should now function as you want it to (keep in mind that you need to run remotes::install_github("pbreheny/visreg") first).

This functionality will appear in visreg 2.7.1.

— Reply to this email directly, view it on GitHub https://github.com/pbreheny/visreg/issues/67#issuecomment-1205533100, or unsubscribe https://github.com/notifications/unsubscribe-auth/A2APW5DBEPNBMGKSR635XTDVXPZ6VANCNFSM4G3TC3YQ . You are receiving this because you commented.Message ID: @.***>

milkweedgirl commented 8 months ago

@pbreheny - I am having this same issue and did the remote github install but noticed my version is showing as 2.7.0.8 instead of 2.7.1 - would this be limiting me being able to use the cond=list function? When I try to specify my weights with this I am still getting the object not found error.

Example code: Visreg(m1, 'suit', by='cage', overlay=FALSE, type="conditional",xlab="suitability", rug=FALSE,partial=TRUE, scale="response", cond=list(total=25))

**total is my weights category

thank-you!

pbreheny commented 8 months ago

2.7.1 doesn't exist yet; 2.7.0.8 is the most current version. Can you provide a reproducible example?

milkweedgirl commented 7 months ago

Hi Patrick- could I email you the dataset I am working with to show you my issue?

thank-you, Emma

my school email: @.***

On Mon, Nov 6, 2023 at 9:46 AM Patrick Breheny @.***> wrote:

2.7.1 doesn't exist yet; 2.7.0.8 is the most current version. Can you provide a reproducible example https://myweb.uiowa.edu/pbreheny/reproducible.html?

— Reply to this email directly, view it on GitHub https://github.com/pbreheny/visreg/issues/67#issuecomment-1795461329, or unsubscribe https://github.com/notifications/unsubscribe-auth/BDY3YPGSLAGAPEM3MS3DCJTYDEH7TAVCNFSM4G3TC3Y2U5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TCNZZGU2DMMJTGI4Q . You are receiving this because you commented.Message ID: @.***>

pbreheny commented 7 months ago

Sure