multiresponse MCMCglmm and location
Rafael Maia
8/21/2017
Simulate data: 3 continuous variables + 1 categorical variable with 2 levels. Categorical variable has effect on 2 of the continuous variables, in opposite directions.
dat <- data.frame(
x = c(rnorm(50,-0.5,1), rnorm(50,0.5,1)),
y = rnorm(100),
z = c(rnorm(50, 0.5, 1), rnorm(50,-0.5,1)),
sex = factor(rep(c('F','M'), each=50))
)
aggregate(dat[,c('x','y','z')], dat[, 'sex', drop=FALSE], mean)## sex x y z
## 1 F -0.7442389 -0.0156031 0.4421233
## 2 M 0.3870102 -0.1305658 -0.4977732
Univariate linear models on each variable individually
summary(lm(cbind(x,y,z)~sex, dat))## Response x :
##
## Call:
## lm(formula = x ~ sex, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.63463 -0.71305 0.05548 0.65474 1.88706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.7442 0.1320 -5.638 1.66e-07 ***
## sexM 1.1312 0.1867 6.059 2.54e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9335 on 98 degrees of freedom
## Multiple R-squared: 0.2725, Adjusted R-squared: 0.2651
## F-statistic: 36.72 on 1 and 98 DF, p-value: 2.541e-08
##
##
## Response y :
##
## Call:
## lm(formula = y ~ sex, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.63787 -0.71225 -0.09715 0.56246 2.81944
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0156 0.1378 -0.113 0.910
## sexM -0.1150 0.1948 -0.590 0.556
##
## Residual standard error: 0.9741 on 98 degrees of freedom
## Multiple R-squared: 0.003541, Adjusted R-squared: -0.006627
## F-statistic: 0.3482 on 1 and 98 DF, p-value: 0.5565
##
##
## Response z :
##
## Call:
## lm(formula = z ~ sex, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1239 -0.6468 0.0244 0.6707 3.0170
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4421 0.1447 3.056 0.00289 **
## sexM -0.9399 0.2046 -4.594 1.29e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.023 on 98 degrees of freedom
## Multiple R-squared: 0.1772, Adjusted R-squared: 0.1688
## F-statistic: 21.11 on 1 and 98 DF, p-value: 1.295e-05Multivariate analysis of variance
summary(manova(lm(cbind(x,y,z)~sex, dat)))## Df Pillai approx F num Df den Df Pr(>F)
## sex 1 0.36524 18.413 3 96 1.624e-09 ***
## Residuals 98
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MCMCglmm: interaction between trait:sex. Estimates match parameters of the data pretty well.
mod1 <- MCMCglmm(
fixed = cbind(x,y,z) ~ trait:sex-1,
rcov = ~us(trait):units,
family = rep('gaussian',3),
data = dat,
verbose=FALSE, nit=110000, burnin=10000, thin=100
)
summary(mod1)##
## Iterations = 10001:109901
## Thinning interval = 100
## Sample size = 1000
##
## DIC: 853.9852
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitx:traitx.units 0.90905 0.6892 1.1979 1123
## traity:traitx.units 0.03727 -0.1431 0.2472 1000
## traitz:traitx.units -0.05197 -0.2458 0.1560 1000
## traitx:traity.units 0.03727 -0.1431 0.2472 1000
## traity:traity.units 0.98986 0.7062 1.2686 1000
## traitz:traity.units -0.06285 -0.2807 0.1274 1000
## traitx:traitz.units -0.05197 -0.2458 0.1560 1000
## traity:traitz.units -0.06285 -0.2807 0.1274 1000
## traitz:traitz.units 1.09162 0.8110 1.4324 1000
##
## Location effects: cbind(x, y, z) ~ trait:sex - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitx:sexF -0.74307 -1.01307 -0.47027 1000.0 <0.001 ***
## traity:sexF -0.01878 -0.28573 0.24697 1236.0 0.900
## traitz:sexF 0.44587 0.14354 0.71499 1000.0 0.004 **
## traitx:sexM 0.38945 0.10639 0.65853 905.4 0.006 **
## traity:sexM -0.13676 -0.40467 0.13872 1000.0 0.346
## traitz:sexM -0.49681 -0.80936 -0.22752 1135.8 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MCMCglmm: without interaction.
mod2 <- MCMCglmm(
fixed = cbind(x,y,z) ~ trait+sex-1,
rcov = ~us(trait):units,
family = rep('gaussian',3),
data = dat,
verbose=FALSE, nit=110000, burnin=10000, thin=100
)
summary(mod2)##
## Iterations = 10001:109901
## Thinning interval = 100
## Sample size = 1000
##
## DIC: 895.2262
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitx:traitx.units 1.2073 0.8582 1.63443 1000
## traity:traitx.units -0.0107 -0.2371 0.22755 1000
## traitz:traitx.units -0.3219 -0.6074 -0.07926 1000
## traitx:traity.units -0.0107 -0.2371 0.22755 1000
## traity:traity.units 0.9903 0.7313 1.25922 1141
## traitz:traity.units -0.0103 -0.2705 0.23004 1000
## traitx:traitz.units -0.3219 -0.6074 -0.07926 1000
## traity:traitz.units -0.0103 -0.2705 0.23004 1000
## traitz:traitz.units 1.3522 0.9498 1.78361 1201
##
## Location effects: cbind(x, y, z) ~ trait + sex - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitx -0.20822 -0.47763 0.03185 1093 0.112
## traity -0.10261 -0.34734 0.14084 867 0.384
## traitz -0.06230 -0.32037 0.19968 1000 0.660
## sexM 0.06354 -0.22760 0.33710 1000 0.648Is there any way to obtain an overall location effect using MCMCglmm?
- model comparison with/without sex using DIC not recommended, as far as I can tell
- bayes factor might have the same issue (plus harmonic mean problem)
second example: no univariate effects but multivariate location difference
set.seed(1328)
dat2 <- data.frame(
x = c(rnorm(50,0,1), rnorm(50,0.2,1)),
y = c(rnorm(50,0,1), rnorm(50,-0.22,1)),
z = c(rnorm(50,0,1), rnorm(50,0.23,1)),
sex = factor(rep(c('F','M'), each=50))
)
aggregate(dat2[,c('x','y','z')], dat2[, 'sex', drop=FALSE], mean)## sex x y z
## 1 F -0.24423886 -0.0156031 -0.05787671
## 2 M 0.08701024 -0.3505658 0.23222679
Univariate linear models on each variable individually
summary(lm(cbind(x,y,z)~sex, dat2))## Response x :
##
## Call:
## lm(formula = x ~ sex, data = dat2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.63463 -0.71305 0.05548 0.65474 1.88706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2442 0.1320 -1.850 0.0673 .
## sexM 0.3312 0.1867 1.774 0.0791 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9335 on 98 degrees of freedom
## Multiple R-squared: 0.03112, Adjusted R-squared: 0.02124
## F-statistic: 3.148 on 1 and 98 DF, p-value: 0.07912
##
##
## Response y :
##
## Call:
## lm(formula = y ~ sex, data = dat2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.63787 -0.71225 -0.09715 0.56246 2.81944
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0156 0.1378 -0.113 0.9101
## sexM -0.3350 0.1948 -1.719 0.0887 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9741 on 98 degrees of freedom
## Multiple R-squared: 0.02928, Adjusted R-squared: 0.01938
## F-statistic: 2.956 on 1 and 98 DF, p-value: 0.0887
##
##
## Response z :
##
## Call:
## lm(formula = z ~ sex, data = dat2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1239 -0.6468 0.0244 0.6707 3.0170
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.05788 0.14467 -0.400 0.690
## sexM 0.29010 0.20459 1.418 0.159
##
## Residual standard error: 1.023 on 98 degrees of freedom
## Multiple R-squared: 0.0201, Adjusted R-squared: 0.01011
## F-statistic: 2.011 on 1 and 98 DF, p-value: 0.1594Multivariate analysis of variance
summary(manova(lm(cbind(x,y,z)~sex, dat2)))## Df Pillai approx F num Df den Df Pr(>F)
## sex 1 0.078261 2.717 3 96 0.04891 *
## Residuals 98
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MCMCglmm: interaction between trait:sex. Estimates match parameters of the data pretty well.
mod1 <- MCMCglmm(
fixed = cbind(x,y,z) ~ trait:sex-1,
rcov = ~us(trait):units,
family = rep('gaussian',3),
data = dat2,
verbose=FALSE, nit=110000, burnin=10000, thin=100
)
summary(mod1)##
## Iterations = 10001:109901
## Thinning interval = 100
## Sample size = 1000
##
## DIC: 853.9852
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitx:traitx.units 0.90905 0.6892 1.1979 1123
## traity:traitx.units 0.03727 -0.1431 0.2472 1000
## traitz:traitx.units -0.05197 -0.2458 0.1560 1000
## traitx:traity.units 0.03727 -0.1431 0.2472 1000
## traity:traity.units 0.98986 0.7062 1.2686 1000
## traitz:traity.units -0.06285 -0.2807 0.1274 1000
## traitx:traitz.units -0.05197 -0.2458 0.1560 1000
## traity:traitz.units -0.06285 -0.2807 0.1274 1000
## traitz:traitz.units 1.09162 0.8110 1.4324 1000
##
## Location effects: cbind(x, y, z) ~ trait:sex - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitx:sexF -0.24307 -0.51307 0.02973 1000.0 0.076 .
## traity:sexF -0.01878 -0.28573 0.24697 1236.0 0.900
## traitz:sexF -0.05413 -0.35646 0.21499 1000.0 0.728
## traitx:sexM 0.08945 -0.19361 0.35853 905.4 0.494
## traity:sexM -0.35676 -0.62467 -0.08128 1000.0 0.010 **
## traitz:sexM 0.23319 -0.07936 0.50248 1135.8 0.110
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MCMCglmm: without interaction.
mod2 <- MCMCglmm(
fixed = cbind(x,y,z) ~ trait+sex-1,
rcov = ~us(trait):units,
family = rep('gaussian',3),
data = dat2,
verbose=FALSE, nit=110000, burnin=10000, thin=100
)
summary(mod2)##
## Iterations = 10001:109901
## Thinning interval = 100
## Sample size = 1000
##
## DIC: 857.3208
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitx:traitx.units 0.91913 0.6655 1.1970 1000
## traity:traitx.units 0.01139 -0.1808 0.2148 1000
## traitz:traitx.units -0.03085 -0.2288 0.1732 1000
## traitx:traity.units 0.01139 -0.1808 0.2148 1000
## traity:traity.units 1.03061 0.7484 1.3088 1136
## traitz:traity.units -0.08078 -0.2749 0.1559 1000
## traitx:traitz.units -0.03085 -0.2288 0.1732 1000
## traity:traitz.units -0.08078 -0.2749 0.1559 1000
## traitz:traitz.units 1.09430 0.8297 1.4658 1320
##
## Location effects: cbind(x, y, z) ~ trait + sex - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitx -0.12932 -0.34587 0.09136 1000 0.280
## traity -0.23289 -0.49557 -0.03634 1000 0.044 *
## traitz 0.03248 -0.18697 0.26328 1000 0.806
## sexM 0.10434 -0.12052 0.34649 1000 0.386
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1