setwd("C:/Users/harari/Desktop/Dropbox/enseignement/M1 EPOG MIF/tobit")# à adapter au chemin de votre dossier
library(tidyverse)
#install.packages("AER")
library(AER)
#install.packages("ggrepel")
library(ggrepel)

db_30_11 <- read.csv("db_30_11.txt")
summary(db_30_11$Classement)

#on indique que le classement est supérieur à 1000 quand une fac n'est pas classée
db_30_11$Classement2<-db_30_11$Classement %>% replace_na(1000)
summary(db_30_11$Classement2)

#on représente avec les noms
ggplot(db_30_11,aes(x=effectif,y=Classement2))+
  geom_point(col='blue',pch='o')+
  geom_text_repel(data=db_30_11[db_30_11$Classement2<1000,],
                  aes(label=Libelle),size=3)

#modèle avec une seule variable explicative
out1<-tobit(data=db_30_11,
      Classement2~effectif,left=-Inf, right = 1000)
summary(out1)

#modèle avec plusieurs variables explicatives
out2<-tobit(data=db_30_11,
Classement2~effectif+prof_sup+bac_S+ouvrier+employe+masters,
left=-Inf ,right = 1000)
summary(out2)


# base de données avec les 2 modèles 
db_30_11<-db_30_11 %>%
  mutate(fit2=out2[["linear.predictors"]]) %>% 
  mutate(fit2_tronc=fit2*(fit2<1000)+1000*(fit2>=1000)) %>% 
  mutate(fit1_tronc=fit1*(fit1<1000)+1000*(fit1>=1000))

# graphique avec les deux modèles
ggplot(db_30_11,aes(x=effectif,y=Classement2, size=3))+
  geom_point(col='blue',pch='o')+
  geom_text_repel(data=db_30_11[db_30_11$Classement2<300,],
                  aes(label=Libelle))+
  geom_line(aes(x=effectif,y=fit1_tronc),col='green',pch='x',size=1)+
  geom_point(aes(x=effectif,y=fit2_tronc),col='red',pch='x')