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)
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