Hi all and especially Jarrod and David,

many thanks for your insights, they were very helpful.

David Duffy looked off-list at the complete data and pointed out a 
potential role of another nuisance parameter "run" characterizing the 
way the data were collected. Basically several sires are almost always 
nested within one value of that parameter. My colleague who did the 
measurements tells me that this parameter should have little or no 
bearing on the results. I have noticed that MCMCglmm with run, sire and 
dam as random factors tends to behave poorly (drifting traces of the 
random effects) even if I increased the burnin and sampling intervals 
and played a bit with the fixed priors. I suppose this arises due to the 
tight correlation between two random effects (in our case run and sire) 
so that the model has difficulties to separate them.

I have also looked at plain lmer models with offspring sex as fixed 
factor and run, sire and dam as random factors - inclusion of the run 
parameter decreased the AIC value but not dramatically so. So I now tend 
to think that, as Jarrod pointed out in another off-list conversation, 
non-standard inheritance of the studied trait seems to be a likely 
explanation of the high heritability estimates in the sire+dam model.

Cheers,
David


On 08/08/2013 10:55, Jarrod Hadfield wrote:
> Hi,
>
> The animal model is only equivalent to the dam/sire model when the 
> heritability estimate is less than one. The animal model bounds the 
> estimate at 1 (or less in the MCMCglmm models because of the priors on 
> the maternal and residual effects). If you run the analysis with 
> asreml you will get the same result - the animal model fails because 
> the animal term explains more than 100% of the variation. In a 
> dam/sire model the constraint is that the sire component cannot exceed 
> 100% of the variance, implying that the heritability is constrained to 
> be less than 4. Estimates of h2>1 can arise through sampling error, 
> but given that the lower 95% credible interval for the heritability 
> from the dam/sire  mode is >> 1, I would say that the inheritance of 
> this trait does not conform to the standard Mendelian genetics (e.g., 
> paternal effects, sex-linkage etc.) Please post these findings back to 
> the list , so others can learn.
>
> Cheers,
>
> Jarrod

On 08/08/2013 11:36, David Duffy wrote:
> Yes, your problem was that _run_ has a significant effect on _val_, 
> and is correlated with _sire_. Either including it as a fixed or 
> random effect does the trick:
>
> V_A=14.9560, V_RUN=6.11257, V_E=2.81424
>
>
> | David Duffy (MBBS PhD) ,-_|\
> | email: davidD@qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101 /     *
> | Epidemiology Unit, Queensland Institute of Medical Research \_,-._/
> | 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v


| David Duffy (MBBS PhD) ,-_|\
| email: davidD@qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101  /     *
| Epidemiology Unit, Queensland Institute of Medical Research \_,-._/
| 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v

> ------------------------------------------------------------------------
> *From:* David Boukal [boukal@entu.cas.cz]
> *Sent:* Wednesday, August 07, 2013 7:35 PM
> *To:* Jarrod Hadfield
> *Subject:* Re: Analysis of unbalanced nested half-sib/full-sib design 
> with MCMCglmm/lmer
>
> Dear Jarrod,
>
> I am attaching the data and the code. It would be great if you could 
> have a look. I have quickly tried the recommended priors and see some 
> differences but the main problem with the seeming discrepancy between 
> the results from the animal model and the sire+dam formulation (at 
> least given the way I deal with them) remains.
>
> Cheers,
> David
>
> On 07/08/2013 18:59, Jarrod Hadfield wrote:
>> Hi,
>>
>> The problem can only be really diagnosed if the data/code are made 
>> available. Is it possible that the data are not very informative with 
>> regards to separating maternal and genetic effects, and the strong 
>> prior you are placing on the genetic variance in the dam/sire model 
>> (2*Vp with 1 df) is causing problems. Those priors are not recommended:
>>
>> http://onlinelibrary.wiley.com/doi/10.1111/j.1365-2656.2011.01896.x/pdf
>>
>> Also, David's point about DIC changes under the two parameterisations 
>> highlights the problem with DIC focused at the highest level. In the 
>> animal model the deviance is conditional on breeding values, in the 
>> dam/sire model the Mendelian sampling deviation bit of the breeding 
>> value is not conditioned on. DIC is essentially focused at two 
>> different levels in the different models and is not comparable.
>>
>>
>> Cheers,
>>
>> Jarrod
>>
>> ------------------------------------------------------------------------
>> *From:* David Boukal [boukal@entu.cas.cz]
>> *Sent:* Wednesday, August 07, 2013 3:18 PM
>> *To:* alastair.wilson@ed.ac.uk; Jarrod Hadfield; Ane Laugen; Bruno 
>> ERNANDE
>> *Cc:* jussi.alho@helsinki.fi; rsardell@bio.fsu.edu; 
>> lara.hallsson@ebc.uu.se; lecastane@gmail.com; davidD@qimr.edu.au
>> *Subject:* Fwd: Analysis of unbalanced nested half-sib/full-sib 
>> design with MCMCglmm/lmer
>>
>> Dear MCMCglmm/animal model gurus and heritability experts,
>>
>> first I would like to apologize for an unsolicited mail. I have 
>> posted a question regarding the analysis of an unbalanced nested 
>> half-sib/full-sib experimental design on the R-sig-mixed-models forum 
>> (copied below). David Duffy (many thanks!) has answered some of my 
>> questions (post also copied below). However, I remain baffled by the 
>> link between the VCV components returned by R and the equations in 
>> Lynch & Walsh. This is most likely trivial once you have done it 
>> before, so I would much appreciate your feedback.
>>
>> My ignorance shows best in the comparison of the following two models 
>> examining the data from the experiment, described in some more detail 
>> below, one using the animal model pedigree and the other one ignoring it:
>>
>> m2 <- MCMCglmm(val ~ sex - 1, random = ~animal + mother, data = df, 
>> pedigree = ped, prior = prior2)
>> m2c <- MCMCglmm(val ~ sex - 1, random = ~father + mother, data = df, 
>> pedigree = NULL, prior = prior2)
>>
>> m2 gives posterior modes for VCV components:
>>     animal    mother     units
>> 17.095817  1.985225  2.546613
>>
>> and m2c gives
>>     father    mother     units
>> 11.475095  2.527233 10.192168
>>
>> So
>> h2.m2 <- m2$VCV[,"animal"] / (m2$VCV[,"animal"] + m2$VCV[,"mother"] + 
>> m2$VCV[,"units"])
>> yields posterior mode of heritability estimate equal to 0.72 
>> (slightly different due to MC error from the value I posted before) and
>>
>> h2.m2c <- m2c$VCV[,"father"] / (m2c$VCV[,"father"] + 
>> m2c$VCV[,"mother"] + m2c$VCV[,"units"])
>> yields posterior mode of 0.46 and 4*h2.m2c, based on Lynch & Walsh 
>> (p. 573) , yields 1.84.
>>
>> Question: If animal+maternal model with pedigree is equivalent to 
>> sire+dam, as suggested by David Duffy, why don't the numbers match? I 
>> suppose something is still wrong with my second heritability formula, 
>> 1.84 is too high.
>>
>> Many thanks for any clarifications - I can post their digest at the 
>> R-sig-ME mailing list for posterity,
>> David Boukal
>>
>>
>> -------- Original Message --------
>> Subject: 	Analysis of unbalanced nested half-sib/full-sib design with 
>> MCMCglmm/lmer
>> Date: 	Tue, 06 Aug 2013 15:34:11 +0200
>> From: 	David Boukal <boukal@entu.cas.cz>
>> Reply-To: 	boukal@entu.cas.cz
>> To: 	r-sig-mixed-models@r-project.org
>> CC: 	David Boukal <boukal@entu.cas.cz>
>>
>>
>>
>> Dear all,
>>
>> I am trying to implement a MCMCglmm/lmer analysis to estimate 
>> narrow-sense heritability h2 and sire/maternal+environmental effects 
>> of a single phenotypic trait in an insect. We used a nested half-sib 
>> full-sib design, treating sex of the offspring as fixed effect. This 
>> analysis is outside my main areas of interest and I have no-one to 
>> discuss it with directly, so I would much appreciate any help to help 
>> me wrap my head around the problem. Please suggest a better place to 
>> migrate the question if it does not fit here.
>>
>> The data (44 males x 3 females nested within each female, total ~430 
>> offspring) are unbalanced due to differences in fertility and 
>> offspring survival - many females yielded no or few offspring. I 
>> understand that the basic nested analysis of variance formulas 
>> provided in the Lynch & Walsh classics should not be trusted there.
>>
>> Our raw data look like this:
>> head(df)
>>    animal sex mother father    val
>> 1    201   m     64      5 -23.03
>> 2    202   m     64      5 -17.43
>> 3    203   m     64      5 -14.39
>> 4    204   f     64      5 -17.39
>> 5    205   f     65      5 -22.04
>> 6    206   f     65      5 -19.84
>>
>> The underlying theory of the linear model for the nested design in 
>> Lynch & Walsh (pp. 571-574), i.e.
>> val_ijk = mu + sex_ijk + father_i + mother_ij + e_ijk,
>> says that 4*Var(s)/Var(z) or alternatively 2*(Var(s)+Var(d))/Var(z) 
>> is the best estimate of h2.
>>
>> Question 1: How do the quantities from Lynch & Walsh relate to the 
>> VCV components provided by MCMCglmm or lmer?
>> Question 2a: Can the animal model of Wilson et al (2010, WAMWiki) be 
>> applied "as-is" with the appropriate pedigree to this nested design 
>> and what are the pros/cons of (not) doing that?
>> Question 2b: Why do the results for heritability differ when I use 
>> both father and mother as explanatory variables as opposed to using 
>> only mother and supplying the pedigree (= essentially the 
>> animal+mother+father columns of my data including the NA rows for 
>> parents)?
>> Question 3: MCMCglmm course notes suggest DIC may not be meaningful 
>> to use when priors differ. How does one select between models that 
>> differ in the number of random effects then?
>>
>> I have tried the following priors and models in MCMCglmm; all models 
>> had nitt = 65000, thin = 50, burnin = 15000, skipped for brevity:
>>
>> p.var <- var(df$val, na.rm=TRUE)
>>
>> prior1 <- list(G = list(G1 = list(V = matrix(p.var/2), n = 1)), R = 
>> list(V=matrix(p.var/2), n = 1))
>> prior2 <- list(R = list(V=matrix(p.var/2), nu = 0.2),
>>                  G = list(G1 = list(V = matrix(p.var/2), n = 1), G2 = 
>> list(V=matrix(p.var/2), n = 1))  )
>>
>> m1 <- MCMCglmm(val ~ sex - 1, random = ~animal, data = df, pedigree = 
>> ped, prior = prior1)
>> [DIC: 1854.0]
>> m2 <- MCMCglmm(val ~ sex - 1, random = ~animal + mother, data = df, 
>> pedigree = ped, prior = prior2)
>> [DIC: 1926.6]
>> m2b <- MCMCglmm(val ~ sex - 1, random = ~father + mother, data = df, 
>> pedigree = ped, prior = prior2)
>> [DIC: 2306.1]
>> m2c <- MCMCglmm(val ~ sex - 1, random = ~father + mother, data = df, 
>> pedigree = NULL, prior = prior2)
>> [DIC: 2305.9]
>>
>> The following heritability formulas:
>> h2.m1 <- m1$VCV[,"animal"] / (m1$VCV[,"animal"] + m1$VCV[,"units"])
>> h2.m2<-m2$VCV[,"animal"] / (m2$VCV[,"animal"] + m2$VCV[,"mother"] + 
>> m2$VCV[,"units"])
>> h2.m2b<-m2b$VCV[,"father"] / (m2b$VCV[,"father"] + m2b$VCV[,"mother"] 
>> + m2b$VCV[,"units"])
>> h2.m2c<-m2c$VCV[,"father"] / (m2c$VCV[,"father"] + m2c$VCV[,"mother"] 
>> + m2c$VCV[,"units"])
>>
>> give the following mode of posterior distribution of heritability 
>> estimate:
>> h2.m1 = 0.92
>> h2.m2 = 0.75
>> h2.m2b = 0.47
>> h2.m2c = 0.45
>>
>> So model m1 has the lowest DIC, but the heritability estimate is kind 
>> of high...
>>
>> Many thanks,
>> David Boukal
>> -------------------------------------------------
>> Department of Ecosystems Biology,
>> Faculty of Science, University of South Bohemia
>> 	&&
>> Department of Biosystematics and Ecology,
>> Biology Centre AS CR, Institute of Entomology
>> Branisovska 31
>> 37005 Ceske Budejovice
>> Czech Republic
>>
>> boukal@entu.cas.cz	www.entu.cas.cz/boukal
>> phone +420 387 772 327    fax   +420 385 310 354
>>
>> -------------------------------------------------
>>
>> ----------------------------------------------------------
>> On 07/08/2013 10:50, Dacvid Duffy wrote:
>>> Question 1: How do the quantities from Lynch & Walsh relate to the VCV
>>> components provided by MCMCglmm or lmer?
>>
>> Same. Note your h2.m2b and h2.m2c as written are inconsistent with this.
>>
>>> Question 2a: Can the animal model of Wilson et al (2010, WAMWiki) be
>>> applied "as-is" with the appropriate pedigree to this nested design and
>>> what are the pros/cons of (not) doing that?
>>
>> Yes. No cons.
>>
>>> Question 2b: Why do the results for heritability differ when I use both
>>> father and mother as explanatory variables as opposed to using only
>>> mother and supplying the pedigree (= essentially the
>>> animal+mother+father columns of my data including the NA rows for 
>>> parents)?
>>
>> If no animal term, ped is not used, I would think. So, m2b == m2c.
>>
>>> Question 3: MCMCglmm course notes suggest DIC may not be meaningful to
>>> use when priors differ. How does one select between models that differ
>>> in the number of random effects then?
>>
>> If you are interested in maternal effects, you should look at the 
>> confidence interval for mother in your m2.
>>
>>> So model m1 has the lowest DIC, but the heritability estimate is 
>>> kind of
>>> high...
>>
>> Eyeball the fullsib and halfsib intraclass correlations.
>>
>> More generally, for the half-sib design, the A+D (or A+maternal) 
>> model is equivalent to your sire+dam model, so that m2b and m2 should 
>> fit the data equally well in terms of likelihood (in a frequentist 
>> LMM). In that case, the DIC differences reflect how model complexity 
>> is being assessed.  I believe there are a few different ways one can 
>> calculate DIC for RE models.
>>
>> For one simulated dataset
>>
>> Paternal half-sib design
>>
>> 100 sires, nested 10 dams per sire, 5 offspring per mating, true h2=0.5
>>
>> Relation r      JSE
>> Halfsib  0.1478 0.0092
>> Fullsib  0.3254 0.0169
>>
>> Program  -2*LL or Dev A        Maternal  Sire Dam     E        DIC
>> Wombat    3136.830    0.89100  -         - -       0.56981  -
>> Wombat    3136.618    0.83105  0.025198  - -       0.59548  -
>> Wombat    3136.618    -        -         0.20777 0.23295 1.01101  -
>> MCMCglmm  5652.37     0.9163   -         - -       0.5547   6987.498
>> MCMCglmm  5733.77     0.8858   0.00596   - -       0.5731   7035.543
>> MCMCglmm  7195.74     -        -         0.2147 0.2311  1.015    
>> 7590.391
>>
>> MCMCglmm, default options.
>>
>>
>>
>> Cheers, David Duffy.
>>
>> | David Duffy (MBBS PhD)                                         ,-_|\
>> | email: davidD@qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101  /     *
>> | Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
>> | 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v
>
> -- 
> -------------------------------------------------
> David Boukal, PhD, Assoc Prof
> Department of Ecosystems Biology,
> Faculty of Science, University of South Bohemia
> 	&&
> Department of Biosystematics and Ecology,
> Biology Centre AS CR, Institute of Entomology
> Branisovska 31
> 37005 Ceske Budejovice
> Czech Republic
>
> boukal@entu.cas.cz	www.entu.cas.cz/boukal
> phone +420 387 772 327    fax   +420 385 310 354
>
> -------------------------------------------------

-- 
-------------------------------------------------
David Boukal, PhD, Assoc Prof
Department of Ecosystems Biology,
Faculty of Science, University of South Bohemia
	&&
Department of Biosystematics and Ecology,
Biology Centre AS CR, Institute of Entomology
Branisovska 31
37005 Ceske Budejovice
Czech Republic

boukal@entu.cas.cz	  www.entu.cas.cz/boukal
phone +420 387 772 327    fax   +420 385 310 354

-------------------------------------------------


	[[alternative HTML version deleted]]

