[R-sig-ME] Coefficients interpretation and plot

[Luciano La Sala]

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