[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