Return to Answer

You were on the right track in your original post Partitioning variance from logistic regression, using glmer() for mixed-effects logistic regression.

You were on the right track in your original post Partitioning variance from logistic regression, using glmer() for mixed-effects logistic regression.

I would recommend against: the advisor's "solution", using lm() for logistic regression, and weighting rows equally (weight by N_indiv).

# Seroprevalance has been rounded, that's not OK
# to do binomial regression, (proportion * weight) must equal an integer
prev$seroexact <- round(prev$Seroprevalence * prev$N_indiv)/prev$N_indiv

# Host.Species is nested within Social.system, but you didn't reuse 
# species letters between Social.systems, so you can specify 
# Host.Species as a random effect without explicitly nesting it

# First random effect model
prev1.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species),
                  family=binomial(link="logit"), weights=N_indiv, data = prev)
summary(prev1.glmer)

## Fixed effects:
# Intercept is pathogen A and social.system A.  
# The z-test of the intercept is testing if the logit=0
# I.e. it's testing whether the combination of
# pathogen A and social.system A has prob=0.5.
# The other z-tests are testing whether other levels of the factors
# yield different probabilities than pathogen A and social.system A

## Random effects:
# This doesn't give you separate Host.Species and residual variances,
# Host.Species is treated as a random effect, so this model is the same as if
# you had summed the results of all studies with identical values of
# Host.Species, Pathogen, and Social.System. I.e. sum the results of the
# first 8 rows and create a single proportion and N_indiv, like so:

prevsum<-aggregate(cbind(N_indiv, prop=(seroexact*N_indiv)) ~ 
                   Social.System+Host.Species+Pathogen, data=prev, sum)
prevsum$prop<-prevsum$prop/prevsum$N_indiv

# which gives the same model:
prevsum.glmer = glmer(prop ~ Pathogen + Social.System + (1|Host.Species),
                      family=binomial(link="logit"), weights=N_indiv, data = prevsum)
summary(prevsum.glmer)

# So why are they broken up into multiple rows?  If each row represents
# one geographic area/time/litter/study/etc. then animals in one row
# might be more similar to eachother than they are to animals in
# another row that has the same values of Social, Species, & Pathogen.
# I think this is what the advisor wants as a "residual".

# To allow a random component for each row:
prev2<-cbind(resid=paste("Row_", row.names(prev), sep=""), prev)

prev2.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species) + (1|resid),
                   family=binomial(link="logit"), weights=N_indiv, data = prev2)
summary(prev2.glmer)

# This isn't a bad start, but I'm not comfortable with it because:
table(prev2[,2:3])

# Social.Sytstem D is only observed in Species F.
# This is called confounding, and it makes it hard to draw conclusions
# about Social Sytstem D.  How do you know what is caused by social
# system D and what is caused by species F?  If your friend really wants to
# make inferences about Social System D, she should collect data from
# another host species that uses Social System D.

# Leave out Soc_D:
prev3.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species) + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])
summary(prev3.glmer)

# Even though Host Species is conceptually a random factor, you really need to observe
# more than 2 species per social system for a mixed model to accurately estimate
# the species variance.  As far as species variance is concerned, each species is a
# single sample (not animals or even litters), and you can't hope to estimate variance
# accurately with only two samples.

# We can fit the model with species as a fixed effect, but we don't have
# enough degrees of freedom to estimate all levels of Species:
prev4.glmer = glmer(seroexact ~ Pathogen + Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])

# Your friend doesn't need to estimate the level of each species in order to test
# whether species has any noticeable effect at all.  Unfortunately, we can't just
# Use the F statistic from anova() because calculating the denominator df for a
# GLMM is not straightforward.
anova(prev4.glmer) #Gives you an F statistic, but no denominator df or p-value

# Instead we fit a simpler model without Species:
prev5.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])

# And we'll compare the two models With a Likelihood-Ratio test using anova()
anova(prev5.glmer,prev4.glmer)

# With a p-value of 0.01331 we can say it's worth keeping Species in the model.

# Now let's check the pathogen * social system interaction:
prev6.glmer = glmer(seroexact ~ Pathogen * Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, nAGQ=2, 
                    data = prev2[prev2$Social.System != "Soc_D",])
summary(prev6.glmer) #Neither interaction term is significant
anova(prev6.glmer)
# We don't need a denominator df to know that the F statistic of 0.0774 for
# the interaction is insignificant.

# Since the interaction between Pathogen and Social System was not significant,
# we don't need to include the interaction term.  Similarly, I don't see a 
# statistical reason to  split the model into two separate 'pathogen specific'
# models, but maybe there's a scientific reason to do so:

# Separate tests for each pathogen:
prev7A.glmer = glmer(seroexact ~ Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D" & prev$Pathogen == "Path_A",])
summary(prev7A.glmer)
# Social System B looks different from Social System A in pathogen A prevalance:

# Calculate the odds of having Pathogen A for Social System A vs B
beta7A<-fixef(prev7A.glmer)
exp(-beta7A[2]) #negative sign means odds of A:B instead of B:A
# So animals with Social System A have about 25 times the odds of
# animals with social system B of having Pathogen A

# Test for Pathogen B:
prev7B.glmer = glmer(seroexact ~ Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D" & prev$Pathogen == "Path_B",])
summary(prev7B.glmer)
# The only significant effects are species-specific, which are not of interest

# Let's return to prev4.glmer, which models both pathogens:
summary(prev4.glmer)

# The only significant fixed effect in prev4.glmer is Pathogen.
beta4<-fixef(prev4.glmer)

# For a randomly selected animal, the odds of having Pathogen B to having Pathogen A are:
exp(beta4[2])

# That's about as much as you can interpret with the data she has.

# To answer the Advisor's request for variance components:
# Residual variance is:
getME(prev4.glmer, "theta")^2

# You can't do a good job of estimating species variance with these data.
# If her advisor won't listen, then you can tell him that your estimate is:
getME(prev3.glmer, "theta")[2]^2
# But it's a really crappy estimate.

# There is no such thing as a variance component for Social System because
# it's a fixed effect.  But you can get its sum of squares:
anova(prev4.glmer)

I would recommend against: the advisor's "solution", using lm() for logistic regression, and weighting rows equally (you should weight by N_indiv).

# Seroprevalance has been rounded, that's not OK
# to do logistic regression, (proportion * weight) must equal an integer
prev$seroexact <- round(prev$Seroprevalence * prev$N_indiv)/prev$N_indiv

# Host.Species is nested within Social.system, but you didn't reuse 
# species letters between Social.systems, so you can specify 
# Host.Species as a random effect without explicitly nesting it

# First random effect model
prev1.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species),
                  family=binomial(link="logit"), weights=N_indiv, data = prev)
summary(prev1.glmer)

## Fixed effects:
# Intercept is pathogen A and social.system A.  
# The z-test of the intercept is testing if the logit=0
# I.e. it's testing whether the combination of
# pathogen A and social.system A has prob=0.5.
# The other z-tests are testing whether other levels of the factors
# yield different probabilities than pathogen A and social.system A

## Random effects:
# This doesn't give you separate Host.Species and residual variances,
# Host.Species is treated as a random effect, so this model is the same as if
# you had summed the results of all studies with identical values of
# Host.Species, Pathogen, and Social.System. I.e. sum the results of the
# first 8 rows and create a single proportion and N_indiv, like so:

prevsum<-aggregate(cbind(N_indiv, prop=(seroexact*N_indiv)) ~ 
                   Social.System+Host.Species+Pathogen, data=prev, sum)
prevsum$prop<-prevsum$prop/prevsum$N_indiv

# which gives the same model:
prevsum.glmer = glmer(prop ~ Pathogen + Social.System + (1|Host.Species),
                      family=binomial(link="logit"), weights=N_indiv, data = prevsum)
summary(prevsum.glmer)

# So why are they broken up into multiple rows?  If each row represents
# one geographic area/time/litter/study/etc. then animals in one row
# might be more similar to eachother than they are to animals in
# another row that has the same values of Social, Species, & Pathogen.
# I think this is what the advisor wants as a "residual".

# To allow a random component for each row:
prev2<-cbind(resid=paste("Row_", row.names(prev), sep=""), prev)

prev2.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species) + (1|resid),
                   family=binomial(link="logit"), weights=N_indiv, data = prev2)
summary(prev2.glmer)

# This isn't a bad start, but I'm not comfortable with it because:
table(prev2[,2:3])

# Social.Sytstem D is only observed in Species F.
# This is called confounding, and it makes it hard to draw conclusions
# about Social Sytstem D.  How do you know what is caused by social
# system D and what is caused by species F?  If your friend really wants to
# make inferences about Social System D, she should collect data from
# another host species that uses Social System D.

# Leave out Soc_D:
prev3.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species) + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])
summary(prev3.glmer)

# Even though Host Species is conceptually a random factor, you really need to observe
# more than 2 species per social system for a mixed model to accurately estimate
# the species variance.  As far as species variance is concerned, each species is a
# single sample (not animals or even litters), and you can't hope to estimate variance
# accurately with only two samples.

# We can fit the model with species as a fixed effect, but we don't have
# enough degrees of freedom to estimate all levels of Species:
prev4.glmer = glmer(seroexact ~ Pathogen + Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])

# Your friend doesn't need to estimate the level of each species in order to test
# whether species has any noticeable effect at all.  Unfortunately, we can't just
# Use the F statistic from anova() because calculating the denominator df for a
# GLMM is not straightforward.
anova(prev4.glmer) #Gives you an F statistic, but no denominator df or p-value

# Instead we fit a simpler model without Species:
prev5.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])

# And we'll compare the two models With a Likelihood-Ratio test using anova()
anova(prev5.glmer,prev4.glmer)

# With a p-value of 0.01331 we can say it's worth keeping Species in the model.

# Now let's check the pathogen * social system interaction:
prev6.glmer = glmer(seroexact ~ Pathogen * Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, nAGQ=2, 
                    data = prev2[prev2$Social.System != "Soc_D",])
summary(prev6.glmer) #Neither interaction term is significant
anova(prev6.glmer)
# We don't need a denominator df to know that the F statistic of 0.0774 for
# the interaction is insignificant.

# Since the interaction between Pathogen and Social System was not significant,
# we don't need to include the interaction term.  Similarly, I don't see a 
# statistical reason to  split the model into two separate 'pathogen specific'
# models, but maybe there's a scientific reason to do so:

# Separate tests for each pathogen:
prev7A.glmer = glmer(seroexact ~ Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D" & prev$Pathogen == "Path_A",])
summary(prev7A.glmer)
# Social System B looks different from Social System A in pathogen A prevalance:

# Calculate the odds of having Pathogen A for Social System A vs B
beta7A<-fixef(prev7A.glmer)
exp(-beta7A[2]) #negative sign means odds of A:B instead of B:A
# So animals with Social System A have about 25 times the odds of
# animals with social system B of having Pathogen A

# Test for Pathogen B:
prev7B.glmer = glmer(seroexact ~ Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D" & prev$Pathogen == "Path_B",])
summary(prev7B.glmer)
# The only significant effects are species-specific, which are not of interest

# Let's return to prev4.glmer, which models both pathogens:
summary(prev4.glmer)

# The only significant fixed effect in prev4.glmer is Pathogen.
beta4<-fixef(prev4.glmer)

# For a randomly selected animal, the odds of having Pathogen B to having Pathogen A are:
exp(beta4[2])

# That's about as much as you can interpret with the data she has.

# To answer the Advisor's request for variance components:
# Residual variance is:
getME(prev4.glmer, "theta")^2

# You can't do a good job of estimating species variance with these data.
# If her advisor won't listen, then you can tell him that your estimate is:
getME(prev3.glmer, "theta")[2]^2
# But it's a really crappy estimate.

# There is no such thing as a variance component for Social System because
# it's a fixed effect.  But you can get its sum of squares:
anova(prev4.glmer)

The advisor's "solution" is not good practice, nor is using lm() for logistic regression and weighting rows equally instead of by N_indiv.

# Seroprevalance has been rounded, that's not OK
# to do binomial regression, (proportion * weight) must equal an integer
prev$seroexact <- round(prev$Seroprevalence * prev$N_indiv)/prev$N_indiv

# Host.Species is nested within Social.system, but you didn't reuse 
# species letters between Social.systems, so you can specify 
# Host.Species as a random effect without explicitly nesting it

# First random effect model
prev1.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species),
                  family=binomial(link="logit"), weights=N_indiv, data = prev)
summary(prev1.glmer)

# This doesn't give you separate Host.Species and residual variances,
# Host.Species treated as a random effect, so this model is the same as if
# you had summed the results of all studies with identical values of
# Host.Species, Pathogen, and Social.System. I.e. sum the results of the
# first 8 rows and create a single proportion and N_indiv, like so:

prevsum<-aggregate(cbind(N_indiv, prop=(seroexact*N_indiv)) ~ 
                   Social.System+Host.Species+Pathogen, data=prev, sum)
prevsum$prop<-prevsum$prop/prevsum$N_indiv

# which gives the same model:
prevsum.glmer = glmer(prop ~ Pathogen + Social.System + (1|Host.Species),
                      family=binomial(link="logit"), weights=N_indiv, data = prevsum)
summary(prevsum.glmer)
 

# So why are they broken up into multiple rows?  If each row represents
# one geographic area/time/litter/study/etc. then animals in one row might
# be more similar to eachother than they are to animals in another row
# that has the same values of Social, Species, & Pathogen.  I think this
# is what the advisor wants as a "residual".

# To allow a random component for each row:
prev2<-cbind(resid=paste("Row_", row.names(prev), sep=""), prev)

prev2.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species) + (1|resid),
                   family=binomial(link="logit"), weights=N_indiv, data = prev2)
summary(prev2.glmer)

# This isn't a bad start, but I'm not comfortable with it because:
table(prev2[,2:3])

# Social.Sytstem D is only observed in Species F.
# This is called confounding, and it makes it hard to draw conclusions
# about Social Sytstem D.  How do you know what is caused by social
# system D and what is caused by species F?  If your friend really wants to
# make inferences about Social System D, she should collect data from
# another host species that uses Social System D.

# Leave out Soc_D:
prev3.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species) + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])
summary(prev3.glmer)
# Intercept is pathogen A and social.system A.  
# The z-test of the intercept is testing if the logit=0
# i.e. it's testing whether the combination of
# pathogen A and social.system A has prob=0.5
# The other z-tests are testing whether other levels of the factor
# yield different probabilities than pathogen A and social.system A

# The only significant effect is Pathogen, but first we need to check whether
# Pathogen interacts with Social System
prev4.glmer = glmer(seroexact ~ Pathogen * Social.System + (1|Host.Species) + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])
summary(prev4.glmer)

# Neither of the interactions are significant, so we can remove the interaction
# and use model prev3.glmer.  For the effect of Pathogen:
betas3<-fixef(prev3.glmer)

# For any given animal, the odds of having Pathogen B to having Pathogen A are:
exp(betas3[2])


# Separate tests for each pathogen:
prev3A.glmer = glmer(seroexact ~ Social.System + (1|Host.Species) + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D" & prev$Pathogen == "Path_A",])
summary(prev3A.glmer)

prev3B.glmer = glmer(seroexact ~ Social.System + (1|Host.Species) + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D" & prev$Pathogen == "Path_B",])
summary(prev3B.glmer)


# Even though Host Species is conceptually a random factor, you really need to observe
# more than 2 species per social system for a mixed model to accurately estimate
# the species variance.  As far as species variance is concerned, each species is a
# single unit (not animals or even litters), and you can't hope to estimate variance
# accurately with only two samples.

# We can fit the model with species as a fixed effect, but we don't have
# enough degrees of freedom to estimate all levels of Species.
prev5.glmer = glmer(seroexact ~ Pathogen + Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])

# Your friend doesn't need to estimate the level of each species in order to test
# whether species has any noticeable effect at all.  Unfortunately, we can't just
# Use the F statistic from anova() because calculating the denominator df for a
# GLMM is not straightforward.
anova(prev5.glmer) #Gives you an F statistic, but no denominator df or p-value

# Instead we fit a simpler model without Species:
prev6.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])

# And we'll compare the two models With a Likelihood-Ratio test using anova()
anova(prev6.glmer,prev5.glmer)

# With a p-value of 0.01331 we can say it's worth keeping Species in the model.

# Now let's go back and make sure that pathogen * social system interaction
# is not important now that we're treating species as fixed

prev7.glmer = glmer(seroexact ~ Pathogen * Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, nAGQ=2, 
                    data = prev2[prev2$Social.System != "Soc_D",])
summary(prev7.glmer)
anova(prev7.glmer)
# The interaction is insignificant

# Since the interaction between Pathogen and Social System was not significant,
# we don't need to include the interaction term.  Similarly, I don't see a 
# statistical reason to  split the model into two separate 'pathogen specific'
# models, but maybe there's a scientific reason to do so:

# Separate tests for each pathogen:
prev5A.glmer = glmer(seroexact ~ Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D" & prev$Pathogen == "Path_A",])
summary(prev5A.glmer)
# Social System B looks different from Social System A in pathogen A prevalance:

# Calculate the odds of having Pathogen A for Social System A vs B
betas5A<-fixef(prev5A.glmer)
exp(-betas5A[2]) #negative sign means odds of A:B instead of B:A
# So animals with Social System A have about 25 times the odds of
# animals with social system B of having Pathogen A

# Test for Pathogen B:
prev5B.glmer = glmer(seroexact ~ Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D" & prev$Pathogen == "Path_B",])
summary(prev5B.glmer)
# The only significant effects are species specific, which are not of interest
 

# Returning to prev5.glmer, which models both pathogens:
summary(prev5.glmer)

# The only significant fixed effect in prev5.glmer is Pathogen.
betas5<-fixef(prev5.glmer)

# For any given animal, the odds of having Pathogen B to having Pathogen A are:
exp(betas5[2])
# This is nearly the same as when we treated species as random


# That's about as much as you can interpret with the data she has.
 

# To answer the Advisor's request for variance components:
# Residual variance is:
getME(prev5.glmer, "theta")^2

# You can't do a good job of estimating species variance with these data.
# If her advisor won't listen, then you can tell him that your estimate is:
getME(prev3.glmer, "theta")[2]^2
# But it's a really crappy estimate.

# There is no such thing as a variance component for Social System because
# it's a fixed effect.  You can get its sum of squares:
anova(prev5.glmer)

I would recommend against: the advisor's "solution", using lm() for logistic regression, and weighting rows equally (weight by N_indiv).

# Seroprevalance has been rounded, that's not OK
# to do binomial regression, (proportion * weight) must equal an integer
prev$seroexact <- round(prev$Seroprevalence * prev$N_indiv)/prev$N_indiv

# Host.Species is nested within Social.system, but you didn't reuse 
# species letters between Social.systems, so you can specify 
# Host.Species as a random effect without explicitly nesting it

# First random effect model
prev1.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species),
                  family=binomial(link="logit"), weights=N_indiv, data = prev)
summary(prev1.glmer)

## Fixed effects:
# Intercept is pathogen A and social.system A.  
# The z-test of the intercept is testing if the logit=0
# I.e. it's testing whether the combination of
# pathogen A and social.system A has prob=0.5.
# The other z-tests are testing whether other levels of the factors
# yield different probabilities than pathogen A and social.system A

## Random effects:
# This doesn't give you separate Host.Species and residual variances,
# Host.Species is treated as a random effect, so this model is the same as if
# you had summed the results of all studies with identical values of
# Host.Species, Pathogen, and Social.System. I.e. sum the results of the
# first 8 rows and create a single proportion and N_indiv, like so:

prevsum<-aggregate(cbind(N_indiv, prop=(seroexact*N_indiv)) ~ 
                   Social.System+Host.Species+Pathogen, data=prev, sum)
prevsum$prop<-prevsum$prop/prevsum$N_indiv

# which gives the same model:
prevsum.glmer = glmer(prop ~ Pathogen + Social.System + (1|Host.Species),
                      family=binomial(link="logit"), weights=N_indiv, data = prevsum)
summary(prevsum.glmer)

# So why are they broken up into multiple rows?  If each row represents
# one geographic area/time/litter/study/etc. then animals in one row
# might be more similar to eachother than they are to animals in
# another row that has the same values of Social, Species, & Pathogen.
# I think this is what the advisor wants as a "residual".

# To allow a random component for each row:
prev2<-cbind(resid=paste("Row_", row.names(prev), sep=""), prev)

prev2.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species) + (1|resid),
                   family=binomial(link="logit"), weights=N_indiv, data = prev2)
summary(prev2.glmer)

# This isn't a bad start, but I'm not comfortable with it because:
table(prev2[,2:3])

# Social.Sytstem D is only observed in Species F.
# This is called confounding, and it makes it hard to draw conclusions
# about Social Sytstem D.  How do you know what is caused by social
# system D and what is caused by species F?  If your friend really wants to
# make inferences about Social System D, she should collect data from
# another host species that uses Social System D.

# Leave out Soc_D:
prev3.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|Host.Species) + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])
summary(prev3.glmer)

# Even though Host Species is conceptually a random factor, you really need to observe
# more than 2 species per social system for a mixed model to accurately estimate
# the species variance.  As far as species variance is concerned, each species is a
# single sample (not animals or even litters), and you can't hope to estimate variance
# accurately with only two samples.

# We can fit the model with species as a fixed effect, but we don't have
# enough degrees of freedom to estimate all levels of Species:
prev4.glmer = glmer(seroexact ~ Pathogen + Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])

# Your friend doesn't need to estimate the level of each species in order to test
# whether species has any noticeable effect at all.  Unfortunately, we can't just
# Use the F statistic from anova() because calculating the denominator df for a
# GLMM is not straightforward.
anova(prev4.glmer) #Gives you an F statistic, but no denominator df or p-value

# Instead we fit a simpler model without Species:
prev5.glmer = glmer(seroexact ~ Pathogen + Social.System + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D",])

# And we'll compare the two models With a Likelihood-Ratio test using anova()
anova(prev5.glmer,prev4.glmer)

# With a p-value of 0.01331 we can say it's worth keeping Species in the model.

# Now let's check the pathogen * social system interaction:
prev6.glmer = glmer(seroexact ~ Pathogen * Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, nAGQ=2, 
                    data = prev2[prev2$Social.System != "Soc_D",])
summary(prev6.glmer) #Neither interaction term is significant
anova(prev6.glmer)
# We don't need a denominator df to know that the F statistic of 0.0774 for
# the interaction is insignificant.

# Since the interaction between Pathogen and Social System was not significant,
# we don't need to include the interaction term.  Similarly, I don't see a 
# statistical reason to  split the model into two separate 'pathogen specific'
# models, but maybe there's a scientific reason to do so:

# Separate tests for each pathogen:
prev7A.glmer = glmer(seroexact ~ Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D" & prev$Pathogen == "Path_A",])
summary(prev7A.glmer)
# Social System B looks different from Social System A in pathogen A prevalance:

# Calculate the odds of having Pathogen A for Social System A vs B
beta7A<-fixef(prev7A.glmer)
exp(-beta7A[2]) #negative sign means odds of A:B instead of B:A
# So animals with Social System A have about 25 times the odds of
# animals with social system B of having Pathogen A

# Test for Pathogen B:
prev7B.glmer = glmer(seroexact ~ Social.System + Host.Species + (1|resid),
                    family=binomial(link="logit"), weights=N_indiv, 
                    data = prev2[prev2$Social.System != "Soc_D" & prev$Pathogen == "Path_B",])
summary(prev7B.glmer)
# The only significant effects are species-specific, which are not of interest

# Let's return to prev4.glmer, which models both pathogens:
summary(prev4.glmer)

# The only significant fixed effect in prev4.glmer is Pathogen.
beta4<-fixef(prev4.glmer)

# For a randomly selected animal, the odds of having Pathogen B to having Pathogen A are:
exp(beta4[2])

# That's about as much as you can interpret with the data she has.

# To answer the Advisor's request for variance components:
# Residual variance is:
getME(prev4.glmer, "theta")^2

# You can't do a good job of estimating species variance with these data.
# If her advisor won't listen, then you can tell him that your estimate is:
getME(prev3.glmer, "theta")[2]^2
# But it's a really crappy estimate.

# There is no such thing as a variance component for Social System because
# it's a fixed effect.  But you can get its sum of squares:
anova(prev4.glmer)