wild animal model implementation with brms
Sasha Mikheyev
9/26/2018
MCMCglmm following tutorial
Data <- as.data.frame(read.table(file = "./gryphon.txt", header = TRUE))
names(Data)[1] <- "animal"
Data$animal <- as.factor(Data$animal)
Data$MOTHER <- as.factor(Data$MOTHER)
Data$BYEAR <- as.factor(Data$BYEAR)
Data$SEX <- as.factor(Data$SEX)
Data$BWT <- as.numeric(Data$BWT)
Data$TARSUS <- as.numeric(Data$TARSUS)
head(Data)## animal MOTHER BYEAR SEX BWT TARSUS
## 1 1029 1145 968 1 10.77 24.77
## 2 1299 811 968 1 9.30 22.46
## 3 643 642 970 2 3.98 12.89
## 4 1183 1186 970 1 5.39 20.47
## 5 1238 1237 970 2 12.12 NA
## 6 891 895 970 1 NA NAPed <- as.data.frame(read.table(file = "./gryphonped.txt", header = TRUE))
for (x in 1:3) Ped[, x] <- as.factor(Ped[, x])
head(Ped)## ID FATHER MOTHER
## 1 1306 <NA> <NA>
## 2 1304 <NA> <NA>
## 3 1298 <NA> <NA>
## 4 1293 <NA> <NA>
## 5 1290 <NA> <NA>
## 6 1288 <NA> <NA>prior1.1 <- list(G = list(G1 = list(V = 1, nu = 0.002)), R = list(V = 1, nu = 0.002))
model1.1 <- MCMCglmm(BWT ~ 1, random = ~animal, pedigree = Ped, data = Data, prior = prior1.1)##
## MCMC iteration = 0
##
## MCMC iteration = 1000
##
## MCMC iteration = 2000
##
## MCMC iteration = 3000
##
## MCMC iteration = 4000
##
## MCMC iteration = 5000
##
## MCMC iteration = 6000
##
## MCMC iteration = 7000
##
## MCMC iteration = 8000
##
## MCMC iteration = 9000
##
## MCMC iteration = 10000
##
## MCMC iteration = 11000
##
## MCMC iteration = 12000
##
## MCMC iteration = 13000##
## Iterations = 3001:12991
## Thinning interval = 10
## Sample size = 1000
##
## DIC: 3912.555
##
## G-structure: ~animal
##
## post.mean l-95% CI u-95% CI eff.samp
## animal 3.402 2.21 4.607 185.8
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 3.847 2.865 4.919 218
##
## Location effects: BWT ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 7.588 7.316 7.858 1000 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## animal units
## 3.335676 3.475024prior2.1 <- list(G = list(G1 = list(V = diag(2), n = 1.002)),
R = list(V = diag(2), n = 1.002))
model2.1 <- MCMCglmm(cbind(BWT, TARSUS) ~ trait - 1, random = ~us(trait):animal,
rcov = ~us(trait):units, family = c("gaussian", "gaussian"),
pedigree = Ped, data = Data, prior = prior2.1)## Warning: 'cBind' is deprecated.
## Since R version 3.2.0, base's cbind() should work fine with S4 objects##
## MCMC iteration = 0
##
## Acceptance ratio for liability set 1 = 0.000871
##
## MCMC iteration = 1000
##
## Acceptance ratio for liability set 1 = 0.322345
##
## MCMC iteration = 2000
##
## Acceptance ratio for liability set 1 = 0.308924
##
## MCMC iteration = 3000
##
## Acceptance ratio for liability set 1 = 0.308784
##
## MCMC iteration = 4000
##
## Acceptance ratio for liability set 1 = 0.349901
##
## MCMC iteration = 5000
##
## Acceptance ratio for liability set 1 = 0.346181
##
## MCMC iteration = 6000
##
## Acceptance ratio for liability set 1 = 0.346789
##
## MCMC iteration = 7000
##
## Acceptance ratio for liability set 1 = 0.347614
##
## MCMC iteration = 8000
##
## Acceptance ratio for liability set 1 = 0.352807
##
## MCMC iteration = 9000
##
## Acceptance ratio for liability set 1 = 0.341111
##
## MCMC iteration = 10000
##
## Acceptance ratio for liability set 1 = 0.353094
##
## MCMC iteration = 11000
##
## Acceptance ratio for liability set 1 = 0.344895
##
## MCMC iteration = 12000
##
## Acceptance ratio for liability set 1 = 0.336404
##
## MCMC iteration = 13000
##
## Acceptance ratio for liability set 1 = 0.345690##
## Iterations = 3001:12991
## Thinning interval = 10
## Sample size = 1000
##
## DIC: 7925.144
##
## G-structure: ~us(trait):animal
##
## post.mean l-95% CI u-95% CI eff.samp
## traitBWT:traitBWT.animal 3.312 2.0861 4.581 144.37
## traitTARSUS:traitBWT.animal 2.430 0.1883 4.540 89.75
## traitBWT:traitTARSUS.animal 2.430 0.1883 4.540 89.75
## traitTARSUS:traitTARSUS.animal 11.972 7.0739 18.930 89.85
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitBWT:traitBWT.units 3.930 2.985 4.942 174.1
## traitTARSUS:traitBWT.units 3.373 1.548 5.269 111.8
## traitBWT:traitTARSUS.units 3.373 1.548 5.269 111.8
## traitTARSUS:traitTARSUS.units 18.213 13.180 23.254 102.5
##
## Location effects: cbind(BWT, TARSUS) ~ trait - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitBWT 7.591 7.318 7.857 623.5 <0.001 ***
## traitTARSUS 20.546 19.997 21.093 872.5 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
brms, following vignette
inv.phylo <- MCMCglmm::inverseA(Ped)
A <- solve(inv.phylo$Ainv)
rownames(A) <- rownames(inv.phylo$Ainv)
model_simple1.1 <- brm(
BWT ~ 1 + (1|animal), data = Data,
family = gaussian(), cov_ranef = list(animal = A),
chains = 2, cores = 2, iter = 1000
)## Warning: Rows containing NAs were excluded from the model.## Compiling the C++ model## Start sampling## Warning: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low## Warning: Examine the pairs() plot to diagnose sampling problems## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: BWT ~ 1 + (1 | animal)
## Data: Data (Number of observations: 854)
## Samples: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
## total post-warmup samples = 1000
##
## Group-Level Effects:
## ~animal (Number of levels: 854)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 1.86 0.17 1.54 2.19 74 1.00
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 7.58 0.14 7.31 7.86 451 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 1.94 0.13 1.67 2.20 74 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Bivariate model
model_simple2.1 <- brm(
cbind(BWT, TARSUS) ~ 1 + (1|animal), data = Data,
family = gaussian(), cov_ranef = list(animal = A),
chains = 2, cores = 2, iter = 1000
)## Setting 'rescor' to TRUE by default for this model## Warning: Rows containing NAs were excluded from the model.## Compiling the C++ model## Start sampling## Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low## Warning: Examine the pairs() plot to diagnose sampling problems## Family: MV(gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: BWT ~ 1 + (1 | animal)
## TARSUS ~ 1 + (1 | animal)
## Data: Data (Number of observations: 683)
## Samples: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
## total post-warmup samples = 1000
##
## Group-Level Effects:
## ~animal (Number of levels: 683)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(BWT_Intercept) 1.60 0.18 1.23 1.94 71 1.00
## sd(TARSUS_Intercept) 2.90 0.43 1.97 3.69 89 1.01
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## BWT_Intercept 7.48 0.15 7.17 7.77 1000 1.00
## TARSUS_Intercept 20.35 0.29 19.75 20.90 1000 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma_BWT 2.13 0.13 1.88 2.36 79 1.00
## sigma_TARSUS 4.59 0.26 4.07 5.10 112 1.01
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## rescor(BWT,TARSUS) 0.54 0.05 0.45 0.63 133 1.01
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## $animal
## $animal$sd
## Estimate Est.Error Q2.5 Q97.5
## BWT_Intercept 1.599120 0.1830747 1.233210 1.942857
## TARSUS_Intercept 2.896177 0.4262761 1.974863 3.691098
##
##
## $residual__
## $residual__$sd
## Estimate Est.Error Q2.5 Q97.5
## BWT 2.125112 0.1283578 1.884407 2.363342
## TARSUS 4.592198 0.2577639 4.070870 5.100268
##
## $residual__$cor
## , , BWT
##
## Estimate Est.Error Q2.5 Q97.5
## BWT 1.0000000 0.00000000 1.0000000 1.0000000
## TARSUS 0.5402788 0.04566532 0.4526675 0.6294334
##
## , , TARSUS
##
## Estimate Est.Error Q2.5 Q97.5
## BWT 0.5402788 0.04566532 0.4526675 0.6294334
## TARSUS 1.0000000 0.00000000 1.0000000 1.0000000
##
##
## $residual__$cov
## , , BWT
##
## Estimate Est.Error Q2.5 Q97.5
## BWT 4.532561 0.5450292 3.550989 5.585384
## TARSUS 5.264724 0.5549383 4.258810 6.402712
##
## , , TARSUS
##
## Estimate Est.Error Q2.5 Q97.5
## BWT 5.264724 0.5549383 4.25881 6.402712
## TARSUS 21.154662 2.3686877 16.57199 26.012740