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

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

# Load the original TREE 1 database with all waves and weights
# D <- load( ... )
# Example for a SPSS dataset:
# D <- read.spss("Data_projet.sav", to.data.frame = T)

n <- nrow(D)

# Création des variables de consommation d'alcool en 3 catégories

D$t1alcool <- NA
D$t1alcool[D$t1drug1=="jamais"] <- 1
D$t1alcool[D$t1drug1=="1-3 fois par mois"] <- 2
D$t1alcool[D$t1drug1=="1-2 fois par semaine" | D$t1drug1=="3-5 fois par semaine" | D$t1drug1=="chaque jour"] <- 3

D$t2alcool <- NA
D$t2alcool[D$t2drug1=="jamais"] <- 1
D$t2alcool[D$t2drug1=="1-3 fois par mois"] <- 2
D$t2alcool[D$t2drug1=="1-2 fois par semaine" | D$t2drug1=="3-5 fois par semaine" | D$t2drug1=="chaque jour"] <- 3

D$t3alcool <- NA
D$t3alcool[D$t3drug1=="jamais"] <- 1
D$t3alcool[D$t3drug1=="1-3 fois par mois"] <- 2
D$t3alcool[D$t3drug1=="1-2 fois par semaine" | D$t3drug1=="3-5 fois par semaine" | D$t3drug1=="chaque jour"] <- 3

D$t4alcool <- NA
D$t4alcool[D$t4drug1=="jamais"] <- 1
D$t4alcool[D$t4drug1=="1-3 fois par mois"] <- 2
D$t4alcool[D$t4drug1=="1-2 fois par semaine" | D$t4drug1=="3-5 fois par semaine" | D$t4drug1=="chaque jour"] <- 3

D$t5alcool <- NA
D$t5alcool[D$t5drug1=="jamais"] <- 1
D$t5alcool[D$t5drug1=="1-3 fois par mois"] <- 2
D$t5alcool[D$t5drug1=="1-2 fois par semaine" | D$t5drug1=="3-5 fois par semaine" | D$t5drug1=="chaque jour"] <- 3

D$t6alcool <- NA
D$t6alcool[D$t6drug1=="jamais"] <- 1
D$t6alcool[D$t6drug1=="1-3 fois par mois"] <- 2
D$t6alcool[D$t6drug1=="1-2 fois par semaine" | D$t6drug1=="3-5 fois par semaine" | D$t6drug1=="chaque jour"] <- 3

D$t7alcool <- NA
D$t7alcool[D$t7drug1=="jamais"] <- 1
D$t7alcool[D$t7drug1=="1-3 fois par mois"] <- 2
D$t7alcool[D$t7drug1=="1-2 fois par semaine" | D$t7drug1=="3-5 fois par semaine" | D$t7drug1=="chaque jour"] <- 3

D$t8alcool <- NA
D$t8alcool[D$t8drug1=="jamais"] <- 1
D$t8alcool[D$t8drug1=="1-3 fois par mois"] <- 2
D$t8alcool[D$t8drug1=="1-2 fois par semaine" | D$t8drug1=="3-5 fois par semaine" | D$t8drug1=="chaque jour"] <- 3

D$t9alcool <- NA
D$t9alcool[D$t9drug1=="jamais"] <- 1
D$t9alcool[D$t9drug1=="1-3 fois par mois"] <- 2
D$t9alcool[D$t9drug1=="1-2 fois par semaine" | D$t9drug1=="3-5 fois par semaine" | D$t9drug1=="chaque jour"] <- 3

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")

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

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

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

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

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

D$t9alcool <- as.factor(D$t9alcool)
levels(D$t9alcool) <- 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)

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

###############################################################
# On garde toutes les séquences avec au moins 6 observations
# et on impute les données manquantes
D2 <- D[D$Nb.alc>=6,]

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

# Pondérations
D2$pond <- NA
selwt <- grep("wt",names(D2))

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

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

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

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

# Imputation des données manquantes
DM <- D2[,sel2]
DMi <- seqimpute(DM, np=1, nf=1, nfi=1, npt=1, pastDistrib=FALSE, futureDistrib=FALSE, mi=1, regr="mlogit",ParExec=F)


# 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 <- seqdef(DMi, weights = D2$pond, states = alcool.states,
                     labels = alcool.labels,alphabet=alcool.alphabet,right=NA)

# plots
seqfplot(sso.alcool, main = "Frequence des sequences (100 plus frequentes)", idxs=1:100, border = NA)
seqdplot(sso.alcool, main = "Distribution des états", border = NA)

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

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

# Classification du resultat
# Attention: necessite la librairie cluster
Classif <- agnes(OM, diss = TRUE, method = "ward")
plot(Classif)

# Le dendrogram semble suggérer une solution en 2 ou 3 groupes

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

# Representer cette solution
seqfplot(sso.alcool,group=Classif.3fac, idxs=1:100, main="Frequence des sequences", border=NA)
seqdplot(sso.alcool,group=Classif.3fac, main="Distribution des etats", border=NA)

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

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

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

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


#################################################################################
# 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 <- 6
D2mi <- D[D$Nb.alc>=nbo,]

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

# Pondérations
D2mi$pond <- NA
selwt <- grep("wt",names(D2mi))

# 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:11], 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", keep.diss=F, keep.data = F)

# 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==D2mi$id[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="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 <- 6
D3 <- D[D$Nb.alc>=nbo & is.na(D$Pond)==F,]

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

# Ajout de l'état 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")
levels(D3$t4alcool) <- c("Abstinent","Light","Heavy","Missing")
levels(D3$t5alcool) <- c("Abstinent","Light","Heavy","Missing")
levels(D3$t6alcool) <- c("Abstinent","Light","Heavy","Missing")
levels(D3$t7alcool) <- c("Abstinent","Light","Heavy","Missing")
levels(D3$t8alcool) <- c("Abstinent","Light","Heavy","Missing")
levels(D3$t9alcool) <- c("Abstinent","Light","Heavy","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"
D3$t4alcool[is.na(D3$t4alcool)==T] <- "Missing"
D3$t5alcool[is.na(D3$t5alcool)==T] <- "Missing"
D3$t6alcool[is.na(D3$t6alcool)==T] <- "Missing"
D3$t7alcool[is.na(D3$t7alcool)==T] <- "Missing"
D3$t8alcool[is.na(D3$t8alcool)==T] <- "Missing"
D3$t9alcool[is.na(D3$t9alcool)==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
CClassif <- D3$Classif.4

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)