rethomics / zeitgebr

Circadian rhythm analysis and visualisation in R
11 stars 7 forks source link

Pipeline to find peaks #5

Closed qgeissmann closed 6 years ago

qgeissmann commented 6 years ago
per <- periodogram(...)
peaks <- find_peaks(per, 
                 n_peaks=3, # how many peaks to return
                 use_raw = FALSE # whethert we use raw power vs p-value (latter is default)
                 ...# other arguments?
)

peaks will have the shape of metadata We could also make a new stat to show the peaks on a ggplot...

qgeissmann commented 6 years ago

progress so far: adapting pracma::findpeaks():

find_peaks <- function(d, n_peaks = 3){
  x <- d[, power - signif_threshold]
  out <- findpeaks(x, npeaks=n_peaks, minpeakheight =  0, sortstr=TRUE)
  peak_idx <- rep(NA_real_,n_peaks)
  if(!is.null(out))
    peak_idx[1:nrow(out)]  <- out[,2]
  peak_period <- as.list(d[, period][peak_idx])
  names(peak_period) <- sprintf("peak_%i", 1:n_peaks)
  peak_pow <- as.list(d[, power][peak_idx])
  names(peak_pow) <- sprintf("peak_pow_%i", 1:n_peaks)
  c(peak_period, peak_pow)
}

data(dams_sample)
dt <- copy(dams_sample)
dt[, moving := activity > 0]
per_dt_xs <- periodogram(moving, dt, FUN = chi_sq_periodogram)

peaks <- per_dt_xs[, find_peaks(.SD), by=id]
ggetho::ggperio(per_dt_xs) + geom_line()  +
  geom_line(aes(y=signif_threshold), colour="blue") +
  geom_point(data = peaks, aes(x=peak_1, y=peak_pow_1), col="red") +
  facet_wrap( ~ id, ncol = 4, labeller = id_labeller)

image

and boxplot of peaks:

ggplot(rejoin(peaks), aes(y=peak_1, x=period_group, fill=period_group)) + 
                geom_boxplot() +
                geom_jitter() + 
                scale_y_time()

image

qgeissmann commented 6 years ago

two changes in approach came up:

  1. find peaks add a peak column in a periodogram behavr table
  2. we made a new ggplot layer (geom_peak), that can use the peak aes

alltogether, the mnew layer can be used to draw peaks and values:


ggetho::ggperio(per_dt_xs_with_peaks, aes(peak=peak)) + geom_line()  +
  geom_line(aes(y=signif_threshold), colour="red") +
  geom_peak(colour="blue") +
  facet_wrap( ~ id, ncol = 8, labeller = id_labeller)

image