Example for a latent class analysis with the poLCA-package in R

When you work with R for some time, you really start to wonder why so many R packages have some kind of pun in their name. Intended or not, the poLCA package is one of them. Today i´ll give a glimpse on this package, which doesn´t have to do anything with dancing or nice dotted dresses.

This article is kind of a draft and will be revised anytime.

The “poLCA”-package has its name from “Polytomous Latent Class Analysis”. Latent class analysis is an awesome and still underused (at least in social sciences) statistical method to identify unobserved groups of cases in your data. Polytomous latent class analysis is applicable with categorical data. The unobserved (latent) variable could be different attitude-sets of people which lead to certain response patterns in a survey. In marketing or market research latent class analysis could be used to identify unobserved target-groups with different attitude structures on the basis of their buy-decisions. The data would be binary (not bought, bought) and depending on the products you could perhaps identify a class which chooses the cheapest, the most durable, the most environment friendly, the most status relevant […] product. The latent classes are assumed to be nominal. The poLCA package is not capable of sampling weights, yet.

By the way: There is also another package for latent class models called “lcmm” and another one named “mclust“.

What does a latent class analysis try to do?
A latent class model uses the different response patterns in the data to find similar groups. It tries to assign groups that are “conditional independent”. That means, that inside of a group the correlations between the variables become zero, because the group membership explains any relationship between the variables.

Latent class analysis is different from latent profile analysis, as the latter uses continous data and the former can be used with categorical data.
Another important aspect of latent class analysis is, that your elements (persons, observations) are not assigned absolutely, but on probability. So you get a probability value for each person to be assigned to group 1, group 2, […], group k.

Before you estimate your LCA model you have to choose how many groups you want to have. You aim for a small number of classes, so that the model is still adequate for the data, but also parsimonious.
If you have a theoretically justified number of groups (k) you expect in your data, you perhaps only model this one solution. A typical assumption would be a group that is pro, one group contra and one group neutral to an object. Another, more exploratory, approach would be to compare multiple models – perhaps one with 2 groups, one with 3 groups, one with 4 groups – and compare these models against each other. If you choose this second way, you can decide to take the model that has the most plausible interpretation. Additionally you could compare the different solutions by BIC or AIC information criteria. BIC is preferred over AIC in latent class models, but usually both are used. A smaller BIC is better than a bigger BIC. Next to AIC and BIC you also get a Chi-Square goodness of fit.
I once asked Drew Linzer, the developer of PoLCA, if there would be some kind of LMR-Test (like in MPLUS) implemented anytime. He said, that he wouldn´t rely on statistical criteria to decide which model is the best, but he would look which model has the most meaningful interpretation and has a better answer to the research question.

Latent class models belong to the family of (finite) mixture models. The parameters are estimated by the EM-Algorithm. It´s called EM, because it has two steps: An “E”stimation step and a “M”aximization step. In the first one, class-membership probabilities are estimated (the first time with some starting values) and in the second step those estimates are altered to maximize the likelihood-function. Both steps are iterative and repeated until the algorithm finds the global maximum. This is the solution with the highest possible likelihood. That´s why starting values in latent class analysis are important. I´m a social scientist who applies statistics, not a statistician, but as far as i understand this, depending on the starting values the algorithm can stop at point where one (!) local maximum is reached, but it might not be the “best” local maximum (the global maximum) and so the algorithm perhaps should´ve run further. If you run the estimation multiple times with different starting values and it always comes to the same solution, you can be pretty sure that you found the global maximum.

data preparation
Latent class models don´t assume the variables to be continous, but (unordered) categorical. The variables are not allowed to contain zeros, negative values or decimals as you can read in the poLCA vignette. If your variables are binary 0/1 you should add 1 to every value, so they become 1/2. If you have NA-values, you have to recode them to a new category. Rating Items with values from 1-5 should could be added a value 6 from the NAs.

Running LCA models
First you should install the package and define a formula for the model to be estimated.

The following code stems from this article. It runs a sequence of models with two to ten groups. With nrep=10 it runs every model 10 times and keeps the model with the lowest BIC.

You´ll get the standard-output for the best model from the poLCA-package:

Generate table showing fitvalues of multiple models
Now i want to build a table for comparison of various model-fit values like this:

This table was build by the following code:

Now i calculate the Entropy (a pseudo-r-squared) for each solution. I took the idea from Daniel Oberski´s Presentation on LCA.

Sometimes, an Elbow-Plot (or Scree-Plot) can be used to see, which solution is parsimonius and has good fit-values. You can get it with this ggplot2-code i wrote:


Inspect population shares of classes
If you are interested in the population-shares of the classes, you can get them like this:

or you inspect the estimated class memberships:

Ordering of latent classes
Latent classes are unordered, so which latent class becomes number one, two, three… is arbitrary. The latent classes are supposed to be nominal, so there is no reason for one class to be the first. You can order the latent classes if you must. There is a function for manually reordering the latent classes: poLCA.reorder()
First, you run a LCA-model, extract the starting values and run the model again, but this time with a manually set order.

Now you have reordered your classes. You could save these starting values, if you want to recreate the model anytime.


This is the poLCA-standard Plot for conditional probabilites, which you get if you add “graph=TRUE” to the poLCA-call.

It´s in a 3D-style which is not really my taste. I found some code at dsparks on github or this Blog, that makes very appealing ggplot2-plots and did some little adjustments.

If you want to compare the items directly:

How to plot correlations of rating items with R

Especially when you´re starting a SEM-analysis you should have a look at your observed variables and the correlations between them. As usual with R, there´s an app a package for that. It´s called corrplot and seems to rely on base R graphics.


Pretty basic, huh? If you would rather have a good looking table, try the mighty sjPlot-package. It uses spearman-correlation as default.


Also thanks to the sjPlot-package from Daniel Lüdecke you only need this code:

to produce this plot:
sjPlot correlationmatrix

As usual, i´ve had my own thoughts on how a perfect correlation matrix for my sociological survey-data should look like. This was my idea: In the lower triangle it contains a jittered scatterplot of the responses. A scatterplot of rating-items with 5 categories does not work without jitter. Diagonally are the univariate distributions of the items as barplot, with one bar for each category. In the upper triangle you find the correlation (spearman) and p-value.

Kudos to Stackoverflow User Sany Muspratt for helping me to figure this out.

And here it is in its full beauty (and grey for cheaper printing):
Correlation plot

Plot with background in ggplot2: Visualising line-ups from Hurricane-festival 1997 – 2015

The Hurricane Festival is taking place again this june. It could be interesting to have a look on its development over the years. In this case, the amount of bands for each year.

First, i gathered some data from Wikipedia and put it in a csv-file. You can access the data here:Hurricane Festival Bands 1997-2015.

A simple barplot is the best way to plot this data. But to make it a little more appealing, i want to use a custom font and a wallpaper from the Hurricane Festival website. But before i start plotting, i need to get the data in shape.

With some Dplyr-magic we aggregate the count of bands per year:

I want to use the font “Open Sans”. You can download the font here Open Sans. Then you have to import it in R using the “extrafont”-package.

As the title of this blog-entry suggests, i also want to use jpg-background for my graph. I downloaded the wallpaper from the Hurricane Festival from here http://www.hurricane.de/de/interaktiv/downloads/.

Now we can plot the data:

Hurricane Bands from 1997-2015

With “alpha=0.85” the bars become a little bit transparent, so you can see a bit more of the background-image.


By the way, this is what the plot would look like with ggplot2-defaults:

ggplot2 default
Not that bad, for just one line of code.

Anyway, here are some further graphics:
Here, i tried to give one point for each time a band has played at hurricane-festival. I tried to use an icon (png-file) of a hand instead of a point, but haven´t figured out how to do it.
Most frequent Bands

As a last plot, i did a wordcloud from all Bands that ever played at Hurricane.
Wordcloud Hurricane-Festival

Recoding all vectors in a dataframe at once (R)

This post is just a little reminder for myself, because i had to look this up a few times.
Version 1

I use “recode” from the car package, to recode the variables. This function will be used with “apply” so it will recode all vectors in the dataframe.

Here´s how it´s done:

Version 2
If you only want to recode some of the variables in you dataframe, you can define them in a list and use this list in a for-loop.

Version 3
If the values just have to be reversed, there´s an even simpler way:

Diverging stacked barchart for plotting likert Items

Questionnaires in the social sciences often include rating items to measure the variability of peoples´ attitudes towards something. Respondents are given a statement and have to report how much they agree or disagree on a 5- or 7-point-scale. A set of rating-items like these can be combined to a likert-scale. It´s also common to build an index-value for the respondents, if the items meet certain criteria of quality. There´s still some controversy, if it´s adequate to use ranked (ordinal) data like likert-items to calculate means. Most researchers think it´s approriate if the scale has at least 5 points and the variable can be considered as an ordinal measure of a continuous attitude.

Anyway. Visualization of the data ist always a good starting point. For this purpose, there are a lot of R-Packages like the HH-Package with its Likert-Function, or the likert-package from Jason Bryer and last, but not least: The sjp.likert-Function from Daniel Lüdecke, which would be my favourite.

All these packages produce sophisticated and very appealing plots. Under its hood, the HH-package uses lattice and the likert and sjPlot package are build on ggplot2. I tried HH-package, but as a ggplot2-user i realized, it would take me too long to figure out the little details. The other two packages could do what i want, but they both need raw-data (SPSS-like) and can´t work with already aggregated data. Both also have distinct kinds of dealing with the “neutral”-category of the items.
Long story short, i decided to use ggplot2 directly instead of using packages build on ggplot2 that have developed a lot of complexity on their own.

The Plot
This plot is a small example. If the code seems too messy to you, or you think the plot can be improved: i´m always interested in how to make things better, please leave a comment.
For example, one could criticize, that the x-axis isn´t meaningful, because of the neutral-category should not be splitted in negative/positive like this. So perhaps, the vertical line and the x-axis-labels should be removed. On the other hand, the HH-Plot likert-function does it the same way. It would be possible to add percentage-values inside the stacked bars, but i think that would be too much. I decided, to make a stacked-frequency table with the sjPlot-Package to complement my likert-plot.

climate change attitudes likert items

And this is the code, i´ve written:

R-Tipps: Dplyr in Funktionen nutzen

Seit es das R-Package dplyr gibt, sind viele Aufgaben bei der Datenbearbeitung einfacher geworden. Soll diesselbe Umformung für viele Variablen durchgeführt werden, kann es eine Menge Tipparbeit sparen, wenn man den Code in eine Funktion packt.

Bei der Arbeit mit dplyr ist es wichtig, dass man plyr – sofern man es benötigt – grundsätzlich als erstes läd. Andernfalls treten mit Sicherheit an irgendeinem Punkt der Arbeit Fehlermeldungen auf. Summarise() gibt es beispielsweise in plyr, aber auch in dplyr, was zu Konflikten wie diesem führt: “Error in n() : This function should not be called directly“.

Was macht dplyr so besonders?
Dplyr-Code ist einfach lesbar, logisch und arbeitet sequenziell, da es den Pipe-Operator %>% aus aus dem MagrittR-Package nutzt.
Durch “%>%” (sprich: then) können Manipulationen in einer intuitiven Reihenfolge ausgeführt werden. Während früher eine Funktion umständlich mit verschachtelten Klammern zu schreiben war, kann man mit dplyr einfach der Intuition folgen und einen Schritt nach dem nächsten machen.

Ein einfaches Beispiel: Ich möchte alle Personen eines Datensatzes entsprechend der Kategorien einer Variable in Gruppen einteilen (z.B. Geschlecht, Einkommensstufe, Postleitzahl) und dann die Anzahl der Personen je Kategorie bestimmen. Ich wähle den dataframe aus und gebe ihn mit dem Pipe-Operator (%>%) an die group_by() Funktion von dplyr weiter. Anschließend wird der gruppierten Datensatz an die summarise()-Funktion weitergegeben, die dann mit n() die Fälle zählt. Das Ergebnis wird als dataframe-Objekt “Beispiel” abgelegt.

Durch group_by() werden Gruppen erzeugt und mittels summarise() werden die Daten je Gruppe so aggregiert, wie ich angebe. Hier wird nur n() je Gruppe gezählt.

Dies ist natürlich nur Minimalbeispiel. Um einen Eindruck davon zu bekommen, was dplyr alles kann, empfehle ich folgenden Überblick:
Data Wrangling with dplyr und tidyr Cheat Sheet.
95% aller Aufgaben, die das Umformen von Daten, erstellen neuer Variablen, Gruppieren von Daten, Auswählen von Fällen, Variablen, Werten, oder Zusammenfügen von Datensätzen erfordern, können damit erledigt werden. Und: Alle diese Schritte können miteinander kombiniert werden.

dplyr in Funktionen
Versucht man den obigen dplyr-code in eine Funktion zu schreiben, stößt man auf eine Fehlermeldung. Warum? Innerhalb von Funktionen muss eine andere Version der dplyr-Funktionen verwendet werden: Standard evaluation (SE)

  • Anstatt summarise() ==> summarise_()
  • Anstatt mutate() ==> mutate_()
  • Anstatt filter() ==> filter_()
  • usw.

Zusätzlich zu diesen SE-Versionen von dplyr-Funktionen ist die Übergabe der Input-Objekte an die Funktion leicht unterschiedlich. Entweder müssen die Objekte

  • als Formel “~ Objekt”
  • als “quote(Objekt)”
  • oder als String mit Anführungszeichen ” ‘Objekt’ “

eingefügt werden. Hadley Wickham, der Programmierer des Pakets, empfiehlt die erste Möglichkeit.

Als Funktion sieht der obige dplyr-code so aus:

Die Funktion “test.function(x,y) kann zwei Input-Objekte annehmen, für die sie den dplyr-code durchführt und das Ergebnis als “test”-dataframe abspeichert.

Komplexeres Beispiel
Es ist möglich, ganze Auswertungsprozeduren als Funktion zusammenzufassen. Der folgende Code führt eine Gewichtung durch, bildet Anteilswerte der Antworten und generiert einen Plot der Daten.

Das Ergebnis:Vergleich Gewichtung
An diesem Plot sieht man, dass vor der Gewichtung der Anteil der Männer etwas höher war, während bei den gewichteten Daten der Frauenanteil etwas erhöht wird. Da es in der Sozialstruktur Deutschlands etwas mehr Frauen als Männer gibt (hauptsächlich wegen der höheren Lebenserwartung), spiegeln die gewichteten Daten die Grundgesamtheit etwas besser wider.

Hier gibt es Weitere Informationen zu standard evaluation in dplyr
Und über die R-Console gelangt man an die vignette: vignette(“tidy-data”).

Ideen für neue Raspberry-Pi Projekte

Bisher habe ich meinen Raspberry-Pi für 3 Zwecke verwendet: Als Pi musicbox, als Mediacenter mit Raspbmc und für mein Webscraping Projekt.

Mit der Pi musicbox kann ich den Raspberry Pi an meine Stereoanlage anschließen und Musik über Spotify-Premium hören. Gesteuert wird das ganze über den Browser eines Smartphones oder PCs im gleichen Wlan-Netzwerk. Diesen Job hat jetzt ein ausgemustertes Smartphone mit installierter Spotify-App übernommen.
Raspbmc macht aus dem Raspberry Pi ein Media-Center, das an den Fernseher angeschlossen werden kann. Es ermöglicht den Zugriff auf alle möglichen Medien (Musik, Filme, Fotos) z.B. von Netzwerkfestplatten und wird ebenfalls per Smartphone oder Tablet PC gesteuert. Um Urlaubsfotos von der Netzwerkfestplatte zu gucken, brauche ich jedoch keinen Raspberry Pi – das kann mein SmartTV schon selbst.
Das dritte Projekt (Radiosong-Webscraping Projekt mit dem Raspberry Pi) ist jetzt beendet. Daher suche ich neue Aufgaben für den kleinen PC.

Momentan reizt mich folgende Idee: Ein Info-Center für das Haus.

Der Raspberry Pi wird an ein Display angeschlossen. Hier werden die Abfahrtszeiten der nächsten Bus- und Straßenbahnverbindungen angezeigt und eine “Ampel”, ob man sie noch erreicht, wenn man jetzt losgeht. Außerdem wird der Raspberry Pi mit einem Temperatur- und einem Luftfreuchtigkeitsfühler ausgestattet. Diese Informationen werden visualisiert. Hier soll ebenfalls eine “Ampel” angezeigt werden, ob die Luftfeuchtigkeit zu hoch ist und gelüftet werden muss. Außerdem könnte noch angezeigt werden, welcher Müll als nächstes abgeholt wird. Eventuell können auch Daten vom Smart-Meter (Stromzähler) abgegriffen werden. Dann wäre der aktuelle Stromverbrauch im Haus immer offensichtlich. Das könnte beim Stromsparen helfen.

Eine weitere Anwendungsmöglichkeit ist Haus-Automatisierung. Steckdosen und evtl. auch Leuchten könnten z.B. via smartphone auf dem Raspberry-Pi (Web-interface) gesteuert werden. Hierfür werden entsprechend steuerbare Steckdosen mit DIP Schalter(!) benötigt, die es aber schon günstig zu kaufen gibt: Brennenstuhl Comfort RCS 1000N Funkschalt Se. Für die Heizkörperregelung wäre die Funksteuerung ebenfalls möglich, allerdings sind entsprechende Thermostate am Heizkörper etwas teurer und zudem ist unsere Zentralheizung selbst schon auf bestimmte Uhrzeiten programmiert. Die Notwendigkeit für eine weitere Steuerungsmöglichkeit ist also fraglich. Man kann mit dem raspberry pi auch eigene Alarmanlagen realisieren, wobei ich darauf nicht allzu sehr vertrauen würde.

Das Info-Center könnte man als “Spiegel” konstruieren: Raspberry Pi Mirror. Diese Lösung ist mir aber etwas zu kostspielig, da hier ein anständiger Monitor für zerlegt wird. Interessanter wäre das Recycling eines alten Kindle-E-Readers von amazon, der per Jailbreak als Monitor umfunktionalisiert werden kann. Der Kindle hat den Vorteil, auch bei Sonneneinstrahlung ablesbar zu sein (vor allem der Paperwhite) und lange Zeit ohne Stromversorgung auszukommen.