creating_datasets.Rmd

---
title: "SDTM and ADaM Dataset Processing"
author: "ZW"
date: "`r Sys.Date()`"
output:
  word_document: default
  html_notebook: default
---
# Introduction

This document describes the process of loading, transforming, and analyzing clinical trial data for SDTM and ADaM compliance. It also includes exploratory data analysis.


```{r setup, include=FALSE}
library(ggplot2)
library(dplyr)
library(knitr)
library(pander)
library(moments) # For skewness and kurtosis
library(psych)
library(pROC)
source("automatyzacja.R",encoding = "UTF-8")
source("tabele_dla_modeli_ang.R", encoding = "UTF-8")
```

## Load and Prepare Data

```{r load-data}
# Set working directory and load data
setwd("C:/Users/zuzan/OneDrive/Pulpit/covid")
raw_data <- read.csv2("C:/Users/zuzan/OneDrive/Pulpit/covid/faktyczne_dane.csv", 
                      stringsAsFactors = FALSE, check.names = FALSE, 
                      fileEncoding = "UTF-8")

# Export column names for reference
write.csv(data.frame(Columns = colnames(raw_data)), "columns_list.csv", row.names = FALSE)

# Adjust column names for processing
colnames(raw_data) <- make.names(colnames(raw_data), unique = TRUE)

```

```{r}
#column_types <- sapply(raw_data, class)
```
## Numeric Values
```{r}
numeric_df <- raw_data[sapply(raw_data, is.numeric)]
selected_df <- numeric_df[1:5, 3:(3+5-1)]  # Ogranicza liczbę wierszy do 5 i wybiera odpowiednie kolumny
print(selected_df)

```
## Descriptive statistics
```{r}
res<- tabela_wynikowa_tylko_ogol(df=numeric_df,colnames(numeric_df), wyjatki = NULL)
#pander(res, caption = "Statystyki opisowe dla badań klinicznych")
```

```{r}
# sprawdzić pusta liczbe kolumn 
# empty_cols <- colSums(is.na(raw_data_clean)) == nrow(raw_data_clean)
# length(empty_cols)
# sum(empty_cols)  # is.na() sprawdza tylko wartości NA, ale kolumny mogą mieć inne „puste” wartości, np. puste stringi ("") lub wartości NULL
```
```{r}
#unique(raw_data)
```





```{r}
empty_name_spaces <- sum(trimws(names(raw_data)) == "")

```
```{r}
empty_cols <- sapply(raw_data, function(col) all(is.na(col) | col == "" | trimws(as.character(col)) == ""))

cat("Liczba całkowicie pustych kolumn:", sum(empty_cols), "\n")
cat("Liczba wszystkich kolumn:", length(empty_cols), "\n")
cat("Liczba kolumn w raw_data:", ncol(raw_data), "\n")

```
```{r}
raw_data_clean <- raw_data[, !empty_cols]
dim(raw_data_clean)  # Sprawdzenie nowego wymiaru danych

```


## Add Missing Columns

```{r add-missing-columns}
if (!"STUDYID" %in% colnames(raw_data)) {
  raw_data$STUDYID <- "Study1"
  cat("Column STUDYID added with default value 'Study1'.\n")
}

if (!"USUBJID" %in% colnames(raw_data)) {
  cat("Column USUBJID is missing. Creating it dynamically using STUDYID and row numbers.\n")
  raw_data$USUBJID <- paste0(raw_data$STUDYID, "_", seq_len(nrow(raw_data)))
}
```

## Define and Map SDTM Domains

```{r define-sdtm}
sdtm_mappings <- list(
  DM = list(
    STUDYID = "StudyID",
    DOMAIN = "DM",
    USUBJID = "USUBJID",
    AGE = "AGE",
    SEX = "SEX",
    RACE = "RACE",
    ETHNIC = "ETHNIC"
  ),
  AE = list(
    STUDYID = "StudyID",
    DOMAIN = "AE",
    USUBJID = "USUBJID",
    AETERM = "AdverseEvent",
    AESTDTC = "AE.AESTDAT.DATE..Data.rozpoczęcia.os1",
    AEENDTC = "AE.AEENDAT..Data.zakończenia.os1",
    AESEV = "Severity",
    AESER = "AE.AESER..Czy.zdarzenie.ciężkie..os1"
  ),
  VS = list(
    STUDYID = "StudyID",
    DOMAIN = "VS",
    USUBJID = "USUBJID",
    VSTESTCD = "VitalSignTest",
    VSORRES = "Result",
    VSSTRESN = "StandardizedResult",
    VSDTC = "V0.VS.VSDAT.DATE..Data.wykonania.pomiarów",
    VSPOS = "Position"
  )
)
```

```{r map-domains}
create_sdtm <- function(raw_data, domain_mapping, domain_name) {
  mapped_columns <- sapply(domain_mapping, function(pattern) {
    match <- agrep(pattern, colnames(raw_data), value = TRUE, ignore.case = TRUE, max.distance = 0.2)
    if (length(match) > 0) {
      return(match[1])  # Use the best match
    } else {
      return(NA)  # No match found
    }
  })
  
  cat("Mapped columns for", domain_name, "domain:\n")
  print(mapped_columns)
  
  domain_data <- raw_data[, na.omit(mapped_columns), drop = FALSE]
  colnames(domain_data) <- names(na.omit(mapped_columns))
  
  if (!"DOMAIN" %in% colnames(domain_data)) {
    domain_data$DOMAIN <- domain_name
  }
  
  if (!"USUBJID" %in% colnames(domain_data)) {
    stop(paste("USUBJID is missing in the", domain_name, "domain. Check the raw data or mapping."))
  }
  
  return(domain_data)
}

sdtm_dm <- create_sdtm(raw_data, sdtm_mappings$DM, "DM")
sdtm_ae <- create_sdtm(raw_data, sdtm_mappings$AE, "AE")
sdtm_vs <- create_sdtm(raw_data, sdtm_mappings$VS, "VS")
```
```{r}
# Ensure AGE is numeric
if ("AGE" %in% colnames(raw_data)) {
  raw_data$AGE <- as.numeric(raw_data$AGE)
  
  # Handle missing values in AGE
  if (any(is.na(raw_data$AGE))) {
    raw_data$AGE[is.na(raw_data$AGE)] <- median(raw_data$AGE, na.rm = TRUE)
    cat("Missing values in AGE were replaced with the median.\n")
  }
}

```

## Save SDTM Datasets

```{r save-sdtm}
write.csv(sdtm_dm, "sdtm_dm.csv", row.names = FALSE)
write.csv(sdtm_ae, "sdtm_ae.csv", row.names = FALSE)
write.csv(sdtm_vs, "sdtm_vs.csv", row.names = FALSE)
```

## Create ADaM Datasets

```{r create-adam}
adsl <- sdtm_dm %>%
  mutate(AGEGR1 = ifelse(as.numeric(AGE) < 18, "<18", ifelse(as.numeric(AGE) >= 65, "65+", "18-64")),
         SEXN = ifelse(SEX == "M", 1, ifelse(SEX == "F", 2, NA))) %>%
  select(STUDYID, USUBJID, AGE, AGEGR1, SEX, SEXN, RACE, ETHNIC)

write.csv(adsl, "adam_adsl.csv", row.names = FALSE)
cat("ADSL dataset created and saved as adam_adsl.csv\n")
```

## Descriptive Statistics for ADSL

```{r adsl-stats}
adsl_stats <- adsl %>% 
  summarise(
    AGE_MEAN = mean(AGE, na.rm = TRUE),
    AGE_MEDIAN = median(AGE, na.rm = TRUE),
    AGE_SD = sd(AGE, na.rm = TRUE),
    AGE_VARIANCE = var(AGE, na.rm = TRUE),
    AGE_SKEWNESS = skewness(AGE, na.rm = TRUE),
    AGE_KURTOSIS = kurtosis(AGE, na.rm = TRUE),
    AGE_MIN = min(AGE, na.rm = TRUE),
    AGE_MAX = max(AGE, na.rm = TRUE),
    AGE_QUANTILES = list(quantile(AGE, probs = c(0.25, 0.5, 0.75), na.rm = TRUE))
  )
print(adsl_stats)
```

## Gender Distribution

```{r gender-distribution}
gender_counts <- adsl %>% 
  mutate(SEX = ifelse(SEX == "", "Unknown", SEX)) %>%  # Add Unknown category for missing values
  group_by(SEX) %>% 
  summarise(COUNT = n())
print(gender_counts)

gender_plot <- ggplot(gender_counts, aes(x = SEX, y = COUNT, fill = SEX)) +
  geom_bar(stat = "identity") +
  labs(title = "Gender Distribution", x = "Gender", y = "Count") +
  theme_minimal()
print(gender_plot)
```

## Save Summary to File

```{r save-summary}
sink("adam_summary.txt")
cat("ADSL Summary:\n")
summary(adsl)
cat("\nADAE Summary:\n")
summary(sdtm_ae)
cat("\nADVS Summary:\n")
summary(sdtm_vs)
sink()
```

## Descriptive Statistics for ADVS - VSORRES

#### VSORRES is a generic column name in the vital signs (VS) domain used to store measurement results for various vital signs such as diastolic blood pressure, heart rate, respiratory rate, and body temperature

#### Temperature

```{r}
class("V0.COVIDOBJ.VSTEST_TEMP.VSORRES..Maksymalna.temperatura.ciała...C.")
as.numeric(raw_data$V0.COVIDOBJ.VSTEST_TEMP.VSORRES..Maksymalna.temperatura.ciała...C.)
```


```{r advs-stats}
temp_stats <- raw_data %>% 
  summarise(
    Mean_temperature =mean(as.numeric(raw_data$V0.COVIDOBJ.VSTEST_TEMP.VSORRES..Maksymalna.temperatura.ciała...C.), na.rm = TRUE), 
    median_temp = median(as.numeric(raw_data$V0.COVIDOBJ.VSTEST_TEMP.VSORRES..Maksymalna.temperatura.ciała...C.), na.rm = TRUE),
    sd_tem = sd(as.numeric(raw_data$V0.COVIDOBJ.VSTEST_TEMP.VSORRES..Maksymalna.temperatura.ciała...C. ), na.rm = TRUE),
    var_temp =var(as.numeric(raw_data$V0.COVIDOBJ.VSTEST_TEMP.VSORRES..Maksymalna.temperatura.ciała...C. ), na.rm = TRUE),
    Quartiles = list(quantile(as.numeric(raw_data$V0.COVIDOBJ.VSTEST_TEMP.VSORRES..Maksymalna.temperatura.ciała...C.), probs = c(0.25, 0.5, 0.75), na.rm = TRUE))
  )

print(temp_stats)
```


## Filter Outliers for Body Temperature

```{r filter-outliers}
# Identify outliers based on predefined thresholds
temperature_column <- "V0.COVIDOBJ.VSTEST_TEMP.VSORRES..Maksymalna.temperatura.ciała...C."
threshold_min <- 30
threshold_max <- 45

raw_data <- raw_data %>% 
  mutate(!!sym(temperature_column) := as.numeric(!!sym(temperature_column)))

# Count and filter out extreme outliers
extreme_outliers <- raw_data %>% 
  filter(!!sym(temperature_column) < threshold_min | !!sym(temperature_column) > threshold_max)

valid_data <- raw_data %>% 
  filter(!!sym(temperature_column) >= threshold_min & !!sym(temperature_column) <= threshold_max)

cat("Number of extreme outliers in temperature:", nrow(extreme_outliers), "\n")
cat("Filtered data dimensions:", dim(valid_data), "\n")
```

## Boxplot for Filtered Body Temperature

```{r filtered-temp-boxplot}
boxplot_filtered_temp <- ggplot(valid_data, aes(y = !!sym(temperature_column))) +
  geom_boxplot() +
  labs(title = "Filtered Body Temperature Distribution", y = "Temperature (°C)") +
  theme_minimal()
print(boxplot_filtered_temp)
```
```{r}
filtered_temp_stats <- valid_data %>%
  summarise(
    Mean_Temperature = mean(!!sym(temperature_column), na.rm = TRUE),
    Median_Temperature = median(!!sym(temperature_column), na.rm = TRUE),
    SD_Temperature = sd(!!sym(temperature_column), na.rm = TRUE),
    Variance_Temperature = var(!!sym(temperature_column), na.rm = TRUE),
    Skewness_Temperature = skewness(!!sym(temperature_column), na.rm = TRUE),
    Kurtosis_Temperature = kurtosis(!!sym(temperature_column), na.rm = TRUE),
    Min_Temperature = min(!!sym(temperature_column), na.rm = TRUE),
    Max_Temperature = max(!!sym(temperature_column), na.rm = TRUE),
    Quartiles_Temperature = list(quantile(!!sym(temperature_column), probs = c(0.25, 0.5, 0.75), na.rm = TRUE))
  )

print("Statistics for Filtered Body Temperature:")
print(filtered_temp_stats)

```



### Type of measurements 
```{r Checking type of measurements}
grep("VSORRES", colnames(raw_data), value = TRUE, ignore.case = TRUE)

print("VSORRES Column Content Preview:")
print(head(raw_data[, grep("VSORRES", colnames(raw_data), value = TRUE, ignore.case = TRUE), drop = FALSE]))
```

## Smoking Status Distribution

```{r smoking-status}
if("V0.DD.SU.SUTRT_CIGARETTES.SUENRF..Palenie.papierosów" %in% colnames(raw_data)) {
  smoking_data <- raw_data %>% 
    group_by(`V0.DD.SU.SUTRT_CIGARETTES.SUENRF..Palenie.papierosów`) %>% 
    summarise(Count = n())
  
  print("Smoking Status Counts:")
  print(smoking_data)
  
  

} else {
  print("Smoking status column not found in the dataset.")
}

```
```{r}
smoking_plot <- ggplot(smoking_data, aes(x = `V0.DD.SU.SUTRT_CIGARETTES.SUENRF..Palenie.papierosów`, y = Count, fill = `V0.DD.SU.SUTRT_CIGARETTES.SUENRF..Palenie.papierosów`)) +
  geom_bar(stat = "identity") +
  labs(
    title = "Rozkład statusu palenia papierosów",
    x = "Status palenia",
    y = "Liczba pacjentów",
    fill = "Status"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold") 
  )

print(smoking_plot)

```
# PORÓWNANIE TESTÓW 
## CZY RÓŻNIĄ SIĘ MIEDZY GRUPAMI WYNIKI 


# REGRESJA LOGISTYCZNA - sprawdzanie celu I rzędowego


##skrypt do formatowania  regresja



### Zmienna zależna to V0.PR.INVAVENT_PRTRT.PROCCUR, czyli informacja o tym, czy pacjent otrzymuje inwazyjne wsparcie oddechowe (wentylacja mechaniczna).

```{r}
# Konwersja zmiennej zależnej na faktor
raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR <- as.factor(raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR)

```


### Interpretacja wartości Estimate (β)

- **β > 0** → zwiększa szanse na wentylację mechaniczną.
- **β < 0** → zmniejsza szanse na wentylację mechaniczną.
- **β = 0** → brak wpływu.

**Przykład:**  
Wiek pacjenta ma **β = 0.02782**, co oznacza, że **każdy dodatkowy rok życia zwiększa logarytm szansy na wentylację o 0.02782**.


---

```{r}
logistic_model_final <- glm(
  V0.PR.INVAVENT_PRTRT.PROCCUR ~ V0.DD.DM.AGE..Wiek + V0.DD.DM.SEX..Płeć + V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel,
  data = raw_data, family = binomial()
)
summary(logistic_model_interaction)
```
```{r}
coef(summary(logistic_model_final))
```

```{r}
wyniki_tabeli <- tabela(logistic_model_final, signif, F, ufn=T)
print(wyniki_tabeli)
```

```{r}
logistic_model_interaction <- glm(V0.PR.INVAVENT_PRTRT.PROCCUR ~ `V0.DD.DM.AGE..Wiek` * `V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel`, 
                                  data = raw_data, family = binomial())
summary(logistic_model_interaction)

```

## Logistic Regression Analysis with Interaction Term

A logistic regression model was used to predict the need for mechanical ventilation (`V0.PR.INVAVENT_PRTRT.PROCCUR`) based on patient age, cough presence, and their interaction.

### **Key Findings**
✔ The interaction between **age and cough is not statistically significant** (**p = 0.926**).  
✔ Adding this interaction **increases model deviance** (**140.93 vs. 140.17**) and **AIC** (**148.93 vs. 148.17**), meaning the model is **not improved**.  
✔ **Age was significant in the previous model (p = 0.0402), but is no longer significant after adding interaction (p = 0.218).**  
✔ **The interaction should not be included in the final model** as it does not contribute to better predictions.  

### **Next Steps**
- Try interactions with other variables (**age × oxygen saturation, age × temperature**).
- Include more clinical markers (**CRP, blood pressure, SpO₂**).
- Compare logistic regression with **random forest or decision trees**.
```{r}
table(raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR)

```
```{r}
unique(raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR)

```



```{r}
table(raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR, useNA = "always")
```


```{r}
str(raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR)
```
```{r}
raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR[raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR == ""] <- NA
```

```{r}
raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR <- droplevels(raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR)

```

```{r}
table(raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR, useNA = "always")
str(raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR)
```


```{r}
raw_data$V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel
```
```{r}
sum(is.na(raw_data$V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel))
```

```{r}
raw_data_clean <- raw_data[!is.na(raw_data$V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel), ]
table(raw_data_clean$V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel, useNA = "always")
```

```{r}
logistic_model_clean <- glm(
  V0.PR.INVAVENT_PRTRT.PROCCUR ~ V0.DD.DM.AGE..Wiek * V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel, 
  data = raw_data_clean, 
  family = binomial()
)

summary(logistic_model_clean)

```
```{r}
table(raw_data_clean$V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel, raw_data_clean$V0.PR.INVAVENT_PRTRT.PROCCUR)

```


```{r}
grep("CHESTPRESSURE", names(raw_data), value = TRUE)

```
```{r}
logistic_model_extended <- glm(V0.PR.INVAVENT_PRTRT.PROCCUR ~ 
                                V0.DD.DM.AGE..Wiek + 
                                V0.DD.DM.SEX..Płeć + 
                                V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel +
                                `V0.COVIDOBJ.FATESTCD_OCCUR.CHESTPRESSURE_FAORRES.07.mc..Ucisk.w.klatce.piersiowej`,
                              data = raw_data, family = binomial())
summary(logistic_model_extended)


```
```{r}
length(raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR)

```

```{r}
data_used_in_model <- raw_data[complete.cases(raw_data[, c(
  "V0.PR.INVAVENT_PRTRT.PROCCUR",
  "V0.DD.DM.AGE..Wiek",
  "V0.DD.DM.SEX..Płeć",
  "V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel",
  "V0.COVIDOBJ.FATESTCD_OCCUR.CHESTPRESSURE_FAORRES.07.mc..Ucisk.w.klatce.piersiowej"
)]), ]

```

```{r}

# Oblicz prawdopodobieństwa przewidywane przez model
predicted_probs <- predict(logistic_model_extended, type = "response")
length(predicted_probs)


roc_curve <- roc(data_used_in_model$V0.PR.INVAVENT_PRTRT.PROCCUR, predicted_probs)
plot(roc_curve, main = "ROC Curve for Logistic Model")


# Rysujemy wykres
plot(roc_curve, col = "blue", main = "Krzywa ROC dla modelu logistycznego")
auc(roc_curve)  # Wyświetlenie wartości AUC
```
```{r}
library(ggplot2)

# Tworzymy zbiór nowych wartości wieku do predykcji
new_data <- data.frame(
  V0.DD.DM.AGE..Wiek = seq(min(raw_data_clean$V0.DD.DM.AGE..Wiek, na.rm = TRUE),
                            max(raw_data_clean$V0.DD.DM.AGE..Wiek, na.rm = TRUE), length.out = 100),
  V0.DD.DM.SEX..Płeć = "1. Kobieta",  # Załóżmy płeć dla predykcji
  V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel = 0,
  V0.COVIDOBJ.FATESTCD_OCCUR.CHESTPRESSURE_FAORRES.07.mc..Ucisk.w.klatce.piersiowej = 0
)

# Predykcja wartości dla nowego zbioru
new_data$predicted_probs <- predict(logistic_model_extended, newdata = new_data, type = "response")

# Tworzymy wykres
ggplot(new_data, aes(x = V0.DD.DM.AGE..Wiek, y = predicted_probs)) +
  geom_line(color = "red") +
  labs(title = "Predykowane prawdopodobieństwo wentylacji mechanicznej",
       x = "Wiek pacjenta", y = "Prawdopodobieństwo") +
  theme_minimal()

```
  


```{r}
# raw_data_clean <- na.omit(raw_data)  # Usunięcie wierszy z NA
# raw_data_clean$V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel
```

```{r}

logistic_model_interaction <- glm(V0.PR.INVAVENT_PRTRT.PROCCUR ~ `V0.DD.DM.AGE..Wiek` * `V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel`, 
                                  data = raw_data_clean, family = binomial())

# Sprawdzenie długości po czyszczeniu
length(predict(logistic_model_interaction, type = "response"))

```

```{r}
# Model regresji logistycznej BEZ interakcji (ponieważ interakcja była nieistotna)
logistic_model_final <- glm(
  V0.PR.INVAVENT_PRTRT.PROCCUR ~ `V0.DD.DM.AGE..Wiek` + `V0.DD.DM.SEX..Płeć` + `V0.COVIDOBJ.FATESTCD_OCCUR.COUGH_FAORRES.02.mc..Kaszel`,
  data = raw_data,
  family = binomial()
)

# Podsumowanie modelu
summary(logistic_model_final)

# Obliczenie pseudo-R²
if (!requireNamespace("pscl", quietly = TRUE)) install.packages("pscl")
library(pscl)
pR2(logistic_model_final)

```

```{r}

# Ekstrakcja współczynników i przedziałów ufności
coefs <- as.data.frame(confint(logistic_model_final))
coefs$Estimate <- coef(logistic_model_final)
coefs$Variable <- rownames(coefs)

# Wykres Forest Plot
ggplot(coefs, aes(x = Variable, y = Estimate, ymin = `2.5 %`, ymax = `97.5 %`)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  coord_flip() +
  theme_minimal() +
  labs(title = "Forest Plot for Logistic Regression", x = "Predictor Variables", y = "Estimate (Log Odds)")

```

```{r}
# Obliczenie i wykres ROC
roc_curve <- roc(raw_data$V0.PR.INVAVENT_PRTRT.PROCCUR, predict(logistic_model_final, type = "response"))

plot(roc_curve, main = "ROC Curve for Logistic Regression", col = "blue", lwd = 2)
auc_value <- auc(roc_curve)
cat("AUC:", auc_value, "\n")

```
```{r}
ggplot(raw_data, aes(x = `V0.DD.DM.AGE..Wiek`, y = predict(logistic_model_final, type = "response"))) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", color = "blue") +
  labs(title = "Probability of Mechanical Ventilation by Age", x = "Age", y = "Predicted Probability") +
  theme_minimal()

```

```{r}
# Wykres reszt deviance
plot(logistic_model_final$fitted.values, residuals(logistic_model_final, type = "deviance"),
     main = "Deviance Residuals vs Fitted Values",
     xlab = "Fitted Values",
     ylab = "Deviance Residuals",
     col = "darkred", pch = 19)
abline(h = 0, lty = 2)

```