Подтвердить что ты не робот

Хорошие способы визуализации продольных категориальных данных в R

[ Обновление:. Хотя я принял ответ, добавьте еще один ответ, если у вас есть дополнительные идеи визуализации (будь то в R или другом языке/программе). Тексты по анализу категориальных данных, похоже, не говорят о визуализации продольных данных, в то время как тексты по продольному анализу данных, похоже, не говорят о визуализации внутри-предметных изменений с течением времени в членстве в категории. Имея больше ответов на этот вопрос, он станет лучшим ресурсом по проблеме, которая не получает большого освещения в стандартных ссылках.]

Коллега просто дал мне продольные категориальные данные, которые нужно посмотреть, и я пытаюсь понять, как фиксировать продольный аспект в визуализации. Я размещаю здесь, потому что я хотел бы сделать это в R, но, пожалуйста, дайте мне знать, имеет ли смысл перекрестную пересылку в Cross-Validated, поскольку перекрестная регистрация обычно не рекомендуется.

Краткая справочная информация. Данные отслеживают академическую значимость от термина к термину для студентов, которые прошли академическую консультационную программу. Данные в длинном формате и имеют пять переменных: "id", "когорт", "термин", "постоянный" и "termGPA". Первые два определяют ученика и срок, в котором они участвовали в консультационной программе. Последние три являются условиями, когда студенческая академическая репутация и ГПД были зарегистрированы. Я добавил несколько примеров данных, используя dput.

Я создал мозаичный сюжет (см. ниже), который группирует студентов по когорте, положению и термину. Это показывает, какая доля студентов была в каждой академической категории в каждом членстве. Но это не отражает продольный аспект - тот факт, что отдельные ученики отслеживаются с течением времени. Я хотел бы проследить путь, которым группы студентов с определенной академической успеваемостью занимают время.

Например: ученики со статусом "AP" (академическое стажирование) осенью 2009 года ( "F09" ), какая доля была еще AP в будущем, и какая доля переместилась в другие категории (например, GS), хорошая репутация ")? Существуют ли различия между когортами с точки зрения перемещения между категориями со временем с момента вступления в консультационную программу?

Я не мог понять, как зафиксировать этот продольный аспект на графике R. Пакет vcd имеет возможности для визуализации категориальных данных, но, похоже, не рассматривает продольные категориальные данные. Существуют ли "стандартные" методы визуализации продольных категориальных данных? Есть ли у R пакеты, предназначенные для этого? Является ли длинный формат подходящим для данных такого типа, или мне будет лучше в широком формате?

Я хотел бы получить предложения по решению этой конкретной проблемы, а также предложения по статьям, книгам и т.д., чтобы больше узнать о визуализации продольных категориальных данных.

Вот код, который я использовал для создания мозаичного сюжета. Код использует данные, перечисленные ниже, с помощью dput.

library(RColorBrewer)

# create a table object for plotting
df1.tab = table(df1$cohort, df1$term, df1$standing,
            dnn=c("Cohort\nAcademic Standing", "Term", "Standing"))

# create a mosaic plot
plot(df1.tab, las=1, dir=c("h","v","h"), 
     col=brewer.pal(8,"Dark2"),
     main="Fall 2009 and Fall 2010 Cohorts")

Здесь мозаичный сюжет (вопрос стороны: есть ли способ сделать столбцы для когорты F10 сидеть прямо под и иметь ту же ширину, что и столбцы для когорты F09, даже если нет данных для некоторых членов в F10 когорта):

enter image description here

И вот данные, используемые для создания таблицы и графика:

df1 =
structure(list(id = c(101L, 102L, 103L, 104L, 105L, 106L, 107L, 
108L, 109L, 110L, 111L, 112L, 113L, 114L, 115L, 116L, 117L, 118L, 
119L, 120L, 121L, 122L, 123L, 124L, 125L, 101L, 102L, 103L, 104L, 
105L, 106L, 107L, 108L, 109L, 110L, 111L, 112L, 113L, 114L, 115L, 
116L, 117L, 118L, 119L, 120L, 121L, 122L, 123L, 124L, 125L, 101L, 
102L, 103L, 104L, 105L, 106L, 107L, 108L, 109L, 110L, 111L, 112L, 
113L, 114L, 115L, 116L, 117L, 118L, 119L, 120L, 121L, 122L, 123L, 
124L, 125L, 101L, 102L, 103L, 104L, 105L, 106L, 107L, 108L, 109L, 
110L, 111L, 112L, 113L, 114L, 115L, 116L, 117L, 118L, 119L, 120L, 
121L, 122L, 123L, 124L, 125L, 101L, 102L, 103L, 104L, 105L, 106L, 
107L, 108L, 109L, 110L, 111L, 112L, 113L, 114L, 115L, 116L, 117L, 
118L, 119L, 120L, 121L, 122L, 123L, 124L, 125L, 101L, 102L, 103L, 
104L, 105L, 106L, 107L, 108L, 109L, 110L, 111L, 112L, 113L, 114L, 
115L, 116L, 117L, 118L, 119L, 120L, 121L, 122L, 123L, 124L, 125L, 
101L, 102L, 103L, 104L, 105L, 106L, 107L, 108L, 109L, 110L, 111L, 
112L, 113L, 114L, 115L, 116L, 117L, 118L, 119L, 120L, 121L, 122L, 
123L, 124L, 125L), cohort = structure(c(1L, 1L, 1L, 1L, 2L, 1L, 
1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 
2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 
1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L), .Label = c("F09", "F10"), class = c("ordered", 
"factor")), term = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L), .Label = c("S09", "F09", "S10", 
"F10", "S11", "F11", "S12"), class = c("ordered", "factor")), 
    standing = structure(c(2L, 4L, 1L, 4L, NA, 4L, 1L, NA, NA, 
    NA, NA, 2L, 2L, 1L, 4L, 4L, 1L, 3L, NA, NA, 4L, 3L, 1L, 4L, 
    NA, 2L, 1L, 3L, 3L, NA, 1L, 2L, NA, NA, NA, NA, 2L, 4L, 3L, 
    4L, 4L, 4L, 2L, NA, NA, 4L, 2L, 4L, 4L, NA, 3L, 4L, 6L, 6L, 
    1L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 4L, 6L, 4L, 4L, 1L, 4L, 1L, 
    2L, 4L, 3L, 1L, 4L, 1L, 6L, 1L, 6L, 6L, 7L, 4L, 4L, 2L, 2L, 
    4L, 2L, 6L, 4L, 6L, 7L, 4L, 2L, 4L, 1L, 2L, 4L, 6L, 6L, 4L, 
    2L, 2L, 3L, 6L, 6L, 7L, 4L, 4L, 3L, 4L, 4L, 6L, 2L, 1L, 6L, 
    6L, 4L, 2L, 1L, 7L, 2L, 4L, 6L, 6L, 4L, 4L, 3L, 6L, 4L, 6L, 
    2L, 4L, 4L, 6L, 4L, 4L, 6L, 3L, 2L, 6L, 6L, 4L, 2L, 6L, 3L, 
    4L, 4L, 6L, 6L, 4L, 4L, 5L, 6L, 4L, 6L, 4L, 4L, 4L, 5L, 4L, 
    4L, 6L, 6L, 2L, 6L, 6L, 4L, 3L, 6L, 6L, 4L, 4L, 6L, 6L, 4L, 
    4L), .Label = c("AP", "CP", "DQ", "GS", "DM", "NE", "WD"), class = "factor"), 
    termGPA = c(1.433, 1.925, 1, 1.68, NA, 1.579, 1.233, NA, 
    NA, NA, NA, 2.009, 1.675, 0, 1.5, 1.86, 0.5, 0.94, NA, NA, 
    1.777, 1.1, 1.133, 1.675, NA, 2, 1.25, 1.66, 0, NA, 1.525, 
    2.25, NA, NA, NA, NA, 1.66, 2.325, 0, 2.308, 1.6, 1.825, 
    2.33, NA, NA, 2.65, 2.65, 2.85, 3.233, NA, 1.25, 1.575, NA, 
    NA, 1, 2.385, 3.133, 0, 0, 1.729, 1.075, 0, 4, NA, 2.74, 
    0, 1.369, 2.53, 0, 2.65, 2.75, 0, 0.333, 3.367, 1, NA, 0.1, 
    NA, NA, 1, 2.2, 2.18, 2.31, 1.75, 3.073, 0.7, NA, 1.425, 
    NA, 2.74, 2.9, 0.692, 2, 0.75, 1.675, 2.4, NA, NA, 3.829, 
    2.33, 2.3, 1.5, NA, NA, NA, 2.69, 1.52, 0.838, 2.35, 1.55, 
    NA, 1.35, 0.66, NA, NA, 1.35, 1.9, 1.04, NA, 1.464, 2.94, 
    NA, NA, 3.72, 2.867, 1.467, NA, 3.133, NA, 1, 2.458, 1.214, 
    NA, 3.325, 2.315, NA, 1, 2.233, NA, NA, 2.567, 1, NA, 0, 
    3.325, 2.077, NA, NA, 3.85, 2.718, 1.385, NA, 2.333, NA, 
    2.675, 1.267, 1.6, 1.388, 3.433, 0.838, NA, NA, 0, NA, NA, 
    2.6, 0, NA, NA, 1, 2.825, NA, NA, 3.838, 2.883)), .Names = c("id", 
"cohort", "term", "standing", "termGPA"), row.names = c("101.F09.s09", 
"102.F09.s09", "103.F09.s09", "104.F09.s09", "105.F10.s09", "106.F09.s09", 
"107.F09.s09", "108.F10.s09", "109.F10.s09", "110.F10.s09", "111.F10.s09", 
"112.F09.s09", "113.F09.s09", "114.F09.s09", "115.F09.s09", "116.F09.s09", 
"117.F09.s09", "118.F09.s09", "119.F10.s09", "120.F10.s09", "121.F09.s09", 
"122.F09.s09", "123.F09.s09", "124.F09.s09", "125.F10.s09", "101.F09.f09", 
"102.F09.f09", "103.F09.f09", "104.F09.f09", "105.F10.f09", "106.F09.f09", 
"107.F09.f09", "108.F10.f09", "109.F10.f09", "110.F10.f09", "111.F10.f09", 
"112.F09.f09", "113.F09.f09", "114.F09.f09", "115.F09.f09", "116.F09.f09", 
"117.F09.f09", "118.F09.f09", "119.F10.f09", "120.F10.f09", "121.F09.f09", 
"122.F09.f09", "123.F09.f09", "124.F09.f09", "125.F10.f09", "101.F09.s10", 
"102.F09.s10", "103.F09.s10", "104.F09.s10", "105.F10.s10", "106.F09.s10", 
"107.F09.s10", "108.F10.s10", "109.F10.s10", "110.F10.s10", "111.F10.s10", 
"112.F09.s10", "113.F09.s10", "114.F09.s10", "115.F09.s10", "116.F09.s10", 
"117.F09.s10", "118.F09.s10", "119.F10.s10", "120.F10.s10", "121.F09.s10", 
"122.F09.s10", "123.F09.s10", "124.F09.s10", "125.F10.s10", "101.F09.f10", 
"102.F09.f10", "103.F09.f10", "104.F09.f10", "105.F10.f10", "106.F09.f10", 
"107.F09.f10", "108.F10.f10", "109.F10.f10", "110.F10.f10", "111.F10.f10", 
"112.F09.f10", "113.F09.f10", "114.F09.f10", "115.F09.f10", "116.F09.f10", 
"117.F09.f10", "118.F09.f10", "119.F10.f10", "120.F10.f10", "121.F09.f10", 
"122.F09.f10", "123.F09.f10", "124.F09.f10", "125.F10.f10", "101.F09.s11", 
"102.F09.s11", "103.F09.s11", "104.F09.s11", "105.F10.s11", "106.F09.s11", 
"107.F09.s11", "108.F10.s11", "109.F10.s11", "110.F10.s11", "111.F10.s11", 
"112.F09.s11", "113.F09.s11", "114.F09.s11", "115.F09.s11", "116.F09.s11", 
"117.F09.s11", "118.F09.s11", "119.F10.s11", "120.F10.s11", "121.F09.s11", 
"122.F09.s11", "123.F09.s11", "124.F09.s11", "125.F10.s11", "101.F09.f11", 
"102.F09.f11", "103.F09.f11", "104.F09.f11", "105.F10.f11", "106.F09.f11", 
"107.F09.f11", "108.F10.f11", "109.F10.f11", "110.F10.f11", "111.F10.f11", 
"112.F09.f11", "113.F09.f11", "114.F09.f11", "115.F09.f11", "116.F09.f11", 
"117.F09.f11", "118.F09.f11", "119.F10.f11", "120.F10.f11", "121.F09.f11", 
"122.F09.f11", "123.F09.f11", "124.F09.f11", "125.F10.f11", "101.F09.s12", 
"102.F09.s12", "103.F09.s12", "104.F09.s12", "105.F10.s12", "106.F09.s12", 
"107.F09.s12", "108.F10.s12", "109.F10.s12", "110.F10.s12", "111.F10.s12", 
"112.F09.s12", "113.F09.s12", "114.F09.s12", "115.F09.s12", "116.F09.s12", 
"117.F09.s12", "118.F09.s12", "119.F10.s12", "120.F10.s12", "121.F09.s12", 
"122.F09.s12", "123.F09.s12", "124.F09.s12", "125.F10.s12"), reshapeLong = structure(list(
    varying = list(c("s09as", "f09as", "s10as", "f10as", "s11as", 
    "f11as", "s12as"), c("s09termGPA", "f09termGPA", "s10termGPA", 
    "f10termGPA", "s11termGPA", "f11termGPA", "s12termGPA")), 
    v.names = c("standing", "termGPA"), idvar = c("id", "cohort"
    ), timevar = "term"), .Names = c("varying", "v.names", "idvar", 
"timevar")), class = "data.frame")
4b9b3361

ОТВЕТЫ

Ответ 1

Вот несколько идей для построения ваших данных. Я использовал ggplot2, и я немного переформатировал данные в местах.

Рисунок 1

enter image description here Я использовал сложный планшет, чтобы подражать вашему мозаичному сюжету и решить проблему выравнивания.

Рисунок 2

enter image description here Точки данных для каждого учащегося связаны серой линией, что напоминает это о параллельном графике. Окраска точек показывает категорическое положение. Использование GPA по оси y помогает распределить точки, чтобы уменьшить перепланировку, и показывает соотношение постоянных и GPA. Основная проблема заключается в том, что многие допустимые standing datapoints выпадают из-за отсутствия соответствующего значения termGPA.

Рисунок 3

enter image description here Здесь я создал новую переменную с именем initial_standing для использования для фасетирования. Каждая панель содержит учеников, которые соответствуют как когортам, так и начальным. Построение идентификатора в виде текста делает эту цифру немного загроможденной, но может быть полезна в некоторых случаях.

Рисунок 4

enter image description here Этот сюжет похож на тепловую карту, где каждый ряд является студентом. Я контролировал порядок оси id, чтобы заставить начальные и групповые группировки оставаться вместе. Если у вас есть еще много строк, вы можете захотеть отсортировать строки по типу кластеризации.

library(ggplot2)

# Create new data frame for determining initial standing.
standing_data = data.frame(id=unique(df1$id), initial_standing=NA, cohort=NA)

for (i in 1:nrow(standing_data)) {
    id = standing_data$id[i]
    subdat = df1[df1$id == id, ]
    subdat = subdat[complete.cases(subdat), ]
    initial_standing = subdat$standing[which.min(subdat$term)]
    standing_data[i, "initial_standing"] = as.character(initial_standing)
    standing_data[i, "cohort"] = as.character(subdat$cohort[1])
}

standing_data$cohort = factor(standing_data$cohort, levels=levels(df1$cohort))
standing_data$initial_standing = factor(standing_data$initial_standing,
                                        levels=levels(df1$standing))

# Add the new column (initial_standing) to df1.
df1 = merge(df1, standing_data[, c("id", "initial_standing")], by="id")

# Remove rows where standing is missing. Make some plots tidier.
df1 = df1[!is.na(df1$standing), ]

# Create id factor, controlling the sort order of the levels.     
id_order = order(standing_data$initial_standing, standing_data$cohort)
df1$id = factor(df1$id, levels=as.character(standing_data$id)[id_order])


p1 = ggplot(df1, aes(x=term, fill=standing)) +
     geom_bar(position="fill", colour="grey20", size=0.5, width=1.0) +
     facet_grid(cohort ~ .) +
     scale_fill_brewer(palette="Set1")

p2 = ggplot(df1, aes(x=term, y=termGPA, group=id)) + 
     geom_line(colour="grey70") + 
     geom_point(aes(colour=standing), size=4) + 
     facet_grid(cohort ~ .) +
     scale_colour_brewer(palette="Set1")

p3 = ggplot(df1, aes(x=term, y=termGPA, group=id)) +
     geom_line(colour="grey70") + 
     geom_point(aes(colour=standing), size=4) + 
     geom_text(aes(label=id), hjust=-0.30, size=3) +
     facet_grid(initial_standing ~ cohort) +
     scale_colour_brewer(palette="Set1")


p4 = ggplot(df1, aes(x=term, y=id, fill=standing)) + 
     geom_tile(colour="grey20") +
     facet_grid(initial_standing ~ ., space="free_y", scales="free_y") +
     scale_fill_brewer(palette="Set1") +
     opts(panel.grid.major=theme_blank()) +
     opts(panel.grid.minor=theme_blank())

ggsave("plot_1.png", p1, width=10, height=6.25, dpi=80)
ggsave("plot_2.png", p2, width=10, height=6.25, dpi=80)
ggsave("plot_3.png", p3, width=10, height=6.25, dpi=80)
ggsave("plot_4.png", p4, width=10, height=6.25, dpi=80)

Ответ 2

Изучая мой вопрос, я нашел несколько других вариантов, которые перечислю здесь.

Ряд относительно новых пакетов R предназначены для визуализации и анализа данных "истории жизни" или "многоступенчатой ​​последовательности". Идея состоит в том, что со временем люди (или объекты) входят и выходят из разных категорий - например, изменения карьеры, брак и развод, здоровье и болезни или, в моем случае, категории академических успехов в колледже.

Пакеты

R для визуализации данных последовательности или жизненного цикла включают biograph, упомянутый @timriffe в комментарии выше, и TraMineR. Автор биографического пакета, Франс Виллекенс, имеет книгу на упаковке "Биограф". Многоэтапный анализ жизненных историй с R, который будет опубликован Springer этой осенью. TraMineR имеет подробное руководство пользователя по ссылке выше, а также более короткую статью JSS. JSS также имеет специальную проблему для моделей с несколькими состояниями в контексте анализа рисков, в которой обсуждаются дополнительные R-пакеты для многоуровневого моделирования.

Я также нашел некоторое специализированное программное обеспечение, предназначенное для визуализации перемещений между категориями с течением времени. Parallel Sets - простая, бесплатная программа для создания основных визуализаций, хотя она обладает ограниченной гибкостью. Lifeflow является более сложным. Он также бесплатный, но вы должны отправить электронное письмо создателю, запрашивающему копию.

Я добавлю более подробные сведения об этом ответе, как только у меня будет возможность попробовать эти инструменты.

Ответ 3

Мне жаль, что я не нашел ответ @bdemarest, прежде чем написал R-пакет для решения этой проблемы, но поскольку OP запросил дополнительные обновления, я поделюсь еще одним решением. На рис. 4 показано, что я называю типом горизонтальной линии.

При разработке пакета longCatEDA R мы обнаружили, что сортировка данных имеет решающее значение для создания полезных графиков (см. example(sorter) и отчет, связанный в комментарии ниже для технических подробностей), тем более, что размер проблемы стал большой. Например, мы начали проблему с ежедневными напитками (абстинент, использование, злоупотребление) для нескольких тысяч участников в течение 3 лет (более 1000 дней).

Ниже приведен код для применения горизонтальной линии к данным @eipi10. Рисунок 1 стратифицирует на term, а рис. 2 стратифицируется по первому состоянию, как показано на рисунке 4 из @bdemarest, хотя результаты не идентичны из-за внутри сортировки страт.

Рисунок 1

Horizontal Line Plot Stratified by Term

Рисунок 2

Horizontal Line Plot Stratified by Initial Status

# libraries
install.packages('longCatEDA')
library(longCatEDA)
library(RColorBrewer)

# transform data long to wide
dfw <- reshape(df1,
           timevar = 'term',
           idvar = c('id', 'cohort'),
           direction = 'wide')

# set up objects required by longCat()
y <- dfw[,seq(3,15,by=2)]
Labels <- levels(df1$standing)
tLabels <- levels(df1$term)
groupLabels <- levels(dfw$cohort)

# use the same colors as bdemarest
cols <- brewer.pal(7, "Set1")

# plot the longCat object
png('plot1.png', width=10, height=6.25, units='in', res=100)
par(bg='cornsilk3', mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
lc <- longCat(y=y, Labels=Labels, tLabels=tLabels, id=dfw$id) 
longCatPlot(lc, cols=cols, xlab='Term', lwd=8, legendBuffer=0)
legend(8.1, 25, legend=Labels, col=cols, lty=1, lwd=4)
dev.off()

# stratify by term
png('plot2.png', width=10, height=6.25, units='in', res=100)
par(bg='cornsilk3', mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
lc.g <- sorter(lc, group=dfw$cohort, groupLabels=groupLabels)
longCatPlot(lc.g, cols=cols, xlab='Term', lwd=8, legendBuffer=0) 
legend(8.1, 25, legend=Labels, col=cols, lty=1, lwd=4)
dev.off()

# stratify by first status, akin to Figure 4 by bdemarest
png('plot2.png', width=10, height=6.25, units='in', res=100)
par(bg='cornsilk3', mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
first <- apply(!is.na(y), 1, function(x) which(x)[1])
first <- y[cbind(seq_along(first), first)]
lc.1 <- sorter(lc, group=factor(first), groupLabels = sort(unique(first)))
longCatPlot(lc.1, cols=cols, xlab='Term', lwd=8, legendBuffer=0) 
legend(8.1, 25, legend=Labels, col=cols, lty=1, lwd=4)
dev.off()