[R-sig-ME] Coefficients interpretation and plot
Luciano La Sala
lucianolasala at yahoo.com.ar
Wed Dec 29 21:17:02 CET 2010
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
More information about the R-sig-mixed-models
mailing list