[R-sig-ME] Coefficients interpretation and plot

[Jarrod Hadfield]

Hi Luciano,

If I understand you correctly, your issue is that the prediction for  
the first egg in the year that is NOT (?) 2007 is greater than second  
eggs in that year, yet the raw data indicate the opposite?

I notice that you have less than 3 eggs for each nest. If there is a  
(positive) relationship between clutch size and egg volume you could  
get such a discrepancy. You could try putting clutch size in the  
model. That being said, the offending coefficient is small with a  
large p-value (0.35) so the discrepancy may not be that surprising.

Also, I'm not sure what the state of play with pvals.func is. mcmcsamp  
used to behave oddly, and from your output the fixed effects look OK,  
but the 95% MCMC CI's for the variance components do not seem to  
overlap the REML estimates. Its possible there on a different scale,  
but I would check.

Cheers,

Jarrod




Quoting Luciano La Sala <lucianolasala at yahoo.com.ar>:

> Hello everyone,
>
> Since I'm not entirely sure this is THE list I should be referring too, feel
> free to blow me off and refer me to another mailing list if necessary.
> I am analyzing a small dataset using lmer from lme4 package. My model has
> "egg volume" as dependent variable and "hatching order" and "year" as
> dependent variables. The best fit model has these two variables plus their
> interaction (hatching order*year). I included Nest_ID as random intercepts.
> Output follows:
>
>> best <- lmer(EggVolume~HatchOrder+Year+HatchOrder*Year+(1|NestID), data =
> Data)
>> summary(best)
>
> Linear mixed model fit by REML
>
> Formula: EggVolume ~ HatchOrder + Year + HatchOrder * Year + (1 | NestID)
>    Data: Data
>    AIC BIC logLik deviance REMLdev
>  736.1 759 -360.1    729.1   720.1
>
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  NestID   (Intercept) 26.2931  5.1277
>  Residual              6.2175  2.4935
>
> Number of obs: 130, groups: NestID, 55
>
> Fixed effects:
>                       Estimate Std. Error t value
> (Intercept)           79.7261     1.1350   70.24
> HatchSecond           -0.7227     0.7758   -0.93
> HatchThird            -4.8455     0.9112   -5.32
> Year2007               3.5548     1.5750    2.26
> HatchSecond:Year2007  -2.6914     1.0752   -2.50
> HatchThird:Year2007   -2.7999     1.2294   -2.28
>
> Correlation of Fixed Effects:
>             (Intr) HtchOS HtchOT Yr2007 HOS:Y2
> HtchOrdrScn -0.277
> HtchOrdrThr -0.229  0.388
> Year2007    -0.721  0.199  0.165
> HtcOS:Y2007  0.200 -0.722 -0.280 -0.299
> HtcOT:Y2007  0.170 -0.287 -0.741 -0.301  0.415
>
> I used the "pvals.fnc" function in the "coda" package to estimate p-values.
> Output follows:
>
>> pvals.fnc(best, nsim = 10000, ndigits = 4, withMCMC = FALSE,
> addPlot=FALSE)
>
> $fixed
>                   Estimate MCMCmean HPD95lower HPD95upper  pMCMC Pr(>|t|)
> (Intercept)       79.7261  79.5990     77.687    81.5521 0.0001   0.0000
> HatchSecond       -0.7227   0.1239     -2.630     2.8256 0.9468   0.3534
> HatchThird        -4.8455  -4.3177     -7.391    -1.0711 0.0086   0.0000
> Year2007           3.5548   3.9605      1.090     6.8664 0.0080   0.0258
> HatchSecond:2007  -2.6914  -3.4393     -7.235     0.3782 0.0772   0.0136
> HatchThird:2007   -2.7999  -3.6649     -7.768     0.5046 0.0830   0.0245
>
> $random
>     Groups        Name Std.Dev. MCMCmedian MCMCmean HPD95lower HPD95upper
> 1   NestID (Intercept)   5.1277     2.3265   2.3166     1.5388     3.1619
> 2 Residual               2.4935     4.5507   4.5744     3.8179     5.4108
>
> I understand that, even in mixed models, one should not interpret main
> effects' coefficients by themselves when a significant interaction is
> present. Instead, coefficients of the main effect and the interaction term
> should be added. In my example:
>
> # Coefficient for HatchSecond:Year2007: -0.7227 + (-2.6914) = -3.4141
>
> Then, in 2007 the volume of HatchSecond eggs was 3.41 units lower than that
> of HathFirst eggs.
>
> # HatchThird:Year2007: -4.8455 + (-2.7999) = -7.6454
>
> Then, in 2007 the volumen of HatchThird eggs was 7.64 units lower than that
> of HathFirst eggs.
>
> 1. Are these interpretations correct in the contest of mixed modeling?
>
> 2. When I plot my raw data (means of egg volume for each year and stratified
> by hatching order), the plot looks good for 2007 (decreasing egg volumes
> along the hatching sequence). However, HatchSecond eggs have a mean volume
> slightly larger than that of HatchFirst eggs (80,02 vs. 79,47) which doesn't
> reconcile with my GLMM: HatchSecond eggs from 2007 are 3.41 units lower than
> that of HathFirst eggs.
>
> That said, I was wondering if this differences is due to the fact that the
> GLMM includes a random effect (Nest) while plotting raw data ignores it.
>
> Thank you very much in advance!
>
> LFLS
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.