inst/doc/circacompare.R

In circacompare: Analyses of Circadian Data

## ---- include=F---------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 5,
  fig.height = 3.5
)
set.seed(42)

## ---- message=F, warning=F----------------------------------------------------
library(circacompare)
library(ggplot2)

## -----------------------------------------------------------------------------
set.seed(42)
data_single <- make_data(k1 = 0, alpha1 = 0, phi1 = 0, noise_sd = 1)[c("time", "measure")]
head(data_single)

## -----------------------------------------------------------------------------
data_grouped <- make_data(phi1 = 6, noise_sd = 1)
head(data_grouped)
tail(data_grouped)

## -----------------------------------------------------------------------------
result <- circa_single(x = data_single, col_time = "time", col_outcome = "measure", period = 24)

result

## -----------------------------------------------------------------------------
result2 <- circacompare(x = data_grouped, col_time = "time", col_group = "group", col_outcome = "measure")

result2

## ---- warning=FALSE, message=FALSE--------------------------------------------
phi1_in <- 3
mixed_data <- function(n) {
  counter <- 1
  for (i in 1:n) {
    x <- make_data(k1 = 0, alpha1 = 5, phi1 = rnorm(1, phi1_in, 1), hours = 72, noise_sd = 2)
    x$id <- counter
    counter <- counter + 1
    if (i == 1) {
      res <- x
    } else {
      res <- rbind(res, x)
    }
  }
  return(res)
}
df <- mixed_data(20)
out <- circacompare_mixed(
  x = df,
  col_time = "time",
  col_group = "group",
  col_outcome = "measure",
  col_id = "id",
  control = list(grouped_params = c("alpha", "phi"), random_params = c("phi1")),
  period = 24,
  suppress_all = TRUE
)

ggplot(data = df[df$id %in% c(1:6), ], aes(time, measure)) +
  geom_point(aes(col = group)) +
  geom_smooth(aes(group = interaction(as.factor(id), group)), span = 0.3) +
  facet_wrap(~id)

## -----------------------------------------------------------------------------
out$plot

summary <- as.data.frame(circacompare:::extract_model_coefs(out$fit))
summary$`95% CI ll` <- summary$estimate - summary$std_error * 1.96
summary$`95% CI ul` <- summary$estimate + summary$std_error * 1.96

summary

## -----------------------------------------------------------------------------
tau_in <- runif(1, 8, 20)
alpha_decay1_in <- runif(1, 0.02, 0.05)

df <- make_data(k1 = 0, alpha1 = 10, phi1 = 0, seed = 42, hours = 120, noise_sd = 2)
df$time <- df$time / 24 * tau_in

# Note that when decay is estimated, it is on a scale of time in hours. There is no relation decay rate and the period.
df$measure[df$group == "g2"] <- df$measure[df$group == "g2"] * exp(-alpha_decay1_in * df$time[df$group == "g2"])

out_alpha_decay <-
  circacompare(
    x = df, "time", "group", "measure", period = NA,
    control = list(
      main_params = c("k", "alpha", "phi", "tau"),
      decay_params = c("alpha"),
      grouped_params = c("alpha", "alpha_decay"),
      period_min = 2, period_max = 24
    )
  )

out_alpha_decay$plot

summary <- as.data.frame(circacompare:::extract_model_coefs(out_alpha_decay$fit))
summary$`95% CI ll` <- summary$estimate - summary$std_error * 1.96
summary$`95% CI ul` <- summary$estimate + summary$std_error * 1.96
summary$p_value <- NULL
summary

## -----------------------------------------------------------------------------
cat("Real period: ", tau_in, "\n",
  "Real alpha_decay: ", alpha_decay1_in,
  sep = ""
)