-
Notifications
You must be signed in to change notification settings - Fork 5
/
nowcasting_age.R
109 lines (95 loc) · 3.99 KB
/
nowcasting_age.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
#' @title nowcasting_age
#'
#' @description Run INLA model on structured data by age-class
#' data has to be in the format of delay-triangle
#'
#' @param dataset data pre formatted in to age classes and delays by week for each cases, delay triangle format
#'
#' @return Trajectories from the inner 'INLA' model
#' @export
nowcasting_age <- function(dataset,
zero_inflated=FALSE){
## Cehcl for zero-inflated
if (zero_inflated){
family <- "zeroinflatednbinomial1"
control.family <- list(
hyper = list("theta1" = list(prior = "loggamma",
param = c(0.01, 0.01)),
"theta2" = list(prior = "gaussian",
param = c(0, 0.4)))
)
} else {
family <- 'nbinomial'
control.family <- list(
hyper = list("theta" = list(prior = "loggamma",
param = c(0.001, 0.001)))
)
}
index.missing <- which(is.na(dataset$Y))
dataset <- dataset |>
dplyr::mutate(
fx_etaria.num = as.numeric(fx_etaria))
## Model equation: intercept + age + f(time random effect | age) + f(Delay random effect | age)
## Y(t,Age) ~ 1 + Age + rw2(t,Age) + rw1(delay, Age),
## prec(rw2) ~ logGamma(10e-3, 10e-3), prec(rw1) ~ logGamma(10e-3, 10e-3)
model <- Y ~ 1 + fx_etaria +
f(Time,
model = "rw2",
hyper = list("prec" = list(prior = "loggamma",
param = c(0.001, 0.001))
),
group = fx_etaria.num, control.group = list(model = "iid")) +
f(delay, model = "rw1",
hyper = list("prec" = list(prior = "loggamma",
param = c(0.001, 0.001))),
group = fx_etaria.num, control.group = list(model = "iid")
)
## Running the Negative Binomial model in INLA
output0 <- INLA::inla(model, family = family,
data = dataset,
control.predictor = list(link = 1, compute = T),
control.compute = list( config = T, waic=F, dic=F),
control.family = control.family
)
## Algorithm to get samples for the predictive distribution for the number of cases
## Step 1: Sampling from the approximate posterior distribution using INLA
srag.samples0.list <- INLA::inla.posterior.sample(n = 1000, output0)
## Give a parameter to trajectories, TO-DO
## Step 2: Sampling the missing triangle from the likelihood using INLA estimates
vector.samples0 <- lapply(X = srag.samples0.list,
FUN = function(x, idx = index.missing){
if(zero_inflated){
unif.log <- as.numeric(runif(idx,0,1) < x$hyperpar[2])
}else{
unif.log = 1
}
stats::rnbinom(n = idx,
mu = exp(x$latent[idx]),
size = x$hyperpar[1]
) * unif.log
} )
## Step 3: Calculate N_{a,t} for each triangle sample {N_{t,a} : t=Tactual-Dmax+1,...Tactual}
gg.age <- function(x, dados, idx){
data.aux <- dados
Tmin <- min(dados$Time[idx])
data.aux$Y[idx] <- x
data.aggregated <- data.aux |>
## Selecionando apenas os dias faltantes a partir
## do domingo da respectiva ultima epiweek
## com dados faltantes
dplyr::filter(Time >= Tmin ) |>
dplyr::group_by(Time, dt_event, fx_etaria, fx_etaria.num) |>
dplyr::summarise(
Y = sum(Y), .groups = "keep"
)
data.aggregated
}
## Step 4: Applying the age aggregation on each posterior
tibble.samples.0 <- lapply( X = vector.samples0,
FUN = gg.age,
dados = dataset,
idx = index.missing)
srag.pred.0 <- dplyr::bind_rows(tibble.samples.0, .id = "sample")
## Returning nowcasting estimation
return(srag.pred.0)
}