-
Notifications
You must be signed in to change notification settings - Fork 0
/
aedseo.R
135 lines (124 loc) · 4.12 KB
/
aedseo.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#' Automated and Early Detection of Seasonal Epidemic Onset
#'
#' @description
#' `r lifecycle::badge("stable")`
#'
#' This function performs automated and early detection of seasonal epidemic
#' onsets (aedseo) on a time series dataset. It estimates growth rates for
#' consecutive time intervals and calculates the sum of cases (sum_of_cases).
#'
#' @param tsd A `aedseo_tsd` object containing time series data with 'time' and
#' 'observed.'
#' @param k An integer specifying the window size for modeling growth rates.
#' @param level The confidence level for parameter estimates, a numeric value
#' between 0 and 1.
#' @param disease_threshold An integer specifying the threshold for considering
#' a disease outbreak.
#' @param family A character string specifying the family for modeling.
#' Choose between "poisson," or "quasipoisson".
#'
#' @return A `aedseo` object containing:
#' - 'reference_time': The time point for which the growth rate is estimated.
#' - 'observed': The observed value in the reference time point.
#' - 'growth_rate': The estimated growth rate.
#' - 'lower_growth_rate': The lower bound of the growth rate's confidence
#' interval.
#' - 'upper_growth_rate': The upper bound of the growth rate's confidence
#' interval.
#' - 'growth_warning': Logical. Is the growth rate significantly higher than
#' zero?
#' - 'sum_of_cases': The sum of cases within the time window.
#' - 'sum_of_cases_warning': Logical. Does the Sum of Cases exceed the
#' disease threshold?
#' - 'seasonal_onset_alarm': Logical. Is there a seasonal onset alarm?
#' - 'converged': Logical. Was the IWLS judged to have converged?
#'
#' @export
#'
#' @examples
#' # Create a tsibble object from sample data
#' tsd_data <- tsd(
#' observed = c(100, 120, 150, 180, 220, 270),
#' time = as.Date(c(
#' "2023-01-01",
#' "2023-01-02",
#' "2023-01-03",
#' "2023-01-04",
#' "2023-01-05",
#' "2023-01-06"
#' )),
#' time_interval = "day"
#' )
#'
#' # Calculate AEDSEO with a 3-day window and a Poisson family model
#' aedseo_results <- aedseo(
#' tsd = tsd_data,
#' k = 3,
#' level = 0.95,
#' disease_threshold = 200,
#' family = "poisson"
#' )
#'
#' # Print the AEDSEO results
#' print(aedseo_results)
aedseo <- function(
tsd,
k = 5,
level = 0.95,
disease_threshold = NA_integer_,
family = c(
"poisson",
"quasipoisson"
# TODO: #10 Include negative.binomial regressions. @telkamp7
)) {
# Throw an error if any of the inputs are not supported
family <- rlang::arg_match(family)
# Extract the length of the series
n <- base::nrow(tsd)
# Allocate space for growth rate estimates
res <- tibble::tibble()
for (i in k:n) {
# Index observations for this iteration
obs_iter <- tsd[(i - k + 1):i, ]
# Calculate growth rates
growth_rates <- fit_growth_rate(
observations = obs_iter$observed,
level = level,
family = family
)
# See if the growth rate is significantly higher than zero
growth_warning <- growth_rates$estimate[2] > 0
# Calculate Sum of Cases (sum_of_cases)
sum_of_cases <- base::sum(obs_iter$observed)
# Evaluate if sum_of_cases exceeds disease_threshold
sum_of_cases_warning <- sum_of_cases > disease_threshold
# Give an seasonal_onset_alarm if both criteria are met
seasonal_onset_alarm <- growth_warning & sum_of_cases_warning
# Collect the results
res <- dplyr::bind_rows(
res,
tibble::tibble(
reference_time = tsd$time[i],
observed = tsd$observed[i],
growth_rate = growth_rates$estimate[1],
lower_growth_rate = growth_rates$estimate[2],
upper_growth_rate = growth_rates$estimate[3],
growth_warning = growth_warning,
sum_of_cases = sum_of_cases,
sum_of_cases_warning = sum_of_cases_warning,
seasonal_onset_alarm = seasonal_onset_alarm,
converged = growth_rates$fit$converged
)
)
}
# Turn the results into an `aedseo` class
ans <- tibble::new_tibble(
x = res,
class = "aedseo",
k = k,
level = level,
disease_threshold = disease_threshold,
family = family
)
return(ans)
}