Hopefully your friend has graduated by now, but if not, the following might help.
You were on the right track in your original post http://stats.stackexchange.com/questions/93450/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 (you should weight by N_indiv).
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 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)