4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ You did change the definition of YvsN as well, right? I can't tell because you still don't show where you create those variables! If YvsN has not been redefined and is instead still using the values you mentioned in your original question, then we would expect for both of these models to be identical. Also, it would be more useful to display the codes by calling table() on them rather than (or in addition to) str(). $\endgroup$ – Jake Westfall Nov 11 '13 at 17:29
  • $\begingroup$ Added values from table() as specified by comment from @JakeWestfall $\endgroup$ – KalahariKev Nov 12 '13 at 1:57
  • $\begingroup$ @JakeWestfall I think this is what I am looking for... Comments? stat.ethz.ch/pipermail/r-sig-mixed-models/2011q4/017352.html $\endgroup$ – KalahariKev Nov 12 '13 at 3:56
  • $\begingroup$ The model you have appended to your question and referred to as your "original model specification" is clearly not your original model specification. The original model, as you wrote in your question, was movedistance ~ Gender+YoungPresent+.... In the other words, it did not contain the YvsN factor that I told you about, but rather contained a YoungPresent factor which "is a binary term, it is always 0 for males, and 0 or 1 for females, 1 when they have young," quite unlike the YvsN factor. The model I suggested to you is clearly not equivalent to the model with YoungPresent! $\endgroup$ – Jake Westfall Nov 12 '13 at 4:19
  • $\begingroup$ The YvsNY term has made a considerable difference to the estimation of the effect of presence of young. This I agree with. Therefore the original model and the new model ARE different. My issue remains that males are still considered within the new term. Previously they were modelled as 0 along with females without young, also 0. Only females with cubs differed. Now I have three levels in this factor, YvsNY. What that does is estimates a single line that still includes males in consideration of the effect of YvsNY just seperately to Females with no young. This is different, but not what I want $\endgroup$ – KalahariKev Nov 12 '13 at 4:38
3
$\begingroup$

You have three groups: males (M), females with young (FY), females without young (FN). Although I understand why it would be conceptually tempting to think of this as two factors (gender and parental status), I think from a statistical point of view it is better to think of these groups as comprising a single factor, "Group", with 3 levels. Under this point of view, it turns out that the two questions you want to ask map perfectly onto a standard set of orthogonal contrast codes you can use to represent this factor. The codes would look like this:

(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

In this new set of codes, the "MvsF" code represents the mean difference between males and females (more specifically, the difference between the Male mean and the mean of the two Female means), and the "YvsN" code represents the mean difference between females with and without young.

$\endgroup$
  • $\begingroup$ Hi Jake. Thanks for the answer. It has some intuitive sense, but I don't know enough to give it the thumbs up yet . I will try it while I wait to see if anyone upvotes your answer, or if another comes up. Cheers. $\endgroup$ – KalahariKev Nov 10 '13 at 3:13
  • $\begingroup$ @KalahariKev No problem. Let me know if you have any questions. $\endgroup$ – Jake Westfall Nov 10 '13 at 3:23
  • $\begingroup$ Sounds like the OP should replace Gender and young present with female+young and female+no young. I disagree about your dummy coding comment though. $\endgroup$ – probabilityislogic Nov 10 '13 at 7:26
  • $\begingroup$ @probabilityislogic Would you care to expand on the dummy coding? Would you code it differently, or do you disagree with the interpretation of the final coefficients? Does your suggestion imply a single variable with three levels would be sufficient? I started to think that there is a way to recode Jake Westfall s variable MvsF to represent all three possibilities, but I personally would be unsure of the interpretation. $\endgroup$ – KalahariKev Nov 11 '13 at 3:23
  • $\begingroup$ @JakeWestfall YvsNY 12.59318 4.54336 2.772 0.005575 ** $\endgroup$ – KalahariKev Nov 11 '13 at 3:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.