[R-sig-ME] MCMCglmm for recip. transplant: 1. Trouble with convergence 2. Error in multivar. model

[Steven Brady]

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