I have an interaction term between season and temperature in my final model (both 4-level categorical), and I want to graph this to show the probability of my outcome variable (a dog testing positive for a disease - dichotomous Y/N) for each category of temperature during each season.

I used melogit for the model with random intercepts for state and county. I then used the margins command and marginsplot to graph the interaction:

melogit Positive i.AgeC GenderCC i.RegionCC i.Year i.DogDensC i.UICC i.PDSIC i.PrecipC i.SeasonC##i.TempC || CustState: || FIPS:

margins TempC, at(Season=(0(1)3))

marginsplot, by(SeasonC)

This produced graphs with the marginal predicted means on the y-axes (see below). I had read in another topic thread within this forum that these marginal predicted means are in fact the probabilities? Is that true? If so, then why are these y-axes starting at -0.05? Is that just because of the large confidence interval on the season=3 graph? (note: my dataset has 40,118 observations, but there are 0 observations when temp=0 and season=3 and only 32 when temp=1 and season=3. That is why the graph for season=3 looks the way it does.)

If the marginal predicted means are not the probabilities, then is it possible to estimate probabilities to graph the interaction given this model?

Array