Libaries

Workflow

This simulation estimates broad sense heritability and assumes the experimental design where clones of the same species are repeatedly and independently sampled. When these clones are measured across an environmental gradient (or any gradient), then the genetic associations of the population of clones can be expressed as the G matrix, or variance-covariance matrix.

Simulate 3 possible outcomes for how performance curves (modeled as a gaussian function) can vary:

Veritical –traits are all correlated across environmental gradient
Horizontal – optimum varies along environmental gradient; so genotypes/clones have different optimum
Generalist-specialist – those with higher max performance, also operate over a smaller window along the environmental gradient

For each possibility:

mermaid("
graph TD
  a{Simulate Data} --20 clones or colonies, 5 reps--> b(Mixed Effect Model) 
  b --> c[Extract G] 
  c --> d[PCA decomposition] 
  d --> e[Visualize]
  a --> e
")

Constructing gaussian shaped performance curves

\[Performance = height * e^{(-0.5(\frac{Temperature - T_{opt}}{width})^2)}\]

gau<-function(a=1,b=30,c=10,t=seq(1,60,1)){
  #z=a*exp(-.5((t-b)/c)^2)
  z=a*exp(-.5*((t-b)/c)^2)
  b<-data.frame(performance=z,t=t)
  return(b)  
}

#gau()
plot(gau()$t,gau()$performance,type="l")

Simulate data with just vertical variation/responses

Simulating the data

#simulate vert values
vval<-rnorm(mean=50,sd=15,n=100)

#pull out 20, randomly

col<-sort(rep(sample(vval,20),5))


#estimate perf curves
random.error<-sample(rnorm(mean=0,sd=.5,n=1000),500)

#generate dataset with 5 reps for 20 colonies for 4 levels of temperature
dat<-rbind(gau(a=col,t=c(20,25,30,35,40)),gau(a=col,t=c(20,25,30,35,40)),gau(a=col,t=c(20,25,30,35,40)),gau(a=col,t=c(20,25,30,35,40)),gau(a=col,t=c(20,25,30,35,40)))
dat$performance<-dat$performance+random.error


#convert colonies and dat into factors
dat$colonies<-as.factor(rep(sort(rep(seq(1,20,1),5)),5))
dat$t<-as.factor(dat$t)

Mixed effects model to estimate G

### building mixed effects model to estimate the variance-covariance structure
# random slope and temperature nested within colonies
vmod<-lmer(formula=performance~1+(0+t|colonies),REML=TRUE,data=dat)

## boundary (singular) fit: see help('isSingular')

## Warning: Model failed to converge with 1 negative eigenvalue: -1.4e-02

#summary(vmod)
Gmatrix.pac<-VarCorr(vmod)$colonies;Gmatrix.pac

##           t20      t25      t30      t35       t40
## t20  849.6332 1225.882 1383.437 1221.864  850.4995
## t25 1225.8815 1768.802 1996.138 1762.978 1227.1459
## t30 1383.4371 1996.138 2252.692 1989.564 1384.8642
## t35 1221.8644 1762.978 1989.564 1757.192 1223.1290
## t40  850.4995 1227.146 1384.864 1223.129  851.3944
## attr(,"stddev")
##      t20      t25      t30      t35      t40 
## 29.14847 42.05713 47.46253 41.91888 29.17866 
## attr(,"correlation")
##           t20       t25       t30       t35       t40
## t20 1.0000000 0.9999842 0.9999833 0.9999945 0.9999837
## t25 0.9999842 1.0000000 1.0000000 0.9999944 0.9999797
## t30 0.9999833 1.0000000 1.0000000 0.9999939 0.9999790
## t35 0.9999945 0.9999944 0.9999939 1.0000000 0.9999936
## t40 0.9999837 0.9999797 0.9999790 0.9999936 1.0000000

Gmatrix.pac[1:5,1:5]

##           t20      t25      t30      t35       t40
## t20  849.6332 1225.882 1383.437 1221.864  850.4995
## t25 1225.8815 1768.802 1996.138 1762.978 1227.1459
## t30 1383.4371 1996.138 2252.692 1989.564 1384.8642
## t35 1221.8644 1762.978 1989.564 1757.192 1223.1290
## t40  850.4995 1227.146 1384.864 1223.129  851.3944

PCA decomposition of G

h<-princomp(Gmatrix.pac[1:5,1:5])
h$loadings[,1:3]

##        Comp.1     Comp.2      Comp.3
## t20 0.3370221  0.2165234  0.84627777
## t25 0.4863022 -0.3928231 -0.12488735
## t30 0.5488050 -0.4713360 -0.17851244
## t35 0.4846879  0.2882604  0.07348448
## t40 0.3373678  0.7025354 -0.48057302

h1<-data.frame(Loadings=h$loadings[,1],t=c(20,25,30,35,40))

aa<-ggplot(h1,aes(x=t,y=Loadings))+geom_line(size=3)+setup+geom_hline(yintercept=0)+ylim(-1,1)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

aa

y axis = loadings
x axis = environmental gradient

Figure

dat$t<-as.numeric(as.character(dat$t))
a<-ggplot(dat,aes(x=t,y=performance,colour=factor(colonies)))+geom_smooth(se=FALSE)+setup+scale_color_hue(l=20, c=50)+ylim(5,100)
a

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Simulate horizontal variation

#simulate Topt values
opt<-rnorm(mean=30,sd=5,n=100)

#pull out 20, randomly
col.opt<-sort(rep(sample(opt,20),5))
col.opt

##   [1] 22.97488 22.97488 22.97488 22.97488 22.97488 24.31647 24.31647 24.31647
##   [9] 24.31647 24.31647 24.33271 24.33271 24.33271 24.33271 24.33271 24.59276
##  [17] 24.59276 24.59276 24.59276 24.59276 26.16661 26.16661 26.16661 26.16661
##  [25] 26.16661 26.71306 26.71306 26.71306 26.71306 26.71306 28.20692 28.20692
##  [33] 28.20692 28.20692 28.20692 28.71378 28.71378 28.71378 28.71378 28.71378
##  [41] 28.84939 28.84939 28.84939 28.84939 28.84939 29.23574 29.23574 29.23574
##  [49] 29.23574 29.23574 29.34596 29.34596 29.34596 29.34596 29.34596 30.09732
##  [57] 30.09732 30.09732 30.09732 30.09732 30.46406 30.46406 30.46406 30.46406
##  [65] 30.46406 30.86803 30.86803 30.86803 30.86803 30.86803 31.12596 31.12596
##  [73] 31.12596 31.12596 31.12596 31.20961 31.20961 31.20961 31.20961 31.20961
##  [81] 31.62681 31.62681 31.62681 31.62681 31.62681 31.67951 31.67951 31.67951
##  [89] 31.67951 31.67951 37.98338 37.98338 37.98338 37.98338 37.98338 38.91675
##  [97] 38.91675 38.91675 38.91675 38.91675

#estimate perf curves
#random.error<-sample(rnorm(mean=0,sd=.5,n=1000),500)

#generate dataset with 5 reps for 20 colonies for 4 levels of temperature
dat2<-rbind(gau(b=col.opt,t=c(20,25,30,35,40)),gau(b=col.opt,t=c(20,25,30,35,40)),gau(b=col.opt,t=c(20,25,30,35,40)),gau(b=col.opt,t=c(20,25,30,35,40)),gau(b=col.opt,t=c(20,25,30,35,40)))
dat2$performance<-dat2$performance#+random.error

###
dat2$colonies<-as.factor(rep(sort(rep(seq(1,20,1),5)),5))
dat2$t<-as.factor(dat2$t)

Mixed effects model to Estimate G

### building mixed effects model to estimate the variance-covariance structure
# random slope and temperature nested within colonies
hmod<-lmer(formula=performance~1+(0+t|colonies),REML=TRUE,data=dat2)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 4 negative eigenvalues

## Warning: Model failed to converge with 2 negative eigenvalues: -2.0e-03
## -6.7e-03

## Warning in as_lmerModLT(model, devfun): Model may not have converged with 4
## eigenvalues close to zero: 5.8e-09 2.1e-09 1.1e-09 6.0e-10

#summary(hmod)
Gmatrix.pac<-VarCorr(hmod)$colonies;Gmatrix.pac

##              t20           t25          t30           t35          t40
## t20  0.005010638  0.0019921445 -0.001217548 -0.0030372003 -0.002433988
## t25  0.001992144  0.0047560858  0.004112264  0.0002479976 -0.004552885
## t30 -0.001217548  0.0041122641  0.005950242  0.0030437430 -0.002780201
## t35 -0.003037200  0.0002479976  0.003043743  0.0035567441  0.001677576
## t40 -0.002433988 -0.0045528854 -0.002780201  0.0016775763  0.006435812
## attr(,"stddev")
##        t20        t25        t30        t35        t40 
## 0.07078586 0.06896438 0.07713781 0.05963844 0.08022352 
## attr(,"correlation")
##            t20         t25        t30         t35        t40
## t20  1.0000000  0.40808391 -0.2229833 -0.71944996 -0.4286178
## t25  0.4080839  1.00000000  0.7730166  0.06029708 -0.8229249
## t30 -0.2229833  0.77301664  1.0000000  0.66162873 -0.4492698
## t35 -0.7194500  0.06029708  0.6616287  1.00000000  0.3506342
## t40 -0.4286178 -0.82292487 -0.4492698  0.35063422  1.0000000

Gmatrix.pac[1:5,1:5]

##              t20           t25          t30           t35          t40
## t20  0.005010638  0.0019921445 -0.001217548 -0.0030372003 -0.002433988
## t25  0.001992144  0.0047560858  0.004112264  0.0002479976 -0.004552885
## t30 -0.001217548  0.0041122641  0.005950242  0.0030437430 -0.002780201
## t35 -0.003037200  0.0002479976  0.003043743  0.0035567441  0.001677576
## t40 -0.002433988 -0.0045528854 -0.002780201  0.0016775763  0.006435812

PCA decomposition of G

h<-princomp(Gmatrix.pac[1:5,1:5])
h$loadings[,1:3]

##         Comp.1      Comp.2    Comp.3
## t20  0.3030303  0.58999029 0.6604798
## t25  0.5556301 -0.10947631 0.2373104
## t30  0.3833898 -0.57731129 0.2300614
## t35 -0.1158491 -0.55051422 0.3573214
## t40 -0.6626003 -0.05976758 0.5717022

h2<-data.frame(Loadings=h$loadings[,1],t=c(20,25,30,35,40))

bb<-ggplot(h2,aes(x=t,y=Loadings))+geom_line(size=3)+setup+geom_hline(yintercept=0)+ylim(-1,1)
bb

y axis = loadings
x axis = environmental gradient

Figure

dat2$t<-as.numeric(as.character(dat2$t))
b<-ggplot(dat2,aes(x=t,y=performance,colour=factor(colonies)))+stat_smooth(se=FALSE)+setup+scale_color_hue(l=20, c=50)+xlim(20,40)+ylim(0.5,1.5)
b

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Simulate generalist-specialist variation

making the data

ones that have small widths, have higher performance

col# height

##   [1] 27.64292 27.64292 27.64292 27.64292 27.64292 32.08113 32.08113 32.08113
##   [9] 32.08113 32.08113 33.50570 33.50570 33.50570 33.50570 33.50570 36.13308
##  [17] 36.13308 36.13308 36.13308 36.13308 39.03085 39.03085 39.03085 39.03085
##  [25] 39.03085 40.89886 40.89886 40.89886 40.89886 40.89886 41.29205 41.29205
##  [33] 41.29205 41.29205 41.29205 41.60609 41.60609 41.60609 41.60609 41.60609
##  [41] 42.50445 42.50445 42.50445 42.50445 42.50445 44.26300 44.26300 44.26300
##  [49] 44.26300 44.26300 44.61538 44.61538 44.61538 44.61538 44.61538 47.09645
##  [57] 47.09645 47.09645 47.09645 47.09645 48.00077 48.00077 48.00077 48.00077
##  [65] 48.00077 49.99801 49.99801 49.99801 49.99801 49.99801 51.43731 51.43731
##  [73] 51.43731 51.43731 51.43731 51.77446 51.77446 51.77446 51.77446 51.77446
##  [81] 56.09633 56.09633 56.09633 56.09633 56.09633 59.95329 59.95329 59.95329
##  [89] 59.95329 59.95329 59.96425 59.96425 59.96425 59.96425 59.96425 65.21479
##  [97] 65.21479 65.21479 65.21479 65.21479

w<-rnorm(mean=5,sd=3,n=100)
col.opt<-sort(rep(sample(w,20),5))
col.opt # width

##   [1] -0.939851 -0.939851 -0.939851 -0.939851 -0.939851  1.929251  1.929251
##   [8]  1.929251  1.929251  1.929251  2.094721  2.094721  2.094721  2.094721
##  [15]  2.094721  2.254068  2.254068  2.254068  2.254068  2.254068  2.462421
##  [22]  2.462421  2.462421  2.462421  2.462421  2.605535  2.605535  2.605535
##  [29]  2.605535  2.605535  2.939121  2.939121  2.939121  2.939121  2.939121
##  [36]  2.953709  2.953709  2.953709  2.953709  2.953709  3.561728  3.561728
##  [43]  3.561728  3.561728  3.561728  4.326515  4.326515  4.326515  4.326515
##  [50]  4.326515  6.129862  6.129862  6.129862  6.129862  6.129862  6.182413
##  [57]  6.182413  6.182413  6.182413  6.182413  6.261775  6.261775  6.261775
##  [64]  6.261775  6.261775  6.460487  6.460487  6.460487  6.460487  6.460487
##  [71]  6.485163  6.485163  6.485163  6.485163  6.485163  6.766178  6.766178
##  [78]  6.766178  6.766178  6.766178  7.031888  7.031888  7.031888  7.031888
##  [85]  7.031888  7.807724  7.807724  7.807724  7.807724  7.807724  8.451728
##  [92]  8.451728  8.451728  8.451728  8.451728  9.004303  9.004303  9.004303
##  [99]  9.004303  9.004303

dat3<-rbind(gau(a=col,c=rev(col.opt),t=c(20,25,30,35,40)),gau(a=col,c=rev(col.opt),t=c(20,25,30,35,40)),gau(a=col,c=rev(col.opt),t=c(20,25,30,35,40)),gau(a=col,c=rev(col.opt),t=c(20,25,30,35,40)),gau(a=col,c=rev(col.opt),t=c(20,25,30,35,40)))
#dat3
dat3$colonies<-as.factor(rep(sort(rep(seq(1,20,1),5)),5))
dat3$t<-as.factor(dat2$t)

Mixed effects model to estimate G

###model
hwmod<-lmer(formula=performance~1+(0+t|colonies),REML=TRUE,data=dat3)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 6 negative eigenvalues

## Warning: Model failed to converge with 2 negative eigenvalues: -2.1e+00
## -6.7e+04

#summary(hmod)
Gmatrix.pac<-VarCorr(hwmod)$colonies;Gmatrix.pac

##             t20        t25        t30        t35        t40
## t20  16.3353944  0.7855368 -25.582030  0.7855365  12.220418
## t25   0.7855368 15.5462116  -2.916343 15.5462118   3.292698
## t30 -25.5820304 -2.9163432  57.650796 -2.9163432 -25.399371
## t35   0.7855365 15.5462118  -2.916343 15.5462120   3.292698
## t40  12.2204181  3.2926979 -25.399371  3.2926981  18.374285
## attr(,"stddev")
##      t20      t25      t30      t35      t40 
## 4.041707 3.942868 7.592812 3.942868 4.286524 
## attr(,"correlation")
##             t20         t25         t30         t35        t40
## t20  1.00000000  0.04929347 -0.83361889  0.04929346  0.7053685
## t25  0.04929347  1.00000000 -0.09741452  1.00000000  0.1948204
## t30 -0.83361889 -0.09741452  1.00000000 -0.09741452 -0.7803961
## t35  0.04929346  1.00000000 -0.09741452  1.00000000  0.1948204
## t40  0.70536845  0.19482038 -0.78039609  0.19482039  1.0000000

Gmatrix.pac[1:5,1:5]

##             t20        t25        t30        t35        t40
## t20  16.3353944  0.7855368 -25.582030  0.7855365  12.220418
## t25   0.7855368 15.5462116  -2.916343 15.5462118   3.292698
## t30 -25.5820304 -2.9163432  57.650796 -2.9163432 -25.399371
## t35   0.7855365 15.5462118  -2.916343 15.5462120   3.292698
## t40  12.2204181  3.2926979 -25.399371  3.2926981  18.374285

PCA decomposition of G

h<-princomp(Gmatrix.pac[1:5,1:5])
h$loadings[,1:3]

##         Comp.1      Comp.2      Comp.3
## t20  0.3908828  0.15080275  0.55119684
## t25  0.0669075 -0.69837362  0.04178372
## t30 -0.8225939 -0.02522854 -0.13279703
## t35  0.0669075 -0.69837363  0.04178365
## t40  0.4019908  0.03421425 -0.82161746

summary(h)

## Importance of components:
##                            Comp.1      Comp.2      Comp.3       Comp.4
## Standard deviation     37.0299793 10.41508906 2.273783877 3.354236e-07
## Proportion of Variance  0.9234649  0.07305323 0.003481864 7.577058e-17
## Cumulative Proportion   0.9234649  0.99651814 1.000000000 1.000000e+00
##                              Comp.5
## Standard deviation     1.299053e-07
## Proportion of Variance 1.136494e-17
## Cumulative Proportion  1.000000e+00

#plot(h$loadings[,1])
h3<-data.frame(Loadings=h$loadings[,1],t=c(20,25,30,35,40))

cc<-ggplot(h3,aes(x=t,y=Loadings))+geom_line(size=3)+setup+geom_hline(yintercept=0)+ylim(-1,1)
cc

y axis = loadings
x axis = environmental gradient

Figure

dat3$t<-as.numeric(as.character(dat3$t))
c<-ggplot(dat3,aes(x=t,y=performance,colour=factor(colonies)))+stat_smooth(se=FALSE)+setup+scale_color_hue(l=20, c=50)
c

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Summary, all figs together

grid.arrange(ncol=2,nrow=3,a,aa,b,bb,c,cc)

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

References:

using mixed effects models, I followed Paccard et al. 2016:

Paccard, A., Van Buskirk, J., & Willi, Y. (2016). Quantitative Genetic Architecture at Latitudinal Range Boundaries: Reduced Variation but Higher Trait Independence. The American Naturalist, 187(5), 667–77. https://doi.org/10.1086/685643

The idea of decomposing the population level variation (with known kinship) of performance curves , I followed Joel Kingsolver’s work:

Kingsolver, J.G., N. Heckman, J. Zhang, P.A. Carter, J.L. Knies, J.R. Stinchcombe, and K. Meyer. 2015. Genetic variation, simplicity, and evolutionary constraints for function-valued traits. The American Naturalist 185: E166-E181.

Stinchcombe, J. R., Function-valued Traits Working Group, M. Kirkpatrick. 2012. Genetics and evolution of function-valued traits: understanding environmentally responsive phenotypes. Trends in Ecology and Evolution.

Kingsolver, J. G., G. J. Ragland, and J. G. Shlichta. 2004. Quantitative genetics of continuous reaction norms: Thermal sensitivity of caterpillar growth rates. Evolution 58:1521-1529.

SessionInfo

sessionInfo()

## R version 4.3.2 (2023-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] gridExtra_2.3     DiagrammeR_1.0.10 lmerTest_3.1-3    lme4_1.1-34      
## [5] Matrix_1.6-1      ggplot2_3.4.3    
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.7          utf8_1.2.3          generics_0.1.3     
##  [4] lattice_0.21-9      digest_0.6.33       magrittr_2.0.3     
##  [7] evaluate_0.21       grid_4.3.2          RColorBrewer_1.1-3 
## [10] fastmap_1.1.1       jsonlite_1.8.7      mgcv_1.9-0         
## [13] fansi_1.0.4         scales_1.2.1        numDeriv_2016.8-1.1
## [16] jquerylib_0.1.4     cli_3.6.1           crayon_1.5.2       
## [19] rlang_1.1.1         visNetwork_2.1.2    ellipsis_0.3.2     
## [22] munsell_0.5.0       splines_4.3.2       withr_2.5.0        
## [25] cachem_1.0.8        yaml_2.3.7          tools_4.3.2        
## [28] nloptr_2.0.3        minqa_1.2.5         dplyr_1.1.2        
## [31] colorspace_2.1-0    boot_1.3-28.1       vctrs_0.6.3        
## [34] R6_2.5.1            lifecycle_1.0.3     htmlwidgets_1.6.2  
## [37] MASS_7.3-60         pkgconfig_2.0.3     pillar_1.9.0       
## [40] bslib_0.5.1         gtable_0.3.4        glue_1.6.2         
## [43] Rcpp_1.0.11         xfun_0.40           tibble_3.2.1       
## [46] tidyselect_1.2.0    highr_0.10          rstudioapi_0.15.0  
## [49] knitr_1.43          farver_2.1.1        htmltools_0.5.6    
## [52] nlme_3.1-163        labeling_0.4.2      rmarkdown_2.24     
## [55] compiler_4.3.2

Simulating performance curves

Andrew D. Nguyen, Quantitative Biologist

2023-12-26

Libaries

Workflow

Constructing gaussian shaped performance curves

Simulate data with just vertical variation/responses

Simulating the data

Mixed effects model to estimate G

PCA decomposition of G

Figure

Simulate horizontal variation

Mixed effects model to Estimate G

PCA decomposition of G

Figure

Simulate generalist-specialist variation

making the data

Mixed effects model to estimate G

PCA decomposition of G

Figure

Summary, all figs together

References:

SessionInfo