#Script for the bivariate animal model analysis #いいいいいいいいいいいいいい ####Start of the script###### #いいいいいいいいいいいいいい Exp_age<-read.table("Pheno_Data.txt",header=T) library(asreml) ped<-read.table("Trimmed_Pedigree.txt", header=T) ainv= asreml.Ainverse(ped[,c(1,3,2)])$ginv pin<-dget("pin.R") Exp_age$Year<-as.factor(Exp_age$Year) Exp_age$Ring<-as.factor(as.character(Exp_age$Ring)) Exp_age$Observer<-as.factor(Exp_age$Observer) Exp_age$BroodID<-as.factor(Exp_age$BroodID) ####################################### #Run univariate animal model for LOCOM ###################################### model.uni<- asreml(LOCOM ~ 1 + June.Days + I(June.Days^2)+ TESTAGE + I(TESTAGE^2) + Sex , random = ~ped(Ring) + Observer + Year, data = Exp_age,na.method.X="include", ginverse = list(Ring=ainv), maxiter=500) summary(model.uni) pin(model.uni, h2~V3/(V4+V3))#h2 wald.asreml(model.uni,ssType="conditional", denDF="numeric")#test significance of fixed effects model.uni$coefficients$fixed#Coefficients of fixed effects sqrt(pin(model.uni, VA~V3)[1])/mean(Exp_age$LOCOM,na.rm=T)*100#CVA model.uni2 <- asreml(LOCOM ~ 1 + June.Days + I(June.Days^2)+ TESTAGE + I(TESTAGE^2) + Sex , random = ~ ped(Ring)+ Observer , data = Exp_age,na.method.X="include", ginverse = list(Ring=ainv), maxiter=500) 1-pchisq(2*(model.uni$loglik-model.uni2$loglik),1)#Test significance of Year variance model.uni3 <- asreml(LOCOM ~ 1 + June.Days + I(June.Days^2)+ TESTAGE + I(TESTAGE^2) + Sex , random = ~ ped(Ring)+ Year , data = Exp_age,na.method.X="include", ginverse = list(Ring=ainv), maxiter=500) 1-pchisq(2*(model.uni$loglik-model.uni3$loglik),1)#Test significance of Observer variance model.uni4 <- asreml(LOCOM ~ 1 + June.Days + I(June.Days^2)+ TESTAGE + I(TESTAGE^2) + Sex , random = ~ Year+ Observer , data = Exp_age,na.method.X="include", maxiter=500) 1-pchisq(2*(model.uni$loglik-model.uni4$loglik),1)#Test significance of VA ################################################# #Run univariate animal models in young and adults ################################################# Exp_age$LOCOMYOUNG<-ifelse(Exp_age$TESTAGE <12& Exp_age$June.Days <283,Exp_age$LOCOM,NA ) Exp_age$LOCOMADULT<-ifelse(Exp_age$TESTAGE %in% c(12:98),Exp_age$LOCOM,NA) #Run the model in young model.1 <- asreml(LOCOMYOUNG ~1 + June.Days + I(June.Days^2) + Sex + TESTAGE + I(TESTAGE^2) , random = ~ped(Ring) + BroodID+ Observer + Year, data = Exp_age,na.method.X="include",na.method.Y="omit", ginverse = list(Ring=ainv), maxiter=500) summary(model.1) pin(model.1,h2~V4/(V4+V5))#h2 wald.asreml(model.1,ssType="conditional", denDF="numeric")#test significance of fixed effects model.1$coefficients$fixed#Coefficients of fixed effects sqrt(pin(model.1, VA~V4)[1])/mean(Exp_age$LOCOMYOUNG,na.rm=T)*100#CVA model.1b <- asreml(LOCOMYOUNG ~1 + June.Days + I(June.Days^2) + Sex + TESTAGE + I(TESTAGE^2) , random = ~ BroodID + Observer + Year, data = Exp_age,na.method.X="include",na.method.Y="omit", maxiter=500) 1-pchisq(2*(model.1$loglik-model.1b$loglik),1)#Test significance of VA in young #Run the model in adults model.2 <- asreml(LOCOMADULT ~1 + June.Days + I(June.Days^2) + Sex + TESTAGE + I(TESTAGE^2) , random = ~ped(Ring) + Observer +Year, data = Exp_age,na.method.X="include",na.method.Y="omit", ginverse = list(Ring=ainv), maxiter=500) summary(model.2) pin(model.2,h2~V3/(V3+V4))#h2 wald.asreml(model.2,ssType="conditional", denDF="numeric")#Test significance of fixed effects model.2$coefficients$fixed#Coefficients of fixed effects sqrt(pin(model.2, VA~V3)[1])/mean(Exp_age$LOCOMADULT,na.rm=T)*100#CVA model.2b <- asreml(LOCOMADULT ~1 + June.Days + I(June.Days^2) + Sex + TESTAGE + I(TESTAGE^2) , random = ~ Observer +Year, data = Exp_age,na.method.X="include",na.method.Y="omit", maxiter=500) 1-pchisq(2*(model.2$loglik-model.2b$loglik),1)#Test significance of VA in adults ################################ #Run the bivariate animal model ################################ model.biv<- asreml(cbind(LOCOMYOUNG,LOCOMADULT) ~ trait + trait:pol(June.Days,2)+ trait:Sex + trait:pol(TESTAGE,2), random = ~us(trait,init=c(14,4,5)):ped(Ring) + diag(trait,init=c(0.1,2)):Observer + at(trait,1):BroodID + diag(trait,init=c(0.01,0.06)):Year, rcov= ~ units:diag(trait,init=c(32,60)), data = Exp_age,na.method.X="include",na.method.Y="include", ginverse = list(Ring=ainv), maxiter=500) summary(model.biv) pin(model.biv, rA~V7/sqrt(V8*V6))#Genetic correlation across ages #Constrain the correlation to 0 model.biv.diag<- asreml(cbind(LOCOMYOUNG,LOCOMADULT) ~ trait + trait:pol(June.Days,2)+ trait:Sex + trait:pol(TESTAGE,2), random = ~diag(trait,init=c(14,10)):ped(Ring) + diag(trait,init=c(0.1,2)):Observer + at(trait,1):BroodID+ diag(trait,init=c(0.01,0.06)):Year, rcov= ~ units:diag(trait,init=c(32,60)), data = Exp_age,na.method.X="include",na.method.Y="include", ginverse = list(Ring=ainv), maxiter=500) 1-pchisq(2*abs(model.biv$loglik - model.biv.diag$loglik),1)#Test if genetic correlation <>0 #Constrain the correlation to 1 using a correlation structure (might take a bit longer) sv<- asreml(cbind(LOCOMYOUNG,LOCOMADULT) ~ trait + trait:pol(June.Days,2)+ trait:Sex + trait:pol(TESTAGE,2), random = ~corgh(trait):ped(Ring)+ at(trait,1):BroodID + diag(trait,init=c(0.1,2)):Observer + diag(trait,init=c(0.01,0.06)):Year, rcov= ~ units:diag(trait,init=c(32,60)), data = Exp_age,na.method.X="include",na.method.Y="include", ginverse = list(Ring=ainv), maxiter=500,start.values=T) gam<-sv$gammas.table gam$Value[1]<-0.999 gam$Constraint[1]<-"F" model.biv2<-asreml(cbind(LOCOMYOUNG,LOCOMADULT) ~ trait + trait:pol(June.Days,2) + trait:Sex +trait:pol(TESTAGE,2), random = ~corgh(trait):ped(Ring)+ at(trait,1):BroodID + diag(trait,init=c(0.1,2)):Observer + diag(trait,init=c(0.01,0.06)):Year, rcov= ~ units:diag(trait,init=c(32,60)), data = Exp_age,na.method.X="include",na.method.Y="include", ginverse = list(Ring=ainv), maxiter=500,G.param=gam, R.param=gam) summary(model.biv2) 1-pchisq(2*abs(model.biv$loglik - model.biv2$loglik),1)# Test if correlation <>1 #Constrain VA(or VR) to be the same in young and adults sv<- asreml(cbind(LOCOMYOUNG,LOCOMADULT) ~ trait + trait:pol(June.Days,2) + trait:Sex + trait:pol(TESTAGE,2), random = ~us(trait):ped(Ring)+ at(trait,1):BroodID + diag(trait,init=c(0.1,2)):Observer + diag(trait,init=c(0.01,0.06)):Year, rcov= ~ units:diag(trait,init=c(32,60)), data = Exp_age,na.method.X="include",na.method.Y="include", ginverse = list(Ring=ainv), maxiter=500,start.values=T) gam<-sv$gammas.table gam$fac<-factor(c(1,2,1,3,4,5,6,7,8,9,10))#option 1: same VA #gam$fac<-factor(c(1,2,3,4,5,6,7,8,9,10,10))#option 2: same VR gam$Value[c(10,11)]<-50 M <- asreml.constraints(~ fac, gammas=gam) model.3<-asreml(cbind(LOCOMYOUNG,LOCOMADULT) ~ trait + trait:pol(June.Days,2) + trait:Sex+trait:pol(TESTAGE,2), random = ~us(trait):ped(Ring)+ at(trait,1):BroodID+ diag(trait,init=c(0.1,2)):Observer + diag(trait,init=c(0.01,0.06)):Year, rcov= ~ units:diag(trait,init=c(50,50)), data = Exp_age,na.method.X="include",na.method.Y="include", ginverse = list(Ring=ainv), maxiter=500,G.param=gam, R.param=gam,constraints=M) summary(model.3) 1-pchisq(2*abs(model.biv$loglik - model.3$loglik),1)# Test if VA(or VR) are statistically different across age groups