In mellux-project/melluxdrc: melluxdrc

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This vignette describes how to simulate dose-response data as parameterised using estimates from Phillips et al. (2017). It also describes how types of experiment can be carried out on such data.

Simulating data

Here, we simulate dose-response data for 41 individuals measured at four lux levels.

library(melluxdrc)
library(dplyr)
library(purrr)
library(ggplot2)

# generate data
lux <- c(1, 10, 100, 1000)
n <- 41
experimental_data <- virtual_experiment(n=n, lux=lux)

# look at it
glimpse(experimental_data)

The package also contains methods for visualising such data

plot_doseresponse(experimental_data)

The estimates from Phillips et al. (2017) exhibited substantial inter-individual variability in individuals' response to light. Built into the package, we allow this individual variation to be reduced.

experimental_data <- virtual_experiment(n=n, lux=lux,
                                        individual_variation_level=0.0)
plot_doseresponse(experimental_data)

Experiments based on virtual population of individuals

The package also contains functionality for performing specific experiments based on a simulated population of dose-response type data.

Within-subject

Here, we first illustrate a within-subject experiment, in which the melatonin values for individuals measured at two lux levels are compared using a paired t-test. Success of these experiments indicates that the difference in melatonin suppression between the two lux levels is of the correct sign and statistically significant (by default at the 5% level).

# first generate large (ideally bigger than this if repeating experiments) virtual population of data
population_df <- virtual_experiment(n=200)

# compare melatonin suppression values at two lux levels using t test for a sample of 5 individuals
# repeat experiment 30 times at calculate fraction of successful results
lux_1 <- 10
lux_2 <- 30
nindiv <- 5
nreps <- 30
is_between <- F
results <- comparison_test(is_between, lux_1, lux_2, nindiv, population_df,
                           nreps=nreps) %>%
    dplyr::mutate(result=dplyr::if_else((result==1) & (p_value < 0.05), 1, 0))

mean(results$result)

If we increase the sample size the fraction of successes should increase.

nindiv <- 20
results <- comparison_test(is_between, lux_1, lux_2, nindiv, population_df,
                           nreps=nreps) %>%
    dplyr::mutate(result=dplyr::if_else((result==1) & (p_value < 0.05), 1, 0))
mean(results$result)

Similarly if we make the two luxes more separated, the fraction of successes increases.

nindiv <- 5
lux_2 <- 100
results <- comparison_test(is_between, lux_1, lux_2, nindiv, population_df,
                           nreps=nreps) %>%
    dplyr::mutate(result=dplyr::if_else((result==1) & (p_value < 0.05), 1, 0))
mean(results$result)

Between-subject

We now illustrate a between-subject experiment, in which two separate samples of individuals have their suppressions measured at two different lux levels. Because of individual variation in dose-response curves, these experiments tend to be less powerful than within-subject experiments.

lux_1 <- 10
lux_2 <- 30
nindiv <- 20
nreps <- 30
results_b <- comparison_test(TRUE, lux_1, lux_2, nindiv, population_df,
                             nreps=nreps) %>%
    dplyr::mutate(result=dplyr::if_else((result==1) & (p_value < 0.05), 1, 0))
mean(results_b$result)

results_w <- comparison_test(FALSE, lux_1, lux_2, nindiv, population_df,
                             nreps=nreps) %>%
    dplyr::mutate(result=dplyr::if_else((result==1) & (p_value < 0.05), 1, 0))
mean(results_w$result)

Experiments involving treatment

We now illustrate how this library can be used to simulate treatments which shift the natural ed50 of individuals. There are two types of experiments that can be done. One involving two samples of individuals: one treated; one not. Another involving the same individuals whom are measured before and after treatment.

We first illustrate the two sample case with a treatment which decreases the ed50 of individuals by a factor of 2.

n <- 200
sample_combined <- virtual_treatment_experiment(n, treated_ed50_multiplier = 0.5)

sample_combined %>% 
  ggplot(aes(x=lux, y=y)) +
  geom_point() +
  geom_line(aes(group=as.factor(id))) +
  facet_wrap(~treated) +
  scale_x_log10() +
  geom_smooth(se=F, method="loess", formula = "y~x") +
  scale_y_continuous(labels=scales::percent) +
  xlab("Lux") +
  ylab("Melatonin suppression")

A note on caution: the model is not able to extrapolate too far away from the measured values of Phillips et al. (2016). If the treatment is too extreme, occasionally, individuals may fall outside our region of confidence. To see this, we include a truncation Boolean which should be checked when doing these treatment-type experiments.

# ok treatment
mean(sample_combined$p1_truncated)

# extreme treatment
sample_treated <- virtual_experiment(n, treated_ed50_multiplier=0.1)
mean(sample_treated$p1_truncated) > 0

We now show simulated data for individuals whom are measured twice: before and after treatment.

n <- 9
sample_treated <- virtual_treatment_experiment(n, treated_ed50_multiplier=0.5)

sample_treated %>% 
  ggplot(aes(x=lux, y=y, colour=treated)) +
  geom_point() +
  geom_line() +
  facet_wrap(~id) +
  scale_x_log10() +
  scale_y_continuous(labels=scales::percent) +
  xlab("Lux") +
  ylab("Melatonin suppression")

# check for truncated p1 values
mean(sample_treated$p1_truncated) > 0

Again it's important to check if treated individuals have parameters resulting in over-extrapolation.

# extreme treatment
sample_treated <- virtual_treatment_experiment(n, treated_ed50_multiplier=0.01)
mean(sample_treated$p1_truncated) > 0

We can also generate treated / untreated data for two separate samples of individuals. I.e. here we are doing a between-subjects type experiment.

n <- 10
sample_treated <- virtual_treatment_experiment(n, treated_ed50_multiplier=0.5,
                                               is_between=TRUE,
                                               individual_variation_level=0.1)

sample_treated %>% 
  ggplot(aes(x=lux, y=y, group=id)) +
  geom_point() +
  geom_line() +
  facet_wrap(~treated) +
  scale_x_log10() +
  scale_y_continuous(labels=scales::percent) +
  xlab("Lux") +
  ylab("Melatonin suppression")