Return to Question

I have several contributing factors in a GLMM; I am using the nlme::lme function. The current form is:

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.

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.

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.

I have several contributing factors in a GLMM; I am using the nlme::lme function. The current form is:

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.

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.

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.

> table(df1$YvsNY)

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

-0.666666667  0.333333333 
    3172         2761 

> 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.

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