library(mvtnorm) evocolor <- c("blue","darkgreen","red","grey90") # setup: create a list of model parameter values # # Arguments: See below # # Value: # A list of model parameters # setup <- function( N=500, z0=c(0,0,0), # breeding values, initial noise reaction norm positive slope Gaa=.1, # G-matrix Gbb=.01, Gcc=.01, Gab=0,Gac=0,Gbc=0, # genetic covariances alpha=.95, # correlation at lag 1 in the O-U process for the environment sigma_env=sqrt(2), # standard deviation of macroenvironmental fluctuations sigma_u=sigma_env, # standard deviation of microenvironmental deviations rho_uv=0, # correlation between microenvironmental deviations at dev. and sel. sigma_e=0, # non-evolving component of the environmental component of the phen. variance tau=0.5, # lag between development and selection omega=1, # width of individual phenotypic Gaussian fitness function omega_b=Inf, # width of Gaussian selection function on plasticity A=0, # optimal elevation B=1 # and slope ) { as.list(environment(), all=TRUE) } # Simulate the process # # Arguments: # par: List of model parameter values (use setup to make this) # tmax: Number of generations to simulate # type: If "recursions" then simulate by iterating eq. (3d) combined with Lande's equation # env.approx: If TRUE, simulate the process replacing G_bb epsilon_{t-tau}^2 by G_bb sigma_epsilon^2, # if FALSE, simulate the process using the phenotypic variance given by in eq. (3b). # # Value: A list with class "mismatch" containing the following components # data: A dataframe containing the following columns # generation: Generation number 1 to tmax # abar: bar a_t # bbar: bar b_t # cbar: bar c_t # varz: var(z-theta) given by eq. (3b) # env: epsilon_t # env.tau: epsilon_{t-tau} # wbar: Mean fitness given by eq. (3c) # par: A list of model parameter values used. # sim <- function( par, tmax=2000, # number of generations to simulate type="recursions", env.approx=FALSE ) { attach(par) abar <- bbar <- cbar <- varz <- wbar <- rep(NA,tmax) # mean additive genetic values among newborns in generation t G <- matrix(c(Gaa,Gab,Gac, Gab,Gbb,Gbc, Gac,Gbc,Gcc),3,3) env <- rep(NA,tmax) # env at time of selection in generation t env.tau <- rep(NA,tmax) # env in generation t at time of development tau before selection for (t in 1:tmax) { env.tau[t] <- alpha^(1-tau)*(if (t==1) rnorm(1,sd=sigma_env) else env[t-1]) + rnorm(1,0,sqrt(1-alpha^(2*(1-tau)))*sigma_env) env[t] <- alpha^tau*env.tau[t] + rnorm(1,0,sqrt(1-alpha^(2*tau))*sigma_env) } switch (type, recursions={ abar[1] <- z0[1] bbar[1] <- z0[2] cbar[1] <- z0[3] for (t in 1:tmax) { if (env.approx) env.tau.2 <- sigma_env^2 else env.tau.2 <- env.tau[t]^2 varz[t] <- Gaa + Gbb*env.tau.2 + (Gbb + bbar[t]^2)*sigma_u^2 + Gcc + cbar[t]^2 + sigma_e^2 + (B^2 - 2*bbar[t]*B*rho_uv)*sigma_u^2 denominator <- omega^2+varz[t] mismatch <- abar[t] + bbar[t]*env.tau[t] - A - B*env[t] commonfactor <- 1 - mismatch^2/denominator wbar[t] <- sqrt(omega^2/denominator)*exp(-mismatch^2/(2*denominator)) if (t==tmax) break() gradient <- -c(mismatch, # from phenotypic selection mismatch*env.tau[t] + sigma_u^2*(bbar[t]-rho_uv*B)*commonfactor, cbar[t]*commonfactor )/denominator - c(0, # from cost of plasticity bbar[t], 0)/(omega_b^2 + Gbb) response <- G%*%gradient abar[t+1] <- abar[t] + response[1] bbar[t+1] <- bbar[t] + response[2] cbar[t+1] <- cbar[t] + response[3] } } ) tmp <- list(data=data.frame(generation=1:tmax, abar=abar, bbar=bbar, cbar=cbar, varz=varz, env=env, env.tau=env.tau, wbar=wbar),par=par) class(tmp) <- "mismatch" detach(par) tmp } weave <- function(x,y) { as.vector(rbind(x,y)) } # validation plot generates Fig 1 and similar plots in Appendix S4. # # Arguments: # param: A list of model parameter values # tshort: Number of generations in middle column # tmax: Number of generation in left column # tlong: Number of generation in right column # scale: Use default value... # seed: seed of the random number generator # validationplot <- function(param,tshort=100,tmax=10000,tlong=1e+5,scale=c(2.5,1.5,2.5),seed=3,...) { longrealisation <- sim(param,tmax=tlong,...) set.seed(seed) realisation <- sim(param,tmax=tmax,...) layout(matrix(1:9,3,3),widths=c(5.2,4,4),heights=c(5,3,4.2)) par(oma=rep(0,4)) plot(realisation,subset=seq(1,tmax,by=round(tmax/1e+3)),legend=FALSE,scale=scale) plot(realisation,subset=(tmax-tshort):tmax,yannot=FALSE,scale=scale) plot(longrealisation,subset=-(1:round(tlong/10)),binormal=TRUE,yannot=FALSE,scale=scale) } plot.mismatch <- function(obj,subset=TRUE,binormal=FALSE,legend=TRUE,yannot=TRUE,mar=.6,scale=c(3,2,2),...) { shade <- function(x,y1,y2) polygon(c(x,x[length(x):1]),c(y1,y2[length(y2):1]),col=evocolor[4],border=NA) sol <- solution(obj$par) with(obj$par,{ with(obj$data[subset,],{ # elevation, optimum etc... par(mar=c(mar,if (yannot) 4.5 else mar,mar,mar)+.2) xlim <- if (binormal) sigma_env*3.2*c(-1,1) else range(generation) x <- env # env.tau changed to env xx <- seq(xlim[1],xlim[2],length=100) thin <- max(1,round(length(abar)/1e+3)) thin <- c(TRUE,rep(FALSE,thin-1)) plot(NA,xlim=xlim,ylim=c(-1,1)*sigma_env*scale[1], axes=FALSE,xlab="",ylab=if (yannot) expression(bar(a[t])) else "",...) box() axis(1,labels=FALSE) axis(2,labels=yannot) if (binormal) { y <- abar #points(x[thin],abar[thin]+bbar[thin]*env.tau[thin],pch=".",col=evocolor[2],cex=.5) points(x[thin],y[thin],pch=".",col=evocolor[1],cex=.5) abline(v=qnorm(c(.025,.975),sd=sigma_env),lty=3) abline(a=0,b=1,lty=3) abline(lm(y ~ -1 + x),lty=2) fit <- lm(y ~ -1 + X1 + X2 + X3 + X4, data=data.frame(y=y,outer(x,c(1,3,5,7),"^"))) lines(xx,predict(fit,newdata=data.frame(outer(xx,c(1,3,5,7),"^"))),col=evocolor[1]) abline(a=0,b=sol["cov.a.env"]/sigma_env^2) } else { shift <- 0 # tau+(1-tau)/2 # shift env such that things align in cronological order shade(generation, abar+bbar*env+1.96*sqrt(varz), abar+bbar*env-1.96*sqrt(varz)) lines(weave(generation-tau,generation)+shift,weave(env.tau,env),col="black",type="l",...) #points(generation-tau+shift,env.tau,col="black",cex=.3,...) lines(generation,abar+bbar*env.tau,col=evocolor[2],...) lines(generation,abar,col=evocolor[1],...) if (legend) legend("topleft",c( expression(A + B*epsilon[t]), expression(bar(a[t])), expression(bar(a[t])+bar(b[t])*epsilon[t-tau]), "95%-int. for z"), lwd=c(1,1,1,9), col=c("black",evocolor[1],evocolor[2],"grey90"),lty=c(1,1,1,1), bty="n",bg="white") } # slope ylim <- scale[2]*sol["bbb"] plot(NA,xlim=xlim,ylim=c(0,ylim), axes=FALSE,xlab="",ylab=if (yannot) expression(bar(b[t])) else "",...) box() axis(1,labels=FALSE) axis(2,labels=yannot) if (binormal) { y <- bbar points(x[thin],y[thin],pch=".",col=evocolor[2],cex=.5) abline(v=qnorm(c(.025,.975),sd=sigma_env),lty=3) abline(lm(y ~ 1),lty=2) fit <- lm(y ~ X1 + X2 + X3 + X4,data=data.frame(y,outer(x,c(2,4,6,8),"^"))) lines(xx,predict(fit,newdata=data.frame(outer(xx,c(2,4,6,8),"^"))),col=evocolor[2]) abline(h=sol["bbb"]) } else { lines(generation,bbar,col=evocolor[2],type="l", ylim=c(0,sol["bbb"]*scale[2]), xlab="",ylab=if (yannot) expression(bar(b[t])) else "",...) } # cbar ylim <- scale[3]*max(.1,sol["bbc2"]) par(mar=c(5,if (yannot) 4.5 else mar,mar,mar)+.2) plot(NA,xlim=xlim,ylim=c(0,ylim), axes=FALSE, xlab=if (binormal) expression(paste("Environment at selection ",epsilon[t])) else "Generation t", ylab=if (yannot) expression(bar(c[t]^2)) else "", ...) box() axis(1,labels=TRUE) axis(2,labels=yannot) if (binormal) { y <- cbar^2 points(x[thin],y[thin],pch=".",col=evocolor[3],cex=.5) abline(v=qnorm(c(.025,.975),sd=sigma_env),lty=3) abline(lm(y ~ 1),lty=2) fit <- lm(y ~ X1 + X2 + X3 + X4,data=data.frame(y,outer(x,c(2,4,6,8),"^"))) lines(xx,predict(fit,newdata=data.frame(outer(xx,c(2,4,6,8),"^"))),col=evocolor[3]) abline(h=sol["bbc2"]) } else { lines(generation,varz,col="black") lines(generation,cbar^2,col=evocolor[3]) } }) }) } # solution: Computes the solution for the joint evolutionary outcome # # Arguments: # # par: a list of parameter values (use setup to create this) # type: If type="complete" then all quantitities characterizing the joint evolutionary outcome # is computed. If type="fig1" and type="fig4", then compute only what's needed for Figs. 2-4 # and Fig 5, respectively. # approximation: If "none", "bethedge" or "smallmismatch", compute the solution based on eq. # (6a), (6b) or (6c) respectively. # constrained: Experimental (use the default) # algorithm: If "fixpoint" solve the equations using fixpoint iterations (works well), if "optim" # then solve the system # f_1(x_1,...,x_n) = 0 # . # . # f_n(x_1,...,x_n) = 0 # by minimizing sum_{i=1}^n f_i(x_1,...,x_n)^2 using the optim R function (appear to be # numerically unstable in some parts of the parameter space) # verbose: If TRUE print number of iterations required for the fixpoint iteration algorithm to # converge # tol: Convergence tolerance for the fixpoint iteration method. # # Value: # # A vector containing quantities of interest at the joint evolutionary equilibrium: # If type="complete" the elements are # # cov.a.env: The covariance between mean reaction norm elevation and the environment at # time of selection, eq. (4c) # var.a: The variance in mean reaction norm elevation, eq. (S1.3) # bbb: Depending on the approximation argument (see above), mean reaction norm slope # based on eqs. (6a), (6b) or (6c) # bbc2: The squared value of barbar c, based on eq. (5c) # varmismatch: The mismatch variance, eq. (S1.4) # meanlogwbar: The expected value of the log of bar w given in eq. (3c) # with respect to the joint distribution of env_t, env_{t-tau} and bar a_t. # exp(meanlogwbar) is thus the long term geometric mean fitness at the joint # evolutionary equilibrium (solid curves in Figs. S4.3C, and S4.6) # solution <- function(par,type="complete",approximation="none",constrained=TRUE,algorithm="fixpoint",verbose=FALSE,tol=1e-12) { attach(par) unknown <- c(Gaa,Gbb,Gcc)!=0 # parameters for which we simulate no evolution are assumed to given if (!unknown[1]) {# by their initial values used in the simulations... # no evolution in elevation cov.a.env <- 0 var.a <- 0 } if (!unknown[2]) bbb <- z0[2] if (!unknown[3]) bbc2 <- z0[3]^2 # sum of squares of all equations newiteration <- function(x,additional=FALSE) { # first extract current values of unknown quantitities from first argument of obj if (unknown[1]) { cov.a.env <- x[1]; x <- x[-1] var.a <- x[1]; x <- x[-1] } if (unknown[2]) { bbb <- x[1]; x <- x[-1] } if (unknown[3]) { bbc2 <- x[1]; x <- x[-1] } # intermediate quantities varz <- Gaa + Gbb*sigma_env^2 + (Gbb+bbb^2)*sigma_u^2 + Gcc + bbc2 + sigma_e^2 + B^2*sigma_u^2 - 2*bbb*B*rho_uv*sigma_u^2 s <- Gaa/(omega^2 + varz) varmismatch <- var.a + (B^2+bbb^2-2*B*bbb*alpha^tau)*sigma_env^2 - 2*(B-bbb*alpha^-tau)*cov.a.env common <- 1 - varmismatch/(omega^2 + varz) if (additional) { meanlogwbar <- .5*log(omega^2/(omega^2 + varz)) - varmismatch/(2*(omega^2+varz)) return(c(varmismatch=unname(varmismatch),meanlogwbar=unname(meanlogwbar))) } else { # create a vector containing current values of functions representing each equation newit <- c( if (unknown[1]) # equations expressing cov(a,env) and var(a) as functions of bbb and bbc c(s*alpha*(B-bbb*alpha^tau)*sigma_env^2/(1-alpha*(1-s)), s*sigma_env^2/(2-s)*(B^2 + bbb^2 - 2*bbb*B*alpha^tau + 2*alpha*(1-s)*(B-bbb*alpha^-tau)*(B-bbb*alpha^tau)/(1-alpha*(1-s))) ), if (unknown[2]) { if (approximation=="none") cm <- common if (approximation=="bethedge") cm <- 0 if (approximation=="smallmismatch") cm <- 1 (B*(alpha^tau*sigma_env^2 + rho_uv*sigma_u^2*cm) - alpha^-tau*cov.a.env)/ (sigma_env^2 + sigma_u^2*cm + (omega^2 + varz)/(omega_b^2 + Gbb)) }, if (unknown[3]) # bbc2 as function of bbb, var(a) and cov(a,env) max(varmismatch - omega^2 - (varz - bbc2),0) ) # return the sum of squares of these functions return(newit) } } obj <- function(x) sum((newiteration(x)-x)^2) # starting values start <- c( if (unknown[1]) c(cov.a.env=.5*sigma_env^2,var.a=.5*sigma_env^2), if (unknown[2]) c(bbb=0.5), if (unknown[3]) c(bbc2=0) ) switch(algorithm, optim={ # the minimum of obj should be the solution ctrl <- list(abstol=1e-10,maxit=500) tmp <- optim(start,obj,method="BFGS",control=ctrl) # first attempt to optimize all three parameters if (tmp$par["bbc2"]<0 & constrained) { # check if solution for environmental variance is negative unknown[3] <- FALSE # if so, treat bbc2 as known #warning("bbc^2 negative, assuming bbc^2=zero") bbc2 <- 0 # and equal to zero start <- start[-length(start)] # remove bbc2 from vector of initial values tmp <- optim(start,obj,method="BFGS",control=ctrl) tmp$par <- c(tmp$par,bbc2=0) } sol <- tmp$par }, fixpoint={ i <- 0 x <- start repeat { prev.x <- x x <- newiteration(x) i <- i + 1 if (abs(sum((x-prev.x)^2))<=tol) break() } if (verbose) cat("Convergence after ",i," iterations.\n") sol <- x } ) additional.quantities <- newiteration(sol,add=TRUE) ssol <- numeric(4) unknown4 <- c(unknown[1],unknown) ssol[unknown4] <- sol ssol[!unknown4] <- 0 names(ssol) <- c("cov.a.env","var.a","bbb","bbc2") tmp <- switch(type, fig1={ c(ssol["cov.a.env"]/sigma_env^2, ssol["bbb"]*alpha^tau, #... ssol["bbc2"]) }, fig4={ c(ssol["cov.a.env"]/sigma_env^2, ssol["bbb"], ssol["bbc2"]) }, complete=c(ssol,additional.quantities)) detach(par) tmp } # Used to make Fig 2 to 5 in paper intersection2d <- function(par, xvar,yvar, xrange,yrange, xtrans=none,ytrans=none, xlab=xvar,ylab=yvar, n=100,nx=n,ny=3, add=TRUE, algorithm="fixpoint", type="fig1", method="edge", evolevels=c(NA,NA,NA), approximations=TRUE, clip=TRUE, tot=FALSE) { k <- length(solution(par,type,algorithm=algorithm)) z <- array(dim=c(nx,ny,k)) if (approximations) z.bh <- z.sm <- z x <- seq(xtrans$fn(xrange[1]),xtrans$fn(xrange[2]),length=nx) # x and y-values on transformed scale y <- seq(ytrans$fn(yrange[1]),ytrans$fn(yrange[2]),length=ny) for (i in 1:nx) for (j in 1:ny) { par[[xvar]] <- xtrans$inv(x[i]) # backtransform x and y and put into par par[[yvar]] <- ytrans$inv(y[j]) z[i,j,] <- solution(par,type=type,algorithm=algorithm) # compute the solution at this point if (approximations) { z.bh[i,j,] <- solution(par,type=type,algorithm=algorithm,approx="bethedge") # compute the solution at this point z.sm[i,j,] <- solution(par,type=type,algorithm=algorithm,approx="smallmismatch") # compute the solution at this point } } layout(matrix(1:3,ny,1),widths=c(6),heights=c(2.3,2.3,3.3)) scale <- c(1,1,2) for (j in ny:1) { par(mar=c(.5, 3.8, .5, 3.8) + 3.8*c((j==1),0,0,0)) plot(NA,xlim=range(x),ylim=c(0,1),xaxt="n",yaxt="n",xlab="",ylab="") axis(1,at=xtrans$fn(xtrans$at),labels=if (j==1) xtrans$labels else FALSE) if (j==1) { axis(1,at=xtrans$fn(xtrans$at),labels=xtrans$labels2,tick=FALSE,line=.8) mtext(xlab,side=1,line=3.5,cex=par("cex")) } axis(4,at=0:10/10,labels=TRUE) cticks <- seq(0,scale[3],by=.5) axis(2,at=cticks/scale[3],labels=cticks) if (j==2) { mtext(expression(paste(cov(bar(a)[t],epsilon[t])/sigma[epsilon]^2*B," ")), side = 4, line = 2.4, cex = par("cex"),col=evocolor[1]) mtext( if (type=="fig1") expression(paste(" ",bar(bar(b))*alpha^tau/B)) else expression(paste(" ",bar(bar(b))/B)) ,side = 4, line = 2.4, cex = par("cex"),col=evocolor[2]) mtext(expression(bar(bar(c))^2),side=2,line=2.2,cex=par("cex"),col=evocolor[3]) } for (i in 1:k) { lines(x,z[,j,i]/scale[i],col=evocolor[i]) } if (approximations) { if (clip) { subset1 <- 1:(nx/2) subset2 <- z[,j,3]<1e-3 } else subset1 <- subset2 <- TRUE lines(x[subset1],z.bh[subset1,j,2]/scale[2],col=evocolor[2],lty=2) lines(x[subset2],z.sm[subset2,j,2]/scale[2],col=evocolor[2],lty=3) } if (type!="fig4") lines(x,z[,j,1]+z[,j,2],col="lightgrey") legend("topleft",legend=substitute(G[aa]==x,list(x=ytrans$inv(y[j]))),bty="n") if (tot) { lines(x,xtrans$inv(x)^(2*.5),lty=2) } } list(z=z,x=xtrans$inv(x),y=ytrans$inv(y)) } # Generates plots similar to Figs. 2 to 5 in the paper but as 3D contour plot ala Botero et. al (2015) in PNAS. intersection <- function(par, xvar,yvar, xrange,yrange, xtrans=none,ytrans=none, xlab=xvar,ylab=yvar, n=100,nx=n,ny=n, add=TRUE, algorithm="optim", type="fig1", method="edge", evolevels=c(NA,NA,NA)) { k <- length(solution(par,type,algorithm=algorithm)) z <- array(dim=c(nx,ny,k)) x <- seq(xtrans$fn(xrange[1]),xtrans$fn(xrange[2]),length=nx) # x and y-values on transformed scale y <- seq(ytrans$fn(yrange[1]),ytrans$fn(yrange[2]),length=ny) for (i in 1:nx) for (j in 1:ny) { par[[xvar]] <- xtrans$inv(x[i]) # backtransform x and y and put into par par[[yvar]] <- ytrans$inv(y[j]) z[i,j,] <- solution(par,type=type,algorithm=algorithm) # compute the solution at this point } plot(NA,xlim=range(x),ylim=range(y),xaxt="n",yaxt="n",xlab=xlab,ylab=ylab) # empty plot of the right size axis(1,at=xtrans$fn(xtrans$at),labels=xtrans$labels) if (!is.null(xtrans$at2)) axis(3,at=xtrans$at2,labels=xtrans$labels2) axis(2,at=ytrans$fn(ytrans$at),labels=ytrans$labels) for (i in 1:k) { if (is.na(evolevels[[i]])) { levels <- if (i==3) c(.02,pretty(range(z[,,i]),10)) else -1:11/10 } else levels <- evolevels[[i]] contour(x,y,z[,,i],add=add,col=evocolor[i],levels=levels,method=method[(i-1)%%length(method)+1]) } z } none <- list( fn=function(x) x,inv=function(x) x,ticks=NULL ) fisher <- list( fn=function(rho) log((1+rho)/(1-rho)), inv=function(x) (exp(x)-1)/(exp(x)+1), at=sort(c(0,c(1:9/10,.9+1:9/100,.99+1:9/1000) -> tmp,-tmp)), labels=sort(c(0,c(1:9/10,.9+1:9/100,.99+1:9/1000) -> tmp,-tmp)) ) loglabs <- c(1,2,NA,NA,5,NA,NA,NA,NA) lablogs <- c(NA,NA,NA,NA,5,NA,NA,NA,9) logtime <- list( fn=function(corr) -log(-log(corr)/.5), inv=function(y) exp(-exp(-y)), at=c(1:9/1000,1:9/100,1:9/10,.9+1:9/100,.99+1:9/1000), labels=c(rep(NA,9),loglabs/100,1:9/10,.9+lablogs/100,.99+lablogs/1000), at2=log(as.vector(outer(1:9,10^(-5:5),"*"))), labels2=as.vector(outer(loglabs,10^(-5:5),"*")) ) corr <- exp(-.5/as.vector(outer(loglabs,10^(-5:5),"*"))) l2 <- rep(NA,length(corr)) for (i in 1:length(corr)) if (!is.na(corr[i])) l2[i] <- prettyNum(corr[i],digits=if (corr[i]>.91) 3 else 2) l2 <- paste("(",l2,")",sep="") l2[l2=="(NA)"] <- NA logtime2 <- list( fn=function(alpha) log(-1/log(alpha)), inv=function(y) exp(-exp(-y)), at=exp(-1/as.vector(outer(1:9,10^(-5:5),"*"))), labels=as.vector(outer(loglabs,10^(-5:5),"*")), labels2=l2 ) log <- list( fn=function(x) log(x), inv=function(x) exp(x), at=as.vector(outer(1:9,10^(-5:5),"*")), labels=as.vector(outer(loglabs,10^(-5:5),"*")) ) logdevsel <- list( fn=function(x) log(x), inv=function(x) exp(x), ticks=sort(as.vector(outer(1:9,10^(-5:5),"*"))) )