Week 8 Exercises: CFA

New packages

Make sure you have these packages installed:

• lavaan
• semPlot

Exercises for the Enthusiastic

Dataset: radakovic_das.csv

Apathy is lack of motivation towards goal-directed behaviours. It is pervasive in a majority of psychiatric and neurological diseases, and impacts everyday life. Traditionally, apathy has been measured as a one-dimensional construct but is in fact composed of different types of demotivation.

The Dimensional Apathy Scale (DAS) is a multidimensional assessment for demotivation, in which 3 subtypes of apathy are assessed:

• Executive: lack of motivation for planning, attention or organisation
• Emotional: lack of emotional motivation (indifference, affective or emotional neutrality, flatness or blunting)
• Initiation: lack of motivation for self-generation of thoughts and/or actions

The DAS measures these subtypes of apathy and allows for quick and easy assessment, through self-assessment, observations by informants/carers or administration by researchers or healthcare professionals.

You can find data for the DAS when administered to 250 healthy adults at https://uoepsy.github.io/data/radakovic_das.csv, and information on the items is below.

DAS Dictionary

All items are measured on a 6-point Likert scale of Always (0), Almost Always (1), Often (2), Occasionally (3), Hardly Ever (4), and Never (5). Certain items (indicated in the table below with a - direction) are reverse scored to ensure that higher scores indicate greater levels of apathy.

item	direction	dimension	question
1	+	Executive	I need a bit of encouragement to get things started
2	-	Initiation	I contact my friends
3	-	Emotional	I express my emotions
4	-	Initiation	I think of new things to do during the day
5	-	Emotional	I am concerned about how my family feel
6	+	Executive	I find myself staring in to space
7	-	Emotional	Before I do something I think about how others would feel about it
8	-	Initiation	I plan my days activities in advance
9	-	Emotional	When I receive bad news I feel bad about it
10	-	Executive	I am unable to focus on a task until it is finished
11	+	Executive	I lack motivation
12	+	Emotional	I struggle to empathise with other people
13	-	Initiation	I set goals for myself
14	-	Initiation	I try new things
15	+	Emotional	I am unconcerned about how others feel about my behaviour
16	-	Initiation	I act on things I have thought about during the day
17	+	Executive	When doing a demanding task, I have difficulty working out what I have to do
18	-	Initiation	I keep myself busy
19	+	Executive	I get easily confused when doing several things at once
20	-	Emotional	I become emotional easily when watching something happy or sad on TV
21	+	Executive	I find it difficult to keep my mind on things
22	-	Initiation	I am spontaneous
23	+	Executive	I am easily distracted
24	+	Emotional	I feel indifferent to what is going on around me

Here are the item numbers that correspond to each dimension.

• Executive: 1, 6, 10, 11, 17, 19, 21, 23
• Emotional: 3, 5, 7, 9, 12, 15, 20, 24
• Initiation: 2, 4, 8, 13, 14, 16, 18, 22

Question 1

Read in the data. It will need a little bit of tidying before we can get to fitting a CFA.

Remember that most of the actions needed for working with those sort of data are described in the Chapter on Data Wrangling for Questionnaires.

Hints

By the looks of things, this is what I would consider doing:

• Rename the variables to easy-to-read strings like q1, q2, q3, etc.
• Set up a data dictionary that records the text of the item q1 corresponds to, the text that q2 corresponds to, etc.
• Recode the Likert scale labels to numbers.
• Reverse-code the questions with a negative direction. Note, you don’t need to this, as they’ll just end up with loadings in the opposite direction, but I would strongly recommend it for interpretation purposes.
  
• Check if there is missing data and if there is, removing those observations.

Solution 1. First let’s just read in the dataset:

library(tidyverse)
library(lavaan)
library(semPlot)

rdas <- read_csv("https://uoepsy.github.io/data/radakovic_das.csv")
head(rdas)

# A tibble: 6 × 24
  I need a bit of encouragement …¹ `I contact my friends` I express my emotion…²
  <chr>                            <chr>                  <chr>                 
1 Often                            Almost Always          Almost Always         
2 Almost Always                    Hardly Ever            Occasionally          
3 Often                            Occasionally           Occasionally          
4 Hardly Ever                      Occasionally           Almost Always         
5 Occasionally                     Hardly Ever            Occasionally          
6 Occasionally                     Occasionally           Almost Always         
# ℹ abbreviated names: ¹​`I need a bit of encouragement to get things started`,
#   ²​`I express my emotions`
# ℹ 21 more variables: `I think of new things to do during the day` <chr>,
#   `I am concerned about how my family feel` <chr>,
#   `I find myself staring in to space` <chr>,
#   `Before I do something I think about how others would feel about it` <chr>,
#   `I plan my days activities in advance` <chr>, …

The names we’re getting are useful in that they show the items, but they’re horrible to have to use in R, so we will ideally replace them with easy to use names. Note also that the data is being read in as the actual response option - e.g., “Almost Always” - and we want to treat these as a numeric scale. So those will have to change too.

Solution 2. I like to make a “data dictionary” whenever I get data like this. While I want to rename the variables to make it easier for me to use, I also want to keep track of what the questions were.
Here I make a “tibble” (the function data.frame() would work too, tibble is just tidyverse version). I indicate what I am going to rename things as (“q1”,“q2”, …, “q24”), and then I have the current names of the variables

rdas_dict <- tibble(
  variable = paste0("q",1:24),
  item = names(rdas)
)

Doing this is really useful because I can’t keep track in my head of what “q5” was.
If I want to know, then I can just do:

rdas_dict[5,]

# A tibble: 1 × 2
  variable item                                   
  <chr>    <chr>                                  
1 q5       I am concerned about how my family feel

Now let’s actually change the names in our data to what we said we would:

names(rdas) <- paste0("q", 1:24)

Solution 3. Okay, so we have all our data in words, not numbers. Views on how to treat Likert data are mixed, but it’s very common to treat it as continuous in Psychology.

Let’s check the response values we have. Just in question 1 for now:

unique(rdas$q1)

[1] "Often"         "Almost Always" "Hardly Ever"   "Occasionally" 
[5] "Always"        NA              "Never"        

A little trick that we can use to find the unique values in an entire dataset is to quickly convert the dataframe into one big long vector. Technically, a dataframe is a “list of vectors”, and the function unlist() will remove this structure.
So we can find all the unique values in all the questions with:

unique(unlist(rdas))

[1] "Often"         "Almost Always" "Hardly Ever"   "Occasionally" 
[5] "Always"        NA              "Never"         "[NO ENTRY]"   

Perfect. So we know we have uniformity of spelling. It happens less often these days as questionnaire software is improving, but you might occasionally encounter typos in some of the questions, or things with and without capital letters (R is a bit thick, and doesn’t recognise that “Often” and “often” are the same thing).
Note that we have the 6 responses that we would expect given the description of the scale, but we also have some NA values, and some [NO ENTRY] values. Not sure how those got there.
We want to turn each “Always” in to 1, each “Almost Always” in to 2, “Often” in to 3, and so on. If we simply leave out the “[NO ENTRY]”, then this will be turned into a missing value NA, which is handy.

rdas <- rdas |> 
  mutate(across(q1:q24, ~case_match(.,
    "Always" ~ 1,
    "Almost Always" ~ 2,
    "Often" ~ 3,
    "Occasionally" ~ 4,
    "Hardly Ever" ~ 5,
    "Never" ~ 6
  )))
head(rdas)

# A tibble: 6 × 24
     q1    q2    q3    q4    q5    q6    q7    q8    q9   q10   q11   q12   q13
  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1     3     2     2     5     2     4     2     5     3     3     4     3     4
2     2     5     4     5    NA     1    NA     3    NA     3     2     3     4
3     3     4     4     5     4     4     3     3     5     2     4     3     4
4     5     4     2     4     2    NA     3     3    NA     3     4     5     3
5     4     5     4     6     3     3     3     4     4     3     3     3     2
6     4     4     2     1     2     3     2     2     1     3     2     5     4
# ℹ 11 more variables: q14 <dbl>, q15 <dbl>, q16 <dbl>, q17 <dbl>, q18 <dbl>,
#   q19 <dbl>, q20 <dbl>, q21 <dbl>, q22 <dbl>, q23 <dbl>, q24 <dbl>

Solution 4. According to the table of items, the ones which need to be reverse scored are:

reversed <- c(2,3,4,5,7,8,9,10,13,14,16,18,20,22)

For these items, we want 6s to become 1s, 5s become 2s, and so on.

The tidyverse solution shown in the readings and solutions to previous labs will work just fine, but if you’re curious, here’s a different way to accomplish the same thing using functions from base R:

rdas[, reversed] <- apply(
  rdas[, reversed], MARGIN = 2, function(x) 7-x)

Note: The above code works nicely because our dataset is currently ordered such that the first column is item 1, 2nd column is item 2, and so on. This means we can use numbers to index the appropriate variables, rather than names. It would need adjusting if, for instance, our first column contained “participant ID”, and our items only began later.

Solution 5. We haven’t learned about more sophisticated methods of handling missing data, so for now we will just remove any rows in which there is missingness - i.e., we’ll do “listwise deletion”:

compl_rdas <- na.omit(rdas)

Question 2

Specify the theoretical model proposed by Radakovic et al. 

For reference, check out the example in the readings.

dasmod <- "





"

Challenge: Before you estimate the model, how many degrees of freedom do you think the model will have? (The readings will help here!)

Hints

You’ll have to use the data dictionary to see which items are associated with which dimensions.

Solution 6. Here is the model structure, according to the list of items in the dictionary.
I’m calling my factors “Em”, “Ex”, and “BCI” for “emotional”, “executive” and “behavioural/cognitive initiation” respectively. The factor correlations will be estimated by default, but I like to write things explicitly.

dasmod <- "
Ex =~ q1 + q6 + q10 + q11 + q17 + q19 + q21 + q23
Em =~ q3 + q5 + q7 + q9 + q12 + q15 + q20 + q24
BCI =~ q2 + q4 + q8 + q13 + q14 + q16 + q18 + q22
Em ~~ Ex
Em ~~ BCI
Ex ~~ BCI
"

Solution 7. Degrees of freedom is computed as the number of “knowns” minus the number of “unknowns”.

Let’s start with figuring out the number of “knowns”: the number of values in the dataset. This number comes from the observed covariance matrix. Let’s imagine a smaller dataset with only five items. It’ll create a covariance matrix like this:

var
covar   var
covar   covar   var
covar   covar   covar   var
covar   covar   covar   covar   var

How many values are in this matrix? In the first row, there’s 1, plus the second row with 2, plus the third row with 3, plus the fourth row with 4, plus the fifth row with 5. In other words, there are

sum(1:5)

[1] 15

values in this covariance matrix.

For the present scenario with 24 items, we will have

sum(1:24)

[1] 300

values in the covariance matrix. (Twenty-four of these will be each item’s own variance, and the other 276 will be covariances between items.)

The alternative formula to calculate this is \(\frac{k \cdot (k+1)}{2}\), and we get the same number by plugging in the number of variables for \(k\): \(\frac{24 \cdot (24+1)}{2} = \frac{600}{2} = 300\)

Now let’s look at the number of “unknowns”: the number of parameters the model has to estimate. This number comes from the number of latent variables and how they relate to each item.

• Each latent variable has its own variance, and there are three latent variables, so the model will have three latent factor variances.
• Each item will load onto one latent variable, and there are 24 items, so the model will have 24 factor loadings.
• Each item will have residual factor variances, and there are 24 items, so the model will have 24 residual factor variances.

Adding these up, we get

3 + 24 + 24

[1] 51

unknown parameters.

Finally, let’s subtract the knowns from the unknowns to get the degrees of freedom:

300 - 51

[1] 249

Question 3

Estimate the model using cfa().
You can choose whether you want to standardise the latent factors or fix the first loading of each factor to be 1 (it’s the same model, just scaled differently).

Examine the model fit - does it fit well?

What modifications do the modification indices suggest? Are the top three suggestions theoretically reasonable, in your opinion?

Remember, we don’t really want to have to make modifications to our models. If you don’t need to (if the model fits well) then don’t bother! (It’s still worth looking at the modification indices though).

Hints

There’s a whole section on “model fit” in the CFA chapter!

And there’s also a whole section on model modifications.

Solution 8. Let’s fit the model!

dasmod.est <- cfa(dasmod, compl_rdas)

The standard test of model fit is a chi-squared test which compares the observed covariance matrix to the model-implied covariance matrix. Ideally, these two matrices will be fairly similar, so we want a non-significant result.

We can get the test statistic and p-value from the model’s chi-squared test as follows:

summary(dasmod.est)$test$standard

$test
[1] "standard"

$stat
[1] 274.8334

$stat.group
[1] 274.8334

$df
[1] 249

$refdistr
[1] "chisq"

$pvalue
[1] 0.1251793

So a chi-squared test with 249 degrees of freedom results in a test statistic of 274.8, associated with a p-value of 0.125. The take-away is that our observed covariance matrix is not significantly different from the model-implied covariance matrix—yay!

Next, let’s check the additional measures of global fit:

fitmeasures(dasmod.est)[c("srmr","rmsea","tli","cfi")]

      srmr      rmsea        tli        cfi 
0.05364507 0.02171599 0.96913532 0.97215469 

All looks pretty good! Cut-offs for SRMR tend to vary, with some using <0.08, or <0.09, and some being stricter with <0.05. Remember, these criteria are somewhat arbitrary.

Modification indices suggest a whole bunch of items that could have some associations beyond that modelled in the factors, but these are all weak correlations at around 0.2.

modindices(dasmod.est, sort=TRUE) |> head()

    lhs op rhs     mi    epc sepc.lv sepc.all sepc.nox
132  q6 ~~  q3 15.896  0.222   0.222    0.284    0.284
113  q1 ~~  q9 11.929 -0.172  -0.172   -0.267   -0.267
216 q19 ~~ q20 10.887  0.168   0.168    0.256    0.256
103  q1 ~~  q6  7.447  0.149   0.149    0.220    0.220
358  q4 ~~  q8  7.417  0.199   0.199    0.207    0.207
105  q1 ~~ q11  6.757 -0.149  -0.149   -0.237   -0.237

These are the top 3 being suggested. I can’t see any obvious link between any of these that would make me think they are related beyond their measuring of ‘apathy’.

rdas_dict[c(3,6),]

# A tibble: 2 × 2
  variable item                             
  <chr>    <chr>                            
1 q3       I express my emotions            
2 q6       I find myself staring in to space

rdas_dict[c(1,9),]

# A tibble: 2 × 2
  variable item                                               
  <chr>    <chr>                                              
1 q1       I need a bit of encouragement to get things started
2 q9       When I receive bad news I feel bad about it        

rdas_dict[c(19,20),]

# A tibble: 2 × 2
  variable item                                                                
  <chr>    <chr>                                                               
1 q19      I get easily confused when doing several things at once             
2 q20      I become emotional easily when watching something happy or sad on TV

Question 4

Are the (standardised) loadings all “big enough”?
There’s no clear threshold that people use here - it depends a lot on the field, and on the wordings of specific items. Ideally, the same value we used in EFA (\(\geq|0.3|\)) would be nice, but not crucial.

Hints

To get out standardised loadings, we can do:

mod.est <- cfa(model_syntax, data = ...)
summary(mod.est, std=TRUE)

And you’ll get out an extra 2 columns in the summary output.

Pay attention to the Estimate, Std.lv and Std.all columns in your output. The way I think of these columns is just to think of how we scale things in regression models:

• Estimate column : item ~ Factor
  
• Std.lv column : item ~ scale(Factor)
  
• Std.all column: scale(item) ~ scale(Factor)

So if, when we fitted the model, we had specified cfa(model, data, std.lv = TRUE), then the factor already has a variance of 1, so Scale(Factor) doesn’t do anything.

See Chapter 4#interpretation.

Solution 9. I’m not going to print all of this right now because there’s so much output, but here’s how we would find standardised loadings. We can find them in the Std.all column.

summary(dasmod.est, std = TRUE)

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  Ex =~                                                                 
    q1                1.000                               0.679    0.694
    q6                0.679    0.121    5.607    0.000    0.461    0.433
    ...
    ...

The standardised loadings are all (just) greater than \(|0.3|\). Questions 13 and 15 are very close…

rdas_dict[c(13,15),]

# A tibble: 2 × 2
  variable item                                                     
  <chr>    <chr>                                                    
1 q13      I set goals for myself                                   
2 q15      I am unconcerned about how others feel about my behaviour

Question 5

Do the factors correlate in the way you would expect?

Is more emotional apathy associated with more executive apathy? and with more initiation apathy?

Hints

If you didn’t reverse code the appropriate items, then this might get confusing, because we’d have to look at factor loadings to know in which direction the factor is going (i.e., are higher numbers “more apathy” or “less apathy”?).

If you did reverse code the appropriate items, then you’re golden, because you made them all point towards “more” apathy.

Solution 10. Here are the correlations we’re interested in. Note that what we are seeing is that the three factors are all positively correlated, but for Em and Ex this is only weak (and not significant).

This isn’t necessary a problem, it just means that these two factors are fairly distinct/orthogonal. We might want to check back in the original paper to see what they proposed!

summary(dasmod.est, std = TRUE)

...
Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  Ex ~~                                                                 
    Em                0.051    0.028    1.807    0.071    0.159    0.159
  Em ~~                                                                 
    BCI               0.043    0.018    2.369    0.018    0.263    0.263
  Ex ~~                                                                 
    BCI               0.151    0.040    3.746    0.000    0.642    0.642

Question 6

Make a diagram of the model.

Hints

For a quick look at the structure of the model, try the semPaths() function from the semPlot package Chapter 4 CFA#making diagrams.

If you were going to use this sort of diagram in a proper write-up, though, it’d be better to make a nicer graphic manually (e.g., in Powerpoint, your favourite graphics software, or semdiag).

Solution 11. Here’s one example:

library(semPlot)
semPaths(dasmod.est, whatLabels = "std", rotation=2)

There are lots of options in semPaths(), so if you can make your graphic more elaborate than this one, then be our guest!

Optional Question 7

Imagine that you’re a clinician administering the DAS to a patient. In clinical settings, it’s common practice to skip the complex factor analysis we’ve been doing here and just create a sum score or a mean score that describe a patient’s responses. Then clinicians can check whether the score is above some threshold to see whether there’s cause for concern.

For each of the dimensions of apathy in the data, calculate sum scores for each of the 250 participants.

Hints

Good ol’ rowSums() to the rescue!

Solution 12. Here are the sets of items associated with each dimension:

Exitems <- c(1,6,10,11,17,19,21,23)
Emitems <- c(3,5,7,9,12,15,20,24)
BCIitems <- c(2,4,8,13,14,16,18,22)

Again, because the item numbers correspond to the column positions in our data, we can just do rowSums indexing on those column numbers to get our scores:

compl_rdas$ExSCORE <- rowSums(compl_rdas[,Exitems])
compl_rdas$EmSCORE <- rowSums(compl_rdas[,Emitems])
compl_rdas$BCIScore <- rowSums(compl_rdas[,BCIitems])

Optional Question 8

How might you think about a sum/mean score in terms of a diagram?

Hints

What does a sum or mean score imply about how each item is weighted compared to the others? How is this different from what a more sophisticated method like EFA or CFA can do?

Solution 13. Computing sum scores can feel like a ‘model free’ calculation, but actually it does pre-suppose a factor structure, and a much more constrained one than those we have been estimating. Specifically, we’re assuming that all items contribute equally to an underlying factor, rather than being weighted differently, and that all items also have the same variance as one another.

For a full explanation of this idea, see “Thinking twice about sum scores”, McNeish & Wolf 2020.

“DOOM” Scrolling

Dataset: doom.csv

The “Domains of Online Obsession Measure” (DOOM) is a fictitious scale that aims to assess the sub types of addictions to online content. It was developed to measure 2 separate domains of online obsession: items 1 to 9 are representative of the “emotional” relationships people have with their internet usage (i.e. how it makes them feel), and items 10 to 15 reflect “practical” relationship (i.e., how it connects or interferes with their day-to-day life). Each item is measured on a 7-point likert scale from “strongly disagree” to “strongly agree”.

We administered this scale to 476 participants in order to assess the validity of the 2 domain structure of the online obsession measure that we obtained during scale development.

The data are available at https://uoepsy.github.io/data/doom.csv, and the table below shows the individual item wordings.

variable	question
item_1	i just can't stop watching videos of animals
item_2	i spend hours scrolling through tutorials but never actually attempt any projects.
item_3	cats are my main source of entertainment.
item_4	life without the internet would be boring, empty, and joyless
item_5	i try to hide how long i’ve been online
item_6	i avoid thinking about things by scrolling on the internet
item_7	everything i see online is either sad or terrifying
item_8	all the negative stuff online makes me feel better about my own life
item_9	i feel better the more 'likes' i receive
item_10	most of my time online is spent communicating with others
item_11	my work suffers because of the amount of time i spend online
item_12	i spend a lot of time online for work
item_13	i check my emails very regularly
item_14	others in my life complain about the amount of time i spend online
item_15	i neglect household chores to spend more time online

Question 9

Assess whether the 2 domain model of online obsession provides a good fit to the validation sample of 476 participants.

Solution 14. From the visual of the correlation matrix, you can see the vague outline of two groups of items correlations. Note there’s a little overlap..

doom <- read_csv("https://uoepsy.github.io/data/doom.csv")
heatmap(cor(doom))

first we write our model:

moddoom <- "
# emotional domain
emot =~ item_1 + item_2 + item_3 + item_4 + item_5 + item_6 + item_7 + item_8 + item_9
# practical domain
pract =~ item_10 + item_11 + item_12 + item_13 + item_14 + item_15
# correlated domains (will be estimated by default)
emot ~~ pract
"

Then we fit it to the data:

moddoom.est <- cfa(moddoom, data = doom)

Then we inspect that fitted model object.
I’m just going to extract the fit indices here first.

fitmeasures(moddoom.est)[c("rmsea","srmr","cfi","tli")]

     rmsea       srmr        cfi        tli 
0.07118809 0.06140483 0.84091826 0.81231930 

Uh-oh.. they don’t look great.

Question 10

Are there any areas of local misfit (certain parameters that are not in the model (and are therefore fixed to zero) but that could improve model fit if they were estimated?).

Solution 15. I’m printing out just the head(), so that I can look at the few parameters with the greatest modification indices.
The top three parameters jump out immediately to me.
item_1 and item_3 have a suggested correlation of c0.5, as do item_7 and item_8. In addition, it’s suggested that including a loading (estimated to be about 0.6) from item_10 to the emot factor would improve the model fit.

modindices(moddoom.est, sort=TRUE) |>
  head()

        lhs op     rhs     mi    epc sepc.lv sepc.all sepc.nox
50   item_1 ~~  item_3 89.051  0.561   0.561    0.460    0.460
34     emot =~ item_10 74.270  1.603   0.772    0.604    0.604
118  item_7 ~~  item_8 61.736  0.369   0.369    0.432    0.432
37     emot =~ item_13 14.708 -0.632  -0.304   -0.264   -0.264
144 item_11 ~~ item_12 11.815 -0.206  -0.206   -0.176   -0.176
133  item_9 ~~ item_10 11.562  0.137   0.137    0.207    0.207

Question 11

Beware: there’s a slightly blurred line here that we’re about to step over, and move from confirmatory back to ‘exploratory’.

Look carefully at the item wordings,do any of the suggested modifications make theoretical sense? Add them to the model and re-fit it. Does this new model fit well?

As a general heuristics:

Less contentious use of modification indices:

• residual covariances for items within a factor (essentially asserts that the two observed variables share some of their specific variance)

More contentious uses of modification indices:

• adding cross-loadings (could argue that an item loading on two factors is not a clean indicator, and so should be removed)
  
• residual covariances for items on different factors - often harder to defend
  
• changing paths between the latent variables - very definitely changing your theory!

Solution 16. There are three main proposed adjustments from our initial model:

1. item_1 ~~ item_3. These questions are both about animals. It would make sense that these are related over and above the underlying “emotional internet usage” factor.
2. item_7 ~~ item_8. These are both about viewing negative content online, so it makes sense here that they would be related beyond the ‘emotional’ factor.
   
3. emot =~ item_10. This item is about communicating with others. It currently loads highly on the pract factor too. It maybe makes sense here that “communicating with others” will capture both a practical element of internet useage and an emotional one.

Putting them all in at once could be a mistake - if we added in emot =~ item_10, then we change slightly the underlying construct of the emot factor, meaning it might make other suggested modifications (item_7 ~~ item_8) less important. It’s a bit like Whac-A-Mole - you make one modification and then a whole new area of misfits appears!

Let’s adjust our model, putting in the covariance between item_1 and item_3 in.

moddoom2 <- "
# emotional domain
emot =~ item_1 + item_2 + item_3 + item_4 + item_5 + item_6 + item_7 + item_8 + item_9
# practical domain
pract =~ item_10 + item_11 + item_12 + item_13 + item_14 + item_15
# correlated domains (will be estimated by default)
emot ~~ pract
# residual covariances
item_1 ~~ item_3
"

Then fit it to the data:

moddoom2.est <- cfa(moddoom2, data = doom)

fitmeasures(moddoom2.est)[c("rmsea","srmr","cfi","tli")]

     rmsea       srmr        cfi        tli 
0.05402943 0.05208128 0.90939355 0.89189003 

The fit is close, but not great, and the same suggested correlations are present in modification indices:

modindices(moddoom2.est, sort=TRUE) |>
  head()

        lhs op     rhs     mi    epc sepc.lv sepc.all sepc.nox
35     emot =~ item_10 72.671  1.844   0.760    0.595    0.595
118  item_7 ~~  item_8 61.401  0.371   0.371    0.437    0.437
38     emot =~ item_13 13.578 -0.708  -0.292   -0.253   -0.253
144 item_11 ~~ item_12 12.189 -0.209  -0.209   -0.179   -0.179
133  item_9 ~~ item_10 11.855  0.139   0.139    0.210    0.210
141 item_10 ~~ item_13 10.150 -0.185  -0.185   -0.255   -0.255

We could try to now put in the covariance between item_7 and item_8 too:

moddoom3 <- "
# emotional domain
emot =~ item_1 + item_2 + item_3 + item_4 + item_5 + item_6 + item_7 + item_8 + item_9
# practical domain
pract =~ item_10 + item_11 + item_12 + item_13 + item_14 + item_15
# correlated domains (will be estimated by default)
emot ~~ pract
# residual covariances
item_1 ~~ item_3
item_7 ~~ item_8
"

moddoom3.est <- cfa(moddoom3, data = doom)

fitmeasures(moddoom3.est)[c("rmsea","srmr","cfi","tli")]

     rmsea       srmr        cfi        tli 
0.03937033 0.04553018 0.95243659 0.94259588 

Whoop! It fits well! It may well be that if we inspect modification indices again, we still see that emot =~ item_10 would improve our model fit. The thing to remember however, is that we could simply keep adding parameters until we run out of degrees of freedom, and our model would “fit better”. But such a model would not be useful. It would not generalise well, because it runs the risk of being overfitted to the nuances of this specific sample.

modindices(moddoom3.est, sort=TRUE) |>
  head()

        lhs op     rhs     mi    epc sepc.lv sepc.all sepc.nox
36     emot =~ item_10 74.498  2.060   0.845    0.662    0.662
39     emot =~ item_13 13.839 -0.781  -0.320   -0.278   -0.278
144 item_11 ~~ item_12 11.589 -0.204  -0.204   -0.173   -0.173
141 item_10 ~~ item_13 10.728 -0.189  -0.189   -0.262   -0.262
143 item_10 ~~ item_15 10.215 -0.212  -0.212   -0.227   -0.227
133  item_9 ~~ item_10  9.674  0.126   0.126    0.199    0.199

Question 12

Based on our analysis of the DOOM measure, which of the following statements accurately reflect our current position? (you can choose multiple!)

• A The theoretical measurement model is now confirmed. Because our final fit indices (CFI/RMSEA) meet the thresholds, the initial theory of DOOM scrolling is validated.
• B We have confirmed our theoretical measurement model, but with some caveats
• C Our analysis has shifted from confirmatory to exploratory; while the model now fits the data, it requires validation on an independent sample.
• D The measure of doom scrolling is fundamentally flawed and should be discarded/updated due to the initial lack of fit.
• E The measure of doom scrolling is likely suffers from poor content validity. The need for correlated errors suggests item wordings are redundant or overlapping, meaning the scale likely needs redesigning.

Solution 17. I would probably opt for C and E here. The suggested factor loading of emot =~ item_10 to me is a bit of a red flag. Ideally we don’t want to have measures that have cross-loadings (items loading on to multiple factors), because it blurs the meaning of the factors.

Similarly, just looking at the wordings of the items - some of them (item_3) don’t actually ask anything specific to being online!

In this case, the likely reason for the initial poor fit of the “DOOM” scale, is that the person who made the items (ahem, me) doesn’t really know anything about the construct they are talking about, and didn’t put much care into constructing the items!

More Conduct Problems

Data: conduct_problems_2.csv

Last week we conducted an exploratory factor analysis of a dataset to try and identify an optimal factor structure for a new measure of conduct (i.e., antisocial behavioural) problems.

This week, we’ll conduct some confirmatory factor analyses (CFA) of the same inventory to assess the extent to which this 2-factor structure fits an independent sample. To do this, we have administered our measure to a new sample of n=600 adolescents.

We have re-ordered the questionnaire items to be grouped into the two types of behaviours:

Non-Aggressive Behaviours

item	behaviour
item 1	Stealing
item 2	Lying
item 3	Skipping school
item 4	Vandalism
item 5	Breaking curfew

Aggressive Behaviours

item	behaviour
item 6	Threatening others
item 7	Bullying
item 8	Spreading malicious rumours
item 9	Using a weapon
item 10	Fighting

The data are available as a .csv at https://uoepsy.github.io/data/conduct_problems_2.csv

Question 13

• Read in the data, and take a quick look around (e.g., cor matrix, quick pairs.panels plots etc).
• Fit the proposed 2 factor model
• Examine the fit of the 2-factor model of conduct problems to this new sample of 600 adolescents.
  
• Evaluate the fit, and make any model modifications if necessary (and only if you feel that there is substantive support for the modification given the items).
  
• Make a diagram of your model, using the standardised factor loadings as labels.
  
• Make a bullet point list of everything you have done so far, and the resulting conclusions. Then, if you feel like it, turn the bulleted list into written paragraphs, and you’ll have a write-up of your analyses!

Solution 18. Here’s the data:

library(tidyverse)
cp2 <- read_csv("https://uoepsy.github.io/data/conduct_problems_2.csv")

cor(cp2)

           item1     item2      item3     item4     item5      item6     item7
item1  1.0000000 0.5254249 0.43498606 0.4831121 0.5644047 0.13450795 0.2806369
item2  0.5254249 1.0000000 0.50099054 0.5119842 0.6740107 0.15239358 0.3031997
item3  0.4349861 0.5009905 1.00000000 0.4897339 0.5990709 0.09463966 0.2557407
item4  0.4831121 0.5119842 0.48973386 1.0000000 0.6199783 0.14847703 0.2492471
item5  0.5644047 0.6740107 0.59907089 0.6199783 1.0000000 0.11353641 0.3112458
item6  0.1345080 0.1523936 0.09463966 0.1484770 0.1135364 1.00000000 0.5461757
item7  0.2806369 0.3031997 0.25574075 0.2492471 0.3112458 0.54617567 1.0000000
item8  0.2576481 0.2748175 0.19355398 0.2553712 0.2463290 0.59441335 0.8009591
item9  0.2444337 0.2759875 0.20798523 0.2260156 0.2333758 0.36915878 0.5988582
item10 0.1783069 0.1618631 0.12924914 0.1433759 0.1394748 0.46316014 0.5316143
           item8     item9    item10
item1  0.2576481 0.2444337 0.1783069
item2  0.2748175 0.2759875 0.1618631
item3  0.1935540 0.2079852 0.1292491
item4  0.2553712 0.2260156 0.1433759
item5  0.2463290 0.2333758 0.1394748
item6  0.5944133 0.3691588 0.4631601
item7  0.8009591 0.5988582 0.5316143
item8  1.0000000 0.6169829 0.5557168
item9  0.6169829 1.0000000 0.3613253
item10 0.5557168 0.3613253 1.0000000

Just from the visual, it looks like the same factor structure is present in this sample.

heatmap(cor(cp2), scale = "none")

This is our proposed model:

library(lavaan)
cpmod <- "
  # the non-aggressive problems factor
  nonagg =~ item1 + item2 + item3 + item4 + item5

  # the aggressive problems factor
  agg =~ item6 + item7 + item8 + item9 + item10

  # covariance between the two factors
  # (this is included by default in cfa)
  agg ~~ nonagg
"

cpmod.est <- cfa(cpmod, data = cp2)

And it appears to fit pretty well!

fitmeasures(cpmod.est)[c("srmr","rmsea","cfi","tli")]

      srmr      rmsea        cfi        tli 
0.03454489 0.03945991 0.98862263 0.98494172 

We can check modification indices anyway, but I don’t plan on making any adjustments given that it already fits well:

modindices(cpmod.est, sort = TRUE) |> head()

      lhs op    rhs     mi    epc sepc.lv sepc.all sepc.nox
72  item6 ~~ item10 10.747  0.082   0.082    0.144    0.144
25 nonagg =~  item7  8.106  0.119   0.080    0.080    0.080
65  item5 ~~  item7  6.720  0.039   0.039    0.160    0.160
71  item6 ~~  item9  6.675 -0.065  -0.065   -0.115   -0.115
24 nonagg =~  item6  6.100 -0.136  -0.092   -0.094   -0.094
33    agg =~  item5  5.273 -0.122  -0.075   -0.072   -0.072

It maybe makes sense that there is some residual covariance between item6 (“threatening others”) and item10 (“fighting”), but it’s only a weak correlation (0.14). Not worth adding.

So let’s get on with making a diagram. We can rotate this however you like. Convention is typically to have it downwards but I like it left to right (not sure why!)

library(semPlot)
semPaths(cpmod.est, 
        whatLabels = "std", 
        rotation = 2)

The lines from agg =~ item6 and nonagg =~ item1 are dotted to indicate that the model was initially fitted with the loading fixed to 1.

Because we’re showing standardised loadings, we could just use the model when fitted with std.lv=TRUE just to stop these dotted lines from appearing:

cpmod.est2 <- cfa(cpmod, data = cp2, std.lv = TRUE)

semPaths(cpmod.est2, 
        whatLabels = "std", 
        rotation = 2)

And let’s give a brief write-up:

A two-factor model was tested. Items 1-5 loaded on a ‘non-aggressive conduct problems’ factor and items 6-10 loaded on an ‘aggression’ factor and these factors were allowed to correlate. Scaling and identification were achieved by fixing the loading of item 1 on the non-aggressive conduct problems factor and item 6 on the aggression factor to 1. The model was estimated using maximum likelihood estimation. The model fit well with CFI=.99, TLI=0.99, RMSEA=.04, and SRMR=.04 (Hu & Bentler, 1999). All loadings were statistically significant and >|.3| on the standardised scale. Overall, therefore, a two-factor oblique model was supported for the conduct problems items. The correlation between the factors was \(r=.38\,\, (p<.001)\).

parameter	est	std.est	se	z	pvalue
nonagg=~item1	1.000	0.664	0.000
nonagg=~item2	1.217	0.765	0.076	15.965	< 0.001
nonagg=~item3	1.012	0.676	0.070	14.424	< 0.001
nonagg=~item4	1.146	0.706	0.077	14.967	< 0.001
nonagg=~item5	1.360	0.872	0.078	17.412	< 0.001
agg=~item6	1.000	0.636	0.000
agg=~item7	1.419	0.878	0.082	17.331	< 0.001
agg=~item8	1.437	0.915	0.081	17.662	< 0.001
agg=~item9	1.093	0.668	0.077	14.125	< 0.001
agg=~item10	0.950	0.608	0.073	13.073	< 0.001
nonagg~~agg	0.156	0.375	0.023	6.856	< 0.001
item1~~item1	0.574	0.559	0.037	15.398	< 0.001
item2~~item2	0.477	0.416	0.035	13.778	< 0.001
item3~~item3	0.551	0.543	0.036	15.262	< 0.001
item4~~item4	0.597	0.501	0.040	14.867	< 0.001
item5~~item5	0.264	0.239	0.028	9.537	< 0.001
item6~~item6	0.561	0.595	0.035	16.193	< 0.001
item7~~item7	0.228	0.229	0.021	10.767	< 0.001
item8~~item8	0.153	0.162	0.019	8.122	< 0.001
item9~~item9	0.566	0.554	0.035	15.978	< 0.001
item10~~item10	0.587	0.630	0.036	16.351	< 0.001
nonagg~~nonagg	0.453	1.000	0.052	8.701	< 0.001
agg~~agg	0.381	1.000	0.045	8.396	< 0.001