Projet pentinat guyot econométrie du risque

#GUYOT EMILIE
#PENTINAT ANTOINE

# On set la seed pour la reproductibilité
set.seed(1234)

#################  Chargement de la base de données  ###########################

mutuelle<-read.csv("C:/Users/Antoine/OneDrive/Documents/Lille/econométrie du risque/Duree_Mutuelle (1).csv", sep=",", dec = ".")

################################################################################


################################################################################

#                     Analyse descriptive des données

################################################################################

##############   STATISTIQUES DESCRIPTIVES UNIVARIEES    ########################

summary(mutuelle)
head(mutuelle)
any(is.na(mutuelle)) #verifier si il y a des valeurs manquantes

############### Analyse pour les variables quantitatives ##############

####### AGE #######
moyenne_age<-mean(mutuelle$AGE)
moyenne_age

e_t_age<-sd(mutuelle$AGE)
e_t_age

min_age<-min(mutuelle$AGE)
min_age
 
max_age<-max(mutuelle$AGE)
max_age

mediane_age<-median(mutuelle$AGE)
mediane_age

#### DUREE #####

#Création de la durée des contrats
#mutuelle$FINOBS <- as.numeric(mutuelle$FINOBS)
mutuelle$PROF <- as.factor(mutuelle$PROF)

mutuelle$DUREE<-mutuelle$FINOBS-mutuelle$DEBUT

#Analyse
moyenne_duree<-mean(mutuelle$DUREE)
moyenne_duree

e_t_duree<-sd(mutuelle$DUREE)
e_t_duree

min_duree<-min(mutuelle$DUREE)
min_duree

max_duree<-max(mutuelle$DUREE)
max_duree

mediane_duree<-median(mutuelle$DUREE)
mediane_duree

############### Analyse pour les variables qualitatives ##############

#### PROF ####
prof<- table(mutuelle$PROF)
prof

prof_prop <- prop.table(prof)*100
prof_prop

#### PROP ####

prop<-table(mutuelle$PROP)
prop

prop_table<-prop.table(prop)*100
prop_table

#### HOMME/FEMME ####

homme_femme<-table(mutuelle$HOMME)
homme_femme

homme_femme_table<-prop.table(homme_femme)*100
homme_femme_table

### GRAPHIQUES ###

#Histogramme de l'âge
hist(mutuelle$AGE, main = "Age des adhérents", xlab = "Age", col = "lightblue", breaks = 5)

#Histogramme de la durée des contrats
hist(mutuelle$DUREE, main = "Durée des Contrats", xlab = "Durée (en années)", col = "lightcoral", breaks=10)

#Diagramme pour les catégories socio-professionnelles
barplot(prof, main = "Répartition des catégories Socio-professionnelles", xlab = "Catégorie", ylab = "Fréquence", col = "lightpink")

#Diagramme en barres pour la variable PROP
barplot(prop_table, main = "Propriétaire ou Non", xlab = "Propriétaire", ylab = "Fréquence", col = "purple", names.arg = c("Non", "Oui"))

#Diagramme en barres pour la variable HOMME
barplot(homme_femme_table, main = "Répartition des genres", xlab = "Genre", ylab = "Fréquence", col = "chocolate", names.arg = c("Femmes", "Hommes"))

##############   STATISTIQUES DESCRIPTIVES BIVARIEES    ########################

########### Analyse pour les variables quantitatives AGE et DUREE ##############

#Calcul de la corrélation
correlation_age_duree<-cor(mutuelle$AGE, mutuelle$DUREE)
correlation_age_duree
 
plot(mutuelle$AGE, mutuelle$DUREE, 
     main = "Relation entre l'age et la durée du contrat",
     xlab = "Age", 
     ylab = "Durée du contrat", 
     col = "black", pch = 1.5)

########### Analyse pour les variables quanti et quali   ##############

### PROF et DUREE ###

#Moyenne et médiane de DUREE par catégorie de PROF
aggregate(mutuelle$DUREE, by = list(mutuelle$PROF), FUN = mean)
aggregate(mutuelle$DUREE, by = list(mutuelle$PROF), FUN = median)

#Graphiques par boxplot
boxplot(mutuelle$DUREE ~ mutuelle$PROF, 
        main = "Durée du contrat par catégorie socio-professionnelle",
        xlab = "Catégorie socio-professionnelle",
        ylab = "Durée du contrat",
        col = "lightblue")

### PROP et DUREE ###

#Moyenne et médiane de DUREE selon le statut de propriété
aggregate(mutuelle$DUREE, by = list(mutuelle$PROP), FUN = mean)
aggregate(mutuelle$DUREE, by = list(mutuelle$PROP), FUN = median)

#Graphique
boxplot(mutuelle$DUREE ~ mutuelle$PROP, 
        main = "Durée du contrat selon le statut de propriétaire",
        xlab = "Propriétaire (0 = Non, 1 = Oui)",
        ylab = "Durée du contrat",
        col = "lightcoral", names = c("Non", "Oui"))

### HOMME et DUREE ###

#Moyenne et médiane de DUREE selon le genre
aggregate(mutuelle$DUREE, by = list(mutuelle$HOMME), FUN = mean)
aggregate(mutuelle$DUREE, by = list(mutuelle$HOMME), FUN = median)

#Graphique
boxplot(mutuelle$DUREE ~ mutuelle$HOMME, 
        main = "Durée du contrat selon le genre",
        xlab = "Genre (0 = Femme, 1 = Homme)",
        ylab = "Durée du contrat",
        col = "purple", names = c("Femme", "Homme"))

### AGE et PROP ###

#Moyenne et médiane de AGE selon le statut de propriété
aggregate(mutuelle$AGE, by = list(mutuelle$PROP), FUN = mean)
aggregate(mutuelle$AGE, by = list(mutuelle$PROP), FUN = median)

#Graphique
boxplot(mutuelle$AGE ~ mutuelle$PROP, 
        main = "Age selon le statut de propriétaire",
        xlab = "Propriétaire (0 = Non, 1 = Oui)",
        ylab = "Age",
        col = "chocolate", names = c("Non", "Oui"))

### PROP et PROF ###

#Tableau croisé
table_prof_prop <- table(mutuelle$PROF,mutuelle$PROP)
table_prof_prop

#Test du Chi 2
chisq.test(table_prof_prop)

### HOMME et PROP ###

#Tableau croisé
table_homme_prop <- table(mutuelle$HOMME,mutuelle$PROP)
table_homme_prop

#Test du Chi 2
chisq.test(table_homme_prop)

################################################################################

#                     Analyse non-paramétriques des données

################################################################################

#Création de la colonne CENSURE
mutuelle$CENSURE<-ifelse(mutuelle$FINOBS == 2015.80, 0, 1)

library(survival)

#Création de SURVIE
surv_object <- Surv(mutuelle$DUREE, mutuelle$CENSURE)

#############################  KAPLAN-MEIER  ##################################

###############Estimation GLOBALE###################

km_all<-survfit(Surv(mutuelle$DUREE, mutuelle$CENSURE) ~ 1)

#Courbe 
plot(km_all, main = "Courbe de Kaplan-Meier sur la durée des contrats",
     xlab = "Durée des contrats (années)", ylab = "Probabilité de survie")

#Estimation du premier quartile
summary(km_all)

quantile(km_all)$quantile

############# Comparaison avec la variable HOMME ##########

km_gender<-survfit(Surv(mutuelle$DUREE, mutuelle$CENSURE) ~ mutuelle$HOMME)

plot(km_gender, col = c("pink", "blue"),
     main = "Courbes de Kaplan-Meier par genre",
     xlab = "Durée des contrats (années)", ylab = "Probabilité de survie")
legend("bottomleft",lty=c(1,1),legend=c('Femmes',"Hommes"),col=c("pink","blue"),inset=0.02)

#Test du LOG-RANG
log_rang_gender<-survdiff(Surv(mutuelle$DUREE, mutuelle$CENSURE) ~ mutuelle$HOMME)
log_rang_gender

#Estimation quantiles
summary(survfit(Surv(mutuelle$DUREE, mutuelle$CENSURE)~mutuelle$HOMME))

quantile(km_gender)$quantile

############# Comparaison avec la variable PROP ##########

km_prop<-survfit(Surv(mutuelle$DUREE, mutuelle$CENSURE) ~ mutuelle$PROP)

plot(km_prop, col = c("chocolate", "lightgreen"),
     main = "Courbes de Kaplan-Meier par statut de propriétaire",
     xlab = "Durée des contrats (années)", ylab = "Probabilité de survie")
legend("bottomleft",lty=c(1,1),legend=c('Non-propriétaires',"Propiétaires"),col=c("chocolate","lightgreen"),inset=0.02)

#Test du LOG-RANG
log_rang_prop<-survdiff(Surv(mutuelle$DUREE, mutuelle$CENSURE) ~ mutuelle$PROP)
log_rang_prop

#Estimation quantiles
summary(survfit(Surv(mutuelle$DUREE, mutuelle$CENSURE)~mutuelle$PROP))

quantile(km_prop)$quantile

############# Comparaison avec la variable PROF ##########

km_prof<-survfit(Surv(mutuelle$DUREE, mutuelle$CENSURE) ~ mutuelle$PROF)

plot(km_prof, col = rainbow(length(unique(mutuelle$PROF))),
     main = "Courbes de Kaplan-Meier par catégorie socio-professionnelle",
     xlab = "Durée des contrats (années)", ylab = "Probabilité de survie")
legend("bottomleft",lty=c(1,1),legend=levels(as.factor(mutuelle$PROF)),col=rainbow(length(unique(mutuelle$PROF))),inset=0.02 , cex=0.6)

#Test du LOG-RANG
log_rang_prof<-survdiff(Surv(mutuelle$DUREE, mutuelle$CENSURE) ~ mutuelle$PROF)
log_rang_prof

#Estimation quantiles
summary(survfit(Surv(mutuelle$DUREE, mutuelle$CENSURE)~mutuelle$PROF))

quantile(km_prof)$quantile

############# Comparaison avec la variable AGE ##########

#Ajouter une variable catégorique pour l'âge
mutuelle$AGE_CAT <- cut(mutuelle$AGE,
                        breaks = c(0, 30, 60, 100),#Catégories: 0-30, 31-60, 61 et plus
                        labels = c("Jeunes", "Adultes", "Seniors"))

#Verification
table(mutuelle$AGE_CAT)

km_age<-survfit(Surv(mutuelle$DUREE, mutuelle$CENSURE) ~ mutuelle$AGE_CAT)

plot(km_age, col = c("lightcoral", "lightblue","purple"),
     main = "Courbes de Kaplan-Meier par catégorie d'age",
     xlab = "Durée des contrats (années)", ylab = "Probabilité de survie")
legend("bottomleft",lty=c(1,1),legend=levels(mutuelle$AGE_CAT),col=c("lightcoral","lightblue","purple"),inset=0.02)

#Test du LOG-RANG
log_rang_age<-survdiff(Surv(mutuelle$DUREE, mutuelle$CENSURE) ~ mutuelle$AGE)
log_rang_age

#Estimation quantiles
summary(survfit(Surv(mutuelle$DUREE, mutuelle$CENSURE)~mutuelle$AGE_CAT))

quantile(km_age)$quantile

################################################################################

#                         Modèles de régressions

################################################################################

#Ici nous allons utiliser les modèles de COX et AFT

##########################  COX  #########################

###### Analyse du modèle avec la variable AGE
cox_age<-coxph(Surv(DUREE, CENSURE) ~ AGE, data = mutuelle)
summary(cox_age)

###### Analyse du modèle avec AGE et PROF
cox_age_prof<-coxph(Surv(DUREE, CENSURE) ~ AGE + PROF, data = mutuelle)
summary(cox_age_prof)

######## Analyse du modèle avec AGE, HOMME et PROF
cox_age_homme_prof<-coxph(Surv(DUREE, CENSURE) ~ AGE + PROF + HOMME, data = mutuelle)
summary(cox_age_homme_prof)
tableau <- as.data.frame(cox_age_homme_prof$coefficients, exp(cox_age_homme_prof$coefficients))

#Analyse du modèle avec toutes les variables
cox_model<-coxph(Surv(DUREE, CENSURE) ~ AGE + PROF + HOMME + PROP , data = mutuelle)
summary(cox_model)

##########  Analyse avec intéractions   
cox_interaction<- coxph(Surv(DUREE, CENSURE) ~ AGE * HOMME + PROF + PROP, data = mutuelle)
summary(cox_interaction)

cox_interaction_age_homme <- coxph(Surv(DUREE, CENSURE) ~ AGE * HOMME + PROF, data = mutuelle)
summary(cox_interaction_age_homme)

cox_interaction_age_prof <- coxph(Surv(DUREE, CENSURE) ~ AGE * PROF + HOMME, data = mutuelle)
summary(cox_interaction_age_prof)

cox_interaction_2<- coxph(Surv(DUREE, CENSURE) ~ AGE*PROF + HOMME + PROP, data = mutuelle)
summary(cox_interaction_2)

#Risques proportionnels

#modèle global
test_risques <- cox.zph(cox_model)
test_risques

#modèle avec AGE, PROF et HOMME
test_risques_2 <- cox.zph(cox_age_homme_prof)
test_risques_2
 
####### Estimation du risque de base cumulé et de la survie conditionnelle 
fit_mutuelle <- coxph(Surv(DUREE, CENSURE) ~ c(AGE - mean(AGE)) + factor(PROF) + 
                        c(PROP - mean(PROP)) + c(HOMME - mean(HOMME)), data = mutuelle)
Hazcum_mutuelle <- basehaz(fit_mutuelle, centered = FALSE)
par(mfrow = c(1, 2))

plot(Hazcum_mutuelle$time, Hazcum_mutuelle$hazard, type = "s", xlab = "Temps (années)", ylab = "Hasard cumulé", main = "Risque cumulé de base")
plot(Hazcum_mutuelle$time, exp(-Hazcum_mutuelle$hazard), type = "s",xlab = "Temps (années)", ylab = "Survie estimée",main = "Courbe de survie estimée")

########   Modèle de Cox avec effet dépendant du temps ####################

#Découper les données en fonction d'une durée seuil ( ici on choisi de prendre 10 ans)
mutuelle1<-survSplit(mutuelle, cut = c(10), end = "DUREE", start = "START", event = "CENSURE")

#Variable PROP

#Créer une nouvelle variable basée
mutuelle1$PROPNEW <- mutuelle1$PROP *(mutuelle1$DUREE > 10)

#Ajuster un modèle de Cox avec différentes variables
fit_prop <- coxph(Surv(START, DUREE, CENSURE) ~ AGE + factor(PROF) + PROP + factor(PROPNEW) + HOMME, data = mutuelle1)
summary(fit_prop)

#Variable HOMME

mutuelle1$HOMME_NEW <- mutuelle1$HOMME * (mutuelle1$DUREE > 10)
fit_homme <- coxph(Surv(START, DUREE, CENSURE) ~ AGE + factor(PROF) + PROP + HOMME + factor(HOMME_NEW), data = mutuelle1)
summary(fit_homme)

#Variable AGE

mutuelle1$AGE_NEW <- mutuelle1$AGE *(mutuelle1$DUREE > 10)
fit_age <- coxph(Surv(START, DUREE, CENSURE) ~ AGE + factor(PROF) + PROP + AGE_NEW + HOMME, data = mutuelle1)
summary(fit_age)

#Variable PROF

mutuelle1$PROF_NEW <- as.numeric(mutuelle1$PROF) * as.numeric(mutuelle1$DUREE > 10)
fit_prof <- coxph(Surv(START, DUREE, CENSURE) ~ AGE + factor(PROF) + PROP + PROF_NEW + HOMME, data = mutuelle1)
summary(fit_prof)

##################################  AFT ####################################

######## Modèle AFT avec distribution Weibull
fit_AFT_weibull <- survreg(Surv(DUREE, CENSURE) ~ AGE + PROF + PROP + HOMME, 
                           data = mutuelle, 
                           dist = "weibull")

summary(fit_AFT_weibull)

newdata <- data.frame(AGE = 45, PROF = factor(3, levels = levels(mutuelle$PROF)), PROP = 1, HOMME = 1)
lin_score <- predict(fit_AFT_weibull, newdata = newdata, type = "lp")
p <- 1 / fit_AFT_weibull$scale
seq_time <- seq(0, 15, by = 0.1)
base_haz_cum <- exp(-p * (lin_score)) * (seq_time)^p
surv_AFT <- exp(-base_haz_cum)

#Tracer la courbe de survie
plot(seq_time, surv_AFT, type = "l", col = "blue", lwd = 2,
     main = "Courbe de survie estimée",
     xlab = "Temps (années)", ylab = "Probabilité de survie")


####### Modèle AFT avec distribution Log-normal
fit_AFT_lognormal <- survreg(Surv(DUREE, CENSURE) ~ AGE + PROF + PROP + HOMME, 
                             data = mutuelle, 
                             dist = "lognormal")

summary(fit_AFT_lognormal)

#Calcul des paramètres pour la courbe de survie
mu <- fit_AFT_lognormal$coefficients[1] + 
  fit_AFT_lognormal$coefficients[2] * 50 + 
  fit_AFT_lognormal$coefficients[3] * 3 + 
  fit_AFT_lognormal$coefficients[4] * 1 + 
  fit_AFT_lognormal$coefficients[5] * 1

Sigma <- fit_AFT_lognormal$scale
seq_time <- seq(0, 15, by = 0.1) 

#Calcul de la probabilité de survie pour chaque temps
surv_lognormal <- plnorm(seq_time, meanlog = mu, sdlog = Sigma, lower.tail = FALSE)

#Tracer la courbe de survie 
plot(seq_time, surv_lognormal, type = "l", col = "red", lwd = 2,
     main = "Courbe de survie estimée",
     xlab = "Temps (années)", ylab = "Probabilité de survie")

#Mettre les 2 sur un meme graph
plot(seq_time, surv_AFT, type = "l", col = "blue", lwd = 2,
     main = "Courbes de survie estimées",
     xlab = "Temps (années)", ylab = "Probabilité de survie")
lines(seq_time, surv_lognormal, col = "red", lwd = 2)

# Ajout de la légende
legend("bottomleft", legend = c("Weibull", "Log-normal"), 
       col = c("blue", "red"), lwd = 2)


##########  Modèle AFT avec approche GEE
require(aftgee)
fit_AFTgee <- aftgee(Surv(DUREE, CENSURE) ~ AGE + PROF + PROP + HOMME, 
                     data = mutuelle)

summary(fit_AFTgee)

#Affichage des coefficients estimés
cat("Coefficients estimés :", fit_AFTgee$coefficients, "\n")

#Tester la significativité globale
cat("Statistiques globales :", fit_AFTgee$test, "\n")

#Tracer les résidus
residuals_gee <- residuals(fit_AFTgee)
plot(residuals_gee, main = "Résidus AFT GEE", xlab = "Index", ylab = "Résidus", pch = 20)
abline(h = 0, col = "red", lwd = 2)