Thursday, April 12, 2018

K is for Cohen's Kappa

Last April, during the A to Z of Statistics, I blogged about Cohen's kappa, a measure of interrater reliability. Cohen's kappa is a way to assess whether two raters or judges are rating something the same way. And thanks to an R package called irr, it's very easy to compute. But first, let's talk about why you would use Cohen's kappa and why it's superior to a more simple measure of interrater reliability, interrater agreement.

I often do research that requires another person to observe the same thing and make their own ratings, using a codebook or similar method. Meta-analysis, in which information from studies on a topic is coded, frequently requires judgment calls. While some things may be very straightforward to code, such pulling out a group sample size that is clearly stated, other things are not; the rater may need to make some decisions about quality of the methods used or exactly what sampling approach was selected, because researchers may use different and/or vague language to describe things. Since the coded data is what ultimately gets analyzed, we need to make sure the coding is done in a way that is systematic and reproducable. Qualitative research, in which the things people say are coded, also requires a codebook that is clear, systematic, and reproducable, and once again, the best way to demonstrate that is to have another person use the same data and codebook and see if they get the same results.

So you want to make sure the degree to which two coders agree on the coded results is high. A simple way of doing that is to look at interrater agreement: the number of times raters agree divided by the number of things being rated. The problem is that, when raters are working with a codebook with a limited number of categories to choose from, they're likely to agree to a certain extent just by chance alone. Even a stopped clock is right twice a day, and even untrained raters coding things willy nilly are going to agree with each other some of the time. In fact, a lot of things we want to happen in research will simply happen by chance alone. Being a good researcher means making certain that the things that happen in our research are unlikely to be due to chance. Cohen's kappa corrects for that, by taking into account how often raters will agree if they were simply to make decisions at random.

You want to set up your data with each coder getting his/her own column. You can put all coded information in a single file, if you'd like, and simply reference the columns you need for your interrater reliability function. For the demonstration with real data (below), I just created two separate files, one for each variable I'm demonstrating.

But first, let's demonstrate with some randomly generated data. Pretend that I have two coins, and I'm going to flip each of them, one and then the other, 20 times. We would expect the resulting pairs of 20 coin flips to be the same at least some of the time. We can easily generate these data using the binomial distribution. I've assigned a theta (probability of a certain outcome) at 0.5, to recreate a "fair" coin. Then I used the cbind (column bind) function to put them together into a data frame.

theta = 0.5
N = 20
flips1 <- rbinom(n = N, size = 1, prob = theta)
flips2 <- rbinom(n = N, size = 1, prob = theta)
coins<-cbind(flips1, flips2)

Now we have a data frame called coins, which contains two columns: flips for coin 1 and flips for coin 2. The irr package will measure simple agreement for us.

install.packages("irr")
library(irr)
## Loading required package: lpSolve
agree(coins, tolerance=0)
##  Percentage agreement (Tolerance=0)
## 
##  Subjects = 20 
##    Raters = 2 
##   %-agree = 40

By setting tolerance to 0, I've forced the agree function to require both columns to have the exact same value for it to be considered agreement. If I was assessing agreement on a rating scale, I might want to allow a small margin of error - perhaps 1 point. As you can see, agreement is 40%, very close to what you would expect by chance alone. And this highlights the issue with using percent agreement: we would expect two raters coding something with 2 categories to agree with each other 50% of the time. This is about how much agreement you would see between two raters who are given a codebook with absolutely no training, though if the categories are even slightly well-defined, you'll see higher agreement just by chance. So training is important, but then, so is using a measure of reliability that takes into account the agreement you would see just by chance.

Now let's run Cohen's kappa on these data.

kappa2(coins)
##  Cohen's Kappa for 2 Raters (Weights: unweighted)
## 
##  Subjects = 20 
##    Raters = 2 
##     Kappa = -0.237 
## 
##         z = -1.08 
##   p-value = 0.279

Kappa is a negative value, showing that they are doing worse than chance - very poor interrater reliability. This would tell me - if I wasn't using randomly generated data - that the codebook is poorly defined and doing little good for my raters, and/or that I may need to retrain my raters.

Now let's demonstrate interrater agreement and Cohen's kappa using some real data. For my meta-analysis, I had a fellow grad student go through and code studies with me. Since I coded many variables for my meta-analysis, and I want to keep this post as short as possible, I've selected 2 to use for this demonstration - 1 that showed poor agreement/kappa initially and 1 that showed high agreement/kappa. I adopted a consensus approach to coding, meaning that when my fellow coder and I disagreed, we met to discuss and come to a decision on how to deal with the discrepant code. Sometimes we changed the codebook as a result, sometimes one or both of us misunderstood the study (and found a better code after reexamining it together), and sometimes we simply had to compromise. We started this process early, getting together after we'd each coded a few studies solo and continuing to meet after coding 3-4 studies each. If we changed the codebook, we'd have to recode earlier studies, and of course, code all new studies with the updated codebook.

The first variable that showed disagreement surprised me: the number of studies in the article that was eligible for the meta-analysis. I was a bit surprised that we disagreed, but I realized, after seeing her coded results that I had not been clear in how I wanted to divide up any subsamples. That discussion led to a better codebook. I've created a tab-delimited file that includes a variable for study ID, then how rater1 and rater2 coded each study on that variable.

numstudies<-read.delim("num_studies.txt", header=TRUE)
agree(numstudies[,2:3], tolerance=0)
##  Percentage agreement (Tolerance=0)
## 
##  Subjects = 62 
##    Raters = 2 
##   %-agree = 79
kappa2(numstudies[,2:3])
##  Cohen's Kappa for 2 Raters (Weights: unweighted)
## 
##  Subjects = 62 
##    Raters = 2 
##     Kappa = 0.521 
## 
##         z = 6.22 
##   p-value = 5.12e-10

Our percent agreement is about 79%, but once you account for chance agreement, our Cohen's kappa is much lower: 0.52. You see why a discussion, and a better codebook, was the right approach here. On the other hand, we showed much better agreement and Cohen's kappa for a variable assessing instructions received by the control group: 0 = Nothing, 1 = A news article not about any kind of crime, 2 = A news article about crime in general, but not the specific case, 3 = A news article about the specific case that contained only neutral information.

CGinstruct<-read.delim("CG_instruct.txt", header=TRUE)
agree(CGinstruct[,2:3], tolerance=0)
##  Percentage agreement (Tolerance=0)
## 
##  Subjects = 97 
##    Raters = 2 
##   %-agree = 96.9
kappa2(CGinstruct[,2:3])
##  Cohen's Kappa for 2 Raters (Weights: unweighted)
## 
##  Subjects = 97 
##    Raters = 2 
##     Kappa = 0.954 
## 
##         z = 14.5 
##   p-value = 0

As you can see, this showed much better results: 97% agreement and a Cohen's kappa of 0.95.

In a publication, you'd want to provide Cohen's kappa for each variable or, if there are a lot of variables, some summary statistics, including range and average. (But this might also be a sign that you have too many variables and should elect only the most important ones for analysis. I ended up dropping some variables, not only because of coding results, but because some key information was missing from most studies.) If you have any Cohen's kappa not in the 0.8 or 0.9 range, you probably want to consider updating your codebook and/or retraining your coders to make sure everyone is on the same page. You also want to come up with a game plan of how to handle disagreements, at the very least because you need to pick a final value to use in your analysis. I prefer the consensus approach myself, but some people will enlist a third coder as a tie breaker.

Wednesday, April 11, 2018

J is for jsonlite Package

J is for jsonlite package Today I'm going to introduce a new method of storing and exchanging data: JSON or JavaScript Object Notation. Up to now, we've been working with delimited text files and R data frames. But JSON (pronounced "Jason") is another way we can store data that can be read by different software packages. In my previous job, some of our tests arrived in Research as JSON files, which we then parsed using Python, SAS, or R. In fact, JSON can be parsed by any programming language. JSON allowed the transfer of large amounts of data, organized into name and value pairs. Each name and value pair is separated by commas, with each object (case) enclosed in curly brackets {}, and also separated by commas. The full dataset is then enclosed in square brackets []. JSON files are human readable; they are self-describing, because you get to pick whatever name to assign to a value, and hierarchical.

Here's what a JSON file might look like for my Blogging A to Z posts:

{"posts": [
  {"postname": "A is for (Cronbach's) Alpha", "date": "20180401", "shorturl": "a-is-for-cronbachs-alpha.html", "posted": true},
  {"postname": "B is for Betas (Standardized Regression Coefficients)", "date": "20180402", "shorturl": "b-is-for-betas-standardized-regression.html", "posted": true},
  {"postname": "C is for Cross Tabs Analysis", "date": "20180403", "shorturl": "c-is-for-cross-tabs-analysis.html", "posted": true},
  {"postname": "D is for Data Frame", "date": "20180404", "shorturl": "d-is-for-data-frame.html", "posted": true},
  {"postname": "E is for Effect Sizes", "date": "20180405", "shorturl": "e-is-for-effect-sizes.html", "posted": true},
  {"postname": "F is for (Confirmatory) Factor Analysis", "date": "20180406", "shorturl": "f-is-for-confirmatory-factor-analysis.html", "posted": true},
  {"postname": "G is for glm Function", "date": "20180407", "shorturl": "g-is-for-glm-function.html", "posted": true},
  {"postname": "H is for Help with R", "date": "20180409", "shorturl": "h-is-for-help-with-r.html", "posted": true},
  {"postname": "I is for (Classical) Item Analysis or I Must Be Flexible", "date": "20180410", "shorturl": "i-is-for-classical-item-analysis-or-i.html", "posted": true},
  {"postname": "J is for jsonlite Package", "date": "20180411", "shorturl": null, "posted": false}
]}

As you can see, the structure is readable and you can make sense out of what information it is communicating. JSON allows many kinds of data, including numeric, string, logical, and null values. This was one reason it was so useful for our test data; because some of our tests were adaptive, each examinee only received certain items, so they would have null values for most of the items in the item bank. We could read in their responses to the items they saw, along with the item ID, and fill in null values for other items they didn't see. We can then put all examinees in a single file, with those who saw an item having values in that column, and those who didn't with null values. Then we can analyze all examinees together and generate item statistics and/or person ability estimates.

JSON does not allow functions or dates, though, so I've create my date variable as a string, enclosed in quotes. To read that information as a date, I would need to do an extra step once I parse it into R, but that's only necessary if you plan on doing any kind of analysis or calculations with dates. For instance, you might have a date of birth variable and want to calculate exact age, using current date, for everyone in your sample. In that case, you'd want to make certain that whatever statistical package you're using knows the variable is a date so it can handle it properly in calculations.

White space is ignored in JSON, with brackets dictating hierarchy and structure, so I could space this file out more if I wanted to, to make it even more readable:

{"posts": [

  {"postname": "A is for (Cronbach's) Alpha",

  "date": "20180401",

  "shorturl": "a-is-for-cronbachs-alpha.html",

  "posted": true}

]}

I saved the object created above, using a simple text editor, as a .json file, which I can then read into R with the jsonlite package. Though the jsonlite package has a way of coercing an object into a data frame, I found it a bit finicky, so I just read the object in then converted it to a data frame.

install.packages("jsonlite")
library(jsonlite)
posts<-fromJSON("posts.json")
posts<-as.data.frame(posts)

Now I have a data frame called "posts", containing all of the information from my JSON file. Let's take a look at how the data were read in, in particular the data types.

str(posts)
## 'data.frame': 10 obs. of  4 variables:
##  $ posts.postname: chr  "A is for (Cronbach's) Alpha" "B is for Betas (Standardized Regression Coefficients)" "C is for Cross Tabs Analysis" "D is for Data Frame" ...
##  $ posts.date    : chr  "20180401" "20180402" "20180403" "20180404" ...
##  $ posts.shorturl: chr  "a-is-for-cronbachs-alpha.html" "b-is-for-betas-standardized-regression.html" "c-is-for-cross-tabs-analysis.html" "d-is-for-data-frame.html" ...
##  $ posts.posted  : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...

If I want to do any kind of date math, I need to convert my post.date column into a date variable. I just need to tell R to turn it into a date and provide the format of the string. (Nerdy note: date variables are actually represented as the number of seconds since January 1, 1970, known as the Unix epoch. This is then converted into a date, formatted in whatever way you specify.)

posts$posts.date <- as.Date(posts$posts.date, "%Y%m%d")
str(posts$posts.date)
##  Date[1:10], format: "2018-04-01" "2018-04-02" "2018-04-03" "2018-04-04" "2018-04-05" ...

Now I can use that variable to compute a new variable - days since posted.

posts$days.since.post <- Sys.Date() - posts$posts.date
str(posts)
## 'data.frame': 10 obs. of  5 variables:
##  $ posts.postname : chr  "A is for (Cronbach's) Alpha" "B is for Betas (Standardized Regression Coefficients)" "C is for Cross Tabs Analysis" "D is for Data Frame" ...
##  $ posts.date     : Date, format: "2018-04-01" "2018-04-02" ...
##  $ posts.shorturl : chr  "a-is-for-cronbachs-alpha.html" "b-is-for-betas-standardized-regression.html" "c-is-for-cross-tabs-analysis.html" "d-is-for-data-frame.html" ...
##  $ posts.posted   : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
##  $ days.since.post:Class 'difftime'  atomic [1:10] 8 7 6 5 4 3 2 0 -1 -2
##   .. ..- attr(*, "units")= chr "days"

But the jsonlite package will not only parse JSON; it can also create a JSON file, for easy sharing. Let's convert the Facebook data file into a JSON file. We'll also add an additional argument, pretty, which adds whitespace to make the file more readable.

Facebook<-read.delim(file="small_facebook_set.txt", header=TRUE)
Facebook_js<-toJSON(Facebook, dataframe=c("rows","columns","values"), pretty=TRUE)
save(Facebook_js, file="FB_JS.JSON")

If you're interested in learning more about JSON files, check out the tutorial on W3Schools.com.

Tuesday, April 10, 2018

I is for (Classical) Item Analysis or I Must Be Flexible

I is for ITEMAN Back when I worked at HMH, I discovered an R package called ITEMAN, which is used for classical item analysis. I've mentioned classical test theory before, which focuses on the overall test or measure, as opposed to individual items. Tests and measures developed with classical test theory can't really be divided up in the way tests and measures developed with item response theory can. But you can still get some useful item statistics when adopting a classical test theory approach, through classical item analysis.

The main item statistic generated in classical item analysis is a P value, not to be confused with the p-value generated in inferential statistical analysis. In this context, P refers to difficulty, and it is abbreviated as P because it is the proportion or percentage of examinees who get the item correct. If almost no one gets the item correct, it is a difficult item. If almost everyone gets the item correct, it is an easy item.

The problem here is that P value is entirely sample-dependent. If you have an exceptionally capable sample take your test, your items will all look easy, even if relatively speaking they are not. Item response theory and Rasch, on the other hand, provide item difficulty that is not sample dependent. If you look at the math behind IRT and Rasch, you can see exactly where sample is being controlled for and therefore partialled out. (Sadly, you'll have to take all of that at face value, because going into how IRT and Rasch is not sample dependent goes beyond the scope of this post/series. But maybe I need to do an A to Z of Rasch next year!)

Basically, when doing classical item analysis, where the capability of your sample can completely change your item statistics, it becomes even more important to validate content and have experts on hand to help determine what items are appropriate for different ability groups. It also highlights the importance of sampling when piloting the test or measure.

In my previous job, I began working on a cognitive ability test developed with classical test theory. This was a fixed form test, or rather, a set of fixed form tests that were written for specific ages. So an 8-year old would take the form created for 8-year olds, which would contain items appropriate for that age level as well as the age levels on either side (7 and 9). Later ages were often combined - for instance, 13/14, 15/16, and 17/18. This was not an adaptive test, but was one way to see if children and adolescents are performing at their age-level with a paper and pencil, multiple-choice test that could be given to a group of pretty much any size. As I've said before, adaptive testing is sometimes the best way, but not always. There are other considerations, such as ease of administration, cost, and what you plan to do with the information.

The other psychometrician on the project was using a piece of software called ITEMAN. But because quality control was one of our top priorities, I was asked to perform my own psychometric analysis with a different piece of software that would give us the same analysis. Fortunately, I found the ITEMAN R package. This allowed us to validate each other, and showed that the ITEMAN package provided nearly exact results to the ITEMAN software.

So my plan was to do a tutorial today on ITEMAN, a rather simple package that comes with two sample datasets and two functions, ITEMAN1 and ITEMAN2. ITEMAN1 is for use with multiple choice tests that have single correct answers; the sample dataset, dichotomous, is for use with ITEMAN1. ITEMAN2 is for polytomous items - items that use rating scales, such as attitude measures; the sample dataset, timms2011_usa works with this function. As I was beginning to write this post, I discovered ITEMAN was no longer available, and I couldn't even find old package files on GitHub. Hence the subtitle for today, I Must Be Flexible. After some searching, I found another package, CTT, that gives similar results to ITEMAN.

So let's demonstrate scoring and classical item analysis with CTT. I pulled an exam key for an old exam I gave in a course I taught about 9 years ago. This was a 30 item test that had 25 multiple choice items, plus 5 short answer and essay items. I generated some random data as responses to the 25 multiple choice items so we can apply the answer key and run item statistics. The dataset contains responses (A, B, C, or D) for all 25 items by 20 students. I'll read in that data, then create an object containing the exam key.

exam_data<-read.delim("madeup_testdata.txt",header=TRUE)
items_only<-exam_data[,2:26]
exam_key<-c("B","C","C","D","A","C","C","B","C","C","A","A","C","C","A","B",
            "B","B","D","A","C","A","A","C","C")

Now I need to have CTT score my test using the key I provided. I requested output.scored, so it gives me a matrix of scored results I can then reference with the name of the score object (exam_score) + $scored. The resulting scored data is then used for item analysis. Be sure to load the CTT package (install first if you haven't yet).

install.packages("CTT")
## Installing package into '\\marge/users$/slocatelli/My Documents/R/win-library/3.4'
## (as 'lib' is unspecified)
library(CTT)
exam_score<-score(items_only, key=exam_key, output.scored=TRUE, rel=TRUE)
## You will find additional options and better formatting using itemAnalysis().
## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero
report<-itemAnalysis(exam_score$scored, itemReport=TRUE)
## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

## Warning in cor(items[, i], Xd): the standard deviation is zero

You notice I'm getting some errors, likely because of some very easy items that everyone got correct. We'll just ignore those for now and move on to examining our results. As part of the itemAnalysis function, CTT creates a data frame called itemReport, which we can access like this:

report[["itemReport"]]
##    itemName  itemMean        pBis         bis alphaIfDeleted
## 1        X1 0.6666667  0.39053161  0.50632146      0.7420752
## 2        X2 0.8571429  0.49165409  0.76245000      0.7361525
## 3        X3 0.4761905  0.02570432  0.03223647      0.7747674
## 4        X4 1.0000000          NA          NA      0.7583576
## 5        X5 0.9047619  0.56092808  0.97239993      0.7346370
## 6        X6 1.0000000          NA          NA      0.7583576
## 7        X7 0.6666667  0.28566255  0.37035946      0.7512439
## 8        X8 1.0000000          NA          NA      0.7583576
## 9        X9 0.4285714  0.63662447  0.80260597      0.7184512
## 10      X10 0.9523810 -0.12122355 -0.26026053      0.7661433
## 11      X11 1.0000000          NA          NA      0.7583576
## 12      X12 0.6666667  0.28566255  0.37035946      0.7512439
## 13      X13 0.9047619  0.56092808  0.97239993      0.7346370
## 14      X14 1.0000000          NA          NA      0.7583576
## 15      X15 0.9523810  0.44556792  0.95661070      0.7445828
## 16      X16 0.4761905  0.73184243  0.91782302      0.7084914
## 17      X17 0.6190476 -0.10360238 -0.13203544      0.7841816
## 18      X18 1.0000000          NA          NA      0.7583576
## 19      X19 1.0000000          NA          NA      0.7583576
## 20      X20 0.4761905  0.73184243  0.91782302      0.7084914
## 21      X21 0.5714286  0.30228059  0.38109155      0.7503664
## 22      X22 0.8571429  0.39819222  0.61751070      0.7423402
## 23      X23 1.0000000          NA          NA      0.7583576
## 24      X24 1.0000000          NA          NA      0.7583576
## 25      X25 0.7142857  0.52730093  0.70081320      0.7301665

itemMean gives us our P value, pBis is our point-biserial (the correlation between score on that item and total score on the other items), biserial (almost the same as point biserial - correlation between item score and total score without that item - but it treats the item differently: as ordinal with an underlying continuity, instead of as a discrete, 0 or 1, value), and alpha if that item were deleted.

This report tells us that there many easy items - p-values of 0.9 and greater - but there are also some moderately difficult to difficult items, which not a bad thing. What is problematic is that some of the items don't seem to relate to overall performance on the test - which can be seen in the low and sometimes negative point-biserial correlations. This gives me some guidance on what items I might want to potentially drop. But deleting any items isn't going to have too much of an impact on reliability of the measure, which sits in the 0.7 range.

CTT has more to offer that I won't go into now, but you can read more about them in the manual available here. Could you do this kind of analysis without a package like CTT? Absolutely. One of the benefits of classical test theory and item analysis approaches is that many of these analyses can be done by hand or simple software, while IRT and Rasch approaches are nearly impossible without computer assistance. If your data is set up as 0=incorrect and 1=correct, you could calculate P values by simply taking the mean. Point-biserial correlations would take a bit more work but are still completely doable on your own. But the nice thing about this package is that it automates the process, so you don't need to write your own functions, and produces a report of all item statistics.

Monday, April 9, 2018

H is for Help with R

My goal this month is to offer some fun tutorials and introductions to the R statistical language and R Studio software. But you may find, at some point, that you need a little...


And that's okay! R can have a very steep learning curve, but even for an expert, you need help figuring things out every once in a while. I constantly look things up while I'm using R - there are very few things I remember exactly how to do off the top of my head. R Studio is great in that it gives you some pop-up guidance when writing code and allows you to search for help in package manuals. But that won't always cut it. To get that help, you need to have some idea of where to begin - which library, which function, and so on. Here are my favorite places to go when I need help doing something in R:
  • Stack Overflow R - This link takes you directly to questions tagged as R; at the top of the screen, you can type in details about your question to narrow it down.
  • Quick-R - This site is completely dedicated to how to do different things in R, in particular statistical analysis. When I can't remember the syntax to create a specific kind of table or need a quick refresher on a plot, I go here.
  • R-bloggers - This page aggregates posts from other bloggers about R; it's a curated feed of some of the best R posts out there. While I tend to go to the first two links on this list for quick help, this is where I go for in-depth tutorials.
  • Variance Explained - One of my favorite blogs - each post features walk-throughs and tons of code by R programmer David Robinson. If you want to behold the awesome things R can do, and maybe get some inspiration for your own awesome things, start here. He also has some great R tutorials.
  • Learning and Using R at Stanford - No, I didn't go to Stanford. Fortunately, you don't have to either to access some of these great resources. (Note, some do require a login, but many are free.) Once again, this is a great site to look for tutorials on specific topics.
  • All else fails, LMGTFY - But in all seriousness, if I can't find the answers I need above, I just do a web search. 
And if you're more of a book learner/lover, here are some of my favorite books on R:
  • R in a Nutshell - I've been carrying this book with me at almost all times for the last month or so, as I've prepared for this blog challenge.
  • A Beginner's Guide to R - Sadly, this book might be outdated at this point, but this was the book I picked up in 2009 when I decided to try this crazy idea of dumping SPSS and going full R (without ever having used R before). Though it sounds like hubris (I was a 5th year doctoral student, so yes), it was also because I was broke and wanted to use something free (I was a 5th year doctoral student, after all). When I found I was in over my head, this book kept me from going back on my idea.
Finally, if you start feeling comfortable with R and want to take things a step farther:
  • Latent Variable Modeling with R - Especially if you'll be doing a lot of psychometric or structural equation modeling work, this is a great book to have. It uses many of the packages I draw upon regularly (and would recommend above similar packages) like lavaan and ltm.
  • Text Mining with R: A Tidy Approach - Julia Silge and David Robinson (above) go through exactly how to use their R package, tidytext, which David frequently uses on his blog.

Sunday, April 8, 2018

Statistics Sunday: Interpreting Confirmatory Factor Analysis Output

It's a Deeply Trivial first - I made my first video! In it, I talk through conducting and interpreting the output of the analyses from the F is for (Confirmatory) Factor Analysis post.



Also referenced in the video:
Hope you enjoy my first ever video! Perhaps I'll make more in the future.

Saturday, April 7, 2018

G is for GLM Function

In the Beta post, I used linear regression and demonstrated how to request standardized regression coefficients (betas). You may remember that I mentioned in that post that there are other types of regression. Linear regression is used with a continouous outcome. But sometimes you want to predict outcomes that aren't continuous - perhaps they're binary outcomes (0 or 1 coded in some meaningful way), counts, or any other types of outcome variables that don't follow the typical normal distribution. Fortunately, there are many types of regression to allow you to work with these different types of variables. You can access many of them with the glm function, a built-in (base R) function for conducting a variety of regression models.

GLM stands for general linear model, which is the basis for many statistical analyses, including regression and structural equation modeling. It provides a mathematical method of relating predictor variables to outcomes in terms of an equation, converting values on the predictor variable(s) to values on the outcome variable.

Using the glm function is very similar to using the lm function for a linear model - you need to symbolically represent your equation:

Y ~ x1 + x2 + ... + lastx

And you need some kind of data object for R to apply the equation to. A key difference in the glm function is that you reference the kind of regression model to use, through the "family" argument. Here are the different family options:
  • Binomial
  • Gaussian
  • Gamma
  • Inverse.Gaussian
  • Poisson
  • Quasi
  • Quasibinomial
  • Quasipoisson
The family refers to the distribution that best represents your outcome data. For instance, Gaussian refers to the normal distribution - unless you specify some additional arguments in your glm function, these results will be essentially the same as using the lm function. Binomial is used for two-level outcomes, which is what I'll demonstrate below. Poisson is used for count data. Quasibinomial is used when you have additional variance in your outcome not explained by the binomial distribution alone (that is, there is more variance than if it followed the binomial distribution perfectly); it provides similar results as binomial, except with an additional error term. Quasipoisson is similar for count data. We won't worry about those for now though - just remember that anytime you have to estimate something extra, you need to provide more data, so each additional term increases how much data you must collect. We'll come back to this topic when we get to power analysis.

So let's have some fun using the GLM function. First up, let's run a binomial regression. But we'll need some binary data first.

Sara's Dissertation

I've been sitting on my dissertation data for several years now and have gotten some fun presentations out it. To briefly summarize, I conducted my dissertation, which if you're so inclined you can read it here, on pretrial publicity: participants were exposed to a certain kind of pretrial publicity at random or to a control condition with no biasing information, then they completed a measure of justice attitudes, read a trial transcript, rendered a verdict, and provided a guilt rating. I included both verdict and guilt rating because of my meta-analysis, which included some studies that used guilt ratings instead of or in addition to verdicts. I was curious how these two measures related to each other, and what impact attitudes might have. So in addition to my planned analyses on pretrial publicity, I conducted some exploratory analyses using the attitude items.

The measure I used provided many attitude items, but based on a literature search, I picked out a handful of specific attitudes that appeared to impact how people make decisions about guilt versus innocence. Though I've chided other researchers for examining individual items instead of overall measure scores, there are some cases where it makes sense to think in terms of specific versus general attitudes. The theoretical background for my dissertation had to do with the tenuous relationship between attitudes and behaviors. Attitudes are poor predictors of behaviors, though more specific attitudes are better predictors than more general attitudes. My hypothesis was that pretrial publicity was not in itself biasing, unless you hold an attitude that a specific piece of information implies a defendant is guilty. And I conducted my exploratory analysis with much the same framework, that specific attitudes relating to the case and factors of the trial would be better predictors of the outcome.

First, I ran a regression with a handful of attitudes: belief that courts should be able to use illegal evidence (e.g., obtained without a search warrant), beliefs about circumstantial evidence, belief that the defendant should be required to testify (even though defendants don't have to do anything, even mount a defense), belief that police abuse power, an attitude that police should be allowed to arrest anyone who seems "suspicious", belief that police are "overzealous", and finally an attitude that upstanding citizens have nothing to fear from police. I conducted this analysis once with verdict as the outcome, and again with guilt rating as the outcome.

Since we want to get started with our glm function, let's run a binomial regression with verdict as the outcome.

dissertation<-read.delim("dissertation_data.txt",
                         header=TRUE)
verdict_binomial<-glm(verdict ~ illev + circumst + deftest + policepower +
                        suspicious + overzealous + upstanding, family=binomial,
                       data=dissertation)
summary(verdict_binomial)
## 
## Call:
## glm(formula = verdict ~ illev + circumst + deftest + policepower + 
##     suspicious + overzealous + upstanding, family = binomial, 
##     data = dissertation)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0267  -1.0100  -0.5999   1.0981   2.0656  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.57140    0.81612  -0.700 0.483834    
## illev        0.13836    0.10659   1.298 0.194280    
## circumst    -0.43640    0.12770  -3.417 0.000632 ***
## deftest      0.09870    0.10859   0.909 0.363398    
## policepower -0.13188    0.11256  -1.172 0.241354    
## suspicious   0.37546    0.12157   3.088 0.002013 ** 
## overzealous -0.01246    0.09991  -0.125 0.900717    
## upstanding   0.22768    0.10820   2.104 0.035362 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 490.08  on 354  degrees of freedom
## Residual deviance: 445.22  on 347  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 461.22
## 
## Number of Fisher Scoring iterations: 4

You would interpret the output in much the same way you would linear regression - significant variables indicate those significantly related to the outcome. I could use this resulting equation to predict a person's verdict based on their attitudes on these items. As I found in my dissertation analysis, beliefs about circumstantial evidence, police ability to arrest "suspicious people", and that upstanding (law-abiding) citizens have no reason to fear the police significantly affected whether the participant convicted in this particular case. The next step in my analysis would be to test how well the equation does at predicting verdict: how often it correctly and incorrectly classified people. (Spoiler alert: I did this as part of my dissertation and found that, though 3 indicators were significant, this equation wasn't great at predicting verdict, doing so correctly only about 55 percent of the time.)

I also conducted a linear regression using guilt rating. Just for fun, let's conduct an lm as well as a glm of the Gaussian family to see how results match up:

guilt_lm<-lm(guilt ~ illev + circumst + deftest + policepower +
                        suspicious + overzealous + upstanding,
                        data=dissertation)
summary(guilt_lm)
## 
## Call:
## lm(formula = guilt ~ illev + circumst + deftest + policepower + 
##     suspicious + overzealous + upstanding, data = dissertation)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0357 -0.7452  0.1828  0.9706  2.5013 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.16081    0.38966  10.678  < 2e-16 ***
## illev        0.11111    0.05816   1.911  0.05689 .  
## circumst    -0.08779    0.06708  -1.309  0.19147    
## deftest     -0.02020    0.05834  -0.346  0.72942    
## policepower  0.02828    0.06058   0.467  0.64090    
## suspicious   0.17286    0.06072   2.847  0.00468 ** 
## overzealous -0.03298    0.04792  -0.688  0.49176    
## upstanding   0.08941    0.05374   1.664  0.09706 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.169 on 347 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.07647, Adjusted R-squared:  0.05784 
## F-statistic: 4.105 on 7 and 347 DF,  p-value: 0.0002387
guilt_gaus<-glm(guilt ~ illev + circumst + deftest + policepower +
                        suspicious + overzealous + upstanding,
                        family="gaussian", data=dissertation)
summary(guilt_gaus)
## 
## Call:
## glm(formula = guilt ~ illev + circumst + deftest + policepower + 
##     suspicious + overzealous + upstanding, family = "gaussian", 
##     data = dissertation)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.0357  -0.7452   0.1828   0.9706   2.5013  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.16081    0.38966  10.678  < 2e-16 ***
## illev        0.11111    0.05816   1.911  0.05689 .  
## circumst    -0.08779    0.06708  -1.309  0.19147    
## deftest     -0.02020    0.05834  -0.346  0.72942    
## policepower  0.02828    0.06058   0.467  0.64090    
## suspicious   0.17286    0.06072   2.847  0.00468 ** 
## overzealous -0.03298    0.04792  -0.688  0.49176    
## upstanding   0.08941    0.05374   1.664  0.09706 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 1.366385)
## 
##     Null deviance: 513.40  on 354  degrees of freedom
## Residual deviance: 474.14  on 347  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 1128.2
## 
## Number of Fisher Scoring iterations: 2

As you can see, the regression results are the same, though the output is slightly different between the two. The lm function gives you your R-squared and F-test for the regression (test that any indicators are significant), while the glm function gives you dispersion parameters and AIC. For this reason, if you're conducting a linear regression, use the lm function. The glm function only makes sense here if you need to specify additional arguments to use a different member of the Gaussian family. The results may be the same, but the applications are slightly different.

I also did a second set of binomial regressions for my dissertation, this one dealing with variable that might affect how a person translates their guilt rating (which was on a scale of 1 to 7, the likelihood that the defendant actually committed the crime) to a verdict. That is, I found that while most people didn't select a verdict of guilty unless they also selected a guilt rating of 6 or 7, some convicted at a much lower rating, while others didn't select guilty even when they gave a guilt rating of 7, which translated to "definitely committed the crime". So I conducted two binomial regressions: one with criminal justice attitudes (e.g., how they think "reasonable doubt" should be defined, belief in the fairness of juries, and so on) as well as guilt rating, and one that also looked at how each attitude item interacts with guilt rating. This lets us see if the impact of guilt rating on verdict depends on a person's attitude. To make things easier to read, I'm going to set up my two equations separately from the glm function.

predictors_only<-'verdict ~ obguilt + reasdoubt + bettertolet + libertyvorder + 
                  jurevidence + guilt'
pred_int<-'verdict ~ obguilt + reasdoubt + bettertolet + libertyvorder + 
                  jurevidence + guilt + obguilt*guilt + reasdoubt*guilt +
                  bettertolet*guilt + libertyvorder*guilt + jurevidence*guilt'
model1<-glm(predictors_only, family="binomial", data=dissertation)
summary(model1)
## 
## Call:
## glm(formula = predictors_only, family = "binomial", data = dissertation)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4814  -0.6050  -0.1296   0.6419   2.7575  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -10.21158    1.33072  -7.674 1.67e-14 ***
## obguilt         0.39024    0.16010   2.437   0.0148 *  
## reasdoubt      -0.12315    0.11591  -1.062   0.2880    
## bettertolet    -0.10054    0.11406  -0.881   0.3781    
## libertyvorder  -0.01425    0.14961  -0.095   0.9241    
## jurevidence     0.08448    0.09897   0.854   0.3933    
## guilt           1.81769    0.20222   8.989  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 490.08  on 354  degrees of freedom
## Residual deviance: 309.66  on 348  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 323.66
## 
## Number of Fisher Scoring iterations: 5
model2<-glm(pred_int, family="binomial", data=dissertation)
summary(model2)
## 
## Call:
## glm(formula = pred_int, family = "binomial", data = dissertation)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6101  -0.5432  -0.1289   0.6422   2.2805  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)  
## (Intercept)         -12.84571    6.10651  -2.104   0.0354 *
## obguilt              -0.34506    1.13742  -0.303   0.7616  
## reasdoubt             1.56658    0.83360   1.879   0.0602 .
## bettertolet          -0.21819    0.86374  -0.253   0.8006  
## libertyvorder         0.88325    0.95459   0.925   0.3548  
## jurevidence          -0.62756    0.94042  -0.667   0.5046  
## guilt                 2.43022    1.19628   2.031   0.0422 *
## obguilt:guilt         0.13062    0.21752   0.600   0.5482  
## reasdoubt:guilt      -0.33930    0.16323  -2.079   0.0376 *
## bettertolet:guilt     0.01426    0.16662   0.086   0.9318  
## libertyvorder:guilt  -0.17482    0.18666  -0.937   0.3490  
## jurevidence:guilt     0.14006    0.18359   0.763   0.4456  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 490.08  on 354  degrees of freedom
## Residual deviance: 300.66  on 343  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 324.66
## 
## Number of Fisher Scoring iterations: 6

This matches what I found initially: that there was a significant interaction between guilt rating and attitudes about reasonable doubt. So depending on how an individual defined reasonable doubt, they would convict at lower or higher guilt ratings.

That's all for now! Tune in tomorrow for my first ever Deeply Trivial video on interpreting CFA output! And back to A to Z posts on Monday!

Friday, April 6, 2018

F is for (Confirmatory) Factor Analysis

Title Back in February and March, I shared three posts on structural equation modeling: one introducing exogenous and endogenous variables, one introducing factor analysis, and one on how factor analysis can be used in psychometric research. Today, I'll demonstrate how to conduct confirmatory factor analysis on measurement (psychometric) data - and look for my next Statistics Sunday post where we dig into understanding and interpreting output from this analysis.

To give a bit of background, I learned structural equation modeling with LISREL. LISREL (LInear Structural RELations) is not only software to conduct SEM; it's also a specific type of notation, developed by Jöreskog, Keesling, and Wiley, that breaks a structural model into matrices. Conducting SEM with LISREL involves typing up a complex control file, defining each matrix by its Greek name from the notation system - Lambda x's and y's to describe factor loadings for exogenous and endogenous variables, respectively; phi to describe the variances and covariances of exogenous latent variables; and so on.

Matrices are an essential part of the math behind statistical analysis, and are also a valuable shorthand to describe complex relationships. When LISREL software was first developed in the 1970s, it made sense that matrices were specified to break down the huge amount of data that would need to be processed for a structural model. Even then, this kind of analysis took time and a lot of computing power.

Today, I'm running on my tiny Macbook Air something that required a mainframe. I don't want to downplay the importance of this matrix approach, because that's still going on on the back end, but when computational power comes cheap and easy, it feels needlessly complicated or, at the very least, outdated to require each matrix be specified to conduct the analysis. And it creates a high intellectual bar for anyone who wants to learn SEM. Is it nice to know that the gamma matrix describes causal paths between exogenous and endogenous variables? Sure. But is it essential to conduct and understand SEM? No, I don't think so.

If you're like me, when you design a structural model to describe your data, you're drawing your variables, some of which are measured directly (observed), some of which are assessed indirectly by combinations of measured variables (latent), and drawing arrows to show how they relate to each other. To use the LISREL approach, I would have to take what I drew and convert it to matrices before I could conduct my analysis. And I would constantly have to double-check which matrix referred to what. Wouldn't it be nice to go straight from model drawing to analysis?

My favorite way to conduct SEM is with an R package called lavaan. Rather than specifying matrices, lavaan has you convert the components of your model into equations, which is a much more direct conversion from a model drawing. A factor is defined followed by the observed variables assessing that factor, as follows:

Factor =~ var1 + var2 + ... + lastvar

It's structured like a linear equation because, at the basic level, that's exactly what it is. You could conduct a factor analysis with 1 or more factors, and factors can be correlated or uncorrelated (also referred to as orthogonal). Hold that thought for now - later on, I'll talk about how orthogonal models can be used for hypothesis testing.

Using my Facebook dataset, I could test many factor models. Three of the measures included in the study assess one topic without subscales, so those could be tested with a single factor model. For simplicity, let's start with one of those. The Diener Satisfaction with Life Scale is a 5-item measure that assesses how satisfied one is with his/her life overall. We would begin by specifying a single factor model. We can name the factor and the overall model whatever we want, but need to use the actual variable names when referring to the observed variables used to assess the factor:

Facebook<-read.delim(file="small_facebook_set.txt", header=TRUE)
SWL_Model<-'SWL =~ LS1 + LS2 + LS3 + LS4 + LS5'

Here's what this model looks like drawn out:


This object we created is then used in lavaan to fit our model. You'll want to load lavaan (install if necessary) before proceeding to the next step:

install.packages("lavaan")
## Installing package into '\\marge/users$/slocatelli/My Documents/R/win-library/3.4'
## (as 'lib' is unspecified)
library(lavaan)
## This is lavaan 0.5-23.1097
## lavaan is BETA software! Please report any bugs.

Now, let's fit our model and check our output:

SWL_Fit<-cfa(SWL_Model, data=Facebook)
summary(SWL_Fit)
## lavaan (0.5-23.1097) converged normally after  24 iterations
## 
##   Number of observations                           257
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic               26.760
##   Degrees of freedom                                 5
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   SWL =~                                              
##     LS1               1.000                           
##     LS2               0.974    0.071   13.625    0.000
##     LS3               0.969    0.065   14.909    0.000
##     LS4               0.855    0.071   12.038    0.000
##     LS5               0.790    0.085    9.290    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .LS1               0.984    0.111    8.869    0.000
##    .LS2               0.890    0.102    8.742    0.000
##    .LS3               0.493    0.073    6.744    0.000
##    .LS4               1.124    0.115    9.793    0.000
##    .LS5               2.092    0.197   10.641    0.000
##     SWL               1.678    0.229    7.332    0.000

Seriously, that's all it takes to conduct SEM with lavaan. We'll dig into interpretation later - for now, note that all observed variables load significantly onto the latent variable SWL - but two things you will want to add when examining the summary output is the standardized solution and fit measures:

summary(SWL_Fit, standardized=TRUE, fit.measures=TRUE)
## lavaan (0.5-23.1097) converged normally after  24 iterations
## 
##   Number of observations                           257
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic               26.760
##   Degrees of freedom                                 5
##   P-value (Chi-square)                           0.000
## 
## Model test baseline model:
## 
##   Minimum Function Test Statistic              635.988
##   Degrees of freedom                                10
##   P-value                                        0.000
## 
## User model versus baseline model:
## 
##   Comparative Fit Index (CFI)                    0.965
##   Tucker-Lewis Index (TLI)                       0.930
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -2111.647
##   Loglikelihood unrestricted model (H1)      -2098.267
## 
##   Number of free parameters                         10
##   Akaike (AIC)                                4243.294
##   Bayesian (BIC)                              4278.785
##   Sample-size adjusted Bayesian (BIC)         4247.082
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.130
##   90 Percent Confidence Interval          0.084  0.181
##   P-value RMSEA <= 0.05                          0.003
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.040
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   SWL =~                                                                
##     LS1               1.000                               1.295    0.794
##     LS2               0.974    0.071   13.625    0.000    1.262    0.801
##     LS3               0.969    0.065   14.909    0.000    1.256    0.873
##     LS4               0.855    0.071   12.038    0.000    1.107    0.722
##     LS5               0.790    0.085    9.290    0.000    1.023    0.578
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .LS1               0.984    0.111    8.869    0.000    0.984    0.370
##    .LS2               0.890    0.102    8.742    0.000    0.890    0.359
##    .LS3               0.493    0.073    6.744    0.000    0.493    0.238
##    .LS4               1.124    0.115    9.793    0.000    1.124    0.478
##    .LS5               2.092    0.197   10.641    0.000    2.092    0.666
##     SWL               1.678    0.229    7.332    0.000    1.000    1.000

Next, let's go through how you could specify a multiple factor model. The Ruminative Response Scale assesses 3 subscales: Depression-Related Rumination, Reflecting, and Brooding. So we can specify a 3 factor model, enclosing all 3 equations, one for each factor, within the quote marks. We'll then fit this model in the same way we did our Satisfaction with Life model above. By default, the 3 factors are allowed to correlate:

RRS_Model<- '
  Depression =~ Rum1 + Rum2 + Rum3 + Rum4 + Rum6 + Rum8 + 
    Rum9 + Rum14 + Rum17 + Rum18 + Rum19 + Rum22
  Reflecting =~ Rum7 + Rum11 + Rum12 + Rum20 + Rum21
  Brooding =~ Rum5 + Rum10 + Rum13 + Rum15 + Rum16
'
RRS_Fit<-cfa(RRS_Model, data=Facebook)
summary(RRS_Fit)
## lavaan (0.5-23.1097) converged normally after  40 iterations
## 
##   Number of observations                           257
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic              600.311
##   Degrees of freedom                               206
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   Depression =~                                       
##     Rum1              1.000                           
##     Rum2              0.867    0.124    6.965    0.000
##     Rum3              0.840    0.124    6.797    0.000
##     Rum4              0.976    0.126    7.732    0.000
##     Rum6              1.167    0.140    8.357    0.000
##     Rum8              1.147    0.141    8.132    0.000
##     Rum9              1.095    0.136    8.061    0.000
##     Rum14             1.191    0.135    8.845    0.000
##     Rum17             1.261    0.141    8.965    0.000
##     Rum18             1.265    0.142    8.887    0.000
##     Rum19             1.216    0.135    8.992    0.000
##     Rum22             1.257    0.142    8.870    0.000
##   Reflecting =~                                       
##     Rum7              1.000                           
##     Rum11             0.906    0.089   10.138    0.000
##     Rum12             0.549    0.083    6.603    0.000
##     Rum20             1.073    0.090   11.862    0.000
##     Rum21             0.871    0.088    9.929    0.000
##   Brooding =~                                         
##     Rum5              1.000                           
##     Rum10             1.092    0.133    8.216    0.000
##     Rum13             0.708    0.104    6.823    0.000
##     Rum15             1.230    0.143    8.617    0.000
##     Rum16             1.338    0.145    9.213    0.000
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   Depression ~~                                       
##     Reflecting        0.400    0.061    6.577    0.000
##     Brooding          0.373    0.060    6.187    0.000
##   Reflecting ~~                                       
##     Brooding          0.419    0.068    6.203    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .Rum1              0.687    0.063   10.828    0.000
##    .Rum2              0.796    0.072   11.007    0.000
##    .Rum3              0.809    0.073   11.033    0.000
##    .Rum4              0.694    0.064   10.857    0.000
##    .Rum6              0.712    0.067   10.668    0.000
##    .Rum8              0.778    0.072   10.746    0.000
##    .Rum9              0.736    0.068   10.768    0.000
##    .Rum14             0.556    0.053   10.442    0.000
##    .Rum17             0.576    0.056   10.370    0.000
##    .Rum18             0.611    0.059   10.418    0.000
##    .Rum19             0.526    0.051   10.352    0.000
##    .Rum22             0.609    0.058   10.428    0.000
##    .Rum7              0.616    0.067    9.200    0.000
##    .Rum11             0.674    0.069    9.746    0.000
##    .Rum12             0.876    0.080   10.894    0.000
##    .Rum20             0.438    0.056    7.861    0.000
##    .Rum21             0.673    0.068    9.867    0.000
##    .Rum5              0.955    0.090   10.657    0.000
##    .Rum10             0.663    0.065   10.154    0.000
##    .Rum13             0.626    0.058   10.819    0.000
##    .Rum15             0.627    0.064    9.731    0.000
##    .Rum16             0.417    0.050    8.368    0.000
##     Depression        0.360    0.072    4.987    0.000
##     Reflecting        0.708    0.111    6.408    0.000
##     Brooding          0.455    0.096    4.715    0.000

If we wanted our 3 factors to be uncorrelated, we would update the RRS_Fit statement as follows:

RRS_Fit<-cfa(RRS_Model, data=Facebook, orthogonal=TRUE)

One more thing, just for fun - instead of having 3 correlated factors, we might instead believe that our 3 factors in turn assess a single, higher-order factor of Rumination. In fact, this type of model makes more sense with this measure. That is, all 22 items can be summed to obtain an overall Rumination score, and combinations of items are summed to get the 3 subscale scores. If we wanted to test this type of model, we would use the following:

RRS_HO_Model<- '
  Depression =~ Rum1 + Rum2 + Rum3 + Rum4 + Rum6 + Rum8 + 
    Rum9 + Rum14 + Rum17 + Rum18 + Rum19 + Rum22
  Reflecting =~ Rum7 + Rum11 + Rum12 + Rum20 + Rum21
  Brooding =~ Rum5 + Rum10 + Rum13 + Rum15 + Rum16
  Rumination =~ Depression + Reflecting + Brooding
'

RRS_HO_Fit<-cfa(RRS_HO_Model, data=Facebook)
summary(RRS_HO_Fit)
## lavaan (0.5-23.1097) converged normally after  33 iterations
## 
##   Number of observations                           257
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic              600.311
##   Degrees of freedom                               206
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   Depression =~                                       
##     Rum1              1.000                           
##     Rum2              0.867    0.124    6.965    0.000
##     Rum3              0.840    0.124    6.797    0.000
##     Rum4              0.976    0.126    7.732    0.000
##     Rum6              1.167    0.140    8.357    0.000
##     Rum8              1.147    0.141    8.132    0.000
##     Rum9              1.095    0.136    8.061    0.000
##     Rum14             1.191    0.135    8.845    0.000
##     Rum17             1.261    0.141    8.965    0.000
##     Rum18             1.265    0.142    8.887    0.000
##     Rum19             1.216    0.135    8.992    0.000
##     Rum22             1.257    0.142    8.870    0.000
##   Reflecting =~                                       
##     Rum7              1.000                           
##     Rum11             0.906    0.089   10.138    0.000
##     Rum12             0.549    0.083    6.603    0.000
##     Rum20             1.073    0.090   11.862    0.000
##     Rum21             0.871    0.088    9.929    0.000
##   Brooding =~                                         
##     Rum5              1.000                           
##     Rum10             1.092    0.133    8.216    0.000
##     Rum13             0.708    0.104    6.823    0.000
##     Rum15             1.230    0.143    8.617    0.000
##     Rum16             1.338    0.145    9.213    0.000
##   Rumination =~                                       
##     Depression        1.000                           
##     Reflecting        1.123    0.144    7.780    0.000
##     Brooding          1.047    0.148    7.067    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .Rum1              0.687    0.063   10.828    0.000
##    .Rum2              0.796    0.072   11.007    0.000
##    .Rum3              0.809    0.073   11.033    0.000
##    .Rum4              0.694    0.064   10.857    0.000
##    .Rum6              0.712    0.067   10.668    0.000
##    .Rum8              0.778    0.072   10.746    0.000
##    .Rum9              0.736    0.068   10.768    0.000
##    .Rum14             0.556    0.053   10.442    0.000
##    .Rum17             0.576    0.056   10.370    0.000
##    .Rum18             0.611    0.059   10.418    0.000
##    .Rum19             0.526    0.051   10.352    0.000
##    .Rum22             0.609    0.058   10.428    0.000
##    .Rum7              0.616    0.067    9.200    0.000
##    .Rum11             0.674    0.069    9.746    0.000
##    .Rum12             0.876    0.080   10.893    0.000
##    .Rum20             0.438    0.056    7.861    0.000
##    .Rum21             0.673    0.068    9.867    0.000
##    .Rum5              0.955    0.090   10.657    0.000
##    .Rum10             0.663    0.065   10.154    0.000
##    .Rum13             0.626    0.058   10.819    0.000
##    .Rum15             0.627    0.064    9.731    0.000
##    .Rum16             0.417    0.050    8.368    0.000
##     Depression        0.003    0.017    0.204    0.838
##     Reflecting        0.259    0.051    5.072    0.000
##     Brooding          0.064    0.026    2.446    0.014
##     Rumination        0.356    0.073    4.858    0.000

That's all for now! Check back for future posts to better understand output and how to assess and improve your models.