Microglia depletion alleviates the disruptions in circadian rhythms and anxiety caused by bmal1 deficiency.

1. #' @title circa_single
2. #' @name circa_single
3. #'
4. #' @description \code{circa_single} performs an analysis on a single rhythmic dataset. It estimates the mesor, amplitude and phase of the data provided.
5. #'
6. #' @param x data.frame. This is the data.frame which contains the rhythmic data in a tidy format.
7. #' @param col_time The name of the column within the data.frame, x, which contains time in hours at which the data were collected.
8. #' @param col_outcome The name of the column within the data.frame, x, which contains outcome measure of interest.
9. #' @param period The period of the rhythm. For circadian rhythms, leave this as the default value, 24.
10. #' @param alpha_threshold The level of alpha for which the presence of rhythmicity is considered. Default is 0.05.
11. #' @param timeout_n The upper limit for the model fitting attempts. Default is 10,000.
12. #' @param return_figure Whether or not to return a ggplot graph of the rhythm and cosine model.
13. #' @param control \code{list}. Used to control the parameterization of the model.
14. #' @param weights An optional numeric vector of (fixed) weights. When present, the objective function is weighted least squares.
15. #' @param suppress_all Logical. Set to \code{TRUE} to avoid seeing errors or messages during model fitting procedure. Default is \code{FALSE}.
16. #'
17. #' @return list
18. #' @export
19. #'
20. #' @examples
21. #' df <- make_data()
22. #' df <- df[df$group == "g1", ]
23. #' out <- circa_single(x = df, col_time = "time", col_outcome = "measure")
24. #' out
25. #'
26. #' # with sample weights (arbitrary weights for demonstration)
27. #' sw <- runif(n = nrow(df))
28. #' out2 <- circa_single(
29. #' x = df,
30. #' col_time = "time",
31. #' col_outcome = "measure",
32. #' weights = sw,
33. #' suppress_all = TRUE
34. #' )
35. #' out2
36. #'
37. circa_single <- function(x,
38. col_time,
39. col_outcome,
40. period = 24,
41. alpha_threshold = 0.05,
42. timeout_n = 10000,
43. return_figure = TRUE,
44. control = list(),
45. weights = NULL,
46. suppress_all = FALSE) {
47. controlVals <- circa_single_control()
48. controlVals[names(control)] <- control
50. if (controlVals$period_param) {
51. controlVals$main_params <- c(controlVals$main_params, "tau")
52. }
53. if ("tau" %in% controlVals$main_params) {
54. controlVals$period_param <- TRUE
55. }
57. x <- x[c(col_time, col_outcome)]
58. colnames(x) <- c("time", "measure")
60. if (!class(x$time) %in% c("numeric", "integer")) {
61. stop(paste("The time variable which you gave was a '",
62. class(x$time), "' \nThis function expects time to be given as hours and be of class 'integer' or 'numeric'.",
63. "\nPlease convert the time variable in your dataframe to be of one of these classes",
64. sep = ""
65. ))
66. }
68. if (!class(x$measure) %in% c("numeric", "integer")) {
69. stop(paste("The measure variable which you gave was a '",
70. class(x$measure), "' \nThis function expects measure to be number and be of class 'integer' or 'numeric'.",
71. "\nPlease convert the measure variable in your dataframe to be of one of these classes",
72. sep = ""
73. ))
74. }
76. if (!class(period) %in% c("numeric", "integer") & !controlVals$period_param) {
77. stop(paste0(
78. "The period argument must be a number representing the period of the rhythm in hours\n",
79. "If you would like the period to be estimated as part of the model, use:\ncontrol=list(period_param=TRUE)"
80. ))
81. }
83. if (controlVals$period_param & !is.na(period)) {
84. message(paste0("control$period_param is TRUE\n'period=", period, "' is being ignored.\nSet 'period=NA' to avoid this message"))
85. }
86. x$time_r <- x$time * 2 * pi
87. if (!controlVals$period_param) {
88. x$period <- period
89. } else {
90. if (is.null(controlVals$period_min) | is.null(controlVals$period_min)) {
91. message(paste0(
92. "If you want the model to estimate the period using a parameter,",
93. "you may get faster convergence if you provide an approximate range using 'period_min' and 'period_max' in control()",
94. "\nCurrently assuming period is between: period_min=", controlVals$period_min,
95. "and period_max=", controlVals$period_max
96. ))
97. }
98. }
100. if (!is.null(weights)) {
101. check_weights(x, weights)
102. x$weights <- weights
103. } else {
104. x$weights <- rep(1, nrow(x))
105. }
107. success <- FALSE
108. n <- 0
109. form <- create_formula(main_params = controlVals$main_params, decay_params = controlVals$decay_params)$formula
111. while (!success) {
112. fit.nls <- try(
113. {
114. stats::nls(
115. formula = form,
116. data = x,
117. start = start_list(outcome = x$measure, controlVals = controlVals),
118. weights = weights
119. )
120. },
121. silent = suppress_all
122. )
123. if (inherits(fit.nls, "try-error")) {
124. n <- n + 1
125. } else {
126. nls_coefs <- extract_model_coefs(fit.nls)
127. V <- nls_coefs[, "estimate"]
128. success <- assess_model_estimates(param_estimates = V, controlVals = controlVals)
129. n <- n + 1
130. }
131. if (n >= timeout_n) {
132. stop("Failed to converge data prior to timeout. \nYou may try to increase the allowed attempts before timeout by increasing the value of the 'timeout_n' argument or setting a new seed before this function.\nIf you have repeated difficulties, please contact me (via github) or Oliver Rawashdeh (contact details in manuscript).")
133. }
134. }
135. data_rhythmic <- nls_coefs["alpha", "p_value"] < alpha_threshold
137. if (!controlVals$period_param) {
138. V["tau"] <- period
139. }
141. eq_expression <- create_formula(main_params = controlVals$main_params, decay_params = controlVals$decay_params)$f_equation
142. eval(parse(text = eq_expression))
145. if (return_figure) {
146. p <- ggplot2::ggplot(x, ggplot2::aes(time, measure)) +
147. ggplot2::geom_point() +
148. ggplot2::xlim(
149. min(floor(x$time / period) * period),
150. max(ceiling(x$time / period) * period)
151. )
152. if (data_rhythmic) {
153. fig_out <- p +
154. ggplot2::stat_function(fun = eq, linewidth = 1) +
155. ggplot2::labs(subtitle = "Data is rhythmic", x = "time (hours)")
156. } else {
157. fig_out <- p +
158. ggplot2::labs(subtitle = "Data is arrhythmic", x = "time (hours)")
159. }
160. }
162. results_summary <-
163. circa_summary(model = fit.nls, period = period, control = controlVals)
165. if (return_figure) {
166. return(list(fit = fit.nls, summary = results_summary, plot = fig_out))
167. } else {
168. return(list(fit = fit.nls, summary = results_summary))
169. }
170. }
173. circa_single_control <- function(period_param = F, period_min = 20, period_max = 28,
174. main_params = c("k", "alpha", "phi"), decay_params = c()) {
175. list(
176. period_param = period_param, period_min = period_min, period_max = period_max,
177. main_params = main_params, decay_params = decay_params
178. )
179. }

circa_single.R at commit 5df1144, under other · at the source