Skip to contents

sees-analysis

Source: vignettes/articles/sees_analysis.Rmd
sees_analysis.Rmd

library(readr)
xs_sample = 
  "extdata/sees_crossSectional_baseline_allCountries.102523.csv" |> 
  path_package(package = "serocalculator") |> 
  read_csv()

curve_params = 
  path_package(
    package = "serocalculator",
    "extdata/sees_dmcmc_09.30.2021.rds") |> 
  readRDS()

noise_params = 
  path_package(
    package = "serocalculator",
    "extdata/SEES_NoiseParam.rds") |> 
  readRDS() |> 
  rename(y.low = llod)

full_data = path_package(
  package = "serocalculator",
  "extdata/SEES_2022-10-24_redacted_2023-10-12.csv") |> 
  read_csv(
    col_types = cols(ae_diagnosis_other = col_character())
  )
library(ggplot2)
library(dplyr)
full_data |> 
  dplyr::filter(
    dag == "bangladesh",
    studyarm %in% c("pros case", "retro case"),
    # timesince0 < 60,
    antigen != "LPS"
  ) |>
  mutate(
    .by = pid,
    n_obs = n()) |> 
  arrange(desc(n_obs), pid, timesince0) |> 
  dplyr::filter(
    pid %in% unique(pid)[1:10]
  ) |> 
  mutate(
    antigen_iso = paste(antigen, isotype, sep = "_")
  ) |> 
  ggplot(
    aes(
      x = timesince0,
      y = result,
      col = antigen_iso
    )
  ) +
  geom_point() + 
  geom_smooth() +
  # geom_line() +
  scale_y_log10() +
  facet_wrap(vars(pid))
library(ggplot2)
xs_sample |>
  dplyr::filter(Country == "Bangladesh") |> 
  ggplot(aes(x = Age, y = result, col = antigen_iso)) +
  geom_point() + 
  geom_smooth() + 
  scale_y_log10()

Here we do one stratum at a time:

library(dplyr)
est1 = est.incidence(
  print.level = 2,
  start = 0.3,
  verbose = TRUE,
  pop_data = 
    xs_sample |> 
    rename(
      a = Age,
      y = result) %>% 
    dplyr::filter(Country == "Bangladesh"),
  
  curve_params = 
    curve_params |> 
    dplyr::filter(Country == "Bangladesh"),
  
  noise_params = 
    noise_params |> 
    dplyr::filter(Country == "Bangladesh"),
  
  c.age = "<5",
  
  antigen_isos = xs_sample$antigen_iso |> unique()
) |> print()

Now we will try the loop function:


est2 = est.incidence.by(
  stepmax = 1,
  verbose = TRUE,
  strata = c("ageCat", "Country"),
  lambda_start = .3,
  print.level = 2,
  num_cores = 1,
  pop_data = 
    xs_sample |> 
    rename(
      a = Age,
      y = result) |> 
    dplyr::filter(
      ageCat != "16+", # missing from curve_params...
      Country != "Nepal"
      ),
  
  curve_params = 
    curve_params |> 
    dplyr::filter(
      Country != "Nepal"
      ),
  
  noise_params = 
    noise_params |> 
    dplyr::filter(
      Country != "Nepal"
      )
)
print(est2)
a = summary(est2) |> print()