-
Notifications
You must be signed in to change notification settings - Fork 2
/
02_run_fit.R
115 lines (90 loc) · 3.39 KB
/
02_run_fit.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
# Notes -------------------------------------------------------------------
# Master file to fit models for different severity items and scores
# - Extent: RW, BinRW, BinMC
# - Intensity signs: MC, OrderedRW
# - Subjective symptoms (and B and C): RW, BinRW
# - SCORAD and oSCORAD: RW, AR1, MixedAR1 and Smoothing
# NB: for itching, sleep and C, the resolution is reso=0.1 (there are decimal values) and so the score are multiplied by 1 / reso ...
# ... to work with integers. Replications (and some parameters) must be scaled by reso.
# Initialisation ----------------------------------------------------------
rm(list = ls()) # Clear Workspace (better to restart the session)
set.seed(1744834965) # Reproducibility (Stan use a different seed)
source(here::here("analysis", "00_init.R"))
#### OPTIONS
mdl_name <- "BinMC"
score <- "extent"
dataset <- "PFDC"
run <- FALSE
n_chains <- 4
n_it <- 2000
####
item_dict <- detail_POSCORAD()
score <- match.arg(score, item_dict[["Name"]])
intensity_signs <- detail_POSCORAD("Intensity signs")$Name
dataset <- match.arg(dataset, c("Derexyl", "PFDC", "Fake"))
mdl_name <- match.arg(mdl_name, available_models(score)$Model)
stopifnot(is_scalar_logical(run),
is_scalar_wholenumber(n_chains),
n_chains > 0,
is_scalar_wholenumber(n_it),
n_it > 0)
is_continuous <- (score %in% c("SCORAD", "oSCORAD"))
item_dict <- item_dict %>% filter(Name == score)
item_lbl <- as.character(item_dict[["Label"]])
max_score <- item_dict[["Maximum"]]
reso <- case_when(is_continuous ~ 1,
TRUE ~ item_dict[["Resolution"]])
M <- round(max_score / reso)
file_dict <- get_results_files(outcome = score, model = mdl_name, dataset = dataset, root_dir = here())
param <- list_parameters(mdl_name)
param[c("Test", "Misc")] <- NULL
if (mdl_name == "BinMC") {
param$PatientTime <- c("ss1", "lambda", "y_rep")
}
# Data --------------------------------------------------------------------
POSCORAD <- load_dataset(dataset)
# Subset dataset
df <- POSCORAD %>%
# filter(Patient %in% 1:30) %>%
rename(Time = Day, Score = all_of(item_lbl)) %>%
select(Patient, Time, Score) %>%
drop_na()
pt <- unique(df[["Patient"]])
id <- get_index(df)
df <- left_join(df, id, by = c("Patient", "Time"))
# Fitting -----------------------------------------------------------------
if (run) {
if (mdl_name == "MC") {
df_MC <- df %>%
rename(y = Score) %>%
group_by(Patient) %>%
mutate(y0 = y + 1, y1 = lead(y0), dt = lead(Time) - Time) %>%
ungroup() %>%
select(y0, y1, dt) %>%
drop_na()
model <- EczemaModel("MC", K = max_score + 1)
fit <- EczemaFit(model,
train = df_MC,
pars = unlist(param),
iter = n_it,
chains = n_chains,
init = 0)
id <- NULL
} else {
if (is_continuous) {
train_tmp <- df
} else {
train_tmp <- df %>% mutate(Score = round(Score / reso))
}
model <- EczemaModel(mdl_name, max_score = M, discrete = !is_continuous)
fit <- EczemaFit(model,
train = train_tmp,
pars = unlist(param),
iter = n_it,
chains = n_chains,
control = list(adapt_delta = .9))
}
saveRDS(fit, file = file_dict$Fit)
par <- extract_parameters(fit, pars = param, id = id)
saveRDS(par, file = file_dict$FitPar)
}