[R-sig-ME] binary trait correlation across environments (experimental trials) using MCMCglmm?
Plough, Louis
lplough @ending from umce@@edu
Sun Jan 6 20:40:23 CET 2019
Hi Jarrod et al.
Thanks for these suggestions. The model coding the interaction as
~us(Trait):animal is working (running), but convergence appears to be
a way off even after 2 or 3 million iterations, especially for the
TrialB variance ( survival in the second environment). Would another
family (ordinal?) be worth considering over the threshold family? Or
changing the prior somehow, e.g. bumping up the alpha.V=diag(2)*1000?
We do expect some positive correlation between the two trials...i.e.
should not be zero.
Coding the model the other way you suggested as animal+animal:Trait
gives me the following error:
> bi_model_trial_pos <- MCMCglmm(phen ~ 1, random = ~animal+animal:Trial, family = "threshold",prior = prior_b, pedigree = pedigree.t, data = trials_up, nitt = 5e+05, burnin = 25000, thin = 20)
Error in buildZ(rmodel.terms[r], data = data, nginverse = names(ginverse)) :
interactions not permitted between standard random effects and those
associated with ginverse
I tried used a similar prior as for the above model, but changed the
number of structures... (not sure if that is the right prior here):
> prior_b <- list(R = list(V = 1, fix=1), G = list(G1 = list(V = 1, nu = 2),G2=list(V=1,alpha.mu=0, nu=2)))
For either model, I have been studying your tutorials for how to tweak
the priors in a sensible way (for a binary trait), but I still
struggle with how the various supplementary parameters alter the prior
distribution or CDF... you had mentioned that my previous priors were
'very informative', but in which way? Towards zero or too high?
LVP
On Thu, Jan 3, 2019 at 4:24 AM HADFIELD Jarrod <j.hadfield using ed.ac.uk> wrote:
>
> Hi,
>
> The prior is not fixed at one. Also, the priors you are using are quite informative. I would use
>
> prior <- list(R = list(V = 1, fix=1), G = list(G1 = list(V = diag(2), nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*100)))
>
> Your equation for the genetic correlation is correct.
>
> You shouldn't expect the correlation in family means to equal the model based estimate of the genetic correlation for a number of reasons; most importantly they are on different scales (data versus latent) and the variances of the family means contains residual variation but the covariance doesn't so the correlation in family means is a (downwardly) biased estimator.
>
> Cheers,
>
> Jarrod
>
>
>
>
>
>
> On 02/01/2019 22:18, Plough, Louis wrote:
>
> Thanks Jarrod for the prompt response!
>
> I think the residual variance was set at 1. My prior looks like this:
>
> prior <- list(R = list(V = 1, nu = 3), G = list(G1 = list(V = diag(2), nu = 2)))
>
> Not sure why the difference for nu in R vs G (3 vs 2)...might have been a typo. Is there guidance on setting the gamma parameter for this kind of binary trait, cross environment correlation model with the 'threshold' family? In the low iteration (toy) run I did for the R-Sig-ME post, the correlation was a bit lower than I would expect based on simply phenotypic correlation of family means, so I think I might need to tweak nu for both?
>
> Can I also confirm that I am estimating the rG correctly with the following code:
>
> corr.gen<-model_trial4$VCV[,'TrialA:TrialB.animal']/sqrt(model_trial4$VCV[,'TrialA:TrialA.animal']*model_trial4$VCV[,'TrialB:TrialB.animal'])
>
> Any reason to use 'TrialA:TrialB.animal' vs 'TrialB:TrialA.animal' in the numerator? Or are they basically equivalent? I wouldn't want add them, would I?
>
> Thanks for your help!!!
>
> LVP
>
>
>
> On Wed, Jan 2, 2019 at 3:17 PM HADFIELD Jarrod <j.hadfield using ed.ac.uk> wrote:
>>
>> Hi,
>>
>> Your model bi_model_trial is the correct one. However you must fix the residual variance at one in the prior. You have estimated the genetic (co)variance matrix for the two trials, from which you can obtain the genetic correlation. Alternatively, you could fit animal+animal:Trial which assumes the genetic variances in the two trials are the same and the correlation is positive. The correlation in this latter model is obtained as VAR(animal)/(VAR(animal)+VAR(animal:Trial)). Also, there seems to be a problem with your Surv data as it has three levels rather than 2 and so a cutpoint is being estimated.
>>
>> Cheers,
>>
>> Jarrod
>>
>>
>>
>>
>>
>>
>> On 2 Jan 2019, at 16:33, Plough, Louis <lplough using umces.edu> wrote:
>>
>> bi_model_trial
>>
>>
>> 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