[R-sig-ME] different results for fixed terms

Cristian Pasquaretta cristian.pasquaretta at unipv.it
Mon Nov 14 23:15:02 CET 2011


Dear List,

I am using MCMCglmm to weight for variance/covariance matrix between 3
correlated responce variables and in order to find wich predictor
variable better explain my three responces.

In paticular I have these three constrained y variable: ad, fid, getin
where ad>fid>getin always

I created a simple prior cause variance of these variable are similar:
prior = list(R = list(V = diag(3), nu = 0.01), G = list(G1 = list(V =
diag(3), nu = 0.01)))

after I pllied the following command adding 3 different explanatory
variable to the model (activity,status,sex) removing the intercept and
inserting interaction with traits.

m1 <- MCMCglmm(
fixed = cbind(ad, fid,sqrt(getin)) ~ trait:(activity+status+sex) - 1,
random = ~ us(trait):ID ,prior=prior,
rcov = ~ us(trait):units,
family = c("gaussian", "gaussian","gaussian"), nitt = 60000, burnin = 10000,
thin = 25, data = dati2)

results for the fixed factors are:
                       post.mean  l-95% CI  u-95% CI eff.samp  pMCMC
traitad:activityFOR     90.93107  78.82041 102.19899     2000 <5e-04 ***
traitfid:activityFOR    80.81136  69.04766  90.83123     2132 <5e-04 ***
traitgetin:activityFOR   6.10295   5.52044   6.68703     2000 <5e-04 ***
traitad:activityLA      90.69215  75.36094 105.44716     2000 <5e-04 ***
traitfid:activityLA     68.69217  55.31163  83.87473     2000 <5e-04 ***
traitgetin:activityLA    6.91179   6.14327   7.64156     2000 <5e-04 ***
traitad:statusS        -12.77827 -28.25378   4.13464     2000  0.131
traitfid:statusS       -16.80222 -32.58547  -2.07351     2231  0.028 *
traitgetin:statusS      -0.87363  -1.76403  -0.06164     2000  0.046 *
traitad:sexM             0.98442 -15.01548  16.52318     2000  0.908
traitfid:sexM            2.18485 -12.58677  16.19221     2000  0.762
traitgetin:sexM         -0.10102  -0.87506   0.68919     2141  0.775

but if I run a similar model just moving the fixed effects
m2 <- MCMCglmm(
fixed = cbind(ad, fid,sqrt(getin)) ~ trait:(sex+status+activity) - 1,
random = ~ us(trait):ID ,prior=prior,
rcov = ~ us(trait):units,
family = c("gaussian", "gaussian","gaussian"), nitt = 60000, burnin = 10000,
thin = 25, data = dati2)
summary(m2)

I get completely different results:

                       post.mean   l-95% CI   u-95% CI eff.samp  pMCMC
traitad:sexF           91.233095  78.717662 101.768045     1985 <5e-04 ***
traitfid:sexF          81.134272  70.233044  92.098223     1847 <5e-04 ***
traitgetin:sexF         6.115784   5.476070   6.687980     2000 <5e-04 ***
traitad:sexM           91.809339  78.188258 105.332424     2000 <5e-04 ***
traitfid:sexM          82.928854  70.979701  96.565077     2000 <5e-04 ***
traitgetin:sexM         6.016756   5.351982   6.750056     2000 <5e-04 ***
traitad:statusS       -12.998756 -29.531207   4.474341     2000  0.121
traitfid:statusS      -17.046561 -33.579056  -2.252378     2000  0.034 *
traitgetin:statusS     -0.883101  -1.721280   0.006665     2000  0.044 *
traitad:activityLA     -0.202191 -16.168013  15.503745     2000  0.988
traitfid:activityLA   -12.191033 -26.588281   3.059936     2176  0.098 .
traitgetin:activityLA   0.805456   0.072213   1.476424     1819  0.029 *

Is it possible that MCMCglmm is influenced by the order of the fixed effects?
Is it possible to have an absolute weight as we had in "lme package"
while digiting anova(mod,type="marginal")?

Maybe my model is wrong and I apologize if there are some mistakes due
to my inexperience using this method.

Thank you

Cristian Pasquaretta




More information about the R-sig-mixed-models mailing list