[R-sig-ME] MCMCglmm for recip. transplant: 1. Trouble with convergence 2. Error in multivar. model
Jarrod Hadfield
j.hadfield at ed.ac.uk
Sun Aug 14 12:26:29 CEST 2011
Hi,
Issue 1: You have only saved 40 samples from the posterior, which may
be part of the reason why the posterior looks "uninspiring".
Nevertheless, reducing the thinning interval shows that posteriors for
the random effect variances are not all mixing well and that some are
getting trapped at zero. Since you used an improper prior (nu=-2) then
this is to some degree to be expected. I often use parameter expanded
priors:
prior1.a<-list(R=list(V=1,nu=0.002),G=list(G1=list(V=1,nu=1,
alpha.mu=0, alpha.V=1000), G2=list(V=1,nu=1, alpha.mu=0,
alpha.V=1000), G3=list(V=1,nu=1, alpha.mu=0, alpha.V=1000)))
The mixing is reasonable (but could be run for longer), and the
posterior modes are close to lmer's ML estimates (especially on the
standard deviations rather than the variances) which may be reassuring.
Issue 2: See my previous post today in reply to "MCMCglmm rcov term",
and of the CourseNotes on multiresponse models.
Cheers,
Jarrod
On 5 Aug 2011, at 19:31, Steven Brady wrote:
> Dear List:
>
> I am analyzing results from a reciprocal transplant experiment
> conducted between two types of environments (r and w). Transplants
> were made in a pairwise fashion: each of five "r" sites (r1:r5) was
> paired with each of five "w" sites (w1:w5). Eggs from each of five
> different clutches were selected from each site and grown out both
> locally and in the transplant environment. A single clutch
> comprises the experimental unit; two clutches (one home; one away)
> comprise each block.
>
> The goal of the experiment is to evaluate the influence of origin
> type, environment type, and their interaction (aka "G X E
> interaction") on the response variables, with the local adaptation
> hypothesis in mind. Thus, the goal of this analysis is to test the
> significance of these two fixed effects and their interaction. I
> would like to use MCMCglmm to do so in order to evaluate a suite of
> response variables using a mixed-MANOVA approach. For this example,
> two response variables are included (1. survival [which is
> binomial], 2. growth). Analogous data are pasted below. I am
> running R 2.11.1 GUI 1.34 Leopard build 32-bit (5589) with MCMCglmm
> v. 2.13.
>
> #begin load data
> data <-
> structure(list(origin = structure(c(1L, 1L, 1L, 1L, 1L, 6L, 6L,
> 6L, 6L, 6L, 1L, 1L, 1L, 1L, 1L, 6L, 6L, 6L, 6L, 6L, 2L, 2L, 2L,
> 2L, 2L, 7L, 7L, 7L, 7L, 7L, 2L, 2L, 2L, 2L, 2L, 7L, 7L, 7L, 7L,
> 7L, 3L, 3L, 3L, 3L, 3L, 8L, 8L, 8L, 8L, 8L, 3L, 3L, 3L, 3L, 3L,
> 8L, 8L, 8L, 8L, 8L, 4L, 4L, 4L, 4L, 4L, 9L, 9L, 9L, 9L, 9L, 4L,
> 4L, 4L, 4L, 4L, 9L, 9L, 9L, 9L, 9L, 5L, 5L, 5L, 5L, 5L, 10L,
> 10L, 10L, 10L, 10L, 5L, 5L, 5L, 5L, 5L, 10L, 10L, 10L, 10L, 10L
> ), .Label = c("r1", "r2", "r3", "r4", "r5", "w1", "w2", "w3",
> "w4", "w5"), class = "factor"), env = structure(c(1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
> 6L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 7L, 7L, 7L, 7L, 7L,
> 7L, 7L, 7L, 7L, 7L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 8L,
> 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
> 4L, 4L, 4L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 5L, 5L, 5L,
> 5L, 5L, 5L, 5L, 5L, 5L, 5L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
> 10L, 10L, 10L), .Label = c("r1", "r2", "r3", "r4", "r5", "w1",
> "w2", "w3", "w4", "w5"), class = "factor"), pair = structure(c(1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
> 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
> 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
> 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L,
> 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L,
> 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
> 5L, 5L, 5L), .Label = c("a", "b", "c", "d", "e"), class = "factor"),
> origin_type = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
> 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
> 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
> 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
> 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
> 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
> 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
> 2L, 2L), .Label = c("r", "w"), class = "factor"), env_type =
> structure(c(1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
> 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
> 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
> 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
> 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("r", "w"), class
> = "factor"),
> clutch = structure(c(1L, 2L, 3L, 4L, 5L, 26L, 27L, 28L, 29L,
> 30L, 1L, 2L, 3L, 4L, 5L, 26L, 27L, 28L, 29L, 30L, 6L, 7L,
> 8L, 9L, 10L, 31L, 32L, 33L, 34L, 35L, 6L, 7L, 8L, 9L, 10L,
> 31L, 32L, 33L, 34L, 35L, 11L, 12L, 13L, 14L, 15L, 36L, 37L,
> 38L, 38L, 39L, 11L, 12L, 13L, 14L, 15L, 36L, 37L, 38L, 38L,
> 39L, 16L, 17L, 18L, 19L, 20L, 40L, 41L, 42L, 43L, 44L, 16L,
> 17L, 18L, 19L, 20L, 40L, 41L, 42L, 43L, 44L, 21L, 22L, 23L,
> 24L, 25L, 45L, 46L, 47L, 48L, 49L, 21L, 22L, 23L, 24L, 25L,
> 45L, 46L, 47L, 48L, 49L), .Label = c("r1_1", "r1_2", "r1_3",
> "r1_4", "r1_5", "r2_1", "r2_2", "r2_3", "r2_4", "r2_5", "r3_1",
> "r3_2", "r3_3", "r3_4", "r3_5", "r4_1", "r4_2", "r4_3", "r4_4",
> "r4_5", "r5_1", "r5_2", "r5_3", "r5_4", "r5_5", "w1_1", "w1_2",
> "w1_3", "w1_4", "w1_5", "w2_1", "w2_2", "w2_3", "w2_4", "w2_5",
> "w3_1", "w3_2", "w3_3", "w3_4", "w4_1", "w4_2", "w4_3", "w4_4",
> "w4_5", "w5_1", "w5_2", "w5_3", "w5_4", "w5_5"), class = "factor"),
> block = structure(c(1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L,
> 26L, 27L, 28L, 29L, 30L, 26L, 27L, 28L, 29L, 30L, 6L, 7L,
> 8L, 9L, 10L, 6L, 7L, 8L, 9L, 10L, 31L, 32L, 33L, 34L, 35L,
> 31L, 32L, 33L, 34L, 35L, 11L, 12L, 13L, 14L, 15L, 11L, 12L,
> 13L, 14L, 15L, 36L, 37L, 38L, 39L, 40L, 36L, 37L, 38L, 39L,
> 40L, 16L, 17L, 18L, 19L, 20L, 16L, 17L, 18L, 19L, 20L, 41L,
> 42L, 43L, 44L, 45L, 41L, 42L, 43L, 44L, 45L, 21L, 22L, 23L,
> 24L, 25L, 21L, 22L, 23L, 24L, 25L, 46L, 47L, 48L, 49L, 50L,
> 46L, 47L, 48L, 49L, 50L), .Label = c("r1_a", "r1_b", "r1_c",
> "r1_d", "r1_e", "r2_a", "r2_b", "r2_c", "r2_d", "r2_e", "r3_a",
> "r3_b", "r3_c", "r3_d", "r3_e", "r4_a", "r4_b", "r4_c", "r4_d",
> "r4_e", "r5_a", "r5_b", "r5_c", "r5_d", "r5_e", "w1_a", "w1_b",
> "w1_c", "w1_d", "w1_e", "w2_a", "w2_b", "w2_c", "w2_d", "w2_e",
> "w3_a", "w3_b", "w3_c", "w3_d", "w3_e", "w4_a", "w4_b", "w4_c",
> "w4_d", "w4_e", "w5_a", "w5_b", "w5_c", "w5_d", "w5_e"), class =
> "factor"),
> surv = c(18L, 19L, 19L, 20L, 19L, 20L, 14L, 10L, 18L, 18L,
> 16L, 19L, 18L, 20L, 18L, 18L, 14L, 19L, 19L, 18L, 12L, 17L,
> 19L, 20L, 17L, 3L, 11L, 12L, 8L, 7L, 14L, 20L, 20L, 20L,
> 16L, 18L, 19L, 17L, 9L, 9L, 3L, 2L, 13L, 3L, 0L, 0L, 0L,
> 0L, 0L, 7L, 14L, 19L, 17L, 10L, 11L, 10L, 7L, 9L, 18L, 9L,
> 9L, 10L, 0L, 2L, 1L, 2L, 1L, 0L, 2L, 12L, 7L, 20L, 19L, 19L,
> 19L, 18L, 17L, 19L, 17L, 20L, 17L, 19L, 18L, 17L, 20L, 17L,
> 17L, 18L, 13L, 19L, 17L, 20L, 19L, 17L, 17L, 18L, 17L, 17L,
> 19L, 17L), stocked = c(20, 20, 20, 20, 20, 20, 20, 20, 20,
> 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
> 20, 20, 20, 20, 10, 10, 20, 20, 20, 20, 20, 20, 20, 20, 10,
> 10, 20, 20, 20, 20, 20, 10, 10, 10, 20, 10, 20, 20, 20, 20,
> 20, 10, 10, 10, 20, 10, 22, 20, 20, 20, 20, 20, 20, 20, 20,
> 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
> 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
> 20)), .Names = c("origin", "env", "pair", "origin_type",
> "env_type", "clutch", "block", "surv", "stocked"), class =
> "data.frame", row.names = c(NA,
> -100L))
> data$growth=runif(100,0.5,0.9) #add growth data; end load data
>
> #Ideally, I would like to test the multivariate response of a suite
> of traits across the GxE interaction (see "issue 2" below). But for
> my first foray in MCMCglmm, I began simply with a univariate model
> of survival containing uninformed priors:
>
> prior1<-list(R=list(V=1e-16,nu=-2),G=list(G1=list(V=1e-16,nu=-2)))
> model1<-MCMCglmm(cbind(surv,(stocked-surv))~origin_type*env_type,
> random = ~pair + clutch + block, family = "multinomial2", data =
> data, thin = 2500, prior = cbind(prior1,prior1), burnin = 30000,
> nitt=130000)
> #response variable is dichotomous: number survived, number died
> plot(model1$VCV)
> plot(model1$Sol)
>
> Issue 1: As you can see, the posterior distributions are not
> terribly convincing. When I attempt to improve this by increasing
> nitt and burnin (each by an order of magnitude), the intercept and
> coefficient posterior distributions improve, but the residual
> variance posterior distributions remain uninspiring. Would the use
> of informed priors improve convergence? If so, are there any
> general rules for determining said priors with these types of data?
>
> Issue 2: When I attempt to fit a multivariate model (the overarching
> goal of this analysis), I receive an error. e.g. given the following
> model with both survival and growth as response data:
>
> model2<-MCMCglmm(cbind(cbind(surv,(stocked-surv)),
> growth)~origin_type*env_type, random = ~pair + clutch + block,
> family = c("multinomial2", "gaussian"), data = data, thin = 2500,
> prior = cbind(prior1,prior1), burnin = 30000, nitt=130000)
>
> I receive this error:
> Error in MCMCglmm(cbind(cbind(surv, (stocked - surv)), growth) ~
> origin_type * :
> R-structure does not define unique residual for each data point
>
> I have no idea how to begin to resolve this error. Any comments
> would be greatly appreciated.
>
> Many thanks in advance,
>
> Steve Brady
>
> Steven P. Brady, Ph.D. Candidate
> School of Forestry & Environmental Studies, Yale University
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list