Latent class mixed models – with graphics

d2 <-lcmm(y~x,random=~x,subject='id',mixture=~x,ng=2,idiag=TRUE,data=dummy,link="linear")
summary(d2) 

Here I am modeling the dependent y with the independent x, allowing x to vary per person (random slope) and using (linear) x as the variable in the mixture term. Below I will look at more appropriate models, am keeping it simple for now. In the output (below) we see the results from the main section of the summary. This highlights, with the ‘x 1′ and ‘x 2′ line, that each class has a different slope, one slightly positive and one slightly negative, with (see ‘intercept 2′) quite different intercepts.

Fixed effects in the class-membership model:
(the class of reference is the last class)
coef           Se       Wald     p-value
intercept    0.4665299   0.1712701     2.723943  0.006450762

Fixed effects in the longitudinal model:
coef           Se       Wald     p-value
intercept 1 (not estimated)                      0.000000000
intercept 2 -4.476219941 0.2091687393 -21.400043 0.000000000
x 1         -0.002339034 0.0008328176  -2.808579 0.004976061
x 2          0.023911923 0.0011908951  20.078950 0.000000000

Goodness of fit statistics are also provide. Of interest is the posterior probabilities of group membership. Each subject is assigned a probability of belong to each (of the 2) class(es).

postprob(d2)
Posterior classification:
class1 class2
N 152.00000 93.00000
%  62.04082 37.95918

Posterior classification table:
--> mean of posterior probabilities in each class
prob1 prob2
class1 0.9330 0.0665
class2 0.0934 0.9070

We see one class has 62% of subjects and the other has 38%. We see that of those in class 1 their mean posterior probability of belonging to class 1 was 0.93 while their mean PP of belonging to class 2 was 0.09. Similarly, of those assigned to class 2, their mean PP of belonging to class 2 was 0.90 and 0.07 for belonging to class 1. We can look further at this with the code below. The minimums are close to 0.5 (somewhat expected), so we know that some subjects trajectories could have really been in either class. I take solace in both lower quartiles being >0.80.

> round(summary(as.numeric(d2$pprob[d2$pprob[,"class"]==1,"prob1"])),2)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.51 0.94 1.00 0.93 1.00 1.00

> round(summary(as.numeric(d2$pprob[d2$pprob[,"class"]==2,"prob2"])),2)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.54 0.83 0.99 0.91 1.00 1.00

Given we are dealing with trends and trajectories, we might like to visualise what this means for our data. Using the code below, we can pull out which class subjects were assigned to, and then plot their data, colouring based on class membership.

dummy$id <- as.character(dummy$id)
people2 <- as.data.frame(d2$pprob[,1:2])
dummy$group2 <- factor(people2$class[sapply(dummy$id, function(x) which(people2$id==x))])
p1 <- ggplot(dummy, aes(x, y, group=id, colour=group2)) + geom_line() + geom_smooth(aes(group=group2), method="loess", size=2, se=F)  + scale_y_continuous(limits = c(13,37)) + labs(x="x",y="y",colour="Latent Class") + opts(title="Raw")
p2 <- ggplot(dummy, aes(x, y, group=id, colour=group2)) + geom_smooth(aes(group=id, colour=group2),size=0.5, se=F) + geom_smooth(aes(group=group2), method="loess", size=2, se=T)  + scale_y_continuous(limits = c(13,37)) + labs(x="x",y="y",colour="Latent Class") + opts(title="Smoothed", legend.position="none")
grid.arrange(p1,p2, ncol=2, main="2 Latent Class")

This produces the following graph.

Here we see what we wanted to see. For starters the lines when smoothing across all data within a class intersect, suggesting we are seeing a real effect of classes. I would not be too interested, myself, in seeing two parallel lines, as that would suggest everyone has the same ‘kind’ of trajectory, and I would suspect the classes to be purely an artifact of the algorithm having to ‘go out and find two classes’. In saying that, we knew from the coefficients that we would see something like this. What I particularly like is that we can see a group (class), class 2, who clearly rocket up at the end. When x hits 200, a vast majority of the subjects in class 2 rapidly increase, where as class 1 generally remain consistent (in y).

If we try 3 classes.

And 4 classes?

As far as choosing how many classes are appropriate, the generally accepted method is to use the BIC (as opposed to the AIC), with the lower BIC being the preferred model. With this data we get a BIC of 32941 for 4 classes, 32972 for 3 classes and 33084 or 2 classes. This suggests 4 classes offers a slightly better fit than 3. Actually truth be told, 5 classes achieves a slightly lower BIC again. However like all modeling, you have to have some human (or really really smart robot circa 2035) input. With 4 classes, the 4th class only holds 10% of subjects, in our reduced modeling subset of the data set, this is 25 people. You have to decide if this is appropriate for your data or not. Here I certainly wouldn’t go any lower than 10%.

Refining the model

As I mentioned earlier, a linear effect may not be the best way to model this data. Here is my preferred number of classes, 3, with a quadratic term introduced. This takes substantially longer to complete (converge). Compare this graph below, to the graph 2 above (for 3 latent classes). There are some subtle differences

As posterior probabilities, and percentage of class membership is concerned, the model with only the linear term and the model that also has the quadratic differ quite a bit. It looks like the quadratic term has a tougher time determining which class each subject falls into. See below.

Just linear term

Posterior classification:
class1 class2 class3
N 139.00000 59.00000 47.00000
%  56.73469 24.08163 19.18367

Posterior classification table:
--> mean of posterior probabilities in each class
prob1 prob2 prob3
class1 0.8800 0.08060 0.039600
class2 0.1390 0.86000 0.000981
class3 0.0887 0.00293 0.908000

> summary(as.numeric(d3$pprob[d3$pprob[,"class"]==1,"prob1"]))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.5024 0.8210 0.9402 0.8798 0.9931 1.0000

> summary(as.numeric(d3$pprob[d3$pprob[,"class"]==2,"prob2"]))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.5136 0.7642 0.9549 0.8597 0.9972 1.0000

> summary(as.numeric(d3$pprob[d3$pprob[,"class"]==3,"prob3"]))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.4962 0.8986 0.9898 0.9083 0.9994 1.0000

Quadratic term included

Posterior classification:
class1 class2 class3
N 110.00000 44.00000 91.00000
%  44.89796 17.95918 37.14286

Posterior classification table:
--> mean of posterior probabilities in each class
prob1 prob2 prob3
class1 0.7990 0.0237 0.1770
class2 0.0465 0.8670 0.0862
class3 0.1630 0.0482 0.7890

> summary(as.numeric(d32$pprob[d32$pprob[,"class"]==1,"prob1"]))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.3598 0.6565 0.8524 0.7992 0.9688 0.9999

> summary(as.numeric(d32$pprob[d32$pprob[,"class"]==2,"prob2"]))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.4284 0.7217 0.9847 0.8673 0.9997 1.0000

> summary(as.numeric(d32$pprob[d32$pprob[,"class"]==3,"prob3"]))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.4159 0.6348 0.7822 0.7890 0.9618 1.0000

However, the BIC actually reduces from 32972 to 32594 which is a much bigger difference than we saw by changing the number of classes. In this instance I would say a) further model exploration is required and b) if I had to, I would chose the model based on the linear term only, giving preference to the better posterior probability performance over the BIC performance.

Hopefully this has been useful as an introduction to latent class modeling and/or and introduction to the lcmm package and/or plotting and visualising longitudinal latent class mixture modeling.

Final note

I have been in contact with the author of this package, who early next year aims to enhance some of the features of the package.

*prior to running the analysis on the data from here on I limited the data set to a subset based on a few criteria, mainly to do with year and initial x value (x0). This was purely to eliminate known calendar and other effects effecting the efficiency of the algorithm. I can provide more information on this if contacted by anyone interested.

**********************************************************************************

Here is the script used in this example. I haven’t included the loading of the data as it was a data set from a project and I have subsequently applied a coefficient to the x and y variable so it is not obvious what they are, for obvious reasons. NOTE: I had to switch between the id column being a factor and character depending on which function I was running at the time. For the lcmm() it needed to be a factor, or the postprob() function would give warnings (not partitioning the data correctly), but then it needed to be a character for the sapply() function to work properly It also needed to be a character for the colour=id in ggplot to avoid the legend ‘blowing out’ with 2000 odd levels available (despite only 3 being plotted). I’m sure there was a better way to do this (probably re write the sapply(), I just ran out of time). Got a question on this, just ask.

Data used is available here. Use `load(“2011-10-04_lcmm_post_mc.rdata”)`.