rethomics / zeitgebr

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

Chi square periodogram has dips dues to folding #3

Closed qgeissmann closed 6 years ago

qgeissmann commented 6 years ago

follow up of https://github.com/rethomics/ggetho/issues/16

qgeissmann commented 6 years ago

I propose a sparse implementation to deal with imcomplete blocks:

#' Calculate Qp
#'
#' @param values activity values (each value represents the measured activity in a minute)
#' @param target_period a period at which the chi-squared statistics is to be calculated
#'
#' @return a numeric of the calculated chi-squared statistics at the given varPer
#' @noRd
calc_Qp <- function(target_period, values, sampling_rate){
  col_num <- round(target_period * sampling_rate)
  #row_num <- ceiling(length(values) / col_num)
  row_num <- length(values) / col_num
  dt <- data.table::data.table( col = (0 : (length(values) -1) %% col_num) + 1,
                                #row = ceiling(1:length(values) / col_num),
                                values = values,
                                key = "col")
  avg_P <- dt[, .(avg_P = mean(values)),by=col]$avg_P
  avg_all <- mean(values)
  numerator <- sum((avg_P - avg_all) ^ 2) *  (nrow(dt) * row_num)

  denom <- sum((values - avg_all) ^ 2)
  numerator / denom
}

corrected_xsp.pdf

in the enclosed pdf, I varied the length of the input time series and applies the naive implementation (blue) vs ours (black). More robustness, also more power as we discard less data! In addition, we get rid of the artefact (see the clear dips in the blue curve when series is 6 days long). Any thoughts @pepelisu ?

qgeissmann commented 6 years ago

btw, we in the original issue: image

now: image

Done with:

library(zeitgebr)
library(ggetho)
data(dams_sample)
summary(dams_sample)
per <- periodogram(activity, dams_sample)
ggetho::ggperio(per, aes(y = power - signif_threshold, colour=period_group)) + stat_pop_etho()