[R-sig-ME] Interpreting fixed effects of mcmcglmm interaction and randomisation
Rob Griffin
robgriffin247 at hotmail.com
Fri May 27 12:44:26 CEST 2016
When you suggest "Mort*SexF" do you mean
mod1 = MCMCglmm(Life ~ Mort*Sex,
... such that the "-1" is removed from my original model specification?
Furthermore, if I added a second categorical fixed effect (e.g. whether the individual had 0, 1, or >1 older siblings [ie. whether it was the first, second, or a later offspring]) would the calculations for predictions be
first born & female: (b1 + b2) * xfirst born & male: (b1 + b3) + (b2 + b6) * xsecond born & female: (b1 + b5) + (b2 + b8) * xsecond born & male: (b1 + b3 + b5 + b10) + (b2 + b6 + b8 + b12) * xlater & female: (b1 + b4) + (b2 + b7) * xlater & male: (b1 + b3 + b4 + b9) + (b2 + b6 + b7 + b11) * x
where
b1 = (Intercept), b2 = Mort, b3 = SexM, b4 = GroupLater, b5 = GroupSecond, b6 = Mort:SexM, b7 = Mort:GroupLater, b8 = Mort:GroupSecond, b9 = SexM:GroupLater, b10 = SexM:GroupSecond, b11 = Mort:SexM:GroupLater, b12 = Mort:SexM:GroupSecond
and should I include all or just the significant effects?
Cheers, Rob
----------------------------------------------------------
Robert M. GriffinPostdoctoral Researcher, University of Turkuwww.griffinevo.wordpress.com
> Subject: Re: [R-sig-ME] Interpreting fixed effects of mcmcglmm interaction and randomisation
> To: robgriffin247 at hotmail.com; r-sig-mixed-models at r-project.org
> From: j.hadfield at ed.ac.uk
> Date: Mon, 16 May 2016 17:26:11 +0100
>
> Hi,
>
> You should use ~Mort*SexF. At the moment you force the male and female
> intercepts through the origin. You will get four terms; an intercept and
> Mort effect (which are the female intercpet and slope) and a SexM and
> Mort:SexM effect (which are the deviations of the male inetcept and
> slope from the female intercept and slope. The two lines will cross at
> some value of Mort. If we denote the four effects as b1 (intercept), b2
> (Mort), b3 (SexM) and b4 (Mort:SexM) and a value of Mort as x then
>
> female prediction = b1+b2*x
> male prediction = b1+b3+(b2+b4)*x
>
> and so the value of x when female prediction = male prediction is -b3/b4
>
> You could test whether the posterior distribution of -b3/b4 lies outside
> the range of observed x in which case sexual-dimorophism always
> increases as x (if -b3/b4 is smaller than the minimum of x) increases or
> always increases as x (if -b3/b4 is larger than the maximum of x)
> decreases. If -b3/b4 lies at intermediate values of x the change in
> sexual dimorphism will flip across the range of x.
>
> Cheers,
>
> Jarrod
>
>
>
>
> On 13/05/2016 13:51, Rob Griffin wrote:
> > Dear list,
> >
> > I am trying learn how to model the relationship between an environmental factor and sexual dimorphism using MCMCglmm, and have based analysis on the data in Garratt et al 2015* - a study looking at the relationship between sexual dimorphism in lifespan, and juvenile mortality. They use two populations measured over a number of years, giving cohorts (for which juvenile mortality can be estimated), and information on the identity, lifespan, and sex of each individual. The script below** creates the dummy data I am using to develop the model which is similar to the data set they would have used (though entirely artificial - I have forced a positive relationship between juvenile mortality and male lifespan, while female lifespan is random). They fit a general linear mixed model to "assess the relationship between survival to the onset of actuarial senescence (dependent variable) and cohort-specific juvenile mortality (independent variable), while including sex as a fixed effect, year and population as random effects, with year nested within population." I have attempted to recreate this analysis (but using MCMCglmm) as the following (where 'Life' is lifespan, 'Mort' is juvenile mortality within the individuals cohort, 'Sex' is the sex of the individual, 'Pop' and 'Year' are the population and year of the individual):
> >
> >
> > mod1 = MCMCglmm(Life ~ Mort:Sex - 1,
> > random = ~ Pop:Year,
> > rcov = ~units,
> > nitt = 13000, burnin = 3000, thin = 10,
> > prior = prior1,
> > pr = T,
> > family = "gaussian",
> > start = list(QUASI = FALSE),
> > data = DF1)
> >
> >> summary(mod1)
> > Iterations = 3001:12991
> > Thinning interval = 10
> > Sample size = 1000
> >
> > DIC: 1437.558
> >
> > G-structure: ~Pop:Year
> >
> > post.mean l-95% CI u-95% CI eff.samp
> > Pop:Year 34.23 19.4 55.2 1000
> >
> > R-structure: ~units
> >
> > post.mean l-95% CI u-95% CI eff.samp
> > units 4.746 3.946 5.487 909.7
> >
> > Location effects: Life ~ Mort:Sex - 1
> >
> > post.mean l-95% CI u-95% CI eff.samp pMCMC
> > Mort:SexF 25.17 20.61 30.46 1000 <0.001 ***
> > Mort:SexM 43.81 38.82 48.71 1000 <0.001 ***
> > ---
> > Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> >
> >
> > Which brings me to my main subject(s) of interest: essentially, how to interpret the fixed effect estimates, given that I am aiming to answer the question, is juvenile mortality related to sexual dimorphism in lifespan?
> >
> > Going from 'Life ~ Mort -1' to 'Life ~ Mort:Sex -1' improves the fit of the model (DIC 1902 -> 1438), so could I conclude that there is a significant sex-specific effect/interaction with juvenile mortality on lifespan? What do the actual values of the posterior mean (post.mean) tell me - is it just that the Mort:SexM interaction effect is significantly larger than the Mort:SexF effect (CI's of the two estimates do not overlap)? Using Mort*Sex gives SexM and SexF estimates which seem to better reflect the true mean values of the population, where the post.means of Mort is ~0, SexF and SexM are ~12, and the Mort:SexM fixed effect is ~18, which makes me more inclined to think that Mort*Sex is a more appropriate model (but brings up the question what does the post.mean of ~18 for Mort:SexM actually mean?).
> >
> > I have also done randomisation where I randomise the mortality rates among cohorts for 50 chains, in the 'Mort:Sex' model this results in the Mort:SexF and Mort:SexM estimates remaining similar to the raw data model - does randomisation offer insight in to the relationship between sexual dimorphism and juvenile mortality? Using a model where 'Mort*Sex' is the fixed effect then the randomisation shifts the estimates of Mort:SexM to ~0, and that difference is absorbed in to the SexM fixed effect (because there is sexual dimorphism, but the relationship to juvenile mortality disappears). Could I effectively use this as a kind of significance test (e.g. I do 1000 randomisations, the mean Mort:SexM estimate is> the actual estimate in 11/1000 cases, thus pseudo-p is ~0.011).
> >
> > Finally, (and on a slightly different topic) I have used 'pr = T' which means I can get posterior distributions for each level of the Pop:Year combination - i.e. a posterior distribution for each cohort. Could this information be used to get an estimate of the sex-specific lifespans for each population? For example, summing the posterior distribution of 'Mort:SexM' with 'Pop:Year.A.1983' would give a distribution of the estimated lifespan for males in population A for the year 1983? I don't think this is the case - as the numbers seem too different from the data - so how could I derive cohort-sex-specific lifespans?
> >
> > Long story short - I'm looking for some insight on how to correctly define and interpret the fixed effects... Overall I'm inclined to think that 'Mort*Sex' gives the correct specification of the model, and the randomisation shows whether the environmental factor has an effect on sexual dimorphism.
> > Thanks - full script below,
> > Rob
> >
> > * Garratt et al 2015, Current Biology "High juvenile mortality is associated with sex-specific adult survival and lifespan in wild roe deer."
> >
> > **
> > #### R-SCRIPT ####
> > rm(list=ls())
> > set.seed(255)
> > library("reshape")
> > library("MCMCglmm")
> >
> > DF1 = data.frame(c(1:320), rep(c("A", "B"), each = 10), rep(1983:1998, each = 20), sample(c("M", "F"), replace = T, size = 320), rnorm(320, 12, 2), rep((1/(1+abs(rnorm(16, 2, 1 )))), each = 20))
> > colnames(DF1) = c("ID", "Pop", "Year", "Sex", "Life", "Mort")
> >
> > # Make SD related to juvenile mortality
> > DF1$Life = ifelse(DF1$Sex == "M", DF1$Life + DF1$Mort*18, DF1$Life)
> > # Mean zero, unit variance
> > DF1$Life0 = (DF1$Life-mean(DF1$Life))/sqrt(var(DF1$Life))
> >
> > prior1 = list( G = list(G1 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000)),
> > R = list(V = 1, nu = 0.002))
> >
> > mod1 = MCMCglmm(Life ~ Mort:Sex - 1,
> > random = ~ Pop:Year,
> > rcov = ~units,
> > nitt = 13000, burnin = 3000, thin = 10,
> > prior = prior1,
> > pr = T,
> > family = "gaussian",
> > start = list(QUASI = FALSE),
> > data = DF1)
> >
> > summary(mod1)
> >
> >
> >
> > ----------------------------------------------------------
> >
> > Robert M. Griffin
> > Postdoctoral Researcher, University of Turku
> > www.griffinevo.wordpress.com
> > _______________________________________________
> > 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.
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list