Stantec Demo code

stantec_2020-03-03_demo.R

#--------------------------------------------#
# ARTEMIS Demo
# Stantec presentation, 2020-03-03
# M. Johnston, Cramer Fish Sciences
# This script presents an introduction to data simulation in {artemis}.
#--------------------------------------------#

# Goal: Simlulate data and visualize relationships

#-------------------------------------------------------#
# Load packages and data

library(artemis)

# model existing data to get a sense of effect size scales
head(eDNA_data)

ggplot(eDNA_data) +
  geom_jitter(aes(x = factor(Distance), y = Cq, 
                  color = factor(Volume)), 
              width = 0.2) +
  theme_minimal()

m1 <- eDNA_lm(Cq ~ Distance + Volume,
              data = eDNA_data,
              std_curve_alpha = 21.2, 
              std_curve_beta = -1.5,
              verbose = TRUE)

summary(m1) # remember that the scale of the betas is the effect on log[eDNA], not Cq.  But they're interpreted just like any other linear model.

#-------------------------------------------------------#
# Simulate data
#-------------------------------------------------------#
# step 1 - create a list of variables of interest and their levels
vars = list(Cq = 1,
         Intercept = 1,
         Distance = c(0, 50, 100),
         Volume = c(50, 100),
         Dead = as.factor(c(0, 1)),
         Biomass = c(1, 22, 100),
         tech_reps = 1:6,
         reps = 1:5
        )

# step 2 - create vector of betas
betas = c(Intercept = -10.96, 
          Distance = -0.04, 
          Volume = 0.001,
          Biomass = 0.02)

# step 3 - simulate data
sims = sim_eDNA_lm(Cq ~ Distance + Volume + Biomass,
                  variable_list = vars,
                  betas = betas,
                  sigma_Cq = 3.5,
                  std_curve_alpha = 21.2, 
                  std_curve_beta = -1.5
                  )

plot(sims)

# convert to a dataframe to make custom plots
ggsims = data.frame(sims)
head(ggsims)

ggplot(ggsims, aes(x = factor(Distance), y = Cq_star)) +
  geom_jitter(aes(color = factor(Volume)), 
              width = 0.1, size = 1.5, alpha = 0.5) +
  theme_minimal()

#-------------------------------------------------------#  # More simulations
# add or subtract variables from the model formula:
vars

betas = c(Intercept = -10.96, 
          Distance = -0.04, 
          Volume = 0.01,
          Biomass = 0.03,
          Dead = 2.9) # binary predictors - effect needs to be large, or data neets to be legion, to observe

sims2 = sim_eDNA_lm(Cq ~ Distance + Volume + Biomass + Dead,
                   variable_list = vars,
                   betas = betas,
                   sigma_Cq = 3.5,
                   std_curve_alpha = 21.2, 
                  std_curve_beta = -1.5
                  )

plot(sims2)