La especificación del modelo GLMM ayuda a los efectos de género + un efecto que se anida sólo dentro de la mujer

Question

Pregunta principal: "¿Cuáles son las variantes que contribuyen a las distancias de desplazamiento diarias?"

Concretamente mi pregunta de hoy se refiere a:

"¿Cuál es la contribución relacionada con el género, y luego dentro de las mujeres cuáles son los efectos relacionados con que esas hembras tengan crías?"

Tengo varias lecturas por animal (cientos por individuo), con 13 individuos (5 hembras y 8 machos). Las hembras a veces tienen crías con ellas, y sé que esto contribuye a la distancia a la que se mueven.

Tengo varios factores contribuyentes en un GLMM; estoy utilizando el nlme::lme función. La forma actual es:

lm1 <- lme(movedistance ~ Gender+YoungPresent+x3+x4+x5+x6, random  = ~1|AnimalID/Month, data = df1)

No hay un efecto significativo de género en este modelo actual; sin embargo, sé que está mal especificado, porque los hombres nunca tienen jóvenes; pero no sé cómo arreglarlo. "YoungPresent" es un término binario, siempre es 0 para los machos, y 0 o 1 para las hembras, 1 cuando tienen crías. Lo que quiero es eliminar de alguna manera la atribución de la variación por "YoungPresent" a los machos en el término "Gender", pero no de las hembras.

Por favor, indíqueme cuál es el término correcto para lo que estoy buscando (¿Cruzado? ¿Anidado?), y cómo puedo especificar correctamente esta estructura en lme.

(EDITADO)

Después de las sugerencias de la primera respuesta, el código tiene ahora el siguiente aspecto

> str(df1$YvsNY
num [1:6308] 0 0 0 0 0 0 0 0 0 0 ...

> str(df1$MvsF)
num [1:6308] -0.667 -0.667 -0.667 -0.667 -0.667 ...
> dmd <- lme(dist~Age+MvsF+TempMax+MeanRain+herb1_dens+herb2_dens+YvsNY+herb3_dens, random  = ~1|ANIMALID/Month, data = df1

> summary(lm1)
Linear mixed-effects model fit by REML
Data: df1 
      AIC      BIC    logLik
 57915.15 57995.23 -28945.58

Random effects:
Formula: ~1 | ANIMALID
   (Intercept)
StdDev:    7.923558
Formula: ~1 | Month %in% ANIMALID
   (Intercept) Residual
StdDev:    7.150394  33.4111

Fixed effects: dist ~ Age + MvsF + TempMax + MeanRain + herb1_dens + herb2_dens +          YvsNY + herb3_dens 
            Value Std.Error   DF   t-value p-value
(Intercept)  86.08050 10.468338 5639  8.222939  0.0000
Age           1.47128  0.967371   10  1.520906  0.1593
MvsF        -10.80126  5.214172   10 -2.071520  0.0651
TempMax      -0.58513  0.136191 5639 -4.296398  0.0000
MeanRain     -0.08233  0.020589  197 -3.998523  0.0001
herb1_dens    0.53651  0.327763  197  1.636886  0.1033
herb2_dens   -0.04928  0.032569  197 -1.513059  0.1319
YvsNY        13.07835  4.435959  197  2.948257  0.0036
herb3_dens    3.51159  1.797992  197  1.953061  0.0522

Sin embargo, cuando comparo esto con la formulación original del modelo, sospecho que las variables ficticias no se abordan adecuadamente en el lme() . Codifiqué "MvsF" como -0,66667 para los hombres y 0,3333 para todas las mujeres y, sin embargo, la estimación, la e.s. y la probalidad son las mismas que utilizando la variante original "género".

> str(df1$Gender)
Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 2 ...

> dmd2 <- lme(dist~Age+Gender+TempMax+MeanRain+herb1_dens+herb2_dens+YvsNY+herb3_dens, random  = ~1|ANIMALID/Month, data = df1

Linear mixed-effects model fit by REML
 Data: df1
   AIC      BIC    logLik
  57915.15 57995.23 -28945.58

Random effects:
 Formula: ~1 | ANIMALID
    (Intercept)
StdDev:    7.924217

 Formula: ~1 | Month %in% ANIMALID
    (Intercept) Residual
StdDev:    7.150576 33.41108

Fixed effects: dist ~ Age + Gender + TempMax + MeanRain + herb1_dens +      herb2_dens + YvsNY + herb3_dens 
           Value Std.Error   DF   t-value p-value
(Intercept) 82.47974 10.290996 5639  8.014748  0.0000
Age          1.47131  0.967439   10  1.520825  0.1593
GenderMale  10.80111  5.214561   10  2.071337  0.0651
TempMax     -0.58512  0.136191 5639 -4.296341  0.0000
MeanRain    -0.08233  0.020590  197 -3.998450  0.0001
herb1_den    0.53652  0.327768  197  1.636892  0.1032
herb2_des   -0.04928  0.032570  197 -1.513020  0.1319
YvsNY       13.07854  4.436096  197  2.948209  0.0036
herb3_dens   3.51152  1.798025  197  1.952988  0.0522

Sospecho que tengo que llamar al MvsF usando algún tipo de señal para lme() para hacerle saber que los valores de MvsF son importantes, de forma similar a como podría usar factor(variate) o I(variate) en línea para denotar la forma lme() debe manejar cada variante. factor(MvsF claramente no tiene efecto (básicamente lo que he mostrado arriba), lme() no trata las variables MvsF y Gender de forma diferente.

Si lo que sugiere mi @JakeWestfall se ha utilizado correctamente, entonces lo único nuevo que he añadido al modelo es el YvsNY donde los "Hombres" se codifican de forma diferente a las "Mujeres sin hijos", mientras que en la variante original se codificaban igual, 0. Esto ha cambiado el modelo sin duda, y parece que está en el camino correcto, pero ¿por qué codifiqué MvsF a ESOS valores, si no cambia nada? Podría fácilmente haber cambiado SOLO YoungPresent (0/1) a YvsNY (0, -.5, .5)....

Uno de los problemas, tal y como yo lo veo, es que los hombres son todavía incluido en la variante YvsN - el parámetro YvsN estima una línea que pasa por tres puntos en el eje x: los tres niveles de esa variante - (-.5,0,.5 = Joven, Hombre, No Joven), y por tanto los Hombres siguen contribuyendo a la estimación de esta variante - cuando creo que no deberían hacerlo. Creo que lo que puedo necesitar es algo similar a una estructura de agrupación (¿en el término aleatorio?) donde YvsN está anidado dentro de Gender (o MvsF Creo que no importa) de manera que los hombres no contribuyen a la estimación de la YvsN parámetro.

--------NEW Valores de NvsNY y MvsF------------

> table(df1$YvsNY)

-0.5    0  0.5 
1180 3172 1581 
> table(df1$MvsF)

-0.666666667  0.333333333 
    3172         2761 

 #Finally to check if this worked, I added a value of 2 to all the male response variates:
  > df2 <- df1
  > df2 <- df2[df2$Gender=="Male",]$dist <- df2[df2$Gender=="Male",]$dist +2

  # and checked that only the males' data was affected:

  > tapply(df1$dist, df1$Gender, mean)
   Female     Male 
  81.01595 92.07785 
  > tapply(df2$dist, df$Gender, mean)
  Female     Male 
  81.01595 94.07785 

  > dmd1 <- nlme(dist~Age+MvsF+TempMax+MeanRain+herb1_dens+herb2_dens+YvsNY+herb3_dens, random  = ~1|AnimalID/Month, data = df1)
  > dmd2 <- nlme(dist~Age+MvsF+TempMax+MeanRain+herb1_dens+herb2_dens+YvsNY+herb3_dens, random  = ~1|AnimalID/Month, data = df2)
  > summary(dmd1)

  #(truncated)

                 Value Std.Error   DF   t-value p-value
  (Intercept)  86.08050 10.468338 5639  8.222939  0.0000
  Age           1.47128  0.967371   10  1.520906  0.1593
  MvsF        -10.80126  5.214172   10 -2.071520  0.0651
  TempMax      -0.58513  0.136191 5639 -4.296398  0.0000
  MeanRain     -0.08233  0.020589  197 -3.998523  0.0001
  herb1_dens    0.53651  0.327763  197  1.636886  0.1033
  herb2_dens   -0.04928  0.032569  197 -1.513059  0.1319
  YvsNY        13.07835  4.435959  197  2.948257  0.0036
  herb3_dens    3.51159  1.797992  197  1.953061  0.0522

  > summary (dmd2)      

  #truncated
                  Value Std.Error   DF   t-value p-value
  (Intercept)  86.74714 10.468406 5639  8.286567  0.0000
  Age           1.47128  0.967379   10  1.520896  0.1593
  MvsF        -12.80125  5.214219   10 -2.455065  0.0340
  TempMax      -0.58513  0.136191 5639 -4.296397  0.0000
  MeanRain     -0.08233  0.020589  197 -3.998520  0.0001
  herb1_dens    0.53651  0.327763  197  1.636889  0.1033
  herb2_dens   -0.04928  0.032569  197 -1.513057  0.1319
  YvsNY        13.07837  4.435970  197  2.948254  0.0036
  herb3_dens    3.51158  1.797993  197  1.953054  0.0522

  #VERY close, but a miniscule difference in coefficient had me a little worried, so I multiplied the response my 2:

  > df3 <- df1
  > df3 <- df1[df1$Gender=="Male",]$dist <- df1[df1$Gender=="Male",]$dist *2

  # and checked that only the males' data was affected:
  > tapply (df3$dist, df3$Gender, mean)
     Female      Male 
    81.01595 184.15570 

  > dmd3 <- nlme(dist~Age+MvsF+TempMax+MeanRain+herb1_dens+herb2_dens+YvsNY+herb3_dens, random  = ~1|AnimalID/Month, data = df3)
  > summary(dmd3)

  #(truncated)
             Value Std.Error   DF    t-value p-value
     (Intercept)  121.22306 17.048079 5639   7.110658  0.0000
     Age            2.75032  1.550867   10   1.773407  0.1066
     MvsF        -101.41686  8.296464   10 -12.224107  0.0000
     TempMax       -1.14168  0.232908 5639  -4.901840  0.0000
     MeanRain      -0.13735  0.035870  197  -3.829013  0.0002
     herb1_dens     0.62596  0.570363  197   1.097478  0.2738
     herb2_dens    -0.12191  0.056353  197  -2.163386  0.0317
     YvsNY         14.71697  7.506579  197   1.960543  0.0513
     herb3_dens     6.15790  3.110133  197   1.979948  0.0491

El efecto parece pequeño, pero aún así parece posible empujar alrededor de la YvsNY estimación, cambiando los valores de respuesta de los machos. Esto es lo que me preocupa.

Answer 1

Tiene tres grupos: machos (M), hembras con crías (FY), hembras sin crías (FN). Aunque entiendo que sería conceptualmente tentador pensar en esto como dos factores (género y situación parental), creo que desde un punto de vista estadístico es mejor pensar que estos grupos comprenden un único factor, "Grupo", con 3 niveles. Bajo este punto de vista, resulta que las dos preguntas que se quieren hacer encajan perfectamente en un conjunto estándar de códigos de contraste ortogonales que se pueden utilizar para representar este factor. Los códigos serían los siguientes

(group <- gl(n=3, k=5, labels=c("M","FY","FN")))
#  [1] M  M  M  M  M  FY FY FY FY FY FN FN FN FN FN
# Levels: M FY FN

# the default contrasts (i.e., dummy codes) suck:
contrasts(group)
#    FY FN
# M   0  0
# FY  1  0
# FN  0  1

# replace with meaningful contrasts:
contrasts(group) <- cbind(MvsF=c(-2/3, 1/3, 1/3), YvsN=c(0, -1/2, 1/2))
contrasts(group)
#          MvsF YvsN
# M  -0.6666667  0.0
# FY  0.3333333 -0.5
# FN  0.3333333  0.5

En este nuevo conjunto de códigos, el código "MvsF" representa la diferencia media entre machos y hembras (más concretamente, la diferencia entre la media de los machos y la media de las dos hembras), y el código "YvsN" representa la diferencia media entre las hembras con y sin cría.