# Sequence_analysis--GenFRee-syntaxe-pour-archivage.r

library(foreign)
library(TraMineR)
library(cluster)
library(seqimpute)

# Load the provided GenFRee dataset into a dataframe named D
# D <- load( ... )
# Example for a Stata dataset:
# D <- read.dta("GenFRee_cohortes_agregees_modele-2_T0T1agg_28062019-Statav12.dta")

n <- nrow(D)

# Creation des variables de consommation d'alcool en 3 categories
# Rem: on utilise 3 variables Q48, Q48_1 et Q48_2, mais la Q48 n'est pas
#      disponible a T4

D$t0alcool <- NA
D$t0alcool[D$T0_Q48=="Non"] <- 1 
D$t0alcool[D$T0_Q48=="Oui" & D$T0_Q48_1=="Jamais" & D$T0_Q48_2=="Jamais"] <- 2 
D$t0alcool[D$T0_Q48=="Oui" & (D$T0_Q48_1=="1 à 2 fois" | D$T0_Q48_1=="3 à 9 fois" | D$T0_Q48_1=="10 fois ou plus") & (D$T0_Q48_2=="1 à 2 fois" | D$T0_Q48_2=="3 à 9 fois" | D$T0_Q48_2=="10 fois ou plus")] <- 3 
table(D$t0alcool)

D$t1alcool <- NA
D$t1alcool[D$T1_Q48=="Non"] <- 1 
D$t1alcool[D$T1_Q48=="Oui" & D$T1_Q48_1=="Jamais" & D$T1_Q48_2=="Jamais"] <- 2 
D$t1alcool[D$T1_Q48=="Oui" & (D$T1_Q48_1=="1 à 2 fois" | D$T1_Q48_1=="3 à 9 fois" | D$T1_Q48_1=="10 fois ou plus") & (D$T1_Q48_2=="1 à 2 fois" | D$T1_Q48_2=="3 à 9 fois" | D$T1_Q48_2=="10 fois ou plus")] <- 3 
table(D$t1alcool)

D$t2alcool <- NA
D$t2alcool[D$T2_Q48=="Non"] <- 1 
D$t2alcool[D$T2_Q48=="Oui" & D$T2_Q48_1=="Jamais" & D$T2_Q48_2=="Jamais"] <- 2 
D$t2alcool[D$T2_Q48=="Oui" & (D$T2_Q48_1=="1 à 2 fois" | D$T2_Q48_1=="3 à 9 fois" | D$T2_Q48_1=="10 fois ou plus") & (D$T2_Q48_2=="1 à 2 fois" | D$T2_Q48_2=="3 à 9 fois" | D$T2_Q48_2=="10 fois ou plus")] <- 3 
table(D$t2alcool)

D$t3alcool <- NA
D$t3alcool[D$T3_Q48=="Non"] <- 1 
D$t3alcool[D$T3_Q48=="Oui" & D$T3_Q48_1=="Jamais" & D$T3_Q48_2=="Jamais"] <- 2 
D$t3alcool[D$T3_Q48=="Oui" & (D$T3_Q48_1=="1 à 2 fois" | D$T3_Q48_1=="3 à 9 fois" | D$T3_Q48_1=="10 fois ou plus") & (D$T3_Q48_2=="1 à 2 fois" | D$T3_Q48_2=="3 à 9 fois" | D$T3_Q48_2=="10 fois ou plus")] <- 3 
table(D$t3alcool)

D$t0alcool <- as.factor(D$t0alcool)
levels(D$t0alcool) <- c("Abstinent","Light","Heavy")

D$t1alcool <- as.factor(D$t1alcool)
levels(D$t1alcool) <- c("Abstinent","Light","Heavy")

D$t2alcool <- as.factor(D$t2alcool)
levels(D$t2alcool) <- c("Abstinent","Light","Heavy")

D$t3alcool <- as.factor(D$t3alcool)
levels(D$t3alcool) <- c("Abstinent","Light","Heavy")


# Nombre de données alcool pour chaque répondant
sel <- grep("alcool",names(D))
D$Nb.alc <- NA
for (i in 1:n){
  D$Nb.alc[i] <- sum(is.na(D[i,sel])==F)
}
table(D$Nb.alc)

# Rem: aucune séquence complète. Cela est du au design de l'enquete et
#      a l'aggregation realisee des deux cohortes


####################################################
# Création d'un sous-ensemble de données utilisables


#################################################################################
# F: Imputation multiple (nb.mi), classification et attribution de chaque personne
#    au groupe le plus probable en fonction du nombre de classifications dans
#    chaque groupe

nb.mi <- 10

# On garde toutes les séquences avec au minimum nbo valeurs observées
nbo <- 2
var.selection <- c(grep("ID_AG",names(D)), grep("alcool",names(D)), grep("pond",names(D)))
D2mi <- D[D$Nb.alc>=nbo,var.selection]

sel2 <- grep("alcool",names(D2mi))

# Pondérations
D2mi$pond <- NA
selwt <- c(9,7,8,16)

# a) On recherche la dernière observation alcool de chaque personne
#    et on prend la pondération correspondante
for (i in 1:nrow(D2mi)){
  qui <- which(is.na(D2mi[i,sel2])==F)
  qui <- qui[length(qui)]
  D2mi$pond[i] <- D2mi[i,selwt[qui]]
}

# Suppression des observations avec NA comme pondération
D2mi <- D2mi[is.na(D2mi$pond)==F,]

# b) Normalisation des pondérations
D2mi$pond <- D2mi$pond/sum(D2mi$pond)*nrow(D2mi)

# On fixe une seed pour pouvoir être reproductible
set.seed(1234)

# Imputation des données manquantes
DMmi <- D2mi[,sel2]
DMmi.i <- seqimpute(DMmi, np=1, nf=1, nfi=1, npt=1, pastDistrib=FALSE, futureDistrib=FALSE, mi=nb.mi, regr="mlogit",ParExec=F)

DMmi.i$pond <- kronecker(rep(1,nb.mi),D2mi$pond)

# Création d'une typologie
alcool.alphabet <- c("Abstinent","Light","Heavy")
alcool.states <- c("Abstinent","Light","Heavy")
alcool.labels <- c("Abstinent","Light","Heavy")

sso.alcool.mi <- seqdef(DMmi.i[,3:6], weights = DMmi.i$pond, states = alcool.states,
                     labels = alcool.labels,alphabet=alcool.alphabet,right=NA)

# Optimal matching
# Matrice des couts de transition
submat.mi <- seqsubm(sso.alcool.mi, method = "TRATE")

# Optimal matching
OMmi <- seqdist(sso.alcool.mi, method = "OM", indel = 1, sm=submat.mi)

# Classification du resultat
# Attention: necessite la librairie cluster
# Classif.mi <- agnes(OMmi, diss = TRUE, method = "ward")
Classif.mi <- agnes(OMmi, diss = TRUE, method = "ward", keep.diss=F, keep.data = F)
#Classif.mi <- hclust(OMmi, method="ward.D2")
#plot(Classif.mi)


# Extraction de la classification en 2 groupes
Classif.mi.2 <- cutree(Classif.mi, k = 2)
Classif.mi.2fac <- factor(Classif.mi.2, labels = paste("Groupe", 1:2))

# Representer cette solution
seqIplot(sso.alcool.mi,group=Classif.mi.2fac, main="Frequence des sequences", border=NA, sortv = "from.start")
seqdplot(sso.alcool.mi,group=Classif.mi.2fac, main="Distribution des etats", border=NA)

# Extraction de la classification en 3 groupes
Classif.mi.3 <- cutree(Classif.mi, k = 3)
Classif.mi.3fac <- factor(Classif.mi.3, labels = paste("Groupe", 1:3))

# Representer cette solution
seqIplot(sso.alcool.mi,group=Classif.mi.3fac, main="Frequence des sequences", border=NA, sortv = "from.start")
seqdplot(sso.alcool.mi,group=Classif.mi.3fac, main="Distribution des etats", border=NA)

# Extraction de la classification en 4 groupes
Classif.mi.4 <- cutree(Classif.mi, k = 4)
Classif.mi.4fac <- factor(Classif.mi.4, labels = paste("Groupe", 1:4))

# Representer cette solution
seqIplot(sso.alcool.mi,group=Classif.mi.4fac, main="Frequence des sequences", border=NA, sortv = "from.start")
seqdplot(sso.alcool.mi,group=Classif.mi.4fac, main="Distribution des etats", border=NA)

# Extraction de la classification en 5 groupes
Classif.mi.5 <- cutree(Classif.mi, k = 5)
Classif.mi.5fac <- factor(Classif.mi.5, labels = paste("Groupe", 1:5))

# Representer cette solution
seqIplot(sso.alcool.mi,group=Classif.mi.5fac, main="Frequence des sequences", border=NA, sortv = "from.start")
seqdplot(sso.alcool.mi,group=Classif.mi.5fac, main="Distribution des états", border=NA)

# Extraction de la classification en 6 groupes
Classif.mi.6 <- cutree(Classif.mi, k = 6)
Classif.mi.6fac <- factor(Classif.mi.6, labels = paste("Groupe", 1:6))

# Representer cette solution
seqIplot(sso.alcool.mi,group=Classif.mi.6fac, main="Frequence des sequences", border=NA, sortv = "from.start")
seqdplot(sso.alcool.mi,group=Classif.mi.6fac, main="Distribution des etats", border=NA)

# Extraction de la classification en 7 groupes
Classif.mi.7 <- cutree(Classif.mi, k = 7)
Classif.mi.7fac <- factor(Classif.mi.7, labels = paste("Groupe", 1:7))

# Representer cette solution
seqIplot(sso.alcool.mi,group=Classif.mi.7fac, main="Frequence des sequences", border=NA, sortv = "from.start")
seqdplot(sso.alcool.mi,group=Classif.mi.7fac, main="Distribution des etats", border=NA)

# Extraction de la classification en 8 groupes
Classif.mi.8 <- cutree(Classif.mi, k = 8)
Classif.mi.8fac <- factor(Classif.mi.8, labels = paste("Groupe", 1:8))

# Representer cette solution
seqIplot(sso.alcool.mi,group=Classif.mi.8fac, main="Frequence des sequences", border=NA, sortv = "from.start")
seqdplot(sso.alcool.mi,group=Classif.mi.8fac, main="Distribution des etats", border=NA)

# Extraction de la classification en 9 groupes
Classif.mi.9 <- cutree(Classif.mi, k = 9)
Classif.mi.9fac <- factor(Classif.mi.9, labels = paste("Groupe", 1:9))

# Representer cette solution
seqIplot(sso.alcool.mi,group=Classif.mi.9fac, main="Frequence des sequences", border=NA, sortv = "from.start")
seqdplot(sso.alcool.mi,group=Classif.mi.9fac, main="Distribution des etats", border=NA)

# Extraction de la classification en 10 groupes
Classif.mi.10 <- cutree(Classif.mi, k = 10)
Classif.mi.10fac <- factor(Classif.mi.10, labels = paste("Groupe", 1:10))

# Representer cette solution
seqIplot(sso.alcool.mi,group=Classif.mi.10fac, main="Frequence des sequences", border=NA, sortv = "from.start")
seqdplot(sso.alcool.mi,group=Classif.mi.10fac, main="Distribution des etats", border=NA)


# Attribution de chaque personne à son groupe final (de 2 à 10 groupes)
FClassif <- matrix(NA,nrow=nrow(D2mi),ncol=10)
FClassif[,1] <- 1

Egal <- rep(0,10)
#
# On attribue chaque réplication d'une même personne à son groupe

for (k in 2:10){
  nb.groups <- k
  if (k==2){
    Classif.k <- Classif.mi.2    
  }
  if (k==3){
    Classif.k <- Classif.mi.3   
  }
  if (k==4){
    Classif.k <- Classif.mi.4    
  }
  if (k==5){
    Classif.k <- Classif.mi.5    
  }
  if (k==6){
    Classif.k <- Classif.mi.6    
  }
  if (k==7){
    Classif.k <- Classif.mi.7    
  }
  if (k==8){
    Classif.k <- Classif.mi.8    
  }
  if (k==9){
    Classif.k <- Classif.mi.9    
  }
  if (k==10){
    Classif.k <- Classif.mi.10    
  }
  #
  GF <- matrix(0,nrow=nrow(D2mi),ncol=nb.groups)
  tt <- 0
  for (i in 1:nb.mi){
    for (j in 1:nrow(D2mi)){
      tt <- tt+1
      GF[j,Classif.k[tt]] <- GF[j,Classif.k[tt]]+1
    }
  }
  
  # On détermine le groupe le plus probable
  for (i in 1:nrow(GF)){
    sel <- which(GF[i,]==max(GF[i,]))
    if (length(sel)>1){
      rnd <- runif(length(sel))
      sel2 <- which(rnd==max(rnd))
      sel2 <- sel2[1]
      sel <- sel[sel2]
      Egal[k] <- Egal[k]+1
    }
    FClassif[i,k] <- sel
  }
}

print(Egal)
table(FClassif)

# Mise de la classification et des pondérations finales dans D
D$Classif.2 <- NA
D$Classif.3 <- NA
D$Classif.4 <- NA
D$Classif.5 <- NA
D$Classif.6 <- NA
D$Classif.7 <- NA
D$Classif.8 <- NA
D$Classif.9 <- NA
D$Classif.10 <- NA
D$Pond <- NA

for (i in 1:nrow(D2mi)){
  sel <- which(D$ID_AG==D2mi$ID_AG[i])
  D$Classif.2[sel] <- FClassif[i,2]
  D$Classif.3[sel] <- FClassif[i,3]
  D$Classif.4[sel] <- FClassif[i,4]
  D$Classif.5[sel] <- FClassif[i,5]
  D$Classif.6[sel] <- FClassif[i,6]
  D$Classif.7[sel] <- FClassif[i,7]
  D$Classif.8[sel] <- FClassif[i,8]
  D$Classif.9[sel] <- FClassif[i,9]
  D$Classif.10[sel] <- FClassif[i,10]
  D$Pond[sel] <- D2mi$pond[i]
}

save(D, file="GenFRee_D10_temp.rda")

# Représentation des séquences à partir de D, sans imputations
# On garde toutes les séquences avec au minimum nbo valeurs observées
nbo <- 2
D3 <- D[D$Nb.alc>=nbo & is.na(D$Pond)==F,]

sel3 <- grep("alcool",names(D3))

# Ajout de l'état Missing
levels(D3$t0alcool) <- c("Abstinent","Light","Heavy","Missing")
levels(D3$t1alcool) <- c("Abstinent","Light","Heavy","Missing")
levels(D3$t2alcool) <- c("Abstinent","Light","Heavy","Missing")
levels(D3$t3alcool) <- c("Abstinent","Light","Heavy","Missing")

D3$t0alcool[is.na(D3$t0alcool)==T] <- "Missing"
D3$t1alcool[is.na(D3$t1alcool)==T] <- "Missing"
D3$t2alcool[is.na(D3$t2alcool)==T] <- "Missing"
D3$t3alcool[is.na(D3$t3alcool)==T] <- "Missing"

# Création d'une typologie
alcool.alphabet.f <- c("Abstinent","Light","Heavy","Missing")
alcool.states.f <- c("Abstinent","Light","Heavy","Missing")
alcool.labels.f <- c("Abstinent","Light","Heavy","Missing")

sso.alcool.f <- seqdef(D3[,sel3], weights = D3$Pond, states = alcool.states.f,
                     labels = alcool.labels.f,alphabet=alcool.alphabet.f,right=NA)

# Representer la solution (changer le .x pour changer le nb de groupes)
CClassif <- D3$Classif.7

seqIplot(sso.alcool.f,group=CClassif, main="Frequence des sequences", border=NA, sortv = "from.start")
seqdplot(sso.alcool.f,group=CClassif, main="Distribution des etats", border=NA)

# Sauvegarde de tout
save.image(file = "Classification_GenFRee_all.RData")