Creating basic psychometric summaries in R

We just finished collecting the first cohort of the WARN-D study, our ambitious effort to create a personalized early warning system for depression. In the 90 minute self report survey, we try to get at many of the most relevant predicators for depression onset. Creating this battery took us many months, and involved a delphi study of experts.

I’m now looking into the data for the first time, and to save myself a bit of time, I created a basic R script to get rough overviews for each scale. This script uses a number of tools (exploratory and confirmatory factor analysis, principal component analysis, network analysis) to get at basic information. I thought I’d share it here if others wanted to use it. It’s set up in a way that all you need to do is select your particular items, and it should all just run automatically.

I’ll walk you through the code briefly, and paste some graphs of how this looks like in example data. In the models, we’re going to ignore distributions (skewed), ordinal vs gaussian data (I just assume data to be normal for this example), and any complications: this is really just a quick setup.

(Note: WordPress keeps screwing with my assign symbol, so I’ll use “←” here instead.)

Step 1: packages and basic data wrangling

library("readr")
library("dplyr")
library("psych")
library("qgraph")
library("bootnet")
library("OpenMx")
library("EGAnet")
library("lavaan")

I show you two ways to select subsets of your data. First, when variables are organized subsequently.

df ← read_csv("data.csv")
start ← which(names(df)=='ItemName1'); start # enter the first questionnaire item
end ← which(names(df)=='ItemName20');  end   # enter the last questionnaire item
df_small ← na.omit(select(df, start:end))

Or second, when variables you want are not subsequent:

selection ← c (
  which(names(df)=='ItemName1'),
  which(names(df)=='ItemName14'),
  which(names(df)=='ItemName27'),
  which(names(df)=='ItemName31'),
  which(names(df)=='ItemName55'),
  which(names(df)=='ItemName102'))
df_small ← na.omit(select(df, c=selection))

I also have a brief section on how to quickly make sum-scores:

df_small ← df %>%
  rowwise() %>%
  mutate(new_sumscore=sum(ItemName1,ItemName2,ItemName3))

Step 2: Descriptives

Now that we have the dataframe we want, we can run a number of psychometric analyses. I usually start with descriptives; I’m sure there are fancier ways to do this, but I really like the dfSummary function from the ‘summarytools’ package. And then I usually visualize my correlation matrix via the ‘qgraph’ package, using a spring layout where more closely related items are visualized more closely together. The vechs function takes the lower triangle of the matrix excluding the diagonal (you don’t want those pesky 1s of the diagonal in your summary).

view(dfSummary(data))
cor_mat ← cor(df_small)
summary(vechs(cor_mat))
qgraph(cor_mat, cut=0, layout="spring", 
       title=paste("Correlation matrix, mean correlation = ", round(means_cor, digits=2), sep=" "))

This results in the following graph (blue = positive correlation, red = negative), hinting at a 5-factor solution with 4 items each:

Step 3: Dimensionality, i.e. how many components, factors, or communities are there

We can do this in three ways. Let’s start with very basic eigenvalue decomposition:

ev ← plot(eigen(cor(df_small))$values, type="b",
        ylab="Value", xlab="Number of Eigenvalues")

If we were to believe in the Kaiser criterion, this would indicate 5 factors to extract:

But there are more sophisticated ways to look at this, such as parallel analysis. The fa.parallel function from the ‘psych’ package does that, with some nice visual output. You can think of this as something similar to eigenvalue decomposition above, but the Kaiser criterion (everything above 1 should be extracted) is based on simulated random data.

fac_comp ← fa.parallel(df_small)

The output of this tells us “Parallel analysis suggests that the number of factors = 5 and the number of components = 5”. The corresponding figure is:

Finally, we can use exploratory graph analysis, i.e. the function EGA from the ‘EGAnet’ package to obtain insights into the dimensionality of the data using network community structures (you can think of communities as statistically equivalent to factors). It’s a convenient tool modeling conditional dependence relations (akin to partial correlations), and the algorithm detects 5 communities (this is based on walktrap community detection).

nw_ega ← EGA(df_small)
summary(nw_ega)

Step 4: Confirmatory factor analysis

As a last step, let’s fit a 5 factor CFA. Before we do so, let’s briefly check what cronbach’s alpha for this scale would be (we wouldn’t expect it to be very high, we already have some indication it is multidimensional). I specify psych::alpha(df_small) in the code because alpha is a common function and leads to conflicts in some cases otherwise. In our case, alpha is around 0.70.

alpha_scale ← psych::alpha(df_small); alpha_scale$total$std.alpha

Next, we can fit a simple 5-factor CFA in the R package ‘lavaan’:

cfa_scale ← 'factor1  =~ ItemName1+ItemName6+ItemName11+ItemName4
              factor2  =~ ItemName2+ItemName7+ItemName12+ItemName8
              factor3  =~ ItemName3+ItemName8+ItemName13+ItemName12
              factor4  =~ ItemName4+ItemName9+ItemName14+ItemName16
              factor5  =~ ItemName5+ItemName10+ItemName15+ItemName20'
fit ← cfa(cfa_scale, data = df_small, std.lv=TRUE)
summary(fit, fit.measures = TRUE)

This results in a 5 factor model that fits quite well.

Number of observations 445 
User Model versus Baseline Model: Comparative Fit Index (CFI) 0.949 
Tucker-Lewis Index (TLI)  0.939 
Root Mean Square Error of Approximation: RMSEA 0.065 
90 Percent confidence interval - lower 0.058 
90 Percent confidence interval - upper 0.072

You can find the R-code here. Enjoy! And please keep in mind that this is a very quick and dirty setup to get a first, rough overview of your data, not a thorough psychometric investigation.