## ====================
## Etapes preliminaires
## ====================

# Installation de deux bibliotheques de fonctions necessaires au programme
# (seulement si elles ne sont pas deja installees)
install.packages(c('TraMineR','cluster'))

# Chargement en memoire de ces bibliotheques
library(TraMineR)
library(cluster)

# Mise en memoire des parametres graphiques par defaut
par.def <- par()

# Importation des trajectoires
donnees <- read.csv("http://nicolas.robette.free.fr/Docs/trajpro.csv",header=T)

# Recodage la variable "generation"
donnees$generation <- factor(donnees$generation,labels=c('1930-38','1939-45','1946-50'))

# Regarder la structure des donnees
str(donnees)

# Definition des labels des etats
labels <- c("agric","acce","cadr","pint","empl","ouvr","etud","inact","smil")

# Mise en forme des trajectoires sous forme d'un "objet" sequence
seq <- seqdef(donnees[,1:37], states=labels)



## =============================
## Optimal matching et typologie
## =============================

# Definition des couts de substitution pour l'optimal matching
couts <- seqsubm(seq,method="CONSTANT", cval=2)

# Calcul des dissimilarites avec l'optimal matching
seq.om <- seqdist(seq, method="OM", indel=1.1, sm=couts)

# Classification ascendante hierarchique
seq.agnes <- agnes(as.dist(seq.om), method="ward", keep.diss=FALSE)

# Dendrogramme de la classification
par(par.def)
plot(as.dendrogram(seq.agnes), leaflab="none")

# Inertie des partitions selon le nombre de classes
par(par.def)
plot(sort(seq.agnes$height, decreasing=TRUE)[1:20], type='s', xlab="nb de classes", ylab="inertie")

# Choix d'une partition en 5 classes
nbcl <- 5
seq.part <- cutree(seq.agnes, nbcl)
seq.part <- factor(seq.part,labels=paste('classe',1:nbcl,sep='.'))



## ==========================================
## Representations graphiques de la typologie
## ==========================================

# Chronogrammes de la typologie (state distribution plots)
seqdplot(seq, group=seq.part, xtlab=14:50, border=NA, withlegend=T)

# Tapis de la typologie (index plots)
seqiplot(seq, group=seq.part, xtlab=14:50, tlim=0, space=0, border=NA, withlegend=T, yaxis=FALSE)

# Idem, tries par multidimensional scaling
ordre <- cmdscale(as.dist(seq.om),k=1)
seqiplot(seq, group=seq.part, sortv=ordre, xtlab=14:50, tlim=0, space=0, border=NA, withlegend=T, yaxis=FALSE)

# Trajectoires les plus frequentes
seqfplot(seq, group=seq.part, withlegend=T)

# Entropie transversale (homogeneite) pour chaque annee
seqHtplot(seq, group=seq.part, xtlab=14:50, withlegend=T)

# Etat modal pour chaque annee
seqmsplot(seq, group=seq.part, xtlab=14:50, withlegend=T)

# Duree moyenne dans chaque etat
seqmtplot(seq, group=seq.part, withlegend=T)



## =========================================
## Description des sequences de la typologie
## =========================================

# Distribution de la partition (effectif et pourcentage)
distri.eff <- table(seq.part)
distri.eff
round(distri.eff/sum(distri.eff)*100,1)

# Distance moyenne au centre de la classe (autre mesure d'homogeneite)
round(aggregate(disscenter(as.dist(seq.om), group=seq.part),list(seq.part),mean)[,-1],1)



## =====================================================
## Poursuivre la description des classes de la typologie
## =====================================================

# Croiser avec autres variables (tris croises et khi-deux)
round(100*prop.table(table(seq.part,donnees$generation),2),1)
chisq.test(table(seq.part,donnees$generation[1:1000]))