Inspect the data structure

## Number of observations and variables
dim(fertility) 

## [1] 8984   22

There are 8,984 respondents and 22 variables in the NLSY97 fertility history data.

## List all the variable names
ls(fertility)

##  [1] "C18"       "C19"       "C20"       "C21"       "C22"       "C23"      
##  [7] "C24"       "C25"       "C26"       "C27"       "C28"       "C29"      
## [13] "C30"       "C31"       "C32"       "C33"       "C34"       "C35"      
## [19] "education" "female"    "pubid"     "race"

The “pubid” is the unique identifier for the respondents. The fertility sequence data is stored in wide format. The variables starting with “C” are the sequence variables ranging from age 18 to 35 with equal interval. There are also 3 sociodemographic variables (i.e., gender, race, and education) that we will use later to conduct analysis to examine sequence-covariate associations.

## Check the unique values of the state variable, e.g., C35
unique(fertility$C35) 

## <labelled<double>[8]>: number of child at age 35
## [1] 8 1 3 2 4 5 6 7
## 
## Labels:
##  value      label
##      1    0 child
##      2    1 child
##      3    2 child
##      4    3 child
##      5    4 child
##      6    5 child
##      7 6 or above
##      8    missing

The state space includes 8 unique values, ranging from having “0 child” to “6+ child”, and the missing value is coded as a unique state.

Extract the 18 sequence variables

## Select the 18 sequence variables: C18-C35
fertility.vars <- fertility %>%
  select(starts_with("C"))

Define labels for the state space

## Create state labels 
shortlabel <- c("C=0", "C=1", "C=2", "C=3", "C=4", "C=5", "C=6/+", "NA")
# The shortlabel will show in plotting as default. 
longlabel  <- c("0 child", "1 child", "2 child", "3 child", 
                "4 child", "5 child", "6 or above child", "attrited") 

Create state sequence object using seqdef() in package “TraMineR”

fertility.seq <- seqdef(fertility.vars,    #the previously extracted state variables
                        states=shortlabel, #the predefined short labels for states
                        labels=longlabel,  #the long labels for states
                        id=fertility$pubid, 
                        xtstep=1)          #the default is 1, the tick-mark is displayed at each position 
## Type ?seqdef() in R console to learn more about the available options

Set the color scheme for the state space

library(RColorBrewer)
display.brewer.all()

nb.cols <- 8
color1<- colorRampPalette(brewer.pal(8, "Spectral"))(nb.cols)
color1<- colorRampPalette(rev(brewer.pal(8, "Spectral")))(nb.cols) #reverse the color scheme
cpal(fertility.seq) <-color1

## the state variables are ordinal except for the attrited, thus use "Spectral" color scheme

## Plot the legend for the states in the defined fertility sequence object 
seqlegend(fertility.seq, cpal=color1, ncol=1, cex=1.5)

Modal state plot: the most frequent state at each position

par(mar=c(4,4,4,4)) #set the margins of plotting area 
seqmsplot(fertility.seq, ylab="Relative Frequency", xlab= "Age", xtlab=18:35, with.legend=F, cex.axis=0.8)

## OR use ggseqmsplot() in package "ggseqplot"
ggseqmsplot(fertility.seq, no.n = TRUE) 

# it allows to show the standard error bar

The modal state remains to be childless from age 18 to 35, reflecting the overall delaying of childbearing in NLSY97 cohort.

State Distribution Plot

State distribution plot (or known as chronogram) is a series of stacked bar charts at each position (i.e., age) of the sequence. It shows the proportions of states at each position (i.e., age for our example).

par(mar=c(4,4,4,4)) #set the margins of plotting area 
seqdplot(fertility.seq, xtlab=18:35, with.legend=F, ylab= "Relative Frequency", xlab= "Age", cex.axis=0.8)

With the increase of age, a higher proportion of respondents start to distribute in the states of having 1 or more children.

## The seqplot() can all be stratified by group of interest. 
# For illustration, we stratify the sequence distribution plot by education level. 
seqdplot(fertility.seq, xtlab=18:35, with.legend=F, 
         xlab= "Age",  group=fertility$education, cex.axis=0.8)

The plotting functions in the main workhorse of sequence analysis “TraMineR” is not too flexible to modify different components of the figure. The new package “ggseqplot” can be used. See an example below:

ggseqdplot(fertility.seq, group=fertility$education, facet_ncol=2, facet_nrow=2) 

From the state distribution plot by education level, we can see that college-educated respondents start to have children at later age. By age 30, over one third of them remain child-free.

Transversal Shannon’s Entropy

Transversal entropy is a measure of diversity of states observed at each position. A value of 0 means that all the cases in are in the same state and its value is maximal when the same proportion of cases are in each state.

## seqstatd() function returns the state relative frequencies, the number of valid states and the entropy of the state distribution at each position (i.e., age) in the sequence. 
sd <- seqstatd(fertility.seq)
plot(sd, type="Ht", main="Age-specific Shannon's Entropy", 
     xlab= "Age", xtlab=18:35)

The age-specific entropy increases over time, indicating that the fertility statuses become more diverse as individuals age.

Sequence Frequency Plot

The sequences are presented as stacks of successive states, ordered by their relative frequency in the dataset.

par(mar=c(4,4,4,4)) #set the margins of plotting area 
seqfplot(fertility.seq, with.legend=F, xlab= "Age", xtlab=18:35, cex.axis=0.8)

Child-free sequence is the most frequent, followed by the sequences with no fertility status observed.

Sequence Index Plot

Sequence index plot shows individual sequences as stacked bars. All the dimensions of temporality, i.e., timing, duration, and order, can be seen in such plot.

par(mar=c(4,4,4,4)) #set the margins of plotting area 
## seqiplot() shows the first 10 sequences
seqiplot(fertility.seq, with.legend=F, xlab= "Age", xtlab=18:35)  

## seqIplot() shows all the sequences
seqIplot(fertility.seq, sortv="from.start", with.legend=F, xlab= "Age", xtlab=18:35, cex.axis=0.8) 

Basic Characteristics within Sequences

## number of transitions 
ntrans <- seqtransn(fertility.seq)
## number of subsequences
subsn  <- seqsubsn(fertility.seq)

## longitudinal entropy 
entropy_long <- seqient(fertility.seq) 
## check the distribution of longitudinal entropy in the sample
hist(entropy_long ) 

By inspecting the distribution of longitudinal entropy, a substantial portion of sequences have value 0. It means that these individuals remain in the same state across the 18 years of observation, e.g., childless or missing all the time.

## turbulence 
turbulence <- seqST(fertility.seq)
## complexity 
complexity <- seqici(fertility.seq)

## Combine all the within-sequence characteristics just calculated
indices <-tibble(ntrans, subsn, turbulence, complexity, entropy_long) 

## correlation between the longitudinal characteristics
cor(indices)  

##                 ntrans     subsn turbulence complexity entropy_long
## ntrans       1.0000000 0.7980644  0.8156056  0.9873341    0.9450613
## subsn        0.7980644 1.0000000  0.5678313  0.7436734    0.6657299
## turbulence   0.8156056 0.5678313  1.0000000  0.8823426    0.9357408
## complexity   0.9873341 0.7436734  0.8823426  1.0000000    0.9847759
## entropy_long 0.9450613 0.6657299  0.9357408  0.9847759    1.0000000

In addition to these within-sequence characteristics, there are also sequence indices reflecting meaningful aspects of sequences, e.g., precarity index. These indices can be used as predictors or outcome variables in further analysis (Ritschard, 2021).

Hierarchical clustering

## Perform hierarchical clustering using "ward.D" method, the input is the OM-based dissimilarity matrix (specified in as.dist)
wardCluster <- hclust(as.dist(OM), method="ward.D")

## Dendrogram represents the agglomerative (bottom up) clustering procedure (see slide 39)
as.dendrogram(wardCluster) %>% plot(leaflab="none")

“One starts out from the observations, each of them being considered as a group. On each iteration, the two closest groups (or initially observations) are grouped, until all the observations form a single group. The agglomeration schedule, i.e. the succession of groupings performed, represents the clustering procedure in the form of a tree which is called a dendrogram” (Studer, 2013:6).

## seq_heatmap(): index plot of sequences ordered according to a dendrogram
seq_heatmap(fertility.seq, wardCluster)

Determine the optimal number of clusters

Various quality measures can be computed to determine the optimal number of clusters. Good clustering needs to maximize the homogeneity within cluster and heterogeneity between group (i.e., whether clusters are distinct from each other). For details, please refer to the “WeightedCluster” manual (Studer, 2013, page 33-34).

## as.clustrange() produces the all the quality indices 
wardRange<-as.clustrange(wardCluster, diss=OM)

## print the quality indices for the two best solutions
summary(wardRange, max.rank=2)

##      1. N groups     1.  stat 2. N groups     2.  stat
## PBC            9 6.669443e-01           8 6.662702e-01
## HG            20 9.443307e-01          19 9.403195e-01
## HGSD          20 9.439064e-01          19 9.399031e-01
## ASW            7 4.556711e-01           8 4.524956e-01
## ASWw           7 4.561074e-01           8 4.530312e-01
## CH             2 3.209675e+03           3 2.726594e+03
## R2            20 7.550712e-01          19 7.483761e-01
## CHsq           7 5.445387e+03           6 5.427821e+03
## R2sq          20 8.988933e-01          19 8.954413e-01
## HC            20 4.379441e-02           9 4.428889e-02

# 7-cluster is the likely best solution 

## Plot the quality indices
plot(wardRange, stat=c("ASWw", "HG", "PBC", "HC"), norm="zscore") 

## cut the dendrogram at 7 clusters 
cluster7 <- cutree(wardCluster, k=7)
table(cluster7) # the distribution of cluster membership 

# Use package "ggseqplot" to produce state distribution plot by cluster membership 
ggseqdplot(fertility.seq, group=cluster7, facet_ncol=3, facet_nrow=3)

Partitioning around medoids (PAM)

PAM aims to “obtain the best partitioning of a data set into a predefined number k of groups.” It aims to identify the k best representatives of groups, called medoids. A medoid is defined as the observation of a group having the smallest weighed sum of distances from the other observations in this group (Studer, 2013:10)

## Considering the result of hierarchical clustering, we set 7 as the number of groups:
pamcluster7 <- wcKMedoids(OM, k=7)

## print the medoids in SPS format
print(fertility.seq[unique(pamcluster7$clustering),], format="SPS")

##      Sequence                               
## 9020 (C=0,18)                               
## 8970 (C=0,5)-(C=1,13)                       
## 1521 (C=0,3)-(C=1,4)-(C=2,11)               
## 634  (C=0,9)-(C=1,3)-(C=2,6)                
## 846  (C=0,13)-(C=1,4)-(C=2,1)               
## 7996 (C=0,2)-(C=1,3)-(C=2,3)-(C=3,8)-(C=4,2)
## 8812 (C=0,5)-(NA,13)

The medoids represented in SPS format can inform the naming of the clusters.

## Calculate quality measures 
pamRange <- wcKMedRange(OM, kvals=2:20)
summary(pamRange, max.rank=2) 

##      1. N groups     1.  stat 2. N groups     2.  stat
## PBC            6 6.690566e-01           8 6.514679e-01
## HG            20 9.695848e-01          16 9.669670e-01
## HGSD          20 9.693359e-01          16 9.667432e-01
## ASW            9 5.063722e-01          10 4.966375e-01
## ASWw           9 5.069428e-01          10 4.972795e-01
## CH             2 3.112973e+03           3 2.893972e+03
## R2            20 7.751902e-01          19 7.714938e-01
## CHsq           6 6.599013e+03           9 6.509341e+03
## R2sq          20 9.217484e-01          19 9.160056e-01
## HC            16 2.411017e-02          20 2.459379e-02

The quality measures show that 6 or 9 clusters may be better than 7 clusters. For illustration purpose, 7 is used. We can always inspect different solutions to make sure that the clustering makes substantive sense while improving quality measures.

Naming the clusters

Labeling the clusters is an important step. Here, each cluster is named according to the medoid of each cluster. The labels have to be succinct and informative to account for the state change and temporality.

## Label the resulting clusters yielded by PAM
type <-factor(pamcluster7$clustering , levels=c(8982, 8933, 1520, 634, 845, 7969, 8778),
              labels=c("Childless",
                       "Early 20, single child",
                       "Early 20, two children",
                       "Mid-20, two children",
                       "Early-30, two children",
                       "Late-teens, multiple children",
                       "incomplete"))
table(type) 

## type
##                     Childless        Early 20, single child 
##                          2675                          1032 
##        Early 20, two children          Mid-20, two children 
##                           984                          1009 
##        Early-30, two children Late-teens, multiple children 
##                          1103                          1300 
##                    incomplete 
##                           881

## Plot the state distribution by cluster membership
ggseqdplot(fertility.seq, group=type, facet_ncol=3, facet_nrow=3)

From the state distribution plots stratified by PAM cluster membership, we can see that the most prevalent fertility pattern is the childless from age 18 to 35 (n=2,675). The other clusters are distinct from each other in terms of timing and number of children. There is also a cluster (i.e., “incomplete”) representing the sequences loaded with missing values.