Using ggplot2 to Graph LCA Parameter Estimates

[Kara Styck]

Hello,

I am new to R and have a question about using ggplot2 and MplusAutomation to generate graphs for LCA parameter estimates (y) by threshold (x) for all latent classes for each of my indicators (separate graphs for each indicator, i.e., estimates for item 1 across all classes, estimates for item 2 across all classes, etc).  How do I specify what variables I want MplusAutomation to select to add to the graph?  I have been successful in getting MplusAutomation to read my Mplus output files, but I am not sure what language and/or commands I am supposed to use to select the "subset" of data I want to be illustrated in the graphs.

Thanks so much in advance for your help! 

Kara

[michael.hallquist]

Hi Kara,

Maybe the most helpful thing I can provide for your circumstance is some example code for an LCA and factor mixture modeling project I did a while back. I hope this will give you a sense of how to approach plotting your data using ggplot2 and MplusAutomation. In general, they play well together and it's straightforward to visualize the results. Let me know if you have remaining questions about this.

To orient you a bit to the dataset, this was a mixed clinical sample and I was looking at items that help to identify/define personality disorders. I turned on PLOT3 so that a gh5 file would be generated. This is useful for plotting the estimated item probabilities in each class. I also used the parameter estimates to plot the thresholds in a separate graph. The graphs are attached as PDFs -- hopefully these are similar to what you're looking for.

Best wishes,

Michael

-----

library(MplusAutomation)

library(plyr)

lcaModels <- readModels("LCA")

#plot estimated probabilities for LCA

vals <- lcaModels[["dsm3_top8_LCA4class.out"]]$gh5$means.and.variances.data$estimated.probs$values

vals <- melt(vals, varnames = c('LC', 'item'), value.name="prob", na.rm = TRUE)

numLC <- max(vals$LC)

vals$LC <- factor(vals$LC)

vals$sev <- factor(rep(0:2, each=numLC))

vals$sevLabel <- factor(rep(0:2, each=numLC), levels=c(2,1,0), labels=c("Strongly Present", "Present", "Absent"), ordered=TRUE)

vals$item <- factor(rep(c("Reacts Criticism", "Seeks Reassurance", "Self-Centered", "Unique Problems", "Anger",

            "Bears Grudges", "Hurt Criticism", "Identity Disturbance"), each=numLC*3)) #x3 because items are trichotomous

vals$LC <- factor(vals$LC, levels=c(3,2,4,1), labels=c("Class 1   (n = 109)", "Class 2   (n = 62)", "Class 3   (n = 48)", "Class 4   (n = 82)"))

lca4estProbs <- vals

pdf("lca4class_estprobs.pdf", width=11, height=8)

g <- ggplot(lca4estProbs, aes(x=item, y=prob, fill=sevLabel)) + geom_bar() + facet_wrap(~LC) + theme_bw(base_size=base_size) + ylab("Criteria Proportions\n") + xlab("GPD Criterion") +

    theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.5)) + scale_fill_grey("Symptom Severity", start=0.2, end=0.8)

print(g)

dev.off()

#And for plotting the threshold estimates (this was for a factor mixture model, but the idea is the same for LCA)

#get just the fmm3 3-class solution thresholds

allThresh <- subset(fmm3Models[["dsm3_top8_FMM3_3class.out"]]$parameters$unstandardized, paramHeader=="Thresholds")

#the basic plot

pdf("fmm3_3class_thresholds.pdf", width=11, height=8)

threshp <- ggplot(allThresh, aes(x=LatentClass, y=est)) + geom_bar() + geom_errorbar(aes(ymin=est-se, ymax=est+se), width=0.5) + facet_wrap(~param, scales="free_y") +

    theme(axis.text.x=theme_text(angle=-90))

print(threshp)

dev.off()

allThresh$item <- factor(sub("^([^\\$]+)\\$.*$", "\\1", allThresh$param, perl=TRUE), 

    levels=c("BORD4", "PAR4", "DEP9", "BORD6", "NAR1", "HIS1", "HIS7", "NAR4"),

    labels=c("Anger", "Bears\nGrudges", "Hurt\nCriticism", "Identity\nDisturb", "Reacts\nCriticism", "Seeks\nReassur", "Self-\nCentered", "Unique\nProblems"))

allThresh$threshold <- factor(sub("^[^\\$]+\\$(.*)$", "\\1", allThresh$param, perl=TRUE), levels=c(1,2), labels=c("Absent vs Present", "Present vs Strong Present"))

#re-sort in terms of severity to match other plots

allThresh$LatentClass <- factor(allThresh$LatentClass, levels=c("3","2","1"), labels=c("1", "2", "3"))

#the pretty plot

pdf("fmm3_3class_thresholds_pretty.pdf", width=14, height=7)

threshp <- ggplot(allThresh, aes(x=LatentClass, y=est)) + geom_point(size=2.5) + geom_errorbar(aes(ymin=est-se, ymax=est+se), width=0.5) + 

    facet_grid(threshold~item, scales="free_y") + theme_bw(base_size=base_size) +

    theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.5)) + ylab("Threshold Estimate\n") + xlab("\nLatent Class")# + ylim(-5,10)

print(threshp)

dev.off()

---

[Kara Styck]

Michael,

Thank-you so much!  I really appreciate your help. 

Kara