Hopefully your friend has graduated by now, but if not, the following might help.
The advisor's "solution" is not good practice, nor is using lm() for logistic regression and weighting rows equally instead of by N_indiv.
You were on the right track in your original post Partitioning variance from logistic regression, using glmer() for mixed-effects logistic regression.
Generalized linear mixed models are tough. http://glmm.wikidot.com/faq has some good information - including the fact that you need many levels of a random factor in order to estimate its variance component.
My code below requires the lme4 package and the data from your link.
# 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)