#Script for the random regression animal model analysis #いいいいいいいいいいいいいい ####Start of the script###### #いいいいいいいいいいいいいい Exp_age<-read.table("Pheno_Data_AGE_ESM.txt",header=T) library(asreml) ped<-read.table("Pedigree_AGE_ESM.txt", header=T) ainv= asreml.Ainverse(ped[,c(1,3,2)])$ginv pin<-dget("pin.R") ######################################## #Prepare the data for RRAM analysis ######################################## #Define age as blocks: 0,1,2,>=3 years old Exp_age[Exp_age$TESTAGE< 12,"block"]<-1 Exp_age[Exp_age$TESTAGE>=12 &Exp_age$TESTAGE<24 ,"block"]<-2 Exp_age[Exp_age$TESTAGE>=24 &Exp_age$TESTAGE<36 ,"block"]<-3 Exp_age[Exp_age$TESTAGE>=36,"block"]<-4 #Mean-center age Exp_age$AGE<-Exp_age$block-mean(Exp_age$block,na.rm=T) #Each individual needs to have one line for each block, even if not measured in each (NAs) blocks<-merge(Exp_age$Ring,unique(Exp_age$block)) names(blocks)<-c("Ring", "block") Exp_age<-merge(Exp_age,blocks,all=T) Exp_age<-Exp_age[order(Exp_age$block),]#Data needs to be sorted by block ##Run the models Exp_age$block<-as.factor(Exp_age$block) 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) ################################################## #Simple model with homegenous residuals and years ################################################## model.LOCOM00 <- asreml(LOCOM ~ 1 + pol(June.Days,2) + as.factor(AGE) + Sex , random = ~ped(Ring)+ Year + Observer , data = Exp_age,na.method.X="include", ginverse = list(Ring=ainv), maxiter=500) summary(model.LOCOM00) ######################################## #Model with heterogeneous residuals ######################################## model.LOCOM0a <- asreml(LOCOM ~ 1 + pol(June.Days,2) + Sex + as.factor(AGE), random = ~ped(Ring)+ Year + Observer , rcov= ~ at(block):units, data = Exp_age,na.method.X="include", ginverse = list(Ring=ainv), maxiter=500) summary(model.LOCOM0a) 1-pchisq(2*abs(model.LOCOM0a$loglik-model.LOCOM00$loglik),3)#Test for heterogeneous vs.homogeneous residuals ######################################## #Model with heterogeneous year variance ######################################## model.LOCOM0c <- asreml(LOCOM ~ 1 + pol(June.Days,2) + Sex + as.factor(AGE), random = ~ped(Ring)+ at(block):Year + Observer , rcov= ~ at(block):units, data = Exp_age,na.method.X="include", ginverse = list(Ring=ainv), maxiter=500) summary(model.LOCOM0c) 1-pchisq(2*abs(model.LOCOM0a$loglik-model.LOCOM0c$loglik),3)#Test for heterogeneous vs.homogeneous year variance ################################################################### #No GxA: Model with heterogeneous residual and year variances #VA writen as a function of age but x=0 ################################################################### model.LOCOM0b <- asreml(LOCOM ~ 1 + pol(June.Days,2)+as.factor(AGE) + Sex, random = ~us(pol(AGE,0)):ped(Ring)+ at(block):Year + Observer , rcov= ~ at(block):units, data = Exp_age,na.method.X="include", ginverse = list(Ring=ainv), maxiter=500) summary(model.LOCOM0b) #hist(model.LOCOM0b$residuals) #randomeffects<-model.LOCOM0b$coefficients$random #blups<-randomeffects[grep("ped",names(randomeffects))] #hist(blups) ################################################################### #GxA: Model where VA varies as a linear function of age (x=1) ################################################################### model.LOCOM1 <- asreml(LOCOM ~ 1 + pol(June.Days,2)+as.factor(AGE)+ Sex, random = ~us(pol(AGE,1)):ped(Ring)+ at(block):Year+ Observer , rcov= ~ at(block):units, data = Exp_age,na.method.X="include", ginverse = list(Ring=ainv), maxiter=100) summary(model.LOCOM1) 1-pchisq(2*abs(model.LOCOM0b$loglik-model.LOCOM1$loglik),2)#Test for GxA ################################################################### ##GxA: Model where VA varies as a quadratic function of age (x=2) ################################################################### #Need to specify starting values, constraints sv <- asreml(LOCOM ~ 1 + pol(June.Days,2) + Sex + as.factor(AGE), random = ~us(pol(AGE,2),init=c(6,3,5,-0.1,-1,0.3)):ped(Ring)+ at(block):Year+ Observer , rcov= ~ at(block):units, data = Exp_age,na.method.X="include", ginverse = list(Ring=ainv),workspace = 32e+6, pworkspace=8e+6, maxiter=50,start.values=T) gam<-sv$gammas.table gam$Value[c(7:15)]<-c(8,9.7,15.7,9.16,2.3,31,57,63,103) gam$Constraint[c(1:6)]<-"U" #Run the model model.LOCOM2<-asreml(LOCOM ~ 1 + pol(June.Days,2) + Sex + as.factor(AGE), random = ~us(pol(AGE,2)):ped(Ring)+ at(block):Year+ Observer , rcov= ~ at(block):units, data = Exp_age,na.method.X="include", ginverse = list(Ring=ainv),workspace = 32e+6, pworkspace=8e+6, maxiter=50,G.param=gam, R.param=gam) summary(model.LOCOM2)#Takes quite a long time 1-pchisq(2*abs(model.LOCOM1$loglik-model.LOCOM2$loglik),3)#Test for quadratic GxA