Progetto Gruppo E.Rmd

---
title: "Public Health Project"
author: "Julia Bui Xuan, Simone Farallo, Michele Salvaterra, Luca Sammarini, Eugenio Tarolli Bramè"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


# Preliminary operation

Caricamento dei pacchetti.

```{r}
library(ggplot2)
library(ggpubr)
library(dplyr)
library(epitools)
library(DescTools)
library(survival)
library(survminer)
library(epiR)
```

Il presente progetto si propone di esaminare quattro dataset relativi alla situazione sanitaria di una piccola cittadina tedesca nel periodo compreso tra il 1984 e il 1988. I dataset che andremo ad analizzare sono:

1. estratto del **German health registry** per l’anno 1984 relativo ad una piccola cittadina (dataset reale openData modificato). Le variabili che lo compongono sono le seguenti (e sono relative alla situazione del cittadino all’inizio del 1984):
- idnum: codice identificativo del soggetto
- smoke (yes or no): se il soggetto fuma
- sex (Female or Male)
- married (yes or no): se è coniugato
- kids (yes or no): se il soggetto ha figli
- work (yes or no): se lavora
- education ( no/low = da nessuna a diploma scuola media inferiore;medium/high = diploma scuola media superiore o oltre)
- age (numerica, in anni)

2. **Registro tumori tedesco** (dataset simulato) relativo al mese di Gennaio 1984. Le variabili che lo compongono sono le seguenti:
- idnum: codice identificativo del paziente
- stadio (I, II, III, IV): stadio del tumore alla diagnosi
- incidenza: data di diagnosi del tumore
- tipo di tumore: seno, polmone, colon, altro
- geneticm: fattore genetico (1=positivo, 0=negativo)

3. **Schede di dimissione ospedaliera** dei soggetti ricoverati in Germania tra gennaio 1984 e ottobre 1984 per trattamenti oncologici (dataset simulato):
- idnum: codice identificativo del soggetto
- Prestazione: tipo di trattamento ricevuto durante il ricovero (chirurgica, chemioterapica o radioterapica)
- data prestazione: data del trattamento
- dimissione: data di dimissione dall’ospedale
- ospedale: codice univoco dell’ospedale

4. Estratto del **Registro di mortalità** della cittadina tedesca che riporta la mortalità dal 1984 al 1988 e lo stato in vita alla fine del 1988 (dataset simulato):
- idnum: codice identificativo del cittadino
- dead: stato in vita alla data enddate
- enddate: data di ultima osservazione (se dead=1 data di morte)

# Loading datasets

```{r}
cancer<-read.csv("https://raw.githubusercontent.com/lucasammarini/public-health-german-cancer/main/dataset/Cancerregister.csv", header = TRUE, sep =";")
death<-read.csv("https://raw.githubusercontent.com/lucasammarini/public-health-german-cancer/main/dataset/Deathregister.csv", header = TRUE, sep =";")
GermanH<-read.csv("https://raw.githubusercontent.com/lucasammarini/public-health-german-cancer/main/dataset/GermanH.csv", header = TRUE, sep =";")
sdo<-read.csv("https://raw.githubusercontent.com/lucasammarini/public-health-german-cancer/main/dataset/SDO.csv", header = TRUE, sep =";")
```


```{r}
head(cancer)
head(death)
head(GermanH)
head(sdo)
```

# Cleaning & Exploration

#### 1. Esaminare i dati e riportare le statistiche descrittive in una tabella per ciascun dataset. Per le date riportare data minima e data massima. Fare attenzione alla possibilità di dati mancanti, incongruenze tra date, record ripetuti che potrebbero creare problemi in fase di linkage e analisi. I record con dati ripetuti o incongruenze tra date (eg. Data trattamento precedente alla data d’incidenza) devono essere segnalati nel report e poi eliminati per le analisi successive.

In questa prima fase l'obiettivo è quello di riportare le statistiche descrittive in una tabella per ciascun dataset, con particolare attenzione alla presenza di dati mancanti, incongruenze tra date e record ripetuti che potrebbero creare problemi in fase di linkage e analisi. Saranno segnalati nel report i record con dati ripetuti o incongruenze tra date e saranno eliminati per le analisi successive.

#### Cancer

Controllo dei duplicati.
```{r}
cancer$idnum[duplicated(cancer$idnum)]
cancer[cancer$idnum== 192,]
cancer[cancer$idnum== 363,]
cancer[cancer$idnum== 1933,]
```

Sono presenti tre duplicati; si prosegue alla rimozione.
```{r}
cancer <- cancer[-158, ]
cancer <- cancer[-305,]
cancer <- cancer[-1637, ]
```
Controllo della rimozione dei duplicati.
```{r}
cancer[cancer$idnum== 363,]
cancer[cancer$idnum== 192,]
cancer[cancer$idnum== 1933,]
```

Si convertono gli spazi vuoti in valori nulli, poiché essi sono dati mancanti.
```{r}
cancer <- data.frame(lapply(cancer, function(x) ifelse(x == "", NA, x)), stringsAsFactors = FALSE)
```

Controllo e rimozione dei valori nulli.
```{r}
sum(is.na(cancer))
cancer[!complete.cases(cancer),]
cancer <- na.omit(cancer)
```

Ordinamento dei dati in base alla data di incidenza del tumore.
```{r}
cancer_clean <- cancer[order(cancer$incidenza),]
head(cancer_clean)
```

```{r}
cancer_clean$Stadio <- as.factor(cancer_clean$Stadio)
cancer_clean$tipotumore <- as.factor(cancer_clean$tipotumore)
cancer_clean$geneticm <- as.factor(cancer_clean$geneticm)
cancer_clean$incidenza <- as.Date(cancer_clean$incidenza, tryFormats = "%d/%m/%Y")
```

Controllo ed esplorazione del dataset pulito.
```{r}
str(cancer_clean)
summary(cancer_clean)
```

Esportazione del dataset pulito.

```{r}
write.csv(cancer_clean, "dataset/Cancerregister_clean.csv", row.names= FALSE)
```

Si procede all'esplorazione del dataset attraverso alcuni grafici.

Grafico a barre della distribuzione dei pazienti per stadio del tumore.

```{r}
ggplot(cancer_clean, aes(x=Stadio, fill=Stadio)) +
  geom_bar() +
  labs(title = "Distribuzione dei decessi per Stadio", 
       y = "Numero di pazienti")
```

Grafico a barre della distribuzione dei pazienti per tipo di tumore.

```{r}
ggplot(cancer_clean, aes(x=tipotumore, fill=tipotumore)) +
  geom_bar()+
  labs(title = "Distribuzione dei decessi per tipo di tumore", 
       y = "Numero di pazienti",
       x = "Tipo di tumore")
```
 
Grafico a barre della distribuzione dei pazienti per fattore genetico.

```{r}
ggplot(cancer_clean, aes(x=geneticm, fill=geneticm)) +
  geom_bar()+
  labs(title = "Distribuzione dei decessi per fattore genetico", 
       y = "Numero di pazienti",
       x = "Presenza di fattore genetico (0 = no, 1 = sì)")
```

#### Death

```{r}
str(death)
```

Correzione del formato delle colonne 'dead' e 'enddate'.
```{r}
death$dead = as.factor(death$dead)
death$enddate = as.Date(death$enddate, tryFormats = "%Y-%m-%d" )
str(death)
summary(death)
dim(death)
names(death)
```

Controllo dei dati mancanti.
```{r}
sum(is.na(death))
```

Non sono presenti dati mancanti.

Controllo dei duplicati.

```{r}
sum(duplicated(death))
length(unique(death$idnum)) == length(death$idnum)
```

Non sono presenti duplicati.

Si procede all'esplorazione del dataset attraverso un grafico.

```{r}
ggplot(death, aes(x = dead, fill = dead)) +
  geom_bar() +
  labs(title = "Distribuzione dei decessi", 
       x = "Stato del paziente (0 = vivo, 1 = morto)", 
       y = "Numero di pazienti")
```

#### GermanH

```{r cars}
str(GermanH)
```

Si analizza quanti valori nulli sono presenti nelle singole colonne.

```{r}
colSums(is.na(GermanH))
```

Eliminazione dei valori nulli.

```{r}
GermanH_ok <- na.omit(GermanH)
```

Controllo della riuscita.

```{r}
colSums(is.na(GermanH_ok))
```

Ricerca dei duplicati.

```{r}
nrow(GermanH_ok[duplicated(GermanH_ok$idnum),])
```

Non sono presenti duplicati.

Correzione del formato delle colonne.
```{r}
GermanH_ok$smoke <- as.factor(GermanH_ok$smoke)
GermanH_ok$sex <- as.factor(GermanH_ok$sex)
GermanH_ok$married <- as.factor(GermanH_ok$married)
GermanH_ok$kids <- as.factor(GermanH_ok$kids)
GermanH_ok$work <- as.factor(GermanH_ok$work)
GermanH_ok$education <- as.factor(GermanH_ok$education)
str(GermanH_ok)
summary(GermanH_ok)
```

Esportazione del dataset pulito.

```{r}
write.csv(GermanH_ok, "dataset/GermanH_clean.csv", row.names= FALSE)
```

Esplorazione del dataset pulito attraverso grafici.

```{r}
colors <- c("#FFA07A", "#ADD8E6")
```

```{r}
barplot(table(GermanH_ok$smoke), main="Numero di pazienti fumatori", xlab="Fumatore", ylab="Numero dei pazienti", col= colors)
```

```{r}
barplot(table(GermanH_ok$sex), main="Numero di pazienti per genere", xlab="Genere", ylab="Numero dei pazienti", col = colors)
```

```{r}
barplot(table(GermanH_ok$married), main="Numero di pazienti per stato sociale", xlab="Coniugato", ylab="Numero dei pazienti", col = colors)
```

```{r}
barplot(table(GermanH_ok$kids), main="Numero di pazienti per presenza di figli", xlab="Presenza di figli", ylab="Numero dei pazienti", col= colors)
```

```{r}
barplot(table(GermanH_ok$work), main="Numero di pazienti per disoccupazione", xlab="Occupazione", ylab="Numero dei pazienti", col = colors)
```

```{r}
ggplot(GermanH_ok, aes(age)) +
  geom_bar(fill = "#0073C2FF")+xlab('Età')+ylab("Numero dei pazienti")+
  labs(title = "Distribuzione dei pazienti per età")
```

#### SDO

```{r}
str(sdo)
```

Ordinamento del dataset.
```{r}
sdo <- sdo[order(sdo$idnum), ]
```

Controllo dei duplicati.
```{r}
nrow(sdo[duplicated(sdo$idnum),])
```

Formattazione delle colonne.
```{r}
sdo$Prestazione <- as.factor(sdo$Prestazione)
sdo$ospedale <- as.factor(sdo$ospedale)
sdo$dataprestazione <- as.Date(sdo$dataprestazione, '%d/%m/%Y')
sdo$dimissione <- as.Date(sdo$dimissione, '%d/%m/%Y')
```


```{r}
summary(sdo)
```

Controllo dei valori nulli.

```{r}
nrow(sdo[is.na(sdo$Prestazione),])
nrow(sdo[is.na(sdo$dataprestazione),])
sdo[is.na(sdo$dataprestazione),]
nrow(sdo[is.na(sdo$dimissione),])
sdo[is.na(sdo$dimissione),]
nrow(sdo[is.na(sdo$ospedale),])
```

Eliminazione dei valori nulli.

```{r}
sdo.no.na <- sdo[!is.na(sdo$dataprestazione) & !is.na(sdo$dimissione), ]
nrow(sdo.no.na[is.na(sdo.no.na$Prestazione),])
nrow(sdo.no.na[is.na(sdo.no.na$dataprestazione),])
nrow(sdo.no.na[is.na(sdo.no.na$dimissione),])
nrow(sdo.no.na[is.na(sdo.no.na$ospedale),])
```

Controllo della presenza di incongruenze tra date e rimozione dei valori incongruenti.

```{r}
nrow(sdo.no.na[sdo.no.na$dataprestazione>sdo.no.na$dimissione, ])
sdo.no.na[sdo.no.na$dataprestazione>sdo.no.na$dimissione, ]
sdo.ok <- sdo.no.na[sdo.no.na$dataprestazione<sdo.no.na$dimissione, ]
```

Statistiche del dataset pulito.

```{r}
summary(sdo.ok)
```

Esportazione del dataset pulito.

```{r}
write.csv(sdo.ok, "dataset/SDO_clean.csv", row.names= FALSE)
```

Esplorazione del dataset pulito attraverso grafici.

```{r}
ggplot(sdo, aes(Prestazione)) +
  geom_bar(fill = "#0073C2FF") +
  theme_pubclean()+
  ylab("Numero dei pazienti")+
  labs(title = "Numero di pazienti per prestazione")
```

```{r}
ggplot(sdo, aes(ospedale)) +
  geom_bar(fill = "#0073C2FF") +
  theme_pubclean()+
  ylab("Numero dei pazienti")+
  xlab("Ospedali")+
  labs(title = "Numero di pazienti per ospedale")
```

I dataset sono stati tutti puliti e controllati ed è stata svolta un'analisi esplorativa iniziale. Ora si può passare alla fase di linkage dei dataset per eseguire gli altri task.

### Prima parte: punti da 2 a 6

#### 2. Effettuare il record-linkage con lo scopo di costruire l’indicatore ‘Intervento chirurgico di asportazione del tumore al seno entro 60 giorni dalla data di diagnosi’ su base mensile per i casi incidenti nel mese di gennaio 1984.

Denominatore: tutte le pazienti di sesso femminile con tumore al seno insorto tra 01/01/1984 e 31/01/1984, in stadio I o II, che hanno subito un intervento chirurgico.

Numeratore: tutte le pazienti al denominatore con intervallo tra la data d'incidenza e la data dell'intervento ≤ 60 giorni.

```{r}
data_sdo = read.csv("dataset/SDO_clean.csv", header = TRUE, sep =",")

data_cancer = read.csv("dataset/Cancerregister_clean.csv", header=TRUE, sep = ",")

dataMerge=merge(data_cancer, data_sdo, by="idnum")

dataMerge$incidenza = as.Date(dataMerge$incidenza, format = "%Y-%m-%d")
dataMerge$dataprestazione= as.Date(dataMerge$dataprestazione)
dataMerge$dimissione= as.Date(dataMerge$dimissione)

data_gerH = read.csv("dataset/GermanH_clean.csv", header=TRUE, sep = ",")
dataMerge2 = merge(dataMerge, data_gerH, by='idnum')

#pulizia
dataMerge2$Prestazione <- as.factor(dataMerge2$Prestazione)
dataMerge2$ospedale <- as.factor(dataMerge2$ospedale)
dataMerge2$Stadio <- as.factor(dataMerge2$Stadio)
dataMerge2$tipotumore <- as.factor(dataMerge2$tipotumore)
dataMerge2$geneticm <- as.factor(dataMerge2$geneticm)
dataMerge2$smoke <- as.factor(dataMerge2$smoke)
dataMerge2$sex <- as.factor(dataMerge2$sex)
dataMerge2$married <- as.factor(dataMerge2$married)
dataMerge2$kids <- as.factor(dataMerge2$kids)
dataMerge2$work <- as.factor(dataMerge2$work)
dataMerge2$education <- as.factor(dataMerge2$education)

summary(dataMerge2)
```

Per costruire l'indicatore, il denominatore deve soddisfare i seguenti criteri di inclusione:

1. sex: female
2. tipotumore: seno
3. Stadio: Stadio I, Stadio II
4. incidenza: dal 1984-01-01 al 1984-01-31
5. prestazione: chirurgica

Al numeratore invece, verranno incluse solamente le pazienti che hanno subito un intervento entro 60 giorni dalla data d'incidenza. Verrà quindi creata una nuova variabile binaria per indicare la presenza dell'evento o meno.

```{r}
#denominatore
data = dataMerge2[dataMerge2$sex == 'Female',]
data = data[data$Stadio == 'Stadio I' | data$Stadio == 'Stadio II',]
data = data[data$tipotumore=='seno', ]
data = data[data$Prestazione=='chirurgica',]
summary(data)
dim(data)
```

La variabile 'incidenza' soddisfa già i criteri di inclusione. Come possiamo notare dalle statistiche descrittive, la data è compresa tra l'11/01/1984 e il 20/01/1984.

Al numeratore verranno incluse tutte le pazienti del denominatore che hanno subito l'intervento entro 60 giorni. 

Attraverso il seguente codice, si aggiunge la variabile "intervallo" per calcolare il numero di giorni fra la data di prestazione e la data di incidenza.

```{r}
data$intervallo = as.numeric(data$dataprestazione - data$incidenza)
summary(data$intervallo)
```
Osservando la nuova variabile 'intervallo' possiamo notare che il tempo minimo di attesa per la prestazione è pari a 9 giorni, invece il tempo massimo di attesa equivale a 193 giorni.

In seguito, viene creata la variabile binaria 'intervento60gg' per individuare le osservazioni che hanno subito o meno l'intervento entro 60 giorni.


```{r}
data$intervento60gg = ifelse(data$intervallo <= 60, 1, 0)
T1 = table(data$intervento60gg); T1 #indicatore tumore seno
prop.table(T1)
indicatore = T1[2]/(T1[2]+T1[1]) #
indicatore
```
 L'indicatore è pari a 0.55. La percentuale di interventi entro 60 giorni dalla data di diagnosi è pari al 55%. 
 
#### 3. Calcolare l’indicatore ‘Intervento chirurgico di asportazione del tumore al seno entro 60 giorni dalla data di diagnosi’ per ospedale e darne rappresentazione grafica, includendo come valore di riferimento nel grafico l’indicatore calcolato sull’intero dataset. Per esempi relativi alla rappresentazione grafica fare riferimento al sito Piano Nazionale Esiti (PNE) o siti analoghi trattati a lezione.

```{r}
table(data$ospedale)
TT<-table(evento = data$intervento60gg, ospedale = data$ospedale);TT 
prop.table(TT,2)


data$intervento60gg = as.factor(data$intervento60gg)
ggplot(data, aes(x = ospedale, fill = intervento60gg)) +
  geom_bar(position = "fill") +
  xlab("Ospedale") +
  ylab("Frazione eventi") +
  ggtitle("Frazione di eventi per ospedale")
```

Come si evince dal grafico sovrastante, gli ospedali 2, 7 e 8 sono quelli con percentuale maggiore relativa alle osservazioni che hanno subito l'evento, rispettivamente pari a 64%, 62% e 63%. L'ospedale 1, invece, presenta la percentuale minore di osservazioni che hanno subito l'evento, pari a 46%.

Nella tabella seguente vengono riportati gli indicatori per ospedale.

```{r}
indice <- data %>%
  group_by(ospedale)%>%
  summarise(indice=sum(intervento60gg == 1)/n())
indice
```

```{r}
bar_plot = ggplot(indice, aes(x = reorder(ospedale, indice), y = indice)) +
  geom_bar(stat = "identity", fill = "darkgreen") +
  xlab("Ospedale") +
  ylab("Indicatore") +
  ggtitle("Indicatore per ospedale") +
  theme(plot.title = element_text(size = 14, hjust = 0.5),
        axis.text.x = element_text(angle = 0, vjust = 0, hjust=0.5)) +
  geom_hline(yintercept = indicatore, linetype='dashed', color='black')
  
bar_plot
```

Considerando il valore dell'indicatore di riferimento, calcolato in precedenza sull'intero dataset, pari a 0.55, gli ospedali 2, 5, 7, 8 presentano un valore superiore. Invece, gli ospedali 1, 3, 4, 6, 9 presentano un valore inferiore.

#### 4. Utilizzare il dataset ottenuto per valutare l’associazione a livello individuale tra il livello di educazione ed il valore dell’indicatore ‘Intervento chirurgico di asportazione del tumore al seno entro 60 giorni dalla data di diagnosi’. Quale misura di effetto è possibile stimare? Calcolate ed interpretate tale misura di effetto grezza. Riportare anche la relativa tabella di contingenza.


Per valutare l'associazione tra il livello di educazione e l'intervento chirurgico di asportazione del tumore al seno entro 60 giorni dalla diagnosi, è possibile utilizzare una tabella di contingenza dove vengono riportati i valori che indicano il numero di casi che rientrano in ogni combinazione di categoria di educazione e di intervento chirurgico.


```{r}
tabella_contingenza <- table(data$education, data$intervento60gg)
tabella_contingenza
OR <- (tabella_contingenza[1,1]/tabella_contingenza[1,2])/(tabella_contingenza[2,1]/tabella_contingenza[2,2])
OR
```

La misura di effetto che è possibile stimare in questo caso è l'ODDS Ratio, esso ci permette di confrontare le probabilità di un evento tra due gruppi. In particolare, in questo caso si vuole confrontare la probabilità che una donna con un basso livello di educazione abbia un intervento chirurgico di asportazione del tumore al seno entro 60 giorni dalla diagnosi con la probabilità che una donna con un medio/alto livello di educazione abbia lo stesso intervento entro lo stesso periodo.

```{r}
epitab(data$education, data$intervento60gg, method = c("oddsratio")) 
```

Come si evince dal risultato sovrastante, l'ODDS ratio è pari a 0.97, con un intervallo di confidenza (0.51, 1.87); poiché l'intervallo di confidenza include il valore 1, non è possibile concludere che l'ODDS ratio sia significativamente diverso da 1.

L'ODDS ratio può essere calcolato anche attraverso un modello logistico:

```{r}
modello_logistico_seno <- glm(intervento60gg ~ education, family = binomial(), data = data)
# Riepilogo del modello
summary(modello_logistico_seno)
#Lettura Coefficienti
exp(cbind("OR" = coef(modello_logistico_seno), confint.default(modello_logistico_seno, level = 0.95)))
```

I risultati coincidono nei tre metodi utilizzati.

#### 5. Calcolate la stessa misura di effetto, questa volta aggiustata per la sola variabile ‘working’, mediante il metodo Mantel Haenszel. Interpretate il risultato.


Per calcolare l'OR tramite il metodo Mantel Haenszel viene costruita una tabella di contingenza per i due gruppi in base allo stato lavorativo. 

```{r}
tabella_contingenza_1 <- table(education = data$education, evento = data$intervento60gg, work = data$work)
tabella_contingenza_1
```

Per verificare la possibilità di utilizzare l'OR di Mantel Haenszel viene utilizzato il test di Breslow-Day, il quale testa l'omogeneità degli ODDS ratio nei due gruppi.

```{r}
BreslowDayTest(tabella_contingenza_1)
```
 
Il test non rifiuta l'ipotesi nulla di omogeneità; si può proseguire nel calcolo dell'OR.

```{r}
OR_adj <- mantelhaen.test(tabella_contingenza_1, conf.level = 0.95)
OR_adj
```
L'OR di Mantel Haenszel non è significativamente diverso da 1, controllando per la variabile 'work'. Inoltre, il rapporto fra l'ODDS ratio crudo e quello di Mantel Haenszel è pari a circa 1.00. Questo significa che lo stato lavorativo non è un confondente nell'associazione tra il livello di educazione e l'intervento chirurgico di asportazione del tumore al seno entro 60 giorni dalla data di diagnosi.

#### 6. Stimate l’associazione a livello individuale tra il livello di educazione ed il valore dell’indicatore, aggiustata per tutte le variabili disponibili che ritenete opportuno inserire come potenziali confondenti, mediante un modello di regressione logistica. Interpretate i risultati. Su quanti soggetti avete effettuato l’analisi? Quali variabili sono associate all’indicatore? In che modo?

```{r}
modello <- glm(intervento60gg ~ education, data = data, family = binomial())
summary(modello)
exp(cbind("OR" = coef(modello), confint.default(modello, level = 0.95)))
```

Una delle variabili che si é ritenuto considerare come confondente è "geneticm".

```{r}
modello_gen<- glm(intervento60gg ~ education + geneticm, data = data, family = binomial())
summary(modello_gen)
exp(cbind("OR" = coef(modello_gen), confint.default(modello_gen, level = 0.95)))
ratio_gen = exp(coef(modello)['educationmedium/high'])/exp(coef(modello_gen)['educationmedium/high'])
ratio_gen
```

Ciò nonostante, il rapporto fra l'OR crudo e quello di MH è pari a circa 1.04. Ciò significa che non è un confondente.

```{r}
modello_stadio <- glm(intervento60gg ~ education + Stadio , data = data, family = binomial())
summary(modello_stadio)
exp(cbind("OR" = coef(modello_stadio), confint.default(modello_stadio, level = 0.95)))
ratio_stadio = exp(coef(modello)['educationmedium/high'])/exp(coef(modello_stadio)['educationmedium/high'])
ratio_stadio
```

Un'altra variabile presa in considerazione come possibile confondente è "Stadio". Si osserva che, anche in questo caso, il rapporto tra fra l'OR crudo e quello di MH è pari a circa 1.05. Ciò significa che la variabile non è un confondente.

Infine,  si va a controllare per la variabile 'age' suddivisa in 5 intervalli. 

```{r}
classe_1 <- 28
classe_2 <- 42
classe_3 <- 56
classe_4 <- 70
classe_5 <- 84
classe_6 <- 98
data$age_class <- as.factor(cut(data$age, breaks = c(classe_1, classe_2, classe_3, classe_4, classe_5, classe_6), labels = FALSE))
modello_age <- glm(intervento60gg ~ education + age_class, data = data, family = binomial())
summary(modello_age)
exp(cbind("OR" = coef(modello_age), confint.default(modello_age, level = 0.95)))
ratio_age = exp(coef(modello)['educationmedium/high'])/exp(coef(modello_age)['educationmedium/high'])
ratio_age
```

Anche in questo caso, la variabile scelta, 'age', non è un confondente.

L'analisi è stata eseguita su tutte le persone del dataset.

```{r}
dim(data)
```

Le persone del dataset sono 311.

Nessuna delle variabili analizzate risulta essere un confondente per la variabile "education".
La variabile "education", inoltre, non è associata all'indicatore.
Tuttavia, la variabile "Stadio" risulta essere positivamente associata all'indicatore: l'OR risulta significativamente superiore a 1.

### Seconda parte: punti da 7 a 12

#### Considerate ora tutti i tumori al colon insorti nel gennaio 1984. Unire i data-set utili per studiare la mortalità del tumore al colon nei soggetti inclusi nell’estrazione del German health Register (dataset 1).

```{r}
data_death = read.csv("dataset/Deathregister.csv", header = TRUE, sep =";")
data_cancer = read.csv("dataset/Cancerregister_clean.csv", header=TRUE, sep = ",")
dataMerge3 = merge(data_cancer, data_death, by="idnum")

data_gerH = read.csv("dataset/GermanH_clean.csv", header=TRUE, sep = ",")
dataMerge4 = merge(dataMerge3, data_gerH, by='idnum')

## Pulizia
dataMerge4$incidenza <- as.Date(dataMerge4$incidenza, format = "%Y-%m-%d")
dataMerge4$enddate <- as.Date(dataMerge4$enddate, format = "%Y-%m-%d")
dataMerge4$Stadio <- as.factor(dataMerge4$Stadio)
dataMerge4$tipotumore <- as.factor(dataMerge4$tipotumore)
dataMerge4$geneticm <- as.factor(dataMerge4$geneticm)
dataMerge4$smoke <- as.factor(dataMerge4$smoke)
dataMerge4$sex <- as.factor(dataMerge4$sex)
dataMerge4$married <- as.factor(dataMerge4$married)
dataMerge4$kids <- as.factor(dataMerge4$kids)
dataMerge4$work <- as.factor(dataMerge4$work)
dataMerge4$education <- as.factor(dataMerge4$education)
dataMerge4$dead <- as.factor(dataMerge4$dead)

summary(dataMerge4)
```

#### 7. Selezionate i record relativi ai tumori al colon e stimate la sopravvivenza a 5 anni. Quanti soggetti sono inclusi nell’analisi? Quanti pazienti sono morti nel periodo di interesse? Riportare graficamente la stima di sopravvivenza nei primi 5 anni dalla diagnosi stimata tramite lo stimatore di Kaplan-Meier. Stimare approssimativamente la sopravvivenza mediana.

```{r}
colon <- dataMerge4[dataMerge4$tipotumore=="colon",]
length(colon[colon$incidenza>colon$enddate])
```

```{r}
dim(colon)
```

Nell'analisi sono inclusi 1304 soggetti.

```{r}
table(colon$dead)
```

Nel periodo di interesse sono morti 707 pazienti.

Si stima la sopravvivenza di Kaplan-Meier, che viene rappresentata graficamente.

```{r}
colon$survtime <- as.numeric(colon$enddate-colon$incidenza)/365.25

fit<-survfit(Surv(survtime, as.numeric(dead)) ~1,data=colon)
summary(fit)

ggsurvplot(fit, data = colon, risk.table = TRUE, conf.int = TRUE, conf.int.fill = "black", conf.int.style = "step", surv.median.line = "h", censor.size=2.5, risk.table.fontsize = 3.5, censor.shape=".", xlab="Time (years)")
```

Nel grafico è riportata la sopravvivenza nei primi 5 anni calcolata con lo stimatore di Kaplan-Meier.

La sopravvivenza mediana, come si osserva nel grafico, è poco inferiore a 4 anni.

#### 8. Stimare la sopravvivenza nei primi 5 anni dalla diagnosi per Stadio e effettuare un test d’ipotesi per verificare se l’azzardo di morte sia diverso per stadio di malattia alla diagnosi.

Si stima la sopravvivenza di Kaplan-Meier per stadio, che viene rappresentata graficamente.

```{r}
fit<-survfit(Surv(survtime, as.numeric(dead)) ~ Stadio,data=colon)
summary(fit)

ggsurvplot(fit, data = colon, risk.table = TRUE, conf.int = TRUE, surv.median.line = "h", pval = TRUE, censor.size=2.5, risk.table.fontsize = 3, risk.table.height = 0.35, censor.shape="|", xlab="Time (years)", legend = "right")
```

Dal grafico della sopravvivenza di Kaplan-Meier risulta piuttosto evidente una differenza sul trend della sopravvivenza per stadio.

Si esegue un test di ipotesi log-rank per valutare se la differenza è significativa.

```{r}
survdiff(Surv(survtime,as.numeric(dead)) ~ Stadio,data=colon)
```

Il test ha un p-value inferiore a 2e-16, dunque si rigetta l'ipotesi nulla: la differenza tra gli azzardi è statisticamente significativa.

#### 9. Applicare un modello per valutare l’associazione tra sesso e mortalità e interpretare la misura di effetto stimata.

Si applica un modello di regressione logistica per valutare l'associazione.

```{r}
modello <- glm(dead ~ relevel(sex, "Male"), data = colon, family = binomial())
summary(modello)
exp(cbind("OR" = coef(modello), confint.default(modello, level = 0.95)))
```

La misura di effetto ricavata dal modello è l'ODDS Ratio, pari a 0.65. Ciò significa che gli uomini hanno un valore di ODDS di morte pari al 65% rispetto a quello delle donne.

#### 10. Quali variabili sono associate alla mortalità? Riportare le relative stime di effetto con gli intervalli di confidenza.

Per studiare l'associazione tra variabili e mortalità si applica un modello di regressione logistica usando come covariate tutte le variabili del dataset.

```{r}
modello10 <- glm(dead ~ relevel(sex, "Male") + Stadio + geneticm + smoke + married + kids + work + education + age, data = colon, family = binomial())
summary(modello10)
exp(cbind("OR" = coef(modello10), confint.default(modello10, level = 0.95)))
```

Tutte le variabili testate, ad eccezione di "married", "kids" e "work", hanno una associazione significativa con la mortalità, e quindi un OR significativamente diverso da 1.

#### 11. Valutare la presenza di confondenti e/o modificatori di effetto tra le variabili disponibili nel German health register e nel registro tumori nella valutazione dell’associazione tra sesso e mortalità. Se identificate un’interazione tra sesso e un’altra variabile riportare le stime di effetto per maschi e femmine separatamente e commentare il tipo di interazione trovato.

Si studia la correlazione tra morte e sesso, valutando la presenza di confondenti o modificatori tramite stratificazione.

La prima variabile studiata come eventuale confondente o modificatore di effetto è il fattore genetico.

```{r}
cont_table_geneticm <- table(sex = colon$sex,
                             dead = relevel(colon$dead, "1"),
                             genetic = colon$geneticm)
strat_geneticm<-epi.2by2(dat=cont_table_geneticm , method="cohort.count")
strat_geneticm$massoc.detail$RR.strata.wald
strat_geneticm$massoc.detail$OR.strata.wald
```

Si osserva una differenza sia tra gli ODDS Ratio che tra i Rischi Relativi. Si verifica se la differenza è significativa tramite il test di omogeneità.

```{r}
strat_geneticm$massoc.detail$wRR.homog
strat_geneticm$massoc.detail$wOR.homog
```

Il test di omogeneità risulta significativo per i Rischi Relativi: il fattore genetico è un modificatore d'effetto per questa misura. Il test di omogeneità sugli OR risulta non significativo.

```{r}
strat_geneticm
```

Il rapporto tra l'OR crudo e l'OR di MH è pari a 1: il fattore genetico non è un confondente per l'OR.

Si valutano i RR separatamente per sesso.

```{r}
cont_table_geneticm_male <- table(genetic = relevel(colon$geneticm, "1")[colon$sex=="Male"],
                                  dead = relevel(colon$dead, "1")[colon$sex=="Male"])
strat_geneticm_male<-epi.2by2(dat=cont_table_geneticm_male, method="cohort.count")
strat_geneticm_male

cont_table_geneticm_female <- table(genetic = relevel(colon$geneticm, "1")[colon$sex=="Female"],
                                    dead = relevel(colon$dead, "1")[colon$sex=="Female"])
strat_geneticm_female<-epi.2by2(dat=cont_table_geneticm_female, method="cohort.count")
strat_geneticm_female
```

L'effetto del fattore genetico negli uomini è minore rispetto alle donne: se si considerano gli uomini come classe di riferimento, l'interazione è positiva.

Si studia ora la variabile "smoke".

```{r}
cont_table_smoke <- table(sex = colon$sex,
                          dead = relevel(colon$dead, "1"),
                          smoke = colon$smoke)
strat_smoke<-epi.2by2(dat=cont_table_smoke , method="cohort.count")
strat_smoke$massoc.detail$RR.strata.wald
strat_smoke$massoc.detail$OR.strata.wald
```

In questo caso entrambe le misure di effetto sono molto simili in entrambi gli strati. Ci si aspetta che i test non siano significativi.

```{r}
strat_smoke$massoc.detail$wRR.homog
strat_smoke$massoc.detail$wOR.homog
```

In entrambi i casi i test non sono significativi: la variabile non è un modificatore di effetto ed è possibile calcolare una misura di effetto comune tramite il metodo di Mantel Haenszel.

```{r}
strat_smoke
```

Dal rapporto tra l'OR crudo e quello di MH, si evince che la variabile "smoke" non è un confondente per questa misura di effetto. Analogamente accade per il RR. 

Si procede ora alla valutazione della variabile "married".

```{r}
cont_table_married <- table(sex = colon$sex,
                            dead = relevel(colon$dead, "1"),
                            married = colon$married)
strat_married<-epi.2by2(dat=cont_table_married , method="cohort.count")
strat_married$massoc.detail$RR.strata.wald
strat_married$massoc.detail$OR.strata.wald
strat_married$massoc.detail$wRR.homog
strat_married$massoc.detail$wOR.homog
strat_married
```

Si osserva che la variabile non è un confondente né un modificatore di effetto per nessuna delle due misure di effetto calcolate.

Si procede ora alla valutazione della variabile "kids".

```{r}
cont_table_kids <- table(sex = colon$sex,
                         dead = relevel(colon$dead, "1"),
                         kids = colon$kids)
strat_kids<-epi.2by2(dat=cont_table_kids , method="cohort.count")
strat_kids$massoc.detail$RR.strata.wald
strat_kids$massoc.detail$OR.strata.wald
strat_kids$massoc.detail$wRR.homog
strat_kids$massoc.detail$wOR.homog
strat_kids
```

La variabile non è né confondente né modificatore di effetto.

Si procede ora alla valutazione della variabile "work".

```{r}
cont_table_work <- table(sex = colon$sex,
                         dead = relevel(colon$dead, "1"),
                         work = colon$work)
strat_work<-epi.2by2(dat=cont_table_work , method="cohort.count")
strat_work$massoc.detail$RR.strata.wald
strat_work$massoc.detail$OR.strata.wald
strat_work$massoc.detail$wRR.homog
strat_work$massoc.detail$wOR.homog
strat_work
```

La variabile non è né confondente né modificatore di effetto.

Si procede ora alla valutazione della variabile "education".

```{r}
cont_table_education <- table(sex = colon$sex,
                              dead = relevel(colon$dead, "1"),
                              education = colon$education)
strat_education<-epi.2by2(dat=cont_table_education , method="cohort.count")
strat_education$massoc.detail$RR.strata.wald
strat_education$massoc.detail$OR.strata.wald
strat_education$massoc.detail$wRR.homog
strat_education$massoc.detail$wOR.homog
```

Il test sull'omogeneità dei Rischi Relativi è significativo: la variabile è un modificatore di effetto per il RR.

```{r}
strat_education
```

"education" non è un confondente per l'OR.

Si procede a stimare il RR per maschi e femmine separatamente.

```{r}
cont_table_education_male <- table(education = relevel(colon$education, "medium/high")[colon$sex=="Male"],
                                   dead = relevel(colon$dead, "1")[colon$sex=="Male"])
strat_education_male<-epi.2by2(dat=cont_table_education_male, method="cohort.count")
strat_education_male

cont_table_education_female <- table(education = relevel(colon$education, "medium/high")[colon$sex=="Female"],
                                     dead = relevel(colon$dead, "1")[colon$sex=="Female"])
strat_education_female<-epi.2by2(dat=cont_table_education_female, method="cohort.count")
strat_education_female
```

L'effetto dell'educazione negli uomini è minore rispetto alle donne: se si considerano gli uomini come classe di riferimento, l'interazione è positiva.

Si procede ora alla valutazione della variabile "Stadio".

```{r}
cont_table_stadio <- table(sex = colon$sex,
                    dead = relevel(colon$dead, "1"),
                    stadio = colon$Stadio)
strat_stadio<-epi.2by2(dat=cont_table_stadio , method="cohort.count")
strat_stadio$massoc.detail$RR.strata.wald
strat_stadio$massoc.detail$OR.strata.wald
strat_stadio$massoc.detail$wRR.homog
strat_stadio$massoc.detail$wOR.homog
strat_stadio
```

I test di omogeneità risultano significativi per entrambe le misure di effetto. Si procede al calcolo del valore dell'effetto separatamente per maschi e femmine.

```{r}
model_bin_male <- glm(dead ~ Stadio, data=colon[colon$sex=="Male",], family="binomial")
summary(model_bin_male)
exp(cbind("ODDS ratio" = coef(model_bin_male), confint.default(model_bin_male, level = 0.95)))

model_bin_female <- glm(dead ~ Stadio, data=colon[colon$sex=="Female",], family="binomial")
summary(model_bin_female)
exp(cbind("ODDS ratio" = coef(model_bin_female), confint.default(model_bin_female, level = 0.95)))
```

Gli uomini hanno un OR minore tra Stadio 1 e 2 e tra Stadio 1 e 3 rispetto alle donne. Hanno invece un OR maggiore tra Stadio 1 e 4 rispetto alle donne.

Si procede alla stima dei RR.

```{r}
model_pois_male <- glm(I(as.numeric(dead)-1) ~ Stadio, data = colon[colon$sex=="Male",], family = poisson)
exp(cbind("Relative risk" = coef(model_pois_male), confint.default(model_pois_male, level = 0.95)))

model_pois_female <- glm(I(as.numeric(dead)-1) ~ Stadio, data = colon[colon$sex=="Female",], family = poisson)
exp(cbind("Relative risk" = coef(model_pois_female), confint.default(model_pois_female, level = 0.95)))
```

Gli uomini hanno in tutti i casi un RR minore rispetto alle donne: prendendo come riferimento il livello "Male", si ha dunque una interazione positiva.

Si procede ora alla valutazione della variabile "age".

```{r}
model <- glm(dead ~ relevel(sex, "Male") + age, data=colon, family="binomial")
summary(model)
model2 <- glm(dead ~ relevel(sex, "Male") * age, data=colon, family="binomial")
summary(model2)
anova(model,model2,test="LRT")
```

Si osserva che l'interazione non è significativa, secondo il test LRT: age non è un modificatore di effetto.

```{r}
model3 <- glm(dead ~ relevel(sex, "Male"), data=colon, family="binomial")
or_crude <- exp(cbind("Crude ODDS ratio" = coef(model3), confint.default(model3, level = 0.95)))
or_mh <- exp(cbind("MH ODDS ratio" = coef(model), confint.default(model, level = 0.95)))
ratio <- or_crude[2]/or_mh[2]
ratio
```

Il rapporto tra i due OR è pari a 1.08. Tenendo come valore limite per il confondimento uno scostamento del 10%, è possibile affermare che non c'è confondimento da parte della variabile "age".

Risultano quindi modificatori di effetto le variabili "geneticm", "education" e "Stadio".

#### 12. A seguito delle considerazioni effettuate nei punti precedenti scegliete un modello finale per valutare i fattori di rischio della mortalità dopo diagnosi di tumore al colon e commentate i risultati.

Poiché lo studio è uno studio di mortalità, si utilizza il modello di Cox per tenere conto anche dei diversi tempi di mortalità. Come covariate si scelgono le variabili risultate significative nel punto 10.

```{r}
model_cox <- coxph(Surv(survtime, I(as.numeric(dead)-1)) ~ relevel(sex, "Male") + geneticm + education + Stadio + smoke + age, data=colon)
summary(model_cox)
```

I test di ipotesi riportano che il modello di Cox è significativo. In particolare, risultano significative le variabili "sex", "geneticm", "Stadio", "age"; "education" è al limite della significatività; "smoke" risulta non significativamente associata.

L'Hazard Ratio risulta significativamente superiore a 1 per le variabili "geneticm", "Stadio III", "Stadio IV" e "age". Queste variabili sono positivamente associate con il rischio di mortalità. L'Hazard Ratio risulta invece significativamente inferiore a 1 per la variabile "sex". Le restanti variabili hanno un Hazard Ratio non significativamente diverso da 1. 

Il modello di Cox assume che l'Hazard Ratio rimanga costante nel tempo. Si verifica questo assunto.

```{r}
ph_test <- cox.zph(model_cox)
print(ph_test)
par(mfrow=c(2,1), mar=c(4,5,1,1))
plot(ph_test)
```

L'assunto è verificato per tutte le variabili prese singolarmente e per il modello nel suo complesso. Infatti, dai dati si può notare come il p-value non risulti significativo.

Si inseriscono ora nel modello le modifiche di effetto studiate nel punto 11.

```{r}
model_cox_2 <- coxph(Surv(survtime, I(as.numeric(dead)-1)) ~ relevel(sex, "Male") * geneticm + relevel(sex, "Male") * education + relevel(sex, "Male") * Stadio + smoke + age, data=colon)
summary(model_cox_2)
```

I test di ipotesi riportano che il modello di Cox è significativo.
In particolare, risultano significative le variabili "sex", "geneticm", "Stadio" per il quarto livello e "age".
"education" e "smoke" risultano non significativamente associate.
In merito alle interazioni inserite nel modello, risultano significative quelle tra "sex" e "Stadio II" e tra "sex" e "Stadio III".

L'Hazard Ratio risulta significativamente superiore a 1 per le variabili "geneticm", "Stadio IV" e "age" e per le interazioni tra "sex" e "Stadio II" e tra "sex" e "Stadio III". Queste variabili sono positivamente associate con il rischio di mortalità. L'Hazard Ratio risulta invece significativamente inferiore a 1 per la variabile "sex". Le restanti variabili hanno un Hazard Ratio non significativamente diverso da 1.

```{r}
anova(model_cox, model_cox_2, test="LRT")
```

I due modelli sono significativamente diversi.

```{r}
ph_test_2 <- cox.zph(model_cox_2)
print(ph_test_2)
```

Anche in questo caso l'assunto di costanza dell'Hazard Ratio risulta verificato.

```{r}
par(mfrow=c(3,1), mar=c(4,5,1,1), cex.lab = 0.8)
plot(ph_test_2)
```