Simulate traits data to show relations with growth

ExampleAGR.R

# Gist to simulate Growth Rate and Trait distribution
# Effects of varying coefficients
# Author: Matthias Grenié

# Packages ####

library(dplyr)  # Manipulate data
library(gtools)  # Combinations
library(ggplot2)  # Plots

# Constants ####
n.spec = 5  # Number of simulated species
n.ind = 100  # Number of individual by species
trait.sd = 3  # Standard deviation to use when simulating trait value
alpha = 0.5  # Coefficient for simple relation with species trait and growth

s = seq(-1.5, 1.5, 0.5)  # Coef sequence

# Creates a data.frame of coefs value to be tested
comb = combinations(length(s), 2, s, repeats.allowed = T) %>% as.data.frame()

colnames(comb) = c("alpha1", "alpha2")

# Function to simulate data ####

simul.data = function(alpha, alpha1, alpha2, seed = 1){
  set.seed(seed)  # Specify seed
  
  # Choose mean trait values
  spec.traits = runif(n.spec, min = 10, max = 50)
  
  # Compute individual trait values
  ind.traits = lapply(spec.traits, function(x) rnorm(n.ind, mean = x, sd = trait.sd))
  
  # Create data frame with individual and species measure with species names
  tot.df = data.frame(meanIndiv = unlist(ind.traits),
                      meanSp = rep(spec.traits, each = n.ind),
                      species = rep(paste("spec", seq(1, n.spec), sep = ""),
                                    each = n.ind))
  
  # Modify data frame to add growth rate with various relations
  tot.df = tot.df %>%
    mutate(lin.agr = alpha * meanSp + rnorm(1, sd = 1.5),
           real.agr = alpha1 * meanSp + alpha2 * (meanIndiv - meanSp) +
             rnorm(1, sd = trait.sd),
           abs.agr = alpha1 * meanSp + alpha2 * abs(meanIndiv - meanSp) +
             rnorm(1, sd = trait.sd),
           alpha1 = alpha1, alpha2 = alpha2)
}

# Generate data ####

df.list = apply(comb, 1, function(x) simul.data(alpha, x[1], x[2]))

# Merge list of dfs from 
# http://stackoverflow.com/questions/8091303/merge-multiple-data-frames-in-a
# -list-simultaneously

merged = Reduce(function(...) merge(..., all=T), df.list)

# Plots ####

# Plot of individual trait vs. AGR computed with real diff with species avg
ggplot(merged, aes(x = meanIndiv, y = real.agr, colour = species)) +
  geom_point() +
  facet_grid(alpha1~alpha2) +
  labs(x = "Individual's mean trait", y = "AGR related to real difference to species mean") +
  theme(panel.margin = unit(0, "lines"),
        # Be sure to put fill = NA for panel.border otherwise panel are blank
        panel.border = element_rect(color = "black", fill = NA))

# AGR with abs diff with species avg
ggplot(merged, aes(x = meanIndiv, y = abs.agr, colour = species)) +
  geom_point() +
  facet_grid(alpha1~alpha2) +
  labs(x = "Individual's mean trait", y = "AGR related to absolute difference to species mean") +
  theme(panel.margin = unit(0, "lines"),
        # Be sure to put fill = NA for panel.border otherwise panel are blank
        panel.border = element_rect(color = "black", fill = NA))