Main question: "What are the contributing variates to daily movement distances?"
Specifically my question today relates to:
"What is the contribution related to gender, and then within female what are the effects related to those females having young?"
I have multiple readings per animal (hundreds per individual), with 13 individuals (5 females and 8 males). The females sometimes have young with them, and I know this contributes to the distance they move.
I have several contributing factors in a GLMM; I am using the `nlme::lme` function. The current form is:
lm1 <- lme(movedistance ~ Gender+YoungPresent+x3+x4+x5+x6, random = ~1|AnimalID/Month, data = df1)
There is no significant gender effect in this current model; however I know that this is mis-specified, because males never have young; but I don't know how to fix it. "YoungPresent" is a binary term, it is always 0 for males, and 0 or 1 for females, 1 when they have young. What I want is to somehow remove the attribution of variation by "YoungPresent" to males in the "Gender" term, but not from females.
Please let me know what is the correct term for what I am looking for (Crossed? Nested?), and how I can correctly specify this structure in lme.
(EDITED)
After suggestions from the first response, the code now looks like this
> 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
However, when I compare this to the original model formulation, I suspect that the dummy variables are not properly addressed in the `lme()`. I coded "MvsF" as -0.66667 for males and 0.3333 for all females, and yet the estimate, s.e. and probaliity is the same as using the original "gender" variate.
> 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
I suspect that I need to call the MvsF using some kind of signal to `lme()` to let it know the values of MvsF are important, similarly to the way I might use factor(variate) or I(variate) inline to denote the way `lme()` should handle each variate. factor(MvsF clearly has no effect (basically what I have shown above), `lme()` does not treat the variables MvsF and Gender any differently.
If what my @JakeWestfall suggests has been used correctly, then the only thing new that I have added to the model is the `YvsNY` where 'Males' are coded differently to 'Females with No Young', where in the original variate they were coded the same, 0. This has changed the model for sure, and looks more like its on the right path, but why did I code `MvsF` to THOSE values, if it changes nothing? I could easily have ONLY changed `YoungPresent` (0/1) to YvsNY (0, -.5, .5)....
One of the problems as I see it, is that males are *still* included in the YvsN variate - the parameter YvsN estimates a line that goes through three points on the x axis: the three levels of that variate - (-.5,0,.5 = Young, Male, No Young), and therefore Males are still contributing to the estimate of this variate - when I think they should not. I believe what I may need is similar to a grouping structure (in the random term?) where `YvsN` is nested within Gender (or `MvsF`, I think it doesn't matter) such that Males do not contribute to estimation of the `YvsN` parameter.
--------NEW Values of NvsNY and 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
The effect seems small, but it still seems possible to push around the `YvsNY` estimate, by changing males response values. This is what worries me.