-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path001_Vedder_et_al_m2_quail_AnalysisCode.r
102 lines (85 loc) · 3.07 KB
/
001_Vedder_et_al_m2_quail_AnalysisCode.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
# Code for "Rapid decline of prenatal maternal effects with age is independent of postnatal environment in a precocial bird"
# Published in Evolution, 2023, doi: 10.1093/evolut/qpad159
# Vedder O, Tschirren B, Postma E, Moiron M
# The code provided here is sufficient to replicate the analyses presented in the above paper
#######################################################################
# DATA ANALYSIS OF BODY MASS GROWTH DATA (Table S2)
#######################################################################
# Load packages
library(nadiv)
library(MCMCglmm)
library(tidyr)
library(dplyr)
library(ggplot2)
# Load phenotypic data
Data <- read.table("growth_model.txt", header=T)
# Response variable
Data$mass = as.numeric(Data$mass)
hist(Data$mass)
# Fixed effects
Data$diet=as.factor(Data$chick.diet)
Data$sex= as.factor(Data$chick.sex)
Data$year = as.factor(Data$year)
Data$f = as.numeric(Data$f)
Data$motherL=as.factor(Data$mother.type)
Data$fatherL=as.factor(Data$father.type)
Data$motherR=as.factor(Data$mother.replicate)
Data$fatherR=as.factor(Data$father.replicate)
Data$ageC=as.factor(Data$age)
# Random effects
Data$ID = as.factor(Data$Offspring.ID)
Data$animal = as.factor(Data$Offspring.ID)
Data$pair=as.factor(paste(Data$father.ID, Data$mother.ID))
Data$mother.ID = as.factor(Data$mother.ID)
Data$father.ID=as.factor(Data$father.ID)
# Load pedigree info
ped<- read.table("pedigree.txt",header=TRUE)
colnames(ped)[1] <- "animal"
ped3=prunePed(ped, Data$animal, make.base=TRUE)
str(ped3)
my_inverse <- inverseA(ped3)$Ainv
# Set number of samples and iterations
nsamp <- 1000
BURN <- 10000; THIN <- 10000
(NITT <- BURN + THIN*nsamp)
# Set prior
prior <- list(R = list(V = diag(1)*100, nu = 1.002),
G = list(G1 = list(V = diag(1)*100, nu = 1.002),
G2 = list(V = diag(1)*100, nu = 1.002),
G3 = list(V = diag(1)*5, nu = 1.002),
G4 = list(V = diag(1)*5, nu = 1.002),
G5 = list(V = diag(1)*5, nu = 1.002)))
# Run model
mod<- MCMCglmm(mass~ageC*sex+f+ageC*diet+year+motherL*fatherL*diet+motherR+fatherR,
random = ~ ID +animal+pair+mother.ID+father.ID,
rcov = ~ units,
data = Data,
prior =prior,
ginverse=list(animal = my_inverse),
family = "gaussian",
nitt = NITT, thin = THIN, burnin = BURN#,
#pr=TRUE
)
summary(mod)
# Asses convergence
effectiveSize(mod$VCV)
heidel.diag(mod$VCV)
autocorr.diag(mod$VCV)
# Fixed effects
posterior.mode(mod$Sol)
HPDinterval(mod$Sol)
#plot(mod$Sol)
# Random effects
round(posterior.mode(mod$VCV),3)
round(HPDinterval(mod$VCV), 3)
plot(mod$VCV)
# Estimate heritability
h2=mod$VCV[,"animal"]/rowSums(mod$VCV)
posterior.mode(h2)
HPDinterval(h2)
plot(h2)
# Estimate repeatability
rpt=(mod$VCV[,"animal"]+mod$VCV[,"ID"])/rowSums(mod$VCV)
posterior.mode(rpt)
HPDinterval(rpt)
plot(rpt)