COVID-19 vaccine effectiveness against in-hospital mortality in immunocompromised older adults: a Brazilian hospital surveillance study, 2024

Preparing data

Load libraries:

library(openxlsx)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(purrr)
library(tidyr)
library(prettyR)
library(scales)

## 
## Attaching package: 'scales'

## The following object is masked from 'package:purrr':
## 
##     discard

library(DescTools)

## 
## Attaching package: 'DescTools'

## The following object is masked from 'package:prettyR':
## 
##     Mode

library(ggplot2)
library(pROC)

## Type 'citation("pROC")' for a citation.

## 
## Attaching package: 'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

library(ResourceSelection)

## ResourceSelection 0.3-6   2023-06-27

compareNA <- function(v1,v2) {
    same <- (v1 == v2) | (is.na(v1) & is.na(v2))
    same[is.na(same)] <- FALSE
    return(same)
}

Load SIVEP-Gripe data:

df <- read.xlsx("SRAG 2024_DOWNLOAD 29_3_25 IMUNODEPRE_RT_PCR OU ANTIGENO POSITIVO.xlsx", detectDates = TRUE)

It has 1436 patients before filtering.

Filtering population

filtered <- df %>%
  filter(
    NOSOCOMIAL == 2,                         # Non-nosocomial
    NU_IDADE_N >= 60,                        # 60 years old or older
    IMUNODEPRE == 1,                         # Immunosuppressed
    HOSPITAL == 1,                           # Hospitalized
    !is.na(DT_INTERNA),                      # With hospitalization date
    CLASSI_FIN == 5,                         # With COVID classification
    PCR_SARS2 == 1 | AN_SARS2 == 1,          # With positive COVID test
    EVOLUCAO %in% c(1, 2),                   # Survival or death (not other causes)
    DT_EVOLUCA >= ymd("2023-12-31"),         # Start of epidemiological year
    DT_EVOLUCA <= ymd("2024-12-28")          # End of epidemiological year
  )

After the population filter, it had 474 patients.

Data cleaning and preparing

Days since the beginning of the infection

# D: DIAS INICIO SIN HOSPITALIZAÇÃO
filtered <- filtered %>%
  mutate(delta_days = as.numeric(DT_INTERNA - DT_SIN_PRI))

#stopifnot(all(compareNA(filtered$"DIAS.INICIO.SIN.HOSPITALIZA<U+00C7><U+00C3>O", filtered$delta_days)))

Additional risk commorbidities

# AR: COMORBIDADES ADICIONAIS
filtered <- filtered %>%
  mutate(additional_risk_comorbidities = 
           as.integer(replace_na(CARDIOPATI, 0) == 1) + 
           as.integer(replace_na(HEMATOLOGI, 0) == 1) + 
           as.integer(replace_na(SIND_DOWN, 0) == 1) + 
           as.integer(replace_na(HEPATICA, 0) == 1) + 
           as.integer((replace_na(ASMA, 0) == 1) | (replace_na(PNEUMOPATI, 0) == 1)) + 
           as.integer(replace_na(DIABETES, 0) == 1) + 
           as.integer(replace_na(NEUROLOGIC, 0) == 1) + 
           as.integer(replace_na(RENAL, 0) == 1) + 
           as.integer(replace_na(OBESIDADE, 0) == 1))


#stopifnot(all(compareNA(filtered$COMORBIDADES.ADICIONAIS, filtered$additional_risk_comorbidities)))

Risk groups

filtered <- filtered %>%
  mutate(
    risk_groups = case_when(
      additional_risk_comorbidities >= 2 ~ '2+',
      additional_risk_comorbidities == 1 ~ '1',
      TRUE ~ 'None'
    ),
    at_least_one_risk = case_when(
      additional_risk_comorbidities >= 1 ~ 1,
      TRUE ~ 0
    )
  )

Dates of vaccines

filtered <- filtered %>%
  mutate(
    dose_dates = pmap(
      list(
        DOSE_1_COV,
        DOSE_2_COV,
        DOSE_ADIC,
        DOSE_REF,
        DOSE_2REF,
        DOS_RE_BI
      ),
      ~ {
        doses <- as.Date(c(...))
        doses[!is.na(doses)]
      }
    )
  )

Column that only consider vaccines that have both the date of administration and the manufacturer.

filtered <- filtered %>%
  mutate(dose_dates_with_fab = pmap(list(FAB_COV_1, DOSE_1_COV, 
                               FAB_COV_2, DOSE_2_COV, 
                               FAB_ADIC, DOSE_ADIC, 
                               FAB_COVREF, DOSE_REF, 
                               FAB_COVRF2, DOSE_2REF,
                               FAB_RE_BI, DOS_RE_BI), 
                           ~ {
                             doses <- c(
                               as.Date(ifelse(!is.na(..1), ..2, NA)),
                               as.Date(ifelse(!is.na(..3), ..4, NA)),
                               as.Date(ifelse(!is.na(..5), ..6, NA)),
                               as.Date(ifelse(!is.na(..7), ..8, NA)),
                               as.Date(ifelse(!is.na(..9), ..10, NA)),
                               as.Date(ifelse(!is.na(..11), ..12, NA))
                             )
                             doses[!is.na(doses)]
                           })
  )

Last vaccine date

# Extract the last dose date (latest non-NA date in the list)
filtered <- filtered %>%
  mutate(last_date = map(dose_dates, ~ if (length(.x) == 0) NA else max(.x, na.rm = TRUE)))

Days between vaccine date and internation

# Calculate the number of days betwenn the vaccine date and internation
filtered <- filtered %>%
  mutate(days_since_last_vaccine = if_else(
      (length(last_date) == 0 || all(is.na(last_date))) | is.na(DT_INTERNA),
      NA,
      as.numeric(as.Date(DT_INTERNA, format = "%Y-%m-%d") - as.Date(map_dbl(filtered$last_date, 1)))
    ))
           
        #   as.numeric( - ))

Intervals between vaccine date and internation

filtered <- filtered %>%
  mutate(
    last_vaccine_interval = case_when(
      days_since_last_vaccine <= 183 ~ "up to 6 months",
      days_since_last_vaccine <= 365 ~ "7–12 months",
      days_since_last_vaccine > 365  ~ "over 1 year",
      TRUE ~ NA_character_
    )
  )

filtered <- filtered %>%
  mutate(
    last_vaccine_interval_broader = case_when(
      days_since_last_vaccine <= 365 ~ "up to 1 year",
      days_since_last_vaccine > 365  ~ "over 1 year",
      TRUE ~ NA_character_
    )
  )

Last vaccine year

# Extract the year from the last_date
filtered <- filtered %>%
  mutate(last_year = map_dbl(last_date, ~ ifelse(is.na(.x), NA, year(.x))))

# Compute the minimum non-NA year
min_last_year <- min(filtered$last_year, na.rm = TRUE)

# Compute year range starting from 1
filtered <- filtered %>%
  mutate(last_year_range = ifelse(is.na(last_year), NA, last_year - min_last_year + 1))

#stopifnot(all(
#  compareNA(filtered$"ANO.DA.<U+00DA>LTIMA.DOSE.VAC_COV;.2021=1;.2022=2;.2023=3;.2024=4", filtered$last_year_range)
#))

Number of doses

# FF: N DOSES VAC_COV   (preenchendo em branco com 0)
filtered <- filtered %>%
  mutate(doses = map_int(dose_dates, length))

#stopifnot(all(compareNA(replace_na(filtered$N.DOSES.VAC_COV, 0), filtered$doses)))

Number of doses groups

filtered <- filtered %>%
  mutate(
    dose_groups = case_when(
      doses >= 3 ~ '3+',
      doses == 2 ~ '2',
      doses == 1 ~ '1',
      TRUE ~ '0'
    ),)

filtered <- filtered %>%
  mutate(
    dose_completion_groups = case_when(
      doses >= 2 ~ '2+',
      TRUE ~ '0-1'
    ),)

Types of vaccines

Define groups and check if the definition is complete

# Define vaccine groups
innactivated_vaccines <- c(
  '86 - COVID-19 SINOVAC/BUTANTAN - CORONAVAC',
  'CORONAVAC',
  'SINOVAC', 
  'SINOVAC/BUTANTAN',
  'SINOVAC/BUTANT', 
  'INSTITUTO BUTANTAN'
)

viral_vector_vaccines <- c(
  '85 - COVID-19 ASTRAZENECA/FIOCRUZ - COVISHIELD',
  '89 - COVID-19 ASTRAZENECA - CHADOX1-S',
  'ASTRAZENECA AB',
  'ASTRAZENICA/FIOCRUZ',
  'FABRICANTE FUNDACAO OSWALDO CRUZ',
  '88 - COVID-19 JANSSEN - AD26.COV2.S',
  'JANSSEN',
  'JANSSEN PHARMACEUTICA NV'
)

mrna_vaccines <- c(
  '87 - COVID-19 PFIZER - COMIRNATY',
  'PFIZER',
  'PFIZER MANUFACTURING BELGIUM NV - BELGICA',
  'PFIZER/ BIVALENTE',
  '103 - COVID-19 PFIZER - COMIRNATY BIVALENTE'
)

mapped_vaccines <- c(innactivated_vaccines, viral_vector_vaccines, mrna_vaccines, '<vazio>')

# Now run assertions similar to Python
stopifnot(
  identical(setdiff(unique(replace_na(filtered$FAB_COV_1, '<vazio>')), mapped_vaccines), character(0))
)
stopifnot(
  identical(setdiff(unique(replace_na(filtered$FAB_COV_2, '<vazio>')), mapped_vaccines), character(0))
)
stopifnot(
  identical(setdiff(unique(replace_na(filtered$FAB_ADIC, '<vazio>')), mapped_vaccines), character(0))
)
stopifnot(
  identical(setdiff(unique(replace_na(filtered$FAB_COVREF, '<vazio>')), mapped_vaccines), character(0))
)
stopifnot(
  identical(setdiff(unique(replace_na(filtered$FAB_COVRF2, '<vazio>')), mapped_vaccines), character(0))
)
stopifnot(
  identical(setdiff(unique(replace_na(as.character(filtered$FAB_RE_BI), '<vazio>')), mapped_vaccines), character(0))
)

Create columns with the count of each valid dose given the type.

filtered <- filtered %>%
  mutate(
    innactivated_doses = pmap_int(
      list(DOSE_1_COV, FAB_COV_1,
           DOSE_2_COV, FAB_COV_2,
           DOSE_ADIC, FAB_ADIC,
           DOSE_REF, FAB_COVREF,
           DOSE_2REF, FAB_COVRF2#,
           # DOS_RE_BI, FAB_RE_BI -- currently, all DOS_RE_BI are mrna vaccines
           ),
      ~ sum(
          (!is.na(..1) & ..2 %in% innactivated_vaccines),
          (!is.na(..3) & ..4 %in% innactivated_vaccines),
          (!is.na(..5) & ..6 %in% innactivated_vaccines),
          (!is.na(..7) & ..8 %in% innactivated_vaccines),
          (!is.na(..9) & ..10 %in% innactivated_vaccines)#,
          #(!is.na(..11) & ..12 %in% innactivated_vaccines)
        )
    ),
    viral_vector_doses = pmap_int(
      list(DOSE_1_COV, FAB_COV_1,
           DOSE_2_COV, FAB_COV_2,
           DOSE_ADIC, FAB_ADIC,
           DOSE_REF, FAB_COVREF,
           DOSE_2REF, FAB_COVRF2#,
           #DOS_RE_BI, FAB_RE_BI -- currently, all DOS_RE_BI are mrna vaccines
           ),
      ~ sum(
          (!is.na(..1) & ..2 %in% viral_vector_vaccines),
          (!is.na(..3) & ..4 %in% viral_vector_vaccines),
          (!is.na(..5) & ..6 %in% viral_vector_vaccines),
          (!is.na(..7) & ..8 %in% viral_vector_vaccines),
          (!is.na(..9) & ..10 %in% viral_vector_vaccines)#,
          #(!is.na(..11) & ..12 %in% viral_vector_vaccines)
        )
    ),
    mrna_doses = pmap_int(
      list(DOSE_1_COV, FAB_COV_1,
           DOSE_2_COV, FAB_COV_2,
           DOSE_ADIC, FAB_ADIC,
           DOSE_REF, FAB_COVREF,
           DOSE_2REF, FAB_COVRF2,
           DOS_RE_BI, FAB_RE_BI),
      ~ sum(
          (!is.na(..1) & ..2 %in% mrna_vaccines),
          (!is.na(..3) & ..4 %in% mrna_vaccines),
          (!is.na(..5) & ..6 %in% mrna_vaccines),
          (!is.na(..7) & ..8 %in% mrna_vaccines),
          (!is.na(..9) & ..10 %in% mrna_vaccines),
          (!is.na(..11)) # & ..12 %in% mrna_vaccines -- all DOS_RE_BI are MRNA vaccines. It does not need to check the FAB
        )
    )
  )

Vaccination schema

inactivated=1; mRNA=2; viral vector=3; inactivated+ mRNA=4; inactivated + viral vector=5; mRNA+ viral vector=6; inactivated+ mRNA + viral vector=7.

# FQ: Esquema vacinal: inativada=1; RNAm=2; vetor viral=3; inativada + RNAm=4; inativada + vetor viral=5; RNAm + vetor viral=6; inativada + RNAm + vetor viral=7
dose_df <- filtered %>%
  transmute(
    i = innactivated_doses,
    v = viral_vector_doses,
    m = mrna_doses
  )

# Apply the logic row-wise
vaccinal_schema <- pmap_int(dose_df, function(i, v, m) {
  if (i > 0 & v > 0 & m > 0) {
    7
  } else if (v > 0 & m > 0) {
    6
  } else if (v > 0 & i > 0) {
    5
  } else if (m > 0 & i > 0) {
    4
  } else if (v > 0) {
    3
  } else if (m > 0) {
    2
  } else if (i > 0) {
    1
  } else {
    0
  }
})

# Attach it back to the `filtered` dataframe
filtered <- filtered %>%
  mutate(
    vaccinal_schema = vaccinal_schema,
    nominal_vaccinal_schema = recode(
      vaccinal_schema, 
      `0` = "None",
      `1` = "Inactivated only",
      `2` = "mRNA only",
      `3` = "Viral vector only",
      `4` = "Inactivated + mRNA",
      `5` = "Inactivated + Viral vector",
      `6` = "mRNA + viral vector",
      `7` = "Inactivated + mRNA + viral vector"
    ),
    vaccinal_schema_count = recode(
      vaccinal_schema, 
      `0` = 0,
      `1` = 1,
      `2` = 1,
      `3` = 1,
      `4` = 2,
      `5` = 2,
      `6` = 2,
      `7` = 3
      )
  )

#stopifnot(all(
#  compareNA(replace_na(filtered$`Esquema.vacinal:.inativada=1;.RNAm=2;.vetor.viral=3;.inativada.+.RNAm=4;.inativada.+.vetor.viral=5;.RNAm.+.vetor.viral=6;.inativada.+.RNAm.+.vetor.viral=7`, 0), vaccinal_schema)
#))

Schema based on inactivated

filtered <- filtered %>%
  mutate(
    vaccinal_schema_based_on_innactivated = pmap_chr(list(innactivated_doses, viral_vector_doses, mrna_doses), function(i, v, m) {
      if (i + v + m <= 1) {
        return('unvaccinated or single-dose')
      } else if (v == 0 & m == 0) {
        return('sequential with inactivated only')
      } else if (i == 0) {
        return('sequential excluding inactivated')
      } else {
        return('sequential with inactivated and others')
      }
    })
  )

Age groups

filtered <- filtered %>%
  mutate(
    age_groups = case_when(
      NU_IDADE_N >= 80 ~ '80+',
      NU_IDADE_N >= 70 ~ '70-79',
      TRUE ~ '60-69'
    )
  )

Race groups

filtered <- filtered %>%
  mutate(
    nominal_race = recode(
      filtered$CS_RACA,
      `1` = "White",
      `2` = "Black",
      `3` = "Yellow",
      `4` = "Brown",
      `5` = "Indigenous",
      `9` = "Unknown"
    )
  )

Educational groups

filtered <- filtered %>%
  mutate(
    educational_attainment = case_when(
      CS_ESCOL_N == 4 ~ 'bachelor',
      CS_ESCOL_N < 4 ~ 'low',
      TRUE ~ 'unknown'
    )
  )

filtered <- filtered %>%
  mutate(
    nominal_education = replace_na(recode(
      filtered$CS_ESCOL_N,
      `0` = "Elementary",
      `1` = "Elementary",
      `2` = "Middle",
      `3` = "Highschool",
      `4` = "Bachelors",
      `5` = "Unknown",
      `9` = "Unknown"
    ), "Unknown")
  )

filtered <- filtered %>%
  mutate(
    education_level = replace_na(recode(
      filtered$CS_ESCOL_N,
      `0` = "Elementary", # Elementary
      `1` = "Elementary", # Elementary
      `2` = "Middle", # Middle
      `3` = "High", # Highschool
      `4` = "High", # Bachelors
      `5` = "Unknown",
      `9` = "Unknown"
    ), "Unknown")
  )

Dose groups

filtered <- filtered %>%
  mutate(
    categorical_doses = case_when(
      doses >= 5 ~ '5+',
      doses == 4 ~ '4',
      doses == 3 ~ '3',
      doses >= 1 ~ '1-2',
      TRUE ~ 'unvaccinated'
    )
  )

ICU admission

filtered <- filtered %>%
  mutate(
    icu_admission = replace_na(recode(UTI, `1` = 1, `2` = 2, `9` = 2, .default = 2), 2),
    nominal_icu_admission = replace_na(recode(UTI, `1` = "Yes", `2` = "No", `9` = "No", .default = "No"), "No")
  )

Invasive support ventilation

filtered <- filtered %>%
  mutate(
    invasive_support_ven = replace_na(recode(SUPORT_VEN, `1` = 1, `3` = 2, `9` = 2, .default = 2), 2),
    nominal_invasive_support_ven = replace_na(recode(SUPORT_VEN, `1` = "Yes", `3` = "No", `9` = "No", .default = "No"), "No")
  )

COVID treatment

filtered <- filtered %>%
  mutate(
    covid_treatment = case_when(
      ((TRAT_COV == 1) & !is.na(TIPO_TRAT)) ~ 1,
      ((TRAT_COV > 1) & !is.na(TIPO_TRAT)) ~ 2,
      TRUE ~ 2
    ),
    nominal_covid_treatment = replace_na(recode(covid_treatment, `1` = "Yes", `9` = "No", .default = "No"), "No")
  )

Split by outcome

filtered <- filtered %>%
  mutate(
    nominal_outcome = recode(filtered$EVOLUCAO, `1` = "Survivor", `2` = "Non-survivor"),
    outcome = ifelse(filtered$EVOLUCAO == 1, 1, 0)
  )

survivors <- filtered %>%
  filter(EVOLUCAO == 1)
nonsurvivors <- filtered %>%
  filter(EVOLUCAO == 2)

Create test spreadsheet

test_dataframe <- filtered %>%
  mutate(
    `Outcome: Survivor = 1; Non-survivor = 2` = EVOLUCAO,
    `Sex: female = 1; male = 2` = recode(CS_SEXO, `F` = 1, `M` = 2),
    `ICU admission: yes = 1; no = 2` = icu_admission,
    `Invasive ventilatory support: yes =1; no =2` = invasive_support_ven,
    `Age group: 60-69 years = 1; 70-79 years =2; 80 years or older = 3` = recode(age_groups, `60-69` = 1, `70-79` = 2, `80+` = 3),
    `Educational attainment: low (up to high school); high (Bachelor) = 2; unknown = 9` = recode(educational_attainment, `low` = 1, `bachelor` = 2, `unknown` = 9),
    `Number of other risk conditions: none = 0; one or more = 1` = if_else(additional_risk_comorbidities >= 1, 1, 0),
    `COVID-19 vaccine doses: Unvaccinated = 0; one to two doses = 1; 3 doses = 2; 4 or more doses = 3` = recode(categorical_doses, `unvaccinated` = 0, `1-2` = 1, `3` = 2, `4` = 3, `5+` = 4),
    ` ` = NA_real_, 
    `Esquemas de vacinação: unvaccinated or single-dose (reference): 1; sequential with inactivated only: 2; sequential excluding inactivated: 3; sequential with inactivated and others: 4` = recode(vaccinal_schema_based_on_innactivated, 
                                                                                                                                          `unvaccinated or single-dose` = 1,
                                                                                                                                          `sequential with inactivated only` = 2,
                                                                                                                                          `sequential excluding inactivated` = 3,
                                                                                                                                          `sequential with inactivated and others` = 4)
  ) %>%
  select(
    `Outcome: Survivor = 1; Non-survivor = 2`,
    `Sex: female = 1; male = 2`,
    `ICU admission: yes = 1; no = 2`,
    `Invasive ventilatory support: yes =1; no =2`,
    `Age group: 60-69 years = 1; 70-79 years =2; 80 years or older = 3`,
    `Educational attainment: low (up to high school); high (Bachelor) = 2; unknown = 9`,
    `Number of other risk conditions: none = 0; one or more = 1`,
    `COVID-19 vaccine doses: Unvaccinated = 0; one to two doses = 1; 3 doses = 2; 4 or more doses = 3`,
    ` `,
    `Esquemas de vacinação: unvaccinated or single-dose (reference): 1; sequential with inactivated only: 2; sequential excluding inactivated: 3; sequential with inactivated and others: 4`
  )

Create test.xlsx file

write.xlsx(test_dataframe, "test.xlsx")

Set references

filtered$outcome <- relevel(factor(filtered$outcome), ref="1")
filtered$age_groups <- relevel(factor(filtered$age_groups), ref = "60-69")
filtered$CS_SEXO <- relevel(factor(filtered$CS_SEXO), ref = "F")
filtered$nominal_race <- relevel(factor(filtered$nominal_race), ref="White")
filtered$nominal_icu_admission <- factor(filtered$nominal_icu_admission)
filtered$nominal_invasive_support_ven <- factor(filtered$nominal_invasive_support_ven)
filtered$educational_attainment <- relevel(factor(filtered$educational_attainment), ref="bachelor")
filtered$nominal_education <- relevel(factor(filtered$nominal_education), ref = "Bachelors")
filtered$education_level <- relevel(factor(filtered$education_level), ref = "High")
filtered$at_least_one_risk <- factor(filtered$at_least_one_risk)
filtered$risk_groups <- relevel(factor(filtered$risk_groups), ref="None")
filtered$categorical_doses <- relevel(factor(filtered$categorical_doses), ref="unvaccinated")
filtered$factor_doses <- relevel(factor(filtered$doses), ref=1)
filtered$nominal_vaccinal_schema <- relevel(factor(filtered$nominal_vaccinal_schema), ref="None")
filtered$nominal_covid_treatment <- relevel(factor(filtered$nominal_covid_treatment), ref="No")
filtered$education_level <- relevel(factor(filtered$education_level), ref = "High")
filtered$ndoses <- relevel(factor(filtered$doses), ref = "0")
filtered$last_year <- relevel(factor(filtered$last_year), ref = "2021")
filtered$last_vaccine_interval <- relevel(factor(filtered$last_vaccine_interval), ref = "over 1 year")
filtered$last_vaccine_interval_broader <- relevel(factor(filtered$last_vaccine_interval_broader), ref = "over 1 year")

Population characterization

In total, 474 immunocompromised older adults hospitalized for COVID-19 in Brazil during the 2024 epidemiologic year were included in the analysis.

Age

Median, min, max

print(median(filtered$NU_IDADE_N))

## [1] 75

print(min(filtered$NU_IDADE_N))

## [1] 60

print(max(filtered$NU_IDADE_N))

## [1] 101

Sex

describe.factor(filtered$CS_SEXO)

##                 
## filtered$CS_SEXO         F         M
##          Count   252.00000 222.00000
##          Percent  53.16456  46.83544

Race

describe.factor(filtered$nominal_race)

##                      
## filtered$nominal_race     White    Brown  Unknown     Black    Yellow
##               Count   311.00000 97.00000 52.00000 10.000000 4.0000000
##               Percent  65.61181 20.46414 10.97046  2.109705 0.8438819

Education

describe.factor(filtered$educational_attainment)

##                                
## filtered$educational_attainment   unknown       low  bachelor
##                         Count   261.00000 187.00000 26.000000
##                         Percent  55.06329  39.45148  5.485232

Percentage without bachelor degree considering only the rows with known data:

low <- nrow(filtered %>% filter(educational_attainment == "low"))
bachelor <- nrow(filtered %>% filter(educational_attainment == "bachelor"))
print(percent(low / (bachelor + low)))

## [1] "88%"

At least one risk condition

describe.factor(filtered$at_least_one_risk)

##                           
## filtered$at_least_one_risk         1         0
##                    Count   314.00000 160.00000
##                    Percent  66.24473  33.75527

Vaccine doses

describe.factor(filtered$dose_groups)

##                     
## filtered$dose_groups        3+        2        0        1
##              Count   318.00000 93.00000 56.00000 7.000000
##              Percent  67.08861 19.62025 11.81435 1.476793

Inactivate

describe.factor(filtered$innactivated_doses >= 1)

##                                 
## filtered$innactivated_doses >= 1      TRUE     FALSE
##                          Count   240.00000 234.00000
##                          Percent  50.63291  49.36709

Viral-vector

describe.factor(filtered$viral_vector_doses >= 1)

##                                 
## filtered$viral_vector_doses >= 1     TRUE    FALSE
##                          Count   247.0000 227.0000
##                          Percent  52.1097  47.8903

mRNA

describe.factor(filtered$mrna_doses >= 1)

##                         
## filtered$mrna_doses >= 1     FALSE      TRUE
##                  Count   275.00000 199.00000
##                  Percent  58.01688  41.98312

Vaccine interval

describe.factor(filtered$last_vaccine_interval)

##                               
## filtered$last_vaccine_interval over 1 year 7–12 months     <NA> up to 6 months
##                        Count     330.00000    79.00000 56.00000       9.000000
##                        Percent    69.62025    16.66667 11.81435       1.898734

describe.factor(filtered$last_vaccine_interval_broader)

##                                       
## filtered$last_vaccine_interval_broader over 1 year up to 1 year     <NA>
##                                Count     330.00000      88.0000 56.00000
##                                Percent    69.62025      18.5654 11.81435

Outcome

table(filtered$nominal_outcome)

## 
## Non-survivor     Survivor 
##          171          303

Bivariate analysis

Total

print(paste("Survivor:", nrow(survivors)))

## [1] "Survivor: 303"

print(paste("Non-Survivor:", nrow(nonsurvivors)))

## [1] "Non-Survivor: 171"

Epidemiologic week of the onset of first symptoms

Survivor

print(paste("Median:", median(survivors$SEM_PRI)))

## [1] "Median: 12"

print(paste("IQR:", IQR(survivors$SEM_PRI)))

## [1] "IQR: 29"

Check normality (Shapiro-Wilk test):

shapiro.test(survivors$SEM_PRI)

## 
##  Shapiro-Wilk normality test
## 
## data:  survivors$SEM_PRI
## W = 0.86589, p-value = 1.466e-15

p-value < 0.05 –> Non-normal

Non-Survivor:

print(paste("Median:", median(nonsurvivors$SEM_PRI)))

## [1] "Median: 14"

print(paste("IQR:", IQR(nonsurvivors$SEM_PRI)))

## [1] "IQR: 30.5"

Check normality (Shapiro-Wilk test):

shapiro.test(nonsurvivors$SEM_PRI)

## 
##  Shapiro-Wilk normality test
## 
## data:  nonsurvivors$SEM_PRI
## W = 0.86711, p-value = 3.812e-11

p-value < 0.05 –> Non-normal

Analysis (Mann-Whitney U test)

wilcox.test(survivors$SEM_PRI, nonsurvivors$SEM_PRI)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  survivors$SEM_PRI and nonsurvivors$SEM_PRI
## W = 24840, p-value = 0.4561
## alternative hypothesis: true location shift is not equal to 0

Raw Odds Ratio

m_raw <- glm(outcome ~ SEM_PRI,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "SEM_PRI" = "Epidemiologic week of the onset of first COVID-19 symptoms"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
epi_week_OR_raw = OR_raw[-1, ]
OR_raw

Days from onset of first symptoms

Survivor

print(paste("Median:", median(survivors$delta_days)))

## [1] "Median: 2"

print(paste("IQR:", IQR(survivors$delta_days)))

## [1] "IQR: 3"

Check normality (Shapiro-Wilk test):

shapiro.test(survivors$delta_days)

## 
##  Shapiro-Wilk normality test
## 
## data:  survivors$delta_days
## W = 0.58882, p-value < 2.2e-16

p-value < 0.05 –> Non-normal

Non-Survivor:

print(paste("Median:", median(nonsurvivors$delta_days)))

## [1] "Median: 2"

print(paste("IQR:", IQR(nonsurvivors$delta_days)))

## [1] "IQR: 5"

Check normality (Shapiro-Wilk test):

shapiro.test(nonsurvivors$delta_days)

## 
##  Shapiro-Wilk normality test
## 
## data:  nonsurvivors$delta_days
## W = 0.70536, p-value < 2.2e-16

p-value < 0.05 –> Non-normal

Analysis (Mann-Whitney U test)

wilcox.test(survivors$delta_days, nonsurvivors$delta_days)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  survivors$delta_days and nonsurvivors$delta_days
## W = 27800, p-value = 0.1799
## alternative hypothesis: true location shift is not equal to 0

Raw Odds Ratio

m_raw <- glm(outcome ~ delta_days,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "delta_days" = "Days from onset of first COVID-19 symptoms to hospitalization"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
delta_days_OR_raw = OR_raw[-1, ]
OR_raw

Age group, years, n

Description

table(filtered$age_groups, filtered$nominal_outcome)

##        
##         Non-survivor Survivor
##   60-69           61      103
##   70-79           53      101
##   80+             57       99

Analysis

chisq.test(filtered$age_groups, filtered$EVOLUCAO)

## 
##  Pearson's Chi-squared test
## 
## data:  filtered$age_groups and filtered$EVOLUCAO
## X-squared = 0.28764, df = 2, p-value = 0.866

Raw Odds Ratio

m_raw <- glm(outcome ~ age_groups,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "age_groups70-79" = "Age group: 70-79",
  "age_groups80+" = "Age group: 80+"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
age_groups_OR_raw = OR_raw[-1, ]
OR_raw

m_raw <- glm(outcome ~ NU_IDADE_N,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "NU_IDADE_N" = "Age (number)"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
age_OR_raw = OR_raw[-1, ]
OR_raw

Gender, n

Description

table(filtered$CS_SEXO, filtered$nominal_outcome)

##    
##     Non-survivor Survivor
##   F           96      156
##   M           75      147

Analysis

chisq.test(filtered$CS_SEXO, filtered$EVOLUCAO)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  filtered$CS_SEXO and filtered$EVOLUCAO
## X-squared = 0.77358, df = 1, p-value = 0.3791

Raw Odds Ratio

m_raw <- glm(outcome ~ CS_SEXO,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "CS_SEXOM" = "Male"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
gender_OR_raw = OR_raw[-1, ]
OR_raw

Race or Ethnicity, n

Description

table(filtered$nominal_race, filtered$nominal_outcome)

##          
##           Non-survivor Survivor
##   White            109      202
##   Black              3        7
##   Brown             43       54
##   Unknown           15       37
##   Yellow             1        3

Analysis

fisher.test(filtered$nominal_race, filtered$EVOLUCAO)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  filtered$nominal_race and filtered$EVOLUCAO
## p-value = 0.3382
## alternative hypothesis: two.sided

Raw Odds Ratio

m_raw <- glm(outcome ~ nominal_race,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "nominal_raceBlack" = "Race: Black",
  "nominal_raceBrown" = "Race: Brown",
  "nominal_raceYellow" = "Race: Yellow",
  "nominal_raceUnknown" = "Race: Unknown"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
race_OR_raw = OR_raw[-1, ]
OR_raw

Educational attainment, n

Description

table(filtered$nominal_education, filtered$nominal_outcome)

##             
##              Non-survivor Survivor
##   Bachelors             5       21
##   Elementary           44       50
##   Highschool           14       41
##   Middle               18       20
##   Unknown              90      171

Analysis

chisq.test(filtered$nominal_education, filtered$EVOLUCAO)

## 
##  Pearson's Chi-squared test
## 
## data:  filtered$nominal_education and filtered$EVOLUCAO
## X-squared = 12.973, df = 4, p-value = 0.01141

Raw Odds Ratio

m_raw <- glm(outcome ~ nominal_education,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "nominal_educationElementary" = "Education: Elementary (Ref: Bachelors)",
  "nominal_educationHighschool" = "Education: High School (Ref: Bachelors)",
  "nominal_educationMiddle" = "Education: Middle School (Ref:Bachelors)",
  "nominal_educationUnknown" = "Education: Unknown (Ref:Bachelors)"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
nominal_education_OR_raw = OR_raw[-1, ]
OR_raw

Educational level, n

Description

table(filtered$education_level, filtered$nominal_outcome)

##             
##              Non-survivor Survivor
##   High                 19       62
##   Elementary           44       50
##   Middle               18       20
##   Unknown              90      171

Analysis

chisq.test(filtered$education_level, filtered$EVOLUCAO)

## 
##  Pearson's Chi-squared test
## 
## data:  filtered$education_level and filtered$EVOLUCAO
## X-squared = 12.677, df = 3, p-value = 0.00539

Raw Odds Ratio

m_raw <- glm(outcome ~ education_level,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "education_levelElementary" = "Education: Elementary (Ref: High school or Higher)",
  "education_levelMiddle" = "Education: Middle (Ref: High school or Higher)",
  "education_levelUnknown" = "Education: Unknown (Ref: High school or Higher)"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
education_level_OR_raw = OR_raw[-1, ]
OR_raw

Other underlying risk conditions, n

Description

table(filtered$at_least_one_risk, filtered$nominal_outcome)

##    
##     Non-survivor Survivor
##   0           64       96
##   1          107      207

Analysis

fisher.test(filtered$at_least_one_risk, filtered$EVOLUCAO)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  filtered$at_least_one_risk and filtered$EVOLUCAO
## p-value = 0.2252
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.5141053 1.1730442
## sample estimates:
## odds ratio 
##  0.7757992

Raw Odds Ratio

m_raw <- glm(outcome ~ at_least_one_risk,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "at_least_one_risk1" = "Other underlying risk conditions"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
at_least_one_risk_OR_raw = OR_raw[-1, ]
OR_raw

Intensive care unit admission, n

Description

table(filtered$nominal_icu_admission, filtered$nominal_outcome)

##      
##       Non-survivor Survivor
##   No            87      221
##   Yes           84       82

Analysis

chisq.test(filtered$nominal_icu_admission, filtered$EVOLUCAO)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  filtered$nominal_icu_admission and filtered$EVOLUCAO
## X-squared = 22.417, df = 1, p-value = 2.194e-06

Raw Odds Ratio

m_raw <- glm(outcome ~ nominal_icu_admission,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "nominal_icu_admissionYes" = "Intensive care unit admission"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
nominal_icu_admission_OR_raw = OR_raw[-1, ]
OR_raw

Invasive ventilatory support, n

Description

table(filtered$nominal_invasive_support_ven, filtered$nominal_outcome)

##      
##       Non-survivor Survivor
##   No           105      290
##   Yes           66       13

Analysis

chisq.test(filtered$nominal_invasive_support_ven, filtered$EVOLUCAO)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  filtered$nominal_invasive_support_ven and filtered$EVOLUCAO
## X-squared = 90.173, df = 1, p-value < 2.2e-16

Raw Odds Ratio

m_raw <- glm(outcome ~ nominal_invasive_support_ven,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "nominal_invasive_support_venYes" = "Invasive ventilatory support"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
nominal_invasive_support_ven_OR_raw = OR_raw[-1, ]
OR_raw

COVID-19 vaccine doses, n

Description

table(filtered$doses, filtered$nominal_outcome)

##    
##     Non-survivor Survivor
##   0           24       32
##   1            1        6
##   2           38       55
##   3           32       63
##   4           40       72
##   5           35       70
##   6            1        5

Analysis (Cochran-Armitage Test)

CochranArmitageTest(table(filtered$nominal_outcome, filtered$doses), alternative="one.sided")

## 
##  Cochran-Armitage test for trend
## 
## data:  table(filtered$nominal_outcome, filtered$doses)
## Z = -1.3483, dim = 7, p-value = 0.08878
## alternative hypothesis: one.sided

Raw Odds Ratio

m_raw <- glm(outcome ~ doses,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "doses" = "Doses (number)"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
doses_OR_raw = OR_raw[-1, ]
OR_raw

m_raw <- glm(outcome ~ ndoses,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "ndoses1" = "1 Dose (Ref: Unvaccinated)",
  "ndoses2" = "2 Doses (Ref: Unvaccinated)",
  "ndoses3" = "3 Doses (Ref: Unvaccinated)",
  "ndoses4" = "4 Doses (Ref: Unvaccinated)",
  "ndoses5" = "5 Doses (Ref: Unvaccinated)",
  "ndoses6" = "6 Doses (Ref: Unvaccinated)"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
ndoses_OR_raw = OR_raw[-1, ]
OR_raw

Categorical

table(filtered$categorical_doses, filtered$nominal_outcome)

##               
##                Non-survivor Survivor
##   unvaccinated           24       32
##   1-2                    39       61
##   3                      32       63
##   4                      40       72
##   5+                     36       75

COVID-19 vaccine types, n

Description

table(filtered$nominal_vaccinal_schema, filtered$nominal_outcome)

##                                    
##                                     Non-survivor Survivor
##   None                                        25       33
##   Inactivated + mRNA                          19       50
##   Inactivated + mRNA + viral vector           19       32
##   Inactivated + Viral vector                  15       18
##   Inactivated only                            27       60
##   mRNA + viral vector                         22       44
##   mRNA only                                    3       10
##   Viral vector only                           41       56

Analysis

fisher.test(filtered$nominal_vaccinal_schema, filtered$nominal_outcome, workspace = 20000000)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  filtered$nominal_vaccinal_schema and filtered$nominal_outcome
## p-value = 0.3018
## alternative hypothesis: two.sided

Raw Odds Ratio

m_raw <- glm(outcome ~ nominal_vaccinal_schema,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "nominal_vaccinal_schemaInactivated + mRNA" = "Inactivated + mRNA vaccines (Ref: Unvaccinated)",
  "nominal_vaccinal_schemaInactivated + mRNA + viral vector" = "Inactivated + mRNA + viral vector (Ref: Unvaccinated)",
  "nominal_vaccinal_schemaInactivated + Viral vector" = "Inactivated + viral vector vaccines (Ref: Unvaccinated)",
  "nominal_vaccinal_schemaInactivated only" = "Inactivated vaccine only (Ref: Unvaccinated)",
  "nominal_vaccinal_schemamRNA + viral vector" = "mRNA + viral vector vaccines (Ref: Unvaccinated)",
  "nominal_vaccinal_schemamRNA only" = "mRNA vaccine only (Ref: Unvaccinated)",
  "nominal_vaccinal_schemaViral vector only" = "Viral vector vaccine only (Ref: Unvaccinated)"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
nominal_vaccinal_schema_OR_raw = OR_raw[-1, ]
OR_raw

Interval since last COVID-19 vaccine dose, n

Description

table(filtered$last_vaccine_interval, filtered$nominal_outcome)

##                 
##                  Non-survivor Survivor
##   over 1 year             121      209
##   7–12 months              25       54
##   up to 6 months            1        8

Broader:

table(filtered$last_vaccine_interval_broader, filtered$nominal_outcome)

##               
##                Non-survivor Survivor
##   over 1 year           121      209
##   up to 1 year           26       62

Analysis

fisher.test(filtered$last_vaccine_interval, filtered$nominal_outcome)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  filtered$last_vaccine_interval and filtered$nominal_outcome
## p-value = 0.2311
## alternative hypothesis: two.sided

m_raw <- glm(outcome ~ last_vaccine_interval,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "last_vaccine_intervalup to 6 months" = "Last vaccine interval: ≤ 6 months (Ref: > 1 year)",
  "last_vaccine_interval7–12 months" = "Last vaccine interval: 7-12 months (Ref: > 1 year)"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
last_vaccine_interval_OR_raw = OR_raw[-1, ]
OR_raw

Broader:

fisher.test(filtered$last_vaccine_interval_broader, filtered$nominal_outcome)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  filtered$last_vaccine_interval_broader and filtered$nominal_outcome
## p-value = 0.2582
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.8099475 2.3996549
## sample estimates:
## odds ratio 
##   1.379529

m_raw <- glm(outcome ~ last_vaccine_interval_broader,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "last_vaccine_interval_broaderup to 1 year" = "Last vaccine interval: ≤ 1 year (Ref: > 1 year)"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
last_vaccine_interval_broader_OR_raw = OR_raw[-1, ]
OR_raw

Year of last COVID-19 vaccine dose, n

Description

table(filtered$last_year, filtered$nominal_outcome)

##       
##        Non-survivor Survivor
##   2021           52      102
##   2022           49       61
##   2023           45      102
##   2024            1        6

Analysis

fisher.test(filtered$last_year, filtered$nominal_outcome)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  filtered$last_year and filtered$nominal_outcome
## p-value = 0.07572
## alternative hypothesis: two.sided

m_raw <- glm(outcome ~ last_year,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "last_year2022" = "Year of last COVID-19 vaccine dose: 2022 (Ref: 2021)",
  "last_year2023" = "Year of last COVID-19 vaccine dose: 2023 (Ref: 2021)",
  "last_year2024" = "Year of last COVID-19 vaccine dose: 2024 (Ref: 2021)"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
last_year_OR_raw = OR_raw[-1, ]
OR_raw

COVID-19 specific treatment, n

Description

table(filtered$nominal_covid_treatment, filtered$nominal_outcome)

##      
##       Non-survivor Survivor
##   No           167      290
##   Yes            4       13

Analysis

chisq.test(filtered$nominal_covid_treatment, filtered$nominal_outcome)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  filtered$nominal_covid_treatment and filtered$nominal_outcome
## X-squared = 0.70543, df = 1, p-value = 0.401

Raw Odds Ratio

m_raw <- glm(outcome ~ nominal_covid_treatment,
               data = filtered,
               family = binomial)

# OR bruto + IC95%
OR_raw <- as.data.frame(confint(m_raw))

## Waiting for profiling to be done...

OR_raw$Estimate <- coef(m_raw)
OR_raw$OR <- exp(OR_raw$Estimate)
OR_raw$LowerCI <- exp(OR_raw[, 1])
OR_raw$UpperCI <- exp(OR_raw[, 2])
label_map <- c(
  "nominal_covid_treatmentYes" = "COVID-19 specific treatment"
)
OR_raw$Variable <- label_map[rownames(OR_raw)]
nominal_covid_treatment_OR_raw = OR_raw[-1, ]
OR_raw

Raw OR Forest Plot

coef_df <- rbind(
  epi_week_OR_raw, delta_days_OR_raw, age_OR_raw, age_groups_OR_raw,
  gender_OR_raw, race_OR_raw, education_level_OR_raw, # nominal_education_OR_raw,
  at_least_one_risk_OR_raw, nominal_icu_admission_OR_raw,
  nominal_invasive_support_ven_OR_raw, doses_OR_raw, ndoses_OR_raw,
  nominal_vaccinal_schema_OR_raw, last_year_OR_raw,
  nominal_covid_treatment_OR_raw
)
coef_df$Variable <- factor(coef_df$Variable, levels=rev(coef_df$Variable))

p = ggplot(coef_df, aes(x = Variable, y = OR)) +  # remove Intercept
  geom_point() +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  coord_flip() +
  labs(#title = "Forest Plot of Logistic Regression",
       x = "Variables",
       y = "Raw Odds Ratio (95% CI)") +
  theme_minimal()

ggsave("raw_forest_plot.pdf", p, width = 8, height = 11)
p

Time betwwen vaccination and outcome

`

min(filtered$days_since_last_vaccine , na.rm = TRUE)

## [1] 63

filtered$days_since_last_vaccine

##   [1]  528  222  774   NA   NA   NA  193 1106  327  854  246   NA  641  556  848
##  [16]  295  317  522  683  768 1056   NA  828  832 1290   NA  871   NA  654   NA
##  [31]  761 1106  126  405  950   NA 1169 1215  609   NA  377 1042  279   NA  740
##  [46]  564  507  542  799 1069  290  323  194 1270  592  704  918   NA   NA  896
##  [61]  906  800 1221  802  668  354  366  628  665  797 1259  789  871  805  875
##  [76]  389 1218  841   NA  301  312  295   NA  623   NA  987  340  606 1084  796
##  [91]  964  951 1140   NA  654  278 1164  861   NA  860  298 1123  872  798  323
## [106]  506   NA  723 1033  886  586 1290   NA   NA  893  255 1016  932   NA  301
## [121]  666  575  482  656   63 1097 1143  244  320  622  854   NA  365  434   NA
## [136]  740  327  726  590  816  614  945 1038 1179  575  491  340 1055 1198  992
## [151]  320   NA  265 1068  337   69   NA 1102 1061 1049  600  536 1296   NA  280
## [166]  550  625  909  850  302  917  630  360  529  901  188 1390  427  428 1058
## [181]  540   NA  296 1231  277  888  319 1056 1106  311  941  241   NA 1051 1091
## [196]   NA  780  877  280  909  449 1203 1128  855  882  791 1108  719   NA  605
## [211]  893 1143 1262  567 1141  589 1222  345  780   NA  500  303  918  820  647
## [226]  680  188  519  373  865  991  686 1012  825  594 1007  347  605  529 1308
## [241] 1027  250  357  358 1131  529  542 1013  548 1150 1271   NA  357 1227  605
## [256]  258  331  805  826  627   NA 1231  256  478  888  513  976  326  540  476
## [271]  331  776  364  918 1054  581   NA  654  319  381  338   NA  155  644  282
## [286] 1067  559  853 1138  144  364  273 1142  628  616  912  514  815 1285  557
## [301]  903  682  751  526  968  414  356  847  600   NA  893  340 1141 1057  881
## [316]  865  113 1089  553  380  552 1148 1074  302 1115  554 1109   NA  501   NA
## [331]   NA  314  766   NA  274  731  812 1036 1063 1031  332  916  714  985  363
## [346]  337  865 1204  822 1254  480  991  871  599  809  702  645   NA  596  418
## [361]  944  683   NA 1075  682 1144  953   NA  619  994  936 1065 1344  213  686
## [376]  865  777  657 1152  654  296  704  835  836  515 1048   NA  923 1018  354
## [391]  837  624  711  798   NA  618  586  730  629  636  476  332  165 1134 1311
## [406]   NA   NA  615  520 1026 1173 1196 1038  599  340  580   NA   89 1047   NA
## [421]   NA  885  414  598 1177  633  697  765  136  873  849  526  366   NA  209
## [436] 1232  843 1194  431  643  581 1253  262  377 1183   NA  834  823 1110  296
## [451]  315  282  994  786  810  585  996  326 1277   NA  497  491  886 1293   NA
## [466] 1036   NA  333  353  964  388  324  584  782

breaks <- seq(
  0,
  max(filtered$days_since_last_vaccine, na.rm = TRUE) + 30,
  by = 30
)

labels <- paste0(
  head(breaks, -1) / 30
)

filtered %>%
  filter(
    !is.na(days_since_last_vaccine),
    !is.na(outcome)
  ) %>%
  ggplot(
    aes(
      x = cut(
        days_since_last_vaccine,
        breaks = breaks,
        labels = labels,
        include.lowest = TRUE
      ),
      fill = factor(outcome)
    )
  ) +
  #geom_bar(position = "fill", na.rm = TRUE) +
  geom_bar(na.rm = TRUE) +
  geom_text(
    stat = "count",
    aes(
      y = after_stat(count), 
      label = after_stat(count),
      group = factor(outcome)
    ),
    #position = position_fill(vjust = 0.5),
    position = position_stack(vjust= 0.5),
    color = "black",
    size = 3,
    na.rm = TRUE
  ) +
  scale_fill_manual(
    values = c("0" = "#faa", "1" = "#afa"),
    labels = c("0" = "Non-survivor", "1" = "Survivor")
  ) +
  labs(
    x = "Months since last vaccine",
    y = "Count",
    fill = "Outcome"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

Model

Main Version - Variables from paper

Category: variable (Ref: value): groups…

• Age: covariate NU_IDADE_N
• Gender: CS_SEXO (Ref: F): M
• Educational Attainment: education_level (Ref: Elementary): Middle, High, Unknown
• ICU: nominal_icu_admission (Ref: No): Yes
• Invasive Support Ven: nominal_invasive_support_ven (Ref: No): Yes
• Doses: categorical_doses (Ref: unvaccinated): 1-2, 3, 4, 5+

Logistic Regression

model <- glm(outcome ~ NU_IDADE_N + CS_SEXO +
               education_level +
               nominal_icu_admission + 
               nominal_invasive_support_ven + categorical_doses,
               family = binomial,
               data = filtered)

summary_model <- summary(model)
print(summary_model)

## 
## Call:
## glm(formula = outcome ~ NU_IDADE_N + CS_SEXO + education_level + 
##     nominal_icu_admission + nominal_invasive_support_ven + categorical_doses, 
##     family = binomial, data = filtered)
## 
## Coefficients:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     -2.040698   0.932136  -2.189  0.02858 *  
## NU_IDADE_N                       0.009214   0.011312   0.815  0.41532    
## CS_SEXOM                         0.069242   0.219842   0.315  0.75279    
## education_levelElementary        1.078364   0.389687   2.767  0.00565 ** 
## education_levelMiddle            1.060408   0.482131   2.199  0.02785 *  
## education_levelUnknown           0.563076   0.340037   1.656  0.09774 .  
## nominal_icu_admissionYes         0.365988   0.254819   1.436  0.15093    
## nominal_invasive_support_venYes  2.529359   0.359145   7.043 1.89e-12 ***
## categorical_doses1-2            -0.407902   0.373804  -1.091  0.27518    
## categorical_doses3              -0.402041   0.383953  -1.047  0.29505    
## categorical_doses4              -0.607060   0.381589  -1.591  0.11164    
## categorical_doses5+             -0.586755   0.378867  -1.549  0.12145    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 619.85  on 473  degrees of freedom
## Residual deviance: 512.16  on 462  degrees of freedom
## AIC: 536.16
## 
## Number of Fisher Scoring iterations: 4

odds_ratios <- exp(coef(model))
conf_int <- exp(confint(model))

## Waiting for profiling to be done...

data.frame(
  variable = names(coef(model)),
  OR = round(odds_ratios, 2),
  IC_95_inf = round(conf_int[,1], 2),
  IC_95_sup = round(conf_int[,2], 2),
  pvalue = round(summary_model$coefficients[,4], 3)
)

Variance inflation factor

vif_results <- car::vif(model)
print(vif_results)

##                                  GVIF Df GVIF^(1/(2*Df))
## NU_IDADE_N                   1.068068  1        1.033474
## CS_SEXO                      1.023759  1        1.011810
## education_level              1.149977  3        1.023564
## nominal_icu_admission        1.254092  1        1.119863
## nominal_invasive_support_ven 1.201453  1        1.096108
## categorical_doses            1.110703  4        1.013211

No multicollinearity detected.

Forest Plot

label_map <- c(
  "NU_IDADE_N" = "Age",
  "CS_SEXOM" = "Male",
  "education_levelElementary" = "Elementary Education",
  "education_levelMiddle" = "Middle Education",
  "education_levelUnknown" = "Unknown Education",
  "nominal_icu_admissionYes" = "ICU Admission",
  "nominal_invasive_support_venYes" = "Invasive Ventilation",
  "categorical_doses1-2" = "1-2 Doses",
  "categorical_doses3" = "3 Doses",
  "categorical_doses4" = "4 Doses",
  "categorical_doses5+" = "5+ Doses"
)


coef_df <- as.data.frame(confint(model))  # 95% CI

## Waiting for profiling to be done...

coef_df$Estimate <- coef(model)
coef_df$OR <- exp(coef_df$Estimate)
coef_df$LowerCI <- exp(coef_df[, 1])
coef_df$UpperCI <- exp(coef_df[, 2])
coef_df$Variable <- label_map[rownames(coef_df)]
coef_df= coef_df[row.names(coef_df) != "education_levelUnknown",]

p = ggplot(coef_df[-1, ], aes(x = reorder(Variable, OR), y = OR)) +  # remove Intercept
  geom_point() +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  coord_flip() +
  labs(title = "Forest Plot of Logistic Regression",
       x = "Variables",
       y = "Odds Ratio (95% CI)") +
  theme_minimal()

ggsave("forest_plot.pdf", p, width = 8, height = nrow(coef_df) * 0.4)
p

Hosmer-Lemeshow Test

hl_test <- hoslem.test(model$y, fitted(model), g = 10)
print(hl_test)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model$y, fitted(model)
## X-squared = 5.0032, df = 8, p-value = 0.7572

p-value > 0.05 -> no evidence of poor fit

ROC and AUC

roc_obj <- roc(response = filtered$outcome,
               predictor = predict(model, type = "response"))

## Setting levels: control = 1, case = 0

## Setting direction: controls < cases

auc(roc_obj)

## Area under the curve: 0.7398

ci.auc(roc_obj)

## 95% CI: 0.6915-0.788 (DeLong)

AUC between 0.7 and 0.8: acceptable ability to discriminate between classes Entire confidence interval above 0.5: model performs better than random

plot(roc_obj, main = "ROC Curve")

Alternative versions

Alternative A - Main version, but it considers the interval since last vaccine

Category: variable (Ref: value): groups…

• Age: covariate NU_IDADE_N
• Gender: CS_SEXO (Ref: F): M
• Educational Attainment: education_level (Ref: Elementary): Middle, High, Unknown
• ICU: nominal_icu_admission (Ref: No): Yes
• Invasive Support Ven: nominal_invasive_support_ven (Ref: No): Yes
• Doses: categorical_doses (Ref: unvaccinated): 1-2, 3, 4, 5+
• Interval: last_vaccine_interval_broader (Ref: >= 1 year): < 1 year

Logistic Regression

model <- glm(outcome ~ NU_IDADE_N + CS_SEXO +
               education_level +
               nominal_icu_admission + 
               nominal_invasive_support_ven +
               categorical_doses +
               last_vaccine_interval_broader,
               family = binomial,
               data = filtered)

summary_model <- summary(model)
print(summary_model)

## 
## Call:
## glm(formula = outcome ~ NU_IDADE_N + CS_SEXO + education_level + 
##     nominal_icu_admission + nominal_invasive_support_ven + categorical_doses + 
##     last_vaccine_interval_broader, family = binomial, data = filtered)
## 
## Coefficients:
##                                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                               -2.680326   0.977639  -2.742  0.00611
## NU_IDADE_N                                 0.009431   0.012039   0.783  0.43340
## CS_SEXOM                                   0.028277   0.237368   0.119  0.90517
## education_levelElementary                  1.346295   0.430485   3.127  0.00176
## education_levelMiddle                      1.274106   0.553342   2.303  0.02130
## education_levelUnknown                     0.844346   0.375416   2.249  0.02451
## nominal_icu_admissionYes                   0.343612   0.268977   1.277  0.20143
## nominal_invasive_support_venYes            2.459607   0.371825   6.615 3.72e-11
## categorical_doses3                         0.039423   0.338643   0.116  0.90733
## categorical_doses4                        -0.145793   0.334708  -0.436  0.66314
## categorical_doses5+                       -0.018990   0.390797  -0.049  0.96124
## last_vaccine_interval_broaderup to 1 year -0.242274   0.378783  -0.640  0.52243
##                                              
## (Intercept)                               ** 
## NU_IDADE_N                                   
## CS_SEXOM                                     
## education_levelElementary                 ** 
## education_levelMiddle                     *  
## education_levelUnknown                    *  
## nominal_icu_admissionYes                     
## nominal_invasive_support_venYes           ***
## categorical_doses3                           
## categorical_doses4                           
## categorical_doses5+                          
## last_vaccine_interval_broaderup to 1 year    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 542.13  on 417  degrees of freedom
## Residual deviance: 443.88  on 406  degrees of freedom
##   (56 observations deleted due to missingness)
## AIC: 467.88
## 
## Number of Fisher Scoring iterations: 4

xtabs(~ categorical_doses + last_vaccine_interval_broader, data = filtered)

##                  last_vaccine_interval_broader
## categorical_doses over 1 year up to 1 year
##      unvaccinated           0            0
##      1-2                   99            1
##      3                     94            1
##      4                     94           18
##      5+                    43           68

odds_ratios <- exp(coef(model))
conf_int <- exp(confint(model))

## Waiting for profiling to be done...

data.frame(
  variable = names(coef(model)),
  OR = round(odds_ratios, 2),
  IC_95_inf = round(conf_int[,1], 2),
  IC_95_sup = round(conf_int[,2], 2),
  pvalue = round(summary_model$coefficients[,4], 3)
)

Variance inflation factor

vif_results <- car::vif(model)
print(vif_results)

##                                   GVIF Df GVIF^(1/(2*Df))
## NU_IDADE_N                    1.079603  1        1.039040
## CS_SEXO                       1.028971  1        1.014382
## education_level               1.158813  3        1.024870
## nominal_icu_admission         1.249966  1        1.118019
## nominal_invasive_support_ven  1.225559  1        1.107050
## categorical_doses             1.649277  3        1.086965
## last_vaccine_interval_broader 1.577699  1        1.256065

No multicollinearity detected.

Forest Plot

label_map <- c(
  "NU_IDADE_N" = "Age",
  "CS_SEXOM" = "Male",
  "education_levelElementary" = "Elementary Education",
  "education_levelMiddle" = "Middle Education",
  "education_levelUnknown" = "Unknown Education",
  "nominal_icu_admissionYes" = "ICU Admission",
  "nominal_invasive_support_venYes" = "Invasive Ventilation",
  "categorical_doses1-2" = "1-2 Doses",
  "categorical_doses3" = "3 Doses",
  "categorical_doses4" = "4 Doses",
  "categorical_doses5+" = "5+ Doses",
  "last_vaccine_interval_broaderup to 1 year" = "Less than 1 year since last"
)


coef_df <- as.data.frame(confint(model))  # 95% CI

## Waiting for profiling to be done...

coef_df$Estimate <- coef(model)
coef_df$OR <- exp(coef_df$Estimate)
coef_df$LowerCI <- exp(coef_df[, 1])
coef_df$UpperCI <- exp(coef_df[, 2])
coef_df$Variable <- label_map[rownames(coef_df)]
coef_df= coef_df[row.names(coef_df) != "education_levelUnknown",]

p = ggplot(coef_df[-1, ], aes(x = reorder(Variable, OR), y = OR)) +  # remove Intercept
  geom_point() +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  coord_flip() +
  labs(title = "Forest Plot of Logistic Regression",
       x = "Variables",
       y = "Odds Ratio (95% CI)") +
  theme_minimal()

#ggsave("forest_plot.pdf", p, width = 8, height = nrow(coef_df) * 0.4)
p

Hosmer-Lemeshow Test

hl_test <- hoslem.test(model$y, fitted(model), g = 10)
print(hl_test)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model$y, fitted(model)
## X-squared = 8.7937, df = 8, p-value = 0.36

p-value > 0.05 -> no evidence of poor fit

ROC and AUC

length(predict(model))

## [1] 418

mf <- model.frame(model)

roc_obj <- roc(
  response  = mf$outcome,
  predictor = predict(model, type = "response")
)

## Setting levels: control = 1, case = 0

## Setting direction: controls < cases

auc(roc_obj)

## Area under the curve: 0.7567

ci.auc(roc_obj)

## 95% CI: 0.7063-0.8071 (DeLong)

AUC between 0.7 and 0.8: acceptable ability to discriminate between classes Entire confidence interval above 0.5: model performs better than random

plot(roc_obj, main = "ROC Curve")

Alternative B - Main version, but splits high education level into Highschool and Bachelors

Category: variable (Ref: value): groups…

• Age: covariate NU_IDADE_N
• Gender: CS_SEXO (Ref: F): M
• Educational Attainment: nominal_education (Ref: Bachelor): Elementary, Middle, Highschool, Unknown
• ICU: nominal_icu_admission (Ref: No): Yes
• Invasive Support Ven: nominal_invasive_support_ven (Ref: No): Yes
• Doses: categorical_doses (Ref: unvaccinated): 1-2, 3, 4, 5+

Logistic Regression

model <- glm(outcome ~ NU_IDADE_N + CS_SEXO +
               nominal_education +
               nominal_icu_admission + 
               nominal_invasive_support_ven + categorical_doses,
               family = binomial,
               data = filtered)

summary_model <- summary(model)
print(summary_model)

## 
## Call:
## glm(formula = outcome ~ NU_IDADE_N + CS_SEXO + nominal_education + 
##     nominal_icu_admission + nominal_invasive_support_ven + categorical_doses, 
##     family = binomial, data = filtered)
## 
## Coefficients:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     -2.590499   1.043181  -2.483   0.0130 *  
## NU_IDADE_N                       0.008608   0.011283   0.763   0.4455    
## CS_SEXOM                         0.058446   0.220281   0.265   0.7908    
## nominal_educationElementary      1.678701   0.653227   2.570   0.0102 *  
## nominal_educationHighschool      0.841233   0.697888   1.205   0.2280    
## nominal_educationMiddle          1.656905   0.709355   2.336   0.0195 *  
## nominal_educationUnknown         1.161896   0.622682   1.866   0.0620 .  
## nominal_icu_admissionYes         0.370002   0.255073   1.451   0.1469    
## nominal_invasive_support_venYes  2.564512   0.363684   7.051 1.77e-12 ***
## categorical_doses1-2            -0.401908   0.373959  -1.075   0.2825    
## categorical_doses3              -0.415682   0.383667  -1.083   0.2786    
## categorical_doses4              -0.613810   0.381905  -1.607   0.1080    
## categorical_doses5+             -0.592463   0.378707  -1.564   0.1177    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 619.85  on 473  degrees of freedom
## Residual deviance: 510.62  on 461  degrees of freedom
## AIC: 536.62
## 
## Number of Fisher Scoring iterations: 4

odds_ratios <- exp(coef(model))
conf_int <- exp(confint(model))

## Waiting for profiling to be done...

data.frame(
  variable = names(coef(model)),
  OR = round(odds_ratios, 2),
  IC_95_inf = round(conf_int[,1], 2),
  IC_95_sup = round(conf_int[,2], 2),
  pvalue = round(summary_model$coefficients[,4], 3)
)

Variance inflation factor

vif_results <- car::vif(model)
print(vif_results)

##                                  GVIF Df GVIF^(1/(2*Df))
## NU_IDADE_N                   1.070363  1        1.034584
## CS_SEXO                      1.025048  1        1.012447
## nominal_education            1.187828  4        1.021749
## nominal_icu_admission        1.253157  1        1.119445
## nominal_invasive_support_ven 1.217552  1        1.103427
## categorical_doses            1.114846  4        1.013682

No multicollinearity detected.

Forest Plot

label_map <- c(
  "NU_IDADE_N" = "Age",
  "CS_SEXOM" = "Male",
  "nominal_educationElementary" = "Elementary Education",
  "nominal_educationHighschool" = "Highschool Education",
  "nominal_educationMiddle" = "Middle Education",
  "nominal_educationUnknown" = "Unknown Education",
  "nominal_icu_admissionYes" = "ICU Admission",
  "nominal_invasive_support_venYes" = "Invasive Ventilation",
  "categorical_doses1-2" = "1-2 Doses",
  "categorical_doses3" = "3 Doses",
  "categorical_doses4" = "4 Doses",
  "categorical_doses5+" = "5+ Doses"
)


coef_df <- as.data.frame(confint(model))  # 95% CI

## Waiting for profiling to be done...

coef_df$Estimate <- coef(model)
coef_df$OR <- exp(coef_df$Estimate)
coef_df$LowerCI <- exp(coef_df[, 1])
coef_df$UpperCI <- exp(coef_df[, 2])
coef_df$Variable <- label_map[rownames(coef_df)]
coef_df= coef_df[row.names(coef_df) != "nominal_educationUnknown",]

ggplot(coef_df[-1, ], aes(x = reorder(Variable, OR), y = OR)) +  # remove Intercept
  geom_point() +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  coord_flip() +
  labs(title = "Forest Plot of Logistic Regression",
       x = "Variables",
       y = "Odds Ratio (95% CI)") +
  theme_minimal()

Hosmer-Lemeshow Test

hl_test <- hoslem.test(model$y, fitted(model), g = 10)
print(hl_test)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model$y, fitted(model)
## X-squared = 7.0241, df = 8, p-value = 0.534

p-value > 0.05 -> no evidence of poor fit

ROC and AUC

roc_obj <- roc(response = filtered$outcome,
               predictor = predict(model, type = "response"))

## Setting levels: control = 1, case = 0

## Setting direction: controls < cases

auc(roc_obj)

## Area under the curve: 0.7431

ci.auc(roc_obj)

## 95% CI: 0.6953-0.7909 (DeLong)

AUC between 0.7 and 0.8: acceptable ability to discriminate between classes Entire confidence interval above 0.5: model performs better than random

plot(roc_obj, main = "ROC Curve")

Alternative C - Alternative B, but using groups from bivariate analyses

Category: variable (Ref: ): groups…

• Age: age_groups (Ref: 60–69): 70-79, 80+
• Gender: CS_SEXO (Ref: F): M
• Educational Attainment: nominal_education (Ref: Bachelor): Elementary, Middle, Highschool, Unknown
• ICU: nominal_icu_admission (Ref: No): Yes
• Invasive Support Ven: nominal_invasive_support_ven (Ref: No): Yes=
• Doses: factor_doses (Ref: 0): 1, 2, 3, 4, 5

Logistic Regression

model <- glm(outcome ~ age_groups + CS_SEXO +
               nominal_education +
               nominal_icu_admission + 
               nominal_invasive_support_ven + factor_doses,
               family = binomial,
               data = filtered)

summary_model <- summary(model)
print(summary_model)

## 
## Call:
## glm(formula = outcome ~ age_groups + CS_SEXO + nominal_education + 
##     nominal_icu_admission + nominal_invasive_support_ven + factor_doses, 
##     family = binomial, data = filtered)
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     -1.94748    0.70280  -2.771  0.00559 ** 
## age_groups70-79                 -0.13463    0.27411  -0.491  0.62333    
## age_groups80+                    0.07620    0.26975   0.282  0.77757    
## CS_SEXOM                         0.07285    0.22132   0.329  0.74202    
## nominal_educationElementary      1.71986    0.66664   2.580  0.00988 ** 
## nominal_educationHighschool      0.87045    0.71065   1.225  0.22063    
## nominal_educationMiddle          1.64160    0.72248   2.272  0.02308 *  
## nominal_educationUnknown         1.18474    0.63559   1.864  0.06232 .  
## nominal_icu_admissionYes         0.37425    0.25539   1.465  0.14280    
## nominal_invasive_support_venYes  2.57636    0.36545   7.050 1.79e-12 ***
## factor_doses1                   -1.82214    1.28285  -1.420  0.15550    
## factor_doses2                   -0.34576    0.37878  -0.913  0.36133    
## factor_doses3                   -0.43667    0.38561  -1.132  0.25746    
## factor_doses4                   -0.63070    0.38391  -1.643  0.10041    
## factor_doses5                   -0.54596    0.38261  -1.427  0.15360    
## factor_doses6                   -1.74952    1.32179  -1.324  0.18564    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 619.85  on 473  degrees of freedom
## Residual deviance: 507.74  on 458  degrees of freedom
## AIC: 539.74
## 
## Number of Fisher Scoring iterations: 4

odds_ratios <- exp(coef(model))
conf_int <- exp(confint(model))

## Waiting for profiling to be done...

data.frame(
  variable = names(coef(model)),
  OR = round(odds_ratios, 2),
  IC_95_inf = round(conf_int[,1], 2),
  IC_95_sup = round(conf_int[,2], 2),
  pvalue = round(summary_model$coefficients[,4], 3)
)

Variance inflation factor

vif_results <- car::vif(model)
print(vif_results)

##                                  GVIF Df GVIF^(1/(2*Df))
## age_groups                   1.084099  2        1.020392
## CS_SEXO                      1.029407  1        1.014597
## nominal_education            1.199530  4        1.023002
## nominal_icu_admission        1.254340  1        1.119973
## nominal_invasive_support_ven 1.218420  1        1.103821
## factor_doses                 1.158127  6        1.012309

No multicollinearity detected.

Forest Plot

label_map <- c(
  "age_groups70-79" = "Age 70-79",
  "age_groups80+" = "Age 80+",
  "CS_SEXOM" = "Male",
  "nominal_educationElementary" = "Elementary Education",
  "nominal_educationHighschool" = "Highschool Education",
  "nominal_educationMiddle" = "Middle Education",
  "nominal_educationUnknown" = "Unknown Education",
  "nominal_icu_admissionYes" = "ICU Admission",
  "nominal_invasive_support_venYes" = "Invasive Ventilation",
  "factor_doses1" = "1 Dose",
  "factor_doses2" = "2 Doses",
  "factor_doses3" = "3 Doses",
  "factor_doses4" = "4 Doses",
  "factor_doses5" = "5 Doses"
)


coef_df <- as.data.frame(confint(model))  # 95% CI

## Waiting for profiling to be done...

coef_df$Estimate <- coef(model)
coef_df$OR <- exp(coef_df$Estimate)
coef_df$LowerCI <- exp(coef_df[, 1])
coef_df$UpperCI <- exp(coef_df[, 2])
coef_df$Variable <- label_map[rownames(coef_df)]
coef_df= coef_df[row.names(coef_df) != "nominal_educationUnknown",]

ggplot(coef_df[-1, ], aes(x = reorder(Variable, OR), y = OR)) +  # remove Intercept
  geom_point() +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  coord_flip() +
  labs(title = "Forest Plot of Logistic Regression",
       x = "Variables",
       y = "Odds Ratio (95% CI)") +
  theme_minimal()

Hosmer-Lemeshow Test

hl_test <- hoslem.test(model$y, fitted(model), g = 10)
print(hl_test)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model$y, fitted(model)
## X-squared = 4.4189, df = 8, p-value = 0.8175

p-value > 0.05 -> no evidence of poor fit

ROC and AUC

roc_obj <- roc(response = filtered$outcome,
               predictor = predict(model, type = "response"))

## Setting levels: control = 1, case = 0

## Setting direction: controls < cases

auc(roc_obj)

## Area under the curve: 0.7508

ci.auc(roc_obj)

## 95% CI: 0.7041-0.7976 (DeLong)

AUC between 0.7 and 0.8: acceptable ability to discriminate between classes Entire confidence interval above 0.5: model performs better than random

plot(roc_obj, main = "ROC Curve")

Alternative D - Alternative C, but age and doses are covariates

Category: variable (Ref: ): groups…

• Age: covariate NU_IDADE_N
• Gender: CS_SEXO (Ref: F): M
• Educational Attainment: nominal_education (Ref: Bachelor): Elementary, Middle, Highschool, Unknown
• ICU: nominal_icu_admission (Ref: No): Yes
• Invasive Support Ven: nominal_invasive_support_ven (Ref: No): Yes=
• Doses: covariate doses

Logistic Regression

model <- glm(outcome ~ NU_IDADE_N + CS_SEXO +
               nominal_education +
               nominal_icu_admission + 
               nominal_invasive_support_ven + doses,
               family = binomial,
               data = filtered)

summary_model <- summary(model)
print(summary_model)

## 
## Call:
## glm(formula = outcome ~ NU_IDADE_N + CS_SEXO + nominal_education + 
##     nominal_icu_admission + nominal_invasive_support_ven + doses, 
##     family = binomial, data = filtered)
## 
## Coefficients:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     -2.687316   1.018435  -2.639  0.00832 ** 
## NU_IDADE_N                       0.008668   0.011257   0.770  0.44127    
## CS_SEXOM                         0.055478   0.219785   0.252  0.80072    
## nominal_educationElementary      1.656088   0.647406   2.558  0.01053 *  
## nominal_educationHighschool      0.846588   0.693633   1.221  0.22227    
## nominal_educationMiddle          1.662971   0.703020   2.365  0.01801 *  
## nominal_educationUnknown         1.157890   0.617952   1.874  0.06096 .  
## nominal_icu_admissionYes         0.358612   0.254405   1.410  0.15865    
## nominal_invasive_support_venYes  2.559883   0.363248   7.047 1.83e-12 ***
## doses                           -0.110484   0.068811  -1.606  0.10836    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 619.85  on 473  degrees of freedom
## Residual deviance: 511.15  on 464  degrees of freedom
## AIC: 531.15
## 
## Number of Fisher Scoring iterations: 4

odds_ratios <- exp(coef(model))
conf_int <- exp(confint(model))

## Waiting for profiling to be done...

data.frame(
  variable = names(coef(model)),
  OR = round(odds_ratios, 2),
  IC_95_inf = round(conf_int[,1], 2),
  IC_95_sup = round(conf_int[,2], 2),
  pvalue = round(summary_model$coefficients[,4], 3)
)

Variance inflation factor

vif_results <- car::vif(model)
print(vif_results)

##                                  GVIF Df GVIF^(1/(2*Df))
## NU_IDADE_N                   1.068426  1        1.033647
## CS_SEXO                      1.021925  1        1.010903
## nominal_education            1.134606  4        1.015911
## nominal_icu_admission        1.247004  1        1.116693
## nominal_invasive_support_ven 1.213735  1        1.101696
## doses                        1.054865  1        1.027066

No multicollinearity detected.

Forest Plot

label_map <- c(
  "NU_IDADE_N" = "Age",
  "CS_SEXOM" = "Male",
  "nominal_educationElementary" = "Elementary Education",
  "nominal_educationHighschool" = "Highschool Education",
  "nominal_educationMiddle" = "Middle Education",
  "nominal_educationUnknown" = "Unknown Education",
  "nominal_icu_admissionYes" = "ICU Admission",
  "nominal_invasive_support_venYes" = "Invasive Ventilation",
  "doses" = "Doses"
)


coef_df <- as.data.frame(confint(model))  # 95% CI

## Waiting for profiling to be done...

coef_df$Estimate <- coef(model)
coef_df$OR <- exp(coef_df$Estimate)
coef_df$LowerCI <- exp(coef_df[, 1])
coef_df$UpperCI <- exp(coef_df[, 2])
coef_df$Variable <- label_map[rownames(coef_df)]
coef_df= coef_df[row.names(coef_df) != "nominal_educationUnknown",]

ggplot(coef_df[-1, ], aes(x = reorder(Variable, OR), y = OR)) +  # remove Intercept
  geom_point() +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  coord_flip() +
  labs(title = "Forest Plot of Logistic Regression",
       x = "Variables",
       y = "Odds Ratio (95% CI)") +
  theme_minimal()

Hosmer-Lemeshow Test

hl_test <- hoslem.test(model$y, fitted(model), g = 10)
print(hl_test)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model$y, fitted(model)
## X-squared = 1.8567, df = 8, p-value = 0.9851

p-value > 0.05 -> no evidence of poor fit

ROC and AUC

roc_obj <- roc(response = filtered$outcome,
               predictor = predict(model, type = "response"))

## Setting levels: control = 1, case = 0

## Setting direction: controls < cases

auc(roc_obj)

## Area under the curve: 0.7419

ci.auc(roc_obj)

## 95% CI: 0.6938-0.79 (DeLong)

AUC between 0.7 and 0.8: acceptable ability to discriminate between classes Entire confidence interval above 0.5: model performs better than random

plot(roc_obj, main = "ROC Curve")

Alternative E - Alternative D, but it includes all categories from bivariate analyses

Category: variable (Ref: ): groups…

• Days from onset of first symtoms to hospitalization: covariate delta_days
• Age: covariate NU_IDADE_N
• Gender: CS_SEXO (Ref: F): M
• Race: nominal_race (Ref: White): Brown, Black, Yellow, Unknown
• Educational Attainment: nominal_education (Ref: Bachelor): Elementary, Middle, Highschool, Unknown
• Risk conditions: covariate additional_risk_comorbidities
• ICU: nominal_icu_admission (Ref: No): Yes
• Invasive Support Ven: nominal_invasive_support_ven (Ref: No): Yes
• Doses: covariate doses
• Vaccinal schema: nominal_vaccinal_schema (Ref: None): Inactivated only, mRNA only, Viral vector only, Inactivated + mRNA, Inactivate + viral vector, mRNA + viral vector, Inactivated + mRNA + viral vector
• Year of last COVID-19 vaccine: covariate last_year
• COVID treatment: nominal_covid_treatment (Ref: No): Yes

Logistic Regression

filtered$nominal_vaccinal_schema <- factor(
  filtered$nominal_vaccinal_schema,
  levels = levels(filtered$nominal_vaccinal_schema)
)

model <- glm(outcome ~ delta_days + NU_IDADE_N + CS_SEXO +
               nominal_race +
               nominal_education +
               additional_risk_comorbidities +
               nominal_icu_admission + 
               nominal_invasive_support_ven +
               doses +
               nominal_vaccinal_schema +
               last_year +
               nominal_covid_treatment,
               family = binomial,
               data = filtered)

summary_model <- summary(model)
print(summary_model)

## 
## Call:
## glm(formula = outcome ~ delta_days + NU_IDADE_N + CS_SEXO + nominal_race + 
##     nominal_education + additional_risk_comorbidities + nominal_icu_admission + 
##     nominal_invasive_support_ven + doses + nominal_vaccinal_schema + 
##     last_year + nominal_covid_treatment, family = binomial, data = filtered)
## 
## Coefficients:
##                                                          Estimate Std. Error
## (Intercept)                                              -3.31409    2.65335
## delta_days                                               -0.04275    0.02848
## NU_IDADE_N                                                0.02482    0.01393
## CS_SEXOM                                                  0.13925    0.24963
## nominal_raceBlack                                         0.35458    1.06268
## nominal_raceBrown                                         0.49023    0.30487
## nominal_raceUnknown                                       0.18784    0.42084
## nominal_raceYellow                                        0.16269    1.23857
## nominal_educationElementary                               1.45709    0.69560
## nominal_educationHighschool                               0.48380    0.75509
## nominal_educationMiddle                                   1.59383    0.78461
## nominal_educationUnknown                                  0.99758    0.66322
## additional_risk_comorbidities                            -0.19310    0.11664
## nominal_icu_admissionYes                                  0.50785    0.29580
## nominal_invasive_support_venYes                           2.63834    0.40998
## doses                                                    -0.12700    0.23107
## nominal_vaccinal_schemaInactivated + mRNA                -0.84109    2.39270
## nominal_vaccinal_schemaInactivated + mRNA + viral vector -0.24897    2.42230
## nominal_vaccinal_schemaInactivated + Viral vector        -0.70372    2.39199
## nominal_vaccinal_schemaInactivated only                  -0.70834    2.39020
## nominal_vaccinal_schemamRNA + viral vector               -0.38586    2.39747
## nominal_vaccinal_schemamRNA only                         -1.31187    2.47551
## nominal_vaccinal_schemaViral vector only                 -0.03772    2.35250
## last_year2022                                             0.55855    0.48196
## last_year2023                                             0.31857    0.70646
## last_year2024                                            -0.95293    1.58518
## nominal_covid_treatmentYes                               -2.17183    1.23411
##                                                          z value Pr(>|z|)    
## (Intercept)                                               -1.249   0.2117    
## delta_days                                                -1.501   0.1334    
## NU_IDADE_N                                                 1.782   0.0748 .  
## CS_SEXOM                                                   0.558   0.5770    
## nominal_raceBlack                                          0.334   0.7386    
## nominal_raceBrown                                          1.608   0.1078    
## nominal_raceUnknown                                        0.446   0.6554    
## nominal_raceYellow                                         0.131   0.8955    
## nominal_educationElementary                                2.095   0.0362 *  
## nominal_educationHighschool                                0.641   0.5217    
## nominal_educationMiddle                                    2.031   0.0422 *  
## nominal_educationUnknown                                   1.504   0.1325    
## additional_risk_comorbidities                             -1.655   0.0978 .  
## nominal_icu_admissionYes                                   1.717   0.0860 .  
## nominal_invasive_support_venYes                            6.435 1.23e-10 ***
## doses                                                     -0.550   0.5826    
## nominal_vaccinal_schemaInactivated + mRNA                 -0.352   0.7252    
## nominal_vaccinal_schemaInactivated + mRNA + viral vector  -0.103   0.9181    
## nominal_vaccinal_schemaInactivated + Viral vector         -0.294   0.7686    
## nominal_vaccinal_schemaInactivated only                   -0.296   0.7670    
## nominal_vaccinal_schemamRNA + viral vector                -0.161   0.8721    
## nominal_vaccinal_schemamRNA only                          -0.530   0.5962    
## nominal_vaccinal_schemaViral vector only                  -0.016   0.9872    
## last_year2022                                              1.159   0.2465    
## last_year2023                                              0.451   0.6520    
## last_year2024                                             -0.601   0.5477    
## nominal_covid_treatmentYes                                -1.760   0.0784 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 542.13  on 417  degrees of freedom
## Residual deviance: 420.12  on 391  degrees of freedom
##   (56 observations deleted due to missingness)
## AIC: 474.12
## 
## Number of Fisher Scoring iterations: 5

odds_ratios <- exp(coef(model))
conf_int <- exp(confint(model))

## Waiting for profiling to be done...

data.frame(
  variable = names(coef(model)),
  OR = round(odds_ratios, 2),
  IC_95_inf = round(conf_int[,1], 2),
  IC_95_sup = round(conf_int[,2], 2),
  pvalue = round(summary_model$coefficients[,4], 3)
)

Variance inflation factor

vif_results <- car::vif(model)
print(vif_results)

##                                    GVIF Df GVIF^(1/(2*Df))
## delta_days                     1.072189  1        1.035466
## NU_IDADE_N                     1.399091  1        1.182832
## CS_SEXO                        1.070896  1        1.034841
## nominal_race                   1.540999  4        1.055541
## nominal_education              1.445236  4        1.047110
## additional_risk_comorbidities  1.150630  1        1.072674
## nominal_icu_admission          1.427255  1        1.194678
## nominal_invasive_support_ven   1.375957  1        1.173012
## doses                          5.056808  1        2.248735
## nominal_vaccinal_schema       11.160259  7        1.188048
## last_year                      9.817672  3        1.463305
## nominal_covid_treatment        1.125500  1        1.060896

Doses and vaccinal schema might have some colinearity.

Forest Plot

label_map <- c(
  "NU_IDADE_N" = "Age",
  "CS_SEXOM" = "Male",
  "nominal_educationElementary" = "Elementary Education",
  "nominal_educationHighschool" = "Highschool Education",
  "nominal_educationMiddle" = "Middle Education",
  "nominal_educationUnknown" = "Unknown Education",
  "nominal_icu_admissionYes" = "ICU Admission",
  "nominal_invasive_support_venYes" = "Invasive Ventilation",
  "additional_risk_comorbidities" = "Comorbidities",
  "doses" = "Doses",
  "nominal_vaccinal_schemaInactivated + mRNA + viral vector" = "Inactivated + mRNA + viral vector",
  "nominal_vaccinal_schemaInactivated + Viral vector" = "Inactivated + Viral vector",
  "nominal_vaccinal_schemaInactivated only" = "Inactivated only",
  "nominal_vaccinal_schemaViral vector only" = "Viral vector only",
  "last_year" = "Last Vaccine Year",
  "nominal_covid_treatmentYes" = "COVID Treatment",
  "nominal_raceBlack" = "Black",
  "nominal_raceBrown" = "Brown",
  "nominal_raceUnknown" = "Unknown Race",
  "nominal_raceYellow" = "Yellow",
  "delta_days" = "Days since the beginning of the infection",
  "nominal_vaccinal_schemamRNA + viral vector" = "mRNA + viral vector",
  "nominal_vaccinal_schemamRNA only" = "mRNA only"
)

coef_df <- as.data.frame(confint(model))  # 95% CI

## Waiting for profiling to be done...

coef_df$Estimate <- coef(model)
coef_df$OR <- exp(coef_df$Estimate)
coef_df$LowerCI <- exp(coef_df[, 1])
coef_df$UpperCI <- exp(coef_df[, 2])
#coef_df$Variable <- rownames(coef_df)

coef_df$Variable <- label_map[rownames(coef_df)]
coef_df= coef_df[row.names(coef_df) != "nominal_educationUnknown" & row.names(coef_df) != "nominal_raceUnknown" ,]

ggplot(coef_df[-1, ], aes(x = reorder(Variable, OR), y = OR)) +  # remove Intercept
  geom_point() +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  coord_flip() +
  labs(title = "Forest Plot of Logistic Regression",
       x = "Variables",
       y = "Odds Ratio (95% CI)") +
  theme_minimal()

Hosmer-Lemeshow Test

hl_test <- hoslem.test(model$y, fitted(model), g = 10)
print(hl_test)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model$y, fitted(model)
## X-squared = 3.2172, df = 8, p-value = 0.92

p-value > 0.05 -> no evidence of poor fit

Alternative F - Alternative D, but doses and are are categorical

Category: variable (Ref: ): groups…

• Age: age_groups (Ref: 60–69): 70-79, 80+
• ICU: nominal_icu_admission (Ref: No): Yes
• Invasive Support Ven: nominal_invasive_support_ven (Ref: No): Yes
• Gender: CS_SEXO (Ref: F): M
• Doses: categorical_doses (Ref: unvaccinated): 1-2, 3, 4, 5+
• Educational Attainment: educational_attainment (Ref: bachelor): low, unknown

Logistic Regression

model <- glm(outcome ~ age_groups + nominal_icu_admission + 
               nominal_invasive_support_ven + CS_SEXO +
               categorical_doses + 
               educational_attainment,
               family = binomial,
               data = filtered)

summary_model <- summary(model)
print(summary_model)

## 
## Call:
## glm(formula = outcome ~ age_groups + nominal_icu_admission + 
##     nominal_invasive_support_ven + CS_SEXO + categorical_doses + 
##     educational_attainment, family = binomial, data = filtered)
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     -1.88617    0.69655  -2.708  0.00677 ** 
## age_groups70-79                 -0.12722    0.27159  -0.468  0.63948    
## age_groups80+                    0.13245    0.26533   0.499  0.61765    
## nominal_icu_admissionYes         0.35219    0.25247   1.395  0.16303    
## nominal_invasive_support_venYes  2.58869    0.36007   7.189  6.5e-13 ***
## CS_SEXOM                         0.04112    0.21918   0.188  0.85119    
## categorical_doses1-2            -0.38303    0.37319  -1.026  0.30471    
## categorical_doses3              -0.53548    0.37957  -1.411  0.15831    
## categorical_doses4              -0.68220    0.37758  -1.807  0.07080 .  
## categorical_doses5+             -0.67934    0.37487  -1.812  0.06996 .  
## educational_attainmentlow        1.45115    0.63799   2.275  0.02293 *  
## educational_attainmentunknown    1.16616    0.63106   1.848  0.06461 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 619.85  on 473  degrees of freedom
## Residual deviance: 514.90  on 462  degrees of freedom
## AIC: 538.9
## 
## Number of Fisher Scoring iterations: 4

odds_ratios <- exp(coef(model))
conf_int <- exp(confint(model))

## Waiting for profiling to be done...

data.frame(
  variable = names(coef(model)),
  OR = round(odds_ratios, 2),
  IC_95_inf = round(conf_int[,1], 2),
  IC_95_sup = round(conf_int[,2], 2),
  pvalue = round(summary_model$coefficients[,4], 3)
)

Variance inflation factor

vif_results <- car::vif(model)
print(vif_results)

##                                  GVIF Df GVIF^(1/(2*Df))
## age_groups                   1.064477  2        1.015744
## nominal_icu_admission        1.235072  1        1.111338
## nominal_invasive_support_ven 1.200929  1        1.095869
## CS_SEXO                      1.025043  1        1.012444
## categorical_doses            1.059239  4        1.007220
## educational_attainment       1.084161  2        1.020407

No multicollinearity detected.

Forest Plot

label_map <- c(
  "age_groups70-79" = "Age 70-79",
  "age_groups80+" = "Age 80+",
  "CS_SEXOM" = "Male",
  "educational_attainmentlow" = "Low Education",
  "educational_attainmentunknown" = "Unknown Education",
  "nominal_educationMiddle" = "Middle Education",
  "nominal_educationUnknown" = "Unknown Education",
  "nominal_icu_admissionYes" = "ICU Admission",
  "nominal_invasive_support_venYes" = "Invasive Ventilation",
  "categorical_doses1-2" = "1-2 Doses",
  "categorical_doses3" = "3 Doses",
  "categorical_doses4" = "4 Doses",
  "categorical_doses5+" = "5+ Doses"
)


coef_df <- as.data.frame(confint(model))  # 95% CI

## Waiting for profiling to be done...

coef_df$Estimate <- coef(model)
coef_df$OR <- exp(coef_df$Estimate)
coef_df$LowerCI <- exp(coef_df[, 1])
coef_df$UpperCI <- exp(coef_df[, 2])
coef_df$Variable <- label_map[rownames(coef_df)]
coef_df= coef_df[row.names(coef_df) != "educational_attainmentunknown",]

ggplot(coef_df[-1, ], aes(x = reorder(Variable, OR), y = OR)) +  # remove Intercept
  geom_point() +
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  coord_flip() +
  labs(title = "Forest Plot of Logistic Regression",
       x = "Variables",
       y = "Odds Ratio (95% CI)") +
  theme_minimal()

Hosmer-Lemeshow Test

hl_test <- hoslem.test(model$y, fitted(model), g = 10)
print(hl_test)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model$y, fitted(model)
## X-squared = 2.2298, df = 8, p-value = 0.9731

p-value > 0.05 -> no evidence of poor fit

ROC and AUC

roc_obj <- roc(response = filtered$outcome,
               predictor = predict(model, type = "response"))

## Setting levels: control = 1, case = 0

## Setting direction: controls < cases

auc(roc_obj)

## Area under the curve: 0.7303

ci.auc(roc_obj)

## 95% CI: 0.681-0.7795 (DeLong)

AUC between 0.7 and 0.8: acceptable ability to discriminate between classes Entire confidence interval above 0.5: model performs better than random

plot(roc_obj, main = "ROC Curve")

Vaccination schema

filtered$vaccinal_schema_based_on_innactivated <- relevel(factor(filtered$vaccinal_schema_based_on_innactivated), ref="unvaccinated or single-dose")
model <- glm(outcome ~ vaccinal_schema_based_on_innactivated,
               family = binomial,
               data = filtered)

summary_model <- summary(model)
print(summary_model)

## 
## Call:
## glm(formula = outcome ~ vaccinal_schema_based_on_innactivated, 
##     family = binomial, data = filtered)
## 
## Coefficients:
##                                                                             Estimate
## (Intercept)                                                                  -0.3677
## vaccinal_schema_based_on_innactivatedsequential excluding inactivated        -0.1368
## vaccinal_schema_based_on_innactivatedsequential with inactivated and others  -0.2672
## vaccinal_schema_based_on_innactivatedsequential with inactivated only        -0.3969
##                                                                             Std. Error
## (Intercept)                                                                     0.2504
## vaccinal_schema_based_on_innactivatedsequential excluding inactivated           0.2962
## vaccinal_schema_based_on_innactivatedsequential with inactivated and others     0.3026
## vaccinal_schema_based_on_innactivatedsequential with inactivated only           0.3420
##                                                                             z value
## (Intercept)                                                                  -1.469
## vaccinal_schema_based_on_innactivatedsequential excluding inactivated        -0.462
## vaccinal_schema_based_on_innactivatedsequential with inactivated and others  -0.883
## vaccinal_schema_based_on_innactivatedsequential with inactivated only        -1.161
##                                                                             Pr(>|z|)
## (Intercept)                                                                    0.142
## vaccinal_schema_based_on_innactivatedsequential excluding inactivated          0.644
## vaccinal_schema_based_on_innactivatedsequential with inactivated and others    0.377
## vaccinal_schema_based_on_innactivatedsequential with inactivated only          0.246
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 619.85  on 473  degrees of freedom
## Residual deviance: 618.18  on 470  degrees of freedom
## AIC: 626.18
## 
## Number of Fisher Scoring iterations: 4

filtered$nominal_vaccinal_schema <- relevel(factor(filtered$nominal_vaccinal_schema), ref="None")
model <- glm(outcome ~ nominal_vaccinal_schema,
               family = binomial,
               data = filtered)

summary_model <- summary(model)
print(summary_model)

## 
## Call:
## glm(formula = outcome ~ nominal_vaccinal_schema, family = binomial, 
##     data = filtered)
## 
## Coefficients:
##                                                          Estimate Std. Error
## (Intercept)                                              -0.27763    0.26515
## nominal_vaccinal_schemaInactivated + mRNA                -0.68995    0.37807
## nominal_vaccinal_schemaInactivated + mRNA + viral vector -0.24367    0.39266
## nominal_vaccinal_schemaInactivated + Viral vector         0.09531    0.43878
## nominal_vaccinal_schemaInactivated only                  -0.52088    0.35215
## nominal_vaccinal_schemamRNA + viral vector               -0.41552    0.37214
## nominal_vaccinal_schemamRNA only                         -0.92634    0.70967
## nominal_vaccinal_schemaViral vector only                 -0.03415    0.33549
##                                                          z value Pr(>|z|)  
## (Intercept)                                               -1.047    0.295  
## nominal_vaccinal_schemaInactivated + mRNA                 -1.825    0.068 .
## nominal_vaccinal_schemaInactivated + mRNA + viral vector  -0.621    0.535  
## nominal_vaccinal_schemaInactivated + Viral vector          0.217    0.828  
## nominal_vaccinal_schemaInactivated only                   -1.479    0.139  
## nominal_vaccinal_schemamRNA + viral vector                -1.117    0.264  
## nominal_vaccinal_schemamRNA only                          -1.305    0.192  
## nominal_vaccinal_schemaViral vector only                  -0.102    0.919  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 619.85  on 473  degrees of freedom
## Residual deviance: 611.32  on 466  degrees of freedom
## AIC: 627.32
## 
## Number of Fisher Scoring iterations: 4

filtered$factor_innactivated_doses <- relevel(factor(filtered$innactivated_doses), ref=1)
filtered$factor_viral_vector_doses<- relevel(factor(filtered$viral_vector_doses ), ref=1)
filtered$factor_mrna_doses <- relevel(factor(filtered$mrna_doses), ref=1)
filtered$nominal_vaccinal_schema <- relevel(factor(filtered$nominal_vaccinal_schema), ref="None")
model <- glm(outcome ~ factor_innactivated_doses + factor_viral_vector_doses + factor_mrna_doses,
               family = binomial,
               data = filtered)

summary_model <- summary(model)
print(summary_model)

## 
## Call:
## glm(formula = outcome ~ factor_innactivated_doses + factor_viral_vector_doses + 
##     factor_mrna_doses, family = binomial, data = filtered)
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                -3.325e-01  2.283e-01  -1.457   0.1453  
## factor_innactivated_doses1 -5.215e-01  6.227e-01  -0.837   0.4023  
## factor_innactivated_doses2 -3.556e-01  2.725e-01  -1.305   0.1919  
## factor_innactivated_doses3 -1.506e+01  5.875e+02  -0.026   0.9795  
## factor_viral_vector_doses1  6.419e-01  2.828e-01   2.270   0.0232 *
## factor_viral_vector_doses2  3.157e-05  2.914e-01   0.000   0.9999  
## factor_viral_vector_doses3 -1.321e-01  4.298e-01  -0.307   0.7586  
## factor_viral_vector_doses4 -1.523e+01  1.455e+03  -0.010   0.9916  
## factor_mrna_doses1         -7.664e-03  3.476e-01  -0.022   0.9824  
## factor_mrna_doses2         -5.404e-01  2.649e-01  -2.040   0.0413 *
## factor_mrna_doses3         -1.956e-01  3.819e-01  -0.512   0.6085  
## factor_mrna_doses4         -9.943e-01  8.609e-01  -1.155   0.2481  
## factor_mrna_doses5         -1.523e+01  1.455e+03  -0.010   0.9916  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 619.85  on 473  degrees of freedom
## Residual deviance: 601.95  on 461  degrees of freedom
## AIC: 627.95
## 
## Number of Fisher Scoring iterations: 14

filtered$received_innactivated_doses <- relevel(factor(filtered$innactivated_doses >= 1), ref=1)
filtered$received_viral_vector_doses<- relevel(factor(filtered$viral_vector_doses >= 1), ref=1)
filtered$received_mrna_doses <- relevel(factor(filtered$mrna_doses >= 1), ref=1)

model <- glm(outcome ~ received_innactivated_doses + received_viral_vector_doses + received_mrna_doses,
               family = binomial,
               data = filtered)

summary_model <- summary(model)
print(summary_model)

## 
## Call:
## glm(formula = outcome ~ received_innactivated_doses + received_viral_vector_doses + 
##     received_mrna_doses, family = binomial, data = filtered)
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)   
## (Intercept)                     -0.54799    0.20156  -2.719  0.00655 **
## received_innactivated_dosesTRUE -0.08046    0.20960  -0.384  0.70109   
## received_viral_vector_dosesTRUE  0.30492    0.20939   1.456  0.14533   
## received_mrna_dosesTRUE         -0.35767    0.20299  -1.762  0.07807 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 619.85  on 473  degrees of freedom
## Residual deviance: 613.81  on 470  degrees of freedom
## AIC: 621.81
## 
## Number of Fisher Scoring iterations: 4