# This is R code for the worked example in: # Henshaw JM, Morrissey MB, Jones AG. Quantifying the causal pathways contributing to natural selection. Evolution. # It was written for R version 1.1.456. # Load the dataset Antechinus.csv before running the code library(mgcv) # Subset of surviving individuals AntechinusSurviving <- Antechinus[which(Antechinus$Survival==1),] # Subset of individuals that mated at least once AntechinusWithMates <- Antechinus[which(Antechinus$Mates>0),] # Generalized additive model (GAM) for survival survival.model <- gam(Survival ~ s(BodySize)+s(TorporDate),family=binomial,data=Antechinus) # GAM for mating success of surviving individuals mating.model <- gam(Mates ~ s(BodySize)+s(TorporDate),data=AntechinusSurviving,family=poisson) # GAM for fecundity of individuals with at least one mate fecund.model <- gam(Fecundity ~ s(BodySize)+s(TorporDate)+s(Mates),data=AntechinusWithMates,family=poisson) # This function calculates the expected fitness of an individual with trait values z (i.e. body size and # torpor date in this example). The parameter h is a 2x3 matrix. It perturbs the trait values for the purposes # of how they affect particlar fitness components (see the main text for more details of the 'twin network' # construction) Specifically, the first, second, and third rows of h perturb z with respect to its effect on # survival, mating success, and fecundity, respectively. expectedFitness<-function(z,h,max.mates=15){ # kludge to remove any dimNames which trip a warning later z<-as.numeric(z) # probability of survival P.s<-predict(survival.model,data.frame(BodySize=z[1]+h[1,1], TorporDate=z[2]+h[1,2]),type="response") # probability distribution for the number of matings, conditioning on survival (note that the density is for # number of mates >=1, because mates == 0 results in no fitness; this is reflected in subsequent calculations) E.m<-predict(mating.model,data.frame(BodySize=z[1]+h[2,1], TorporDate=z[2]+h[2,2]),type="response") D.m<-dpois(1:max.mates,E.m) # expected fecundity given number of mates, conditioning on survival E.f<-predict(fecund.model,data.frame(BodySize=z[1]+h[3,1], TorporDate=z[2]+h[3,2],Mates=1:max.mates), type="response") # expected fitness return(P.s*sum(D.m*E.f)) } # This function calculates the mean fitness across the population, perturbed by h as described above. meanFitness<-function(z,h=matrix(0,3,2)){ return(mean(apply(z,1,expectedFitness,h=h))) } # z<-Antechinus[,c("BodySize","TorporDate")] # store value of mean fitness Wbar<-meanFitness(z) # calculate total effects of body size and torpor date on fitness beta<-array(dim=2) h<-cbind(rep(0.01,3),0) beta[1]<-(meanFitness(z,h)-Wbar)/0.01/Wbar h<-cbind(0,rep(0.01,3)) beta[2]<-(meanFitness(z,h)-Wbar)/0.01/Wbar beta # calculate path-specific effects of body size and torpor date on fitness via survival, mating success and # fecundity (here using delta=0.01) beta.part<-array(dim=c(3,2)) h<-matrix(0,3,2); h[1,1]<-0.01 beta.part[1,1]<-(meanFitness(z,h)-Wbar)/0.01/Wbar h<-matrix(0,3,2); h[2,1]<-0.01 beta.part[2,1]<-(meanFitness(z,h)-Wbar)/0.01/Wbar h<-matrix(0,3,2); h[3,1]<-0.01 beta.part[3,1]<-(meanFitness(z,h)-Wbar)/0.01/Wbar h<-matrix(0,3,2); h[1,2]<-0.01 beta.part[1,2]<-(meanFitness(z,h)-Wbar)/0.01/Wbar h<-matrix(0,3,2); h[2,2]<-0.01 beta.part[2,2]<-(meanFitness(z,h)-Wbar)/0.01/Wbar h<-matrix(0,3,2); h[3,2]<-0.01 beta.part[3,2]<-(meanFitness(z,h)-Wbar)/0.01/Wbar # matrix of path-specific effect (columns represent traits body size and torpor date, rows represent # path-specific effects via survival, mating success, fecundity) beta.part # linear regression selection gradients (after Lande & Arnold 1983) lm(Fecundity/Wbar ~ BodySize + TorporDate, data=Antechinus)