Visualising data from glm/ glmer

Question

Sorry if this question sounds very amateurish, but I am completely new to R and this website. So, I have a data set with counts of infected organisms as a response variable (for three different species) and 5 different predictor variables. The infected organism count data is from over 200 surveys where they did not collect the same number of organisms each time. Some variables were measured for the entire site (constant across all surveys), some are unique to each survey.

An example of how the data is structured:

So, I used the cbind() when creating the glm to account for the fact that this is proportion data.

Disease_glm <- glm(cbind(n_infected, not_infected) ~ 
    var1 + var2 + var3 + var4 + var5, family = "binomial", 
    data = Disease)

Then I used DHARMa to test for dispersion

testDispersion(simulateResiduals(Disease_glm))

There was significant overdispersion so I created a unique ID for each survey

Disease$survey_ID = c(1:dim(Disease)[1])

So the data now has a survey ID column tucked in at the end.

And then I used it as a random effect in my model:

Disease_glmer <- glmer(cbind(n_infected, not_infected) ~ 
    var1 + var2 + var3 + var4 + var5 + (1|survey_ID), 
    family = "binomial", data = Disease)

This solved the overdispersion issue.

I just wanted to visualise the results for each species. How much impact does var1, var2 ... have on infection numbers in Species 1, for instance.

This is where major issues started cropping up. At first, I tried

ggplot(data =Disease D) + 
    geom_point(mapping = aes(x = var1+var2+var3+var4+var5, 
    y = (cbind(n_infected, not_infected)),  color = species))

But I got the error message "Error in check_aesthetics(): ! Aesthetics must be either length 1 or the same as the data (225): y"

Since I am so new that I do not know how to interpret this at all, I just abandoned it and tried a different route:

par(mfrow=c(1, 2))
with(Disease, plot((cbind(n_infected, 
    not_infected)) ~ var1 + var2 + var3 + var4 + var5, pch=16, 
    color = MergedChytrid$species))
curve(exp(cbind(1,x)%*%coef(chytrid_glmer)), add=T, col = "goldenrod", 
    lwd=3)

The first bit gives me the error:

Error in xy.coords(x, y, xlabel, ylabel, log) : 'x' and 'y' 
lengths differ

Trying to run the curve function gave me:

Error in cbind(1, x) %*% coef(chytrid_glmer) : requires 
numeric/complex matrix/vector arguments

Could you please help me understand where I am going wrong? It seems like such a basic thing. Any help would be appreciated

Answer 1

I'm not sure about your use of a random effect model. You create a group for each observation and use that as a random effect... I don't see the reason for doing that.

Regarding the overdispersion, you could account for it by using the quasibinomial distribution family. Maybe an alternative could be to use the negative binomial distribution (glm.nb in MASS package) and set n_surveyed as offset term.

How much impact does var1,var2... have on infection numbers in Species 1, for instance.

It seems to me that you should have Species in your model interacting with each of the 5 variables (perhaps also Area should be in the model). Like:

Disease_glm <- glm(cbind(n_infected, not_infected) ~ var1*Species + var2*Species + var3*Species + var4*Species + var5*Species + Area, family = "quasibinomial", data = Disease)

Then the table of coefficents (from summary(Disease_glm)) will give you the estimates for each species and var relative to the intercept term. It may be difficult to read and the package emmeans could be handy.