ANOVA anidado: Tamaños de muestra desiguales? ¿Componentes de varianza?

Question

Estoy completamente fuera de mi profundidad sobre esto, y de toda la lectura que intente hacer me confunde. Estoy esperando que usted puede explicarle las cosas a mí de una manera que tenga sentido. (Como siempre parece ser el caso, "no debería ser tan difícil!")

Estoy tratando de ayudar a un alumno que está estudiando el efecto de los sistemas sociales en la prevalencia de las enfermedades en diversos cánido especies de acogida. Queremos considerar el sistema social (por ejemplo, del grupo de-vida vs solitario) como un efecto fijo, y albergan especies como efecto aleatorio anidado dentro de un sistema social (es decir, cada una de las especies sólo se dispone de un sistema social de tipo).

Mi entendimiento es que la mejor manera de hacer esto sería hacer una mezcla de efectos de regresión logística. Hemos hecho esto, y funciona, y éramos felices. Por desgracia, su asesor es insistir en que ella calcular la cantidad de variación debido a que el sistema social del huésped contra las especies vs residual. Yo no puedo averiguar cómo hacerlo a través de efectos mixtos de regresión logística, y mi pregunta anterior sobre este tema quedó sin respuesta.

Su asesor sugirió hacer ANOVA lugar, logit-la transformación de la prevalencia de la enfermedad de valores (la fracción de la población que está infectado). Esto presenta un problema debido a que algunos de la prevalencia de los valores son 0 o 1, lo que resultaría en $-\infty$ o $\infty$ una vez logit-transformado. Su asesor de la "solución" fue sustituto $-5$ $5$ $-\infty$ o $\infty$, respectivamente. Esto se siente realmente kludgey y me hace temblar bastante duro. Pero él es el uno de calificaciones de ella, y en este momento solo quiero hacer con esto, así que si él está bien con él, luego lo que sea.

Estamos usando R para este análisis. El código se puede descargar aquí, y los datos de entrada aquí. El archivo de datos incluye datos sobre dos diferentes patógenos (a y B), de la cual estamos analizando por separado (como se muestra en el código).

Aquí está el ANOVA de configuración que hemos hecho para el Patógeno B:

mod1.lm <- lm(Seroprevalence_logit ~ Social.System + Social.System/Host.Species,
              data = prev_B)
print(mod1.anova <- anova(mod1.lm))

Esto lleva a mi primera pregunta: ¿Es correcto y apropiado? Factores a considerar:

• Queremos tener un Modelo II (efecto aleatorio) variable anidada dentro de un Modelo I (efecto fijo) de la variable.
• No todo sistema social tiene el mismo número de especies de acogida anidada dentro de ella.
• No todas las especies de acogida tiene el mismo número de poblaciones examinadas.
• No toda la población examinados tenían el mismo número de individuos (columna N_indiv en mydata.csv). Esto es más de una ponderación problema de algo más fundamental, creo yo.

Mi siguiente pregunta, y la principal de este post, es: ¿Cómo puedo partición de la varianza? He aquí lo que estaban pensando:

MS_A <- mod1.anova$"Mean Sq"[1]
MS_BinA <- mod1.anova$"Mean Sq"[2]
MS_resid <- mod1.anova$"Mean Sq"[3]
n <- length(unique(prev_A$Social.System))
r <- length(unique(prev_A$Host.Species))
VC_A <- (MS_A - MS_BinA)/(n*r)
VC_BinA <- (MS_BinA - MS_resid)/n
VC_resid <- MS_resid

Por desgracia, esto se traduce en tristeza utilizando el análisis de VARIANZA de la especificación he detallado anteriormente. Aquí están los resultados para el Patógeno B:

• VC_A (es decir, Social.Sistema): $-1.48$
• VC_BinA (es decir, en el de Acogida.Las especies): $13.8$
• VC_resid: $5.57$

La investigación me lleva a creer que esto debe resultar en componentes de varianza de los porcentajes de 0%, el 71,3% y 28,7%, respectivamente. Sin embargo, esto es insuficiente por dos razones:

• El p-valor Social.Sistema de la ANOVA se ~$0.025$, lo que sugiere que se deben tener en cuenta al menos algunos de la varianza observada. (Anfitrión.Especies tenían un valor de p ~$3*10^{-5}$.)
• Me preocupa que una desviación negativa componente podría ser una bandera roja para algo.

Por favor, cualquier ayuda que puede representar a cualquiera de estas preguntas, sería muy apreciado. I TA d a un curso de licenciatura en bioestadística, así que tengo un poco de contexto, pero me parece que no puede averiguar estos problemas específicos. Gracias de antemano.

Answer 1

Es de esperar que su amigo se ha graduado por ahora, pero si no, el siguiente podría ser de ayuda.

Que estábamos en el camino correcto en tu post original de Partición de la varianza de la regresión logística, utilizando glmer() de efectos mixtos de regresión logística.

Yo recomendaría contra: el asesor de la "solución", el uso de lm() para la regresión logística, y la ponderación de filas igual (debe de peso por N_indiv).

Los modelos mixtos lineales generalizados son difíciles. http://glmm.wikidot.com/faq tiene buena información - incluyendo el hecho de que usted necesita muchos niveles de un factor aleatorio con el fin de estimar las componentes de varianza.

Mi siguiente código requiere que el paquete lme4 y los datos de su enlace.

# 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)