AAQ-2 Psychometric Analysis with Rasch Measurement Theory

Code

# one package below requires that you use devtools to install it manually:
# first install devtools by
# install.packages('devtools')

library(RISEkbmRasch) # devtools::install_github("pgmj/RISEkbmRasch")
library(grateful)
library(ggrepel)
library(car)
library(kableExtra)
library(readxl)
library(tidyverse)
library(eRm)
library(mirt)
library(psych)
library(psychotree)
library(matrixStats)
library(reshape)
library(knitr)
library(patchwork)
library(formattable) 
library(glue)
library(haven)
library(labelled)

### optional libraries
#library(TAM)
#library(skimr)

### some commands exist in multiple packages, here we define preferred ones that are frequently used
select <- dplyr::select
count <- dplyr::count
recode <- car::recode
rename <- dplyr::rename

1 Background

Data from (Langer et al. 2024). More comments on that paper will be added later. The only other Item Response Theory-based paper found analyzing the AAQ-2 is (Ong et al. 2019) but they do a rather limited analysis using the Graded Response Model, for example omitting test of unidimensionality, local independence of items, and ordering of response categories.

Code

### import data
df.all <- read_spss("data/DataBase AAQ-II_Criteria_rev.sav")

df <- df.all %>% 
  select(starts_with("AAQ_II"),GenderBirth)

dif.sex <- df$GenderBirth
df$GenderBirth <- NULL
names(df) <- paste0("aaq_",c(1:7))

#val_labels(dif.sex)

itemlabels <- data.frame(
  stringsAsFactors = FALSE,
  itemnr = paste0("aaq_",c(1:7)),
  item = c("My painful experiences and memories make it difficult for me to live a life that I would value",
    "I’m afraid of my feelings",
    "I worry about not being able to control my worries and feelings",
    "My painful memories prevent me from having a fulfilling life",
    "Emotions cause problems in my life",
    "It seems like most people are handling their lives better than I am",
    "Worries get in the way of my success"
  )
)

source("RISE_theme.R")

2 All items in the analysis

Code

RIlistitems(df)

itemnr	item
aaq_1	My painful experiences and memories make it difficult for me to live a life that I would value
aaq_2	I’m afraid of my feelings
aaq_3	I worry about not being able to control my worries and feelings
aaq_4	My painful memories prevent me from having a fulfilling life
aaq_5	Emotions cause problems in my life
aaq_6	It seems like most people are handling their lives better than I am
aaq_7	Worries get in the way of my success

3 Demographics

Code

RIdemographics(dif.sex, "Gender")

Gender	n	Percent
1	571	33.5
2	1134	66.5

Code

hist(df.all$Age, title = "Age distribution", xlab = "Age", main = "Age distribution", col = "lightblue")

Code

#abline(v = mean(df.all$Age), col = "red", lwd = 2)

Very young sample, and not a lot of variation in age either.

4 Descriptives of raw data

Response distribution for all items together.

Code

RIallresp(df)

Response category	Number of responses	Percent
1	2541	21.3
2	2263	19.0
3	1443	12.1
4	2241	18.8
5	1525	12.8
6	1046	8.8
7	857	7.2
NA	19	0.2

We need the response scale to start at 0 instead of 1 for the Rasch analysis to work correctly.

Code

df <- df %>% 
  mutate(across(starts_with("aaq_"), ~ .x - 1))

Code

allresp <- data.frame(
    Response.category = c(0L, 1L, 2L, 3L, 4L, 5L, 6L, NA),
  Number.of.responses = c(2541L, 2263L, 1443L, 2241L, 1525L, 1046L, 857L, 19L),
              Percent = c(21.3, 19, 12.1, 18.8, 12.8, 8.8, 7.2, 0.2)
)

ggplot(allresp, aes(x = Response.category, y = Percent)) +
  geom_bar(stat = "identity", fill = "#009ca6") +
  geom_text(aes(label = Number.of.responses), vjust = -0.5, size = 3) +
  labs(title = "Response distribution for all items",
       x = "Response category",
       y = "% of responses") +
  scale_x_continuous(breaks = c(0:6)) +
  theme_rise()

4.1 Descriptives - item level

itemnr	item
aaq_1	My painful experiences and memories make it difficult for me to live a life that I would value
aaq_2	I’m afraid of my feelings
aaq_3	I worry about not being able to control my worries and feelings
aaq_4	My painful memories prevent me from having a fulfilling life
aaq_5	Emotions cause problems in my life
aaq_6	It seems like most people are handling their lives better than I am
aaq_7	Worries get in the way of my success

Code

RItileplot(df)

Code

RIbarstack(df)

Code

RIbarplot(df)

Data is skewed towards the lower end of the scale. Most item responses are a bit oddly distributed. Response category 2 is consistently deviating from the expected (normal-ish) pattern. Makes one wonder about the data collection and the response category wording used. Let us check.

Code

val_labels(df.all$AAQ_II_1)

         Nunca es verdad para mi  Muy raramente es verdad para mi 
                               1                                2 
     Raramente es verdad para mi        A veces es verdad para mi 
                               3                                4 
Frecuentemente es verdad para mi   Casi siempre es verdad para mi 
                               5                                6 
       Siempre es verdad para mi 
                               7

My Spanish is not great, so here is the Google Translated version:

It’s never true for me
It’s very rarely true for me.
It’s rarely true for me
Sometimes it’s true for me
It’s often true for me.
It’s almost always true for me.
It’s always true for me

4.2 Missing data

Code

RImissing(df)

Code

RImissingP(df, n = 20)

One respondent with 100% missing, and 12 with 1 item missing. We’ll remove all respondents with missing data since we have a large dataset.

Code

df <- na.omit(df)
dif.sex <- df.all %>% 
  select(GenderBirth,starts_with("AAQ")) %>%
  na.omit() %>% 
  pull(GenderBirth)

dif.sex <- factor(dif.sex, levels = c(1,2),
                  labels = c("Male","Female"))

5 Rasch analysis 1

The eRm package, which uses Conditional Maximum Likelihood (CML) estimation, will be used primarily. For this analysis, the Partial Credit Model will be used.

itemnr	item
aaq_1	My painful experiences and memories make it difficult for me to live a life that I would value
aaq_2	I’m afraid of my feelings
aaq_3	I worry about not being able to control my worries and feelings
aaq_4	My painful memories prevent me from having a fulfilling life
aaq_5	Emotions cause problems in my life
aaq_6	It seems like most people are handling their lives better than I am
aaq_7	Worries get in the way of my success

Code

RIitemfitPCM2(df, 250, 16
              )

	OutfitMSQ	InfitMSQ	OutfitZSTD	InfitZSTD
aaq_1	0.934	0.888	-0.645	-1.417
aaq_2	0.966	0.948	-0.51	-0.154
aaq_3	0.867	0.862	-1.526	-1.031
aaq_4	0.864	0.847	-1.57	-1.676
aaq_5	0.679	0.672	-3.877	-4.095
aaq_6	1.238	1.133	2.687	0.681
aaq_7	0.794	0.803	-2.377	-2.274

Code

RIpcmPCA(df)

PCA of Rasch model residuals
Eigenvalues
1.88
1.57
1.21
0.93
0.79

Code

RIresidcorr(df, cutoff = 0.2)

	aaq_1	aaq_2	aaq_3	aaq_4	aaq_5	aaq_6
aaq_1
aaq_2	-0.17
aaq_3	-0.31	0.14
aaq_4	0.44	-0.23	-0.3
aaq_5	-0.28	-0.18	-0.11	-0.14
aaq_6	-0.27	-0.23	-0.26	-0.34	-0.15
aaq_7	-0.32	-0.31	-0.1	-0.31	0.09	0.04
Note:
Relative cut-off value (highlighted in red) is 0.042, which is 0.2 above the average correlation.

Code

RIloadLoc(df)

Code

RIitemCats(df)

Code

# increase fig-height above as needed, if you have many items
RItargeting(df)

Code

RIitemHierarchy(df)

Code

erm.df <- PCM(df)
LRtest(erm.df, dif.sex)


Andersen LR-test: 
LR-value: 62.299 
Chi-square df: 41 
p-value:  0.018

Code

RIdifTableLR(df, dif.sex)

	Item locations				Standard errors
Item	Male	Female	MaxDiff	All	SE_Male	SE_Female	SE_All
aaq_1	0.675	0.741	0.066	0.721	0.223	0.130	0.111
aaq_2	0.484	0.507	0.023	0.499	0.195	0.125	0.104
aaq_3	0.13	-0.046	0.176	0.023	0.178	0.121	0.098
aaq_4	0.847	0.877	0.03	0.868	0.226	0.135	0.115
aaq_5	0.482	0.426	0.056	0.434	0.199	0.122	0.102
aaq_6	0.199	0.209	0.01	0.197	0.177	0.121	0.099
aaq_7	0.276	0.318	0.042	0.306	0.184	0.122	0.101
Note:
Values highlighted in red are above the chosen cutoff 0.5 logits. Background color brown and blue indicate the lowest and highest values among the DIF groups.

Item 5 has somewhat low item fit. PCA of residuals is below 2, but residual correlations show several issues. The strongest correlation is between items 1 and 4, followed by 2 and 3. Items 5 and 7 are also above the cutoff.

The tests for DIF indicates possible issues (p < .05), but the DIF table shows that the differences are small, at least on the item average level. The p-value is small due to the large sample.

All items show disordered response categories. The pattern is most clearly seen in the targeting figure. The consistency observed, with threshold 3 below threhold 2 for all items, makes one wonder if there was a mistake in the coding of responses to numbers, but we don’t know if the supplied datafile is the original.

We’ll merge categories 1 and 2:

It’s very rarely true for me.
It’s rarely true for me

Code

df2 <- df %>% 
  mutate(across(starts_with("aaq_"), ~ recode(., "2=1;3=2;4=3;5=4;6=5")))

5.1 1 and 2 merged

Code

mirt(df2, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Let’s also try merging 2 and 3 for comparison.

Code

df2 <- df %>% 
  mutate(across(starts_with("aaq_"), ~ recode(., "3=2;4=3;5=4;6=5")))

5.2 2 and 3 merged instead

Code

mirt(df2, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

It seems 1 and 2 works better, as expected.

Code

df2 <- df %>% 
  mutate(across(starts_with("aaq_"), ~ recode(., "2=1;3=2;4=3;5=4;6=5")))

5.3 Residual correlations

We have several item pairs that correlate above the relative cutoff of 0.2 and will deal with them one at a time. The strongest correlation is between items 1 and 4.

item 1: My painful experiences and memories make it difficult for me to live a life that I would value
item 4: My painful memories prevent me from having a fulfilling life

These items are very much alike, apart from “fulfilling life” vs “a life that I would value” and item 1 adding “experiences” to “memories”, but even these differences are very similar. It’s not surprising that they correlate strongly.

Item 1 has better separation of item response thresholds. We’ll remove item 4.

Code

df2 <- df2 %>% 
  select(!aaq_4)

6 Rasch analysis 2

itemnr	item
aaq_1	My painful experiences and memories make it difficult for me to live a life that I would value
aaq_2	I’m afraid of my feelings
aaq_3	I worry about not being able to control my worries and feelings
aaq_4	My painful memories prevent me from having a fulfilling life
aaq_5	Emotions cause problems in my life
aaq_6	It seems like most people are handling their lives better than I am
aaq_7	Worries get in the way of my success

Code

RIitemfitPCM2(df2, 250, 16)

	OutfitMSQ	InfitMSQ	OutfitZSTD	InfitZSTD
aaq_1	1.068	1.015	0.504	0.28
aaq_2	0.874	0.883	-1.013	-1.258
aaq_3	0.805	0.785	-2.119	-2.567
aaq_5	0.654	0.655	-4.019	-4.227
aaq_6	1.122	1.067	0.904	0.957
aaq_7	0.729	0.743	-2.828	-3.173

Code

RIpcmPCA(df2)

PCA of Rasch model residuals
Eigenvalues
1.66
1.31
1.21
0.97
0.85

Code

RIresidcorr(df2, cutoff = 0.2)

	aaq_1	aaq_2	aaq_3	aaq_5	aaq_6
aaq_1
aaq_2	-0.12
aaq_3	-0.26	0.1
aaq_5	-0.17	-0.21	-0.16
aaq_6	-0.22	-0.31	-0.36	-0.22
aaq_7	-0.23	-0.37	-0.19	0.03	-0.07
Note:
Relative cut-off value (highlighted in red) is 0.015, which is 0.2 above the average correlation.

Code

RIloadLoc(df2)

Code

# increase fig-height above as needed, if you have many items
RItargeting(df2)

Code

RIitemHierarchy(df2)

Items 2 and 3 still correlate quite strongly.

item 2: I’m afraid of my feelings
item 3: I worry about not being able to control my worries and feelings

Item 3 has worse fit, and we’ll remove it. It is also a dual question, “worries and feelings”, which can explain the item fit.

Code

df2 <- df2 %>% 
  select(!aaq_3)

7 Rasch analysis 3

itemnr	item
aaq_1	My painful experiences and memories make it difficult for me to live a life that I would value
aaq_2	I’m afraid of my feelings
aaq_3	I worry about not being able to control my worries and feelings
aaq_4	My painful memories prevent me from having a fulfilling life
aaq_5	Emotions cause problems in my life
aaq_6	It seems like most people are handling their lives better than I am
aaq_7	Worries get in the way of my success

Code

RIitemfitPCM2(df2, 250, 16)

	OutfitMSQ	InfitMSQ	OutfitZSTD	InfitZSTD
aaq_1	0.979	0.936	0.025	-0.694
aaq_2	0.953	0.933	-0.733	-0.716
aaq_5	0.629	0.624	-4.639	-4.585
aaq_6	0.962	0.946	-0.617	-0.567
aaq_7	0.683	0.7	-3.399	-3.553

Code

RIpcmPCA(df2)

PCA of Rasch model residuals
Eigenvalues
1.52
1.35
1.16
0.98
0.01

Code

RIresidcorr(df2, cutoff = 0.2)

	aaq_1	aaq_2	aaq_5	aaq_6
aaq_1
aaq_2	-0.12
aaq_5	-0.22	-0.17
aaq_6	-0.32	-0.33	-0.3
aaq_7	-0.3	-0.34	0	-0.15
Note:
Relative cut-off value (highlighted in red) is -0.026, which is 0.2 above the average correlation.

Code

RIloadLoc(df2)

Code

# increase fig-height above as needed, if you have many items
RItargeting(df2)

Code

RIitemHierarchy(df2)

There is a large gap in targeting due to the dysfunctional response categories that needed merging. Item 5 is low in item fit, “Emotions cause problems in my life”. It is a very general item, which may explain the low fit. We’ll keep it for now.

8 DIF-analysis

8.1 Gender

itemnr	item
aaq_1	My painful experiences and memories make it difficult for me to live a life that I would value
aaq_2	I’m afraid of my feelings
aaq_5	Emotions cause problems in my life
aaq_6	It seems like most people are handling their lives better than I am
aaq_7	Worries get in the way of my success

Code

RIdifTableLR(df2, dif.sex)

	Item locations				Standard errors
Item	Male	Female	MaxDiff	All	SE_Male	SE_Female	SE_All
aaq_1	0.836	0.877	0.041	0.871	0.234	0.131	0.112
aaq_2	0.637	0.628	0.009	0.635	0.201	0.123	0.104
aaq_5	0.639	0.529	0.11	0.558	0.206	0.122	0.103
aaq_6	0.301	0.271	0.03	0.276	0.177	0.116	0.097
aaq_7	0.393	0.399	0.006	0.404	0.188	0.120	0.100
Note:
Values highlighted in red are above the chosen cutoff 0.5 logits. Background color brown and blue indicate the lowest and highest values among the DIF groups.

Code

RIdifFigureLR(df2, dif.sex)

Code

RIdifThreshTblLR(df2, dif.sex)

	Threshold locations				Standard errors
Item threshold	Male	Female	MaxDiff	All	SE_Male	SE_Female	SE_All
aaq_1
c1	-1.606	-1.394	0.212	-1.469	0.129	0.097	0.078
c2	0.693	0.656	0.037	0.673	0.129	0.091	0.074
c3	1.207	1.265	0.058	1.255	0.181	0.118	0.099
c4	2.296	1.555	0.741	1.709	0.313	0.148	0.132
c5	1.592	2.301	0.709	2.185	0.417	0.199	0.179
aaq_2
c1	-1.873	-1.678	0.195	-1.749	0.140	0.110	0.086
c2	0.373	0.138	0.235	0.223	0.128	0.092	0.074
c3	0.984	0.884	0.1	0.912	0.170	0.105	0.089
c4	1.289	1.466	0.177	1.428	0.221	0.128	0.111
c5	2.412	2.329	0.083	2.361	0.344	0.181	0.160
aaq_5
c1	-1.77	-1.763	0.007	-1.762	0.138	0.108	0.085
c2	0.338	0.41	0.072	0.389	0.126	0.093	0.075
c3	1.053	0.923	0.13	0.969	0.168	0.112	0.093
c4	1.754	1.377	0.377	1.469	0.250	0.135	0.118
c5	1.82	1.696	0.124	1.724	0.347	0.161	0.146
aaq_6
c1	-1.829	-1.658	0.171	-1.719	0.145	0.112	0.088
c2	0.239	0.18	0.059	0.203	0.138	0.100	0.081
c3	0.421	0.617	0.196	0.558	0.164	0.113	0.093
c4	1.014	0.948	0.066	0.977	0.191	0.125	0.104
c5	1.659	1.266	0.393	1.360	0.248	0.133	0.116
aaq_7
c1	-2.002	-1.979	0.023	-1.983	0.146	0.116	0.090
c2	0.332	0.129	0.203	0.198	0.132	0.092	0.075
c3	0.597	0.933	0.336	0.834	0.162	0.108	0.090
c4	1.491	1.246	0.245	1.320	0.216	0.130	0.111
c5	1.548	1.665	0.117	1.649	0.283	0.153	0.135
Note:
Values highlighted in red are above the chosen cutoff 0.5 logits. Background color brown and blue indicate the lowest and highest values among the DIF groups.

Code

RIdifThreshFigLR(df2, dif.sex)

Item 1 shows a surprising pattern for the two highest threholds.

9 Other metrics

9.1 Reliability

Code

RItif(df2)

9.2 Person location & infit ZSTD

Code

RIpfit(df)

9.3 Item parameters

Code

RIitemparams(df2)

	Threshold 1	Threshold 2	Threshold 3	Threshold 4	Threshold 5	Item location
aaq_1	-1.47	0.67	1.25	1.71	2.18	0.87
aaq_2	-1.75	0.22	0.91	1.43	2.36	0.64
aaq_5	-1.76	0.39	0.97	1.47	1.72	0.56
aaq_6	-1.72	0.20	0.56	0.98	1.36	0.28
aaq_7	-1.98	0.20	0.83	1.32	1.65	0.4
Note:
Item location is the average of the thresholds for each item.

9.4 Transformation table

Code

RIscoreSE(df2)

Ordinal sum score	Logit score	Logit std.error
0	-4.000	1.487
1	-2.862	0.963
2	-2.172	0.793
3	-1.639	0.705
4	-1.182	0.643
5	-0.789	0.587
6	-0.460	0.537
7	-0.189	0.495
8	0.038	0.461
9	0.234	0.435
10	0.407	0.415
11	0.566	0.401
12	0.714	0.391
13	0.855	0.384
14	0.992	0.381
15	1.128	0.382
16	1.265	0.385
17	1.407	0.392
18	1.555	0.404
19	1.714	0.422
20	1.890	0.446
21	2.093	0.482
22	2.338	0.537
23	2.661	0.627
24	3.148	0.803
25	4.000	1.245

9.5 Ordinal/interval figure

Code

RIscoreSE(df2, output = "figure")

10 Summary of Rasch analysis

Two item pairs had strongly correlated residuals, and all items had issues with disordered thresholds related to the second lowest response category. There was no DIF for gender at birth.

Targeting is not great. Items have rather similar locations, and there is a large gap where the average person location is. If the intended use is clinical, there is decent reliability for those with above average locations.

11 Comparison of latent scores

The “standard” way to use the AAQ-2 is to sum/average the items, after recoding item categories to numerics. We will compare this to the Rasch estimated latent scores.

For ordinal sum scores, we will use the original 7 items with their original 7 response categories, even though the analysis does not support this. We do this to illustrate the difference that can be hidden behind data that is unjustifiedly sum scored and data based on psychometric analysis.

Code

df.viz <- df %>% 
  na.omit() %>% 
  mutate(sumscore = rowSums(across(starts_with("aaq_")))) %>% 
  mutate(sumscore = as.numeric(sumscore)) %>% 
  rownames_to_column("id") %>% 
  mutate(id = as.numeric(id))

df.viz$latentscore <- RIestThetas2(df2)

Code

ggplot(df.viz, aes(x = latentscore, y = sumscore)) +
  geom_point() +
  labs(x = "Rasch estimated latent score", y = "Ordinal sum score") +
  theme_rise()

Code

pfit <- RIpfit(df2, output = "dataframe") %>% 
  janitor::clean_names() %>% 
  rename(id = rownumber)

df.viz <- left_join(df.viz, pfit, by = "id")

df.viz <- df.viz %>% 
  mutate(person_fit = ifelse(person_infit_zstd > 2 | person_infit_zstd < -2, "outlier", "normal"))

ggplot(df.viz, aes(x = latentscore, y = sumscore, color = person_fit)) +
  geom_point() +
  labs(x = "Rasch estimated latent score", y = "Ordinal sum score") +
  theme_rise() +
  scale_color_viridis_d() +
  geom_smooth(method = "loess", color = "darkgrey")

Code

hist(df.viz$latentscore, breaks = 20, main = "Histogram of Rasch estimated (WL) latent scores", col = "lightpink")

Code

hist(df.viz$sumscore, breaks = 20, main = "Histogram of ordinal sum scores", col = "lightblue")
mtext(text = "Note: this is based on the original 7 items and 7 response categories.", 
      side = 3, line = 0.4)

12 CFA

Let’s make a brief comparison.

Code

library(lavaan)
library(lavaanExtra)

Here is some info about the CFA methodology: https://pgmj.github.io/ki_irt_mhcsf/cfa.html It will be copied to this document later.

We will use both WLSMV and ML/MLR estimators.

Code

# Define latent variables
latent <- list(
  aaq2 = names(df)
)

# Write the model
cfa.model <- write_lavaan(latent = latent)

fit.wlsmv <- cfa(model = cfa.model, 
               data = df,
               estimator = "WLSMV",
               ordered = TRUE,
               rotation = "oblimin")

fit.ml <- cfa(model = cfa.model, 
               data = df,
               estimator = "ML",
               rotation = "oblimin")

fit.mlr <- cfa(model = cfa.model, 
               data = df,
               estimator = "MLR",
               rotation = "oblimin")

Code

fit_metrics_robust <- c("chisq.scaled", "df", "pvalue.scaled", 
                         "cfi.robust", "tli.robust", "rmsea.robust", 
                        "rmsea.ci.lower.robust","rmsea.ci.upper.robust",
                        "srmr")

mfit1 <- 
  fitmeasures(fit.wlsmv, fit_metrics_robust) %>% 
  rbind() %>% 
  as.data.frame() %>% 
  mutate(across(where(is.numeric),~ round(.x, 3))) %>%
  rename(Chi2 = chisq.scaled,
         p = pvalue.scaled,
         CFI = cfi.robust,
         TLI = tli.robust,
         RMSEA = rmsea.robust,
         CI_low = rmsea.ci.lower.robust,
         CI_high = rmsea.ci.upper.robust,
         SRMR = srmr) %>% 
  add_column(Model = "WLSMV", .before = "Chi2")

Code

fit_metrics.ml <- c("chisq", "df", "pvalue", 
                  "cfi", "tli", "rmsea", 
                  "rmsea.ci.lower","rmsea.ci.upper",
                  "srmr")

mfit2 <- 
  fitmeasures(fit.ml, fit_metrics.ml) %>% 
  rbind() %>% 
  as.data.frame() %>% 
  mutate(across(where(is.numeric),~ round(.x, 3))) %>% 
    rename(Chi2 = chisq,
         p = pvalue,
         CFI = cfi,
         TLI = tli,
         RMSEA = rmsea,
         CI_low = rmsea.ci.lower,
         CI_high = rmsea.ci.upper,
         SRMR = srmr) %>% 
  add_column(Model = "ML", .before = "Chi2")

Code

fit_metrics.mlr <- c("chisq.scaled", "df.scaled", "pvalue", 
                  "cfi.scaled", "tli.scaled", "rmsea.scaled", 
                  "rmsea.ci.lower.scaled","rmsea.ci.upper.scaled",
                  "srmr")

mfit3 <- 
  fitmeasures(fit.mlr, fit_metrics.mlr) %>% 
  rbind() %>% 
  as.data.frame() %>% 
  mutate(across(where(is.numeric),~ round(.x, 3))) %>% 
    rename(Chi2 = chisq.scaled,
         p = pvalue,
         df = df.scaled,
         CFI = cfi.scaled,
         TLI = tli.scaled,
         RMSEA = rmsea.scaled,
         CI_low = rmsea.ci.lower.scaled,
         CI_high = rmsea.ci.upper.scaled,
         SRMR = srmr) %>% 
  add_column(Model = "MLR", .before = "Chi2")

12.1 Model fit comparison

Code

mfit1 %>% 
  bind_rows(mfit2) %>% 
  bind_rows(mfit3) %>% 
  rownames_to_column("remove") %>%
  select(!remove) %>%
  kbl_rise()

Model	Chi2	df	CFI	TLI	RMSEA	CI_low	CI_high	SRMR
WLSMV	1349.580	14	0.853	0.780	0.251	0.238	0.264	0.065
ML	1033.195	14	0.882	0.822	0.207	0.197	0.218	0.055
MLR	739.565	14	0.876	0.814	0.175	0.166	0.184	0.055

12.2 Modification indices

We’ll look at the MLR output, since that model had the least bad fit.

Code

modificationIndices(fit.mlr,
                    standardized = T) %>% 
  as.data.frame(row.names = NULL) %>% 
  filter(mi > 30,
         op == "~~") %>% 
  arrange(desc(mi)) %>% 
  mutate(across(where(is.numeric),~ round(.x, 3))) %>%
  kbl_rise(fontsize = 14, tbl_width = 75)

lhs	op	rhs	mi	epc	sepc.lv	sepc.all	sepc.nox
aaq_1	~~	aaq_4	694.861	0.948	0.948	0.754	0.754
aaq_2	~~	aaq_3	180.075	0.501	0.501	0.390	0.390
aaq_6	~~	aaq_7	95.556	0.403	0.403	0.284	0.284
aaq_5	~~	aaq_7	71.931	0.276	0.276	0.290	0.290
aaq_4	~~	aaq_7	58.763	-0.273	-0.273	-0.230	-0.230
aaq_1	~~	aaq_7	56.133	-0.259	-0.259	-0.223	-0.223
aaq_3	~~	aaq_4	54.749	-0.277	-0.277	-0.217	-0.217
aaq_1	~~	aaq_5	46.764	-0.221	-0.221	-0.218	-0.218
aaq_2	~~	aaq_7	44.276	-0.236	-0.236	-0.198	-0.198
aaq_1	~~	aaq_3	43.648	-0.241	-0.241	-0.192	-0.192
aaq_4	~~	aaq_6	42.108	-0.279	-0.279	-0.182	-0.182

Huge residual correlations between items 1 and 4, and 2 and 3. Amusingly, GitHub Copilot now suggests to me that we add these to the model, which is a really bad practice that too often is used. Unidimensionality is a key assumption of this CFA, and we should not add items to the model just because they have high residual correlations.

12.3 CFA based on Rasch model

Code

# Define latent variables
latent <- list(
  aaq2 = names(df2)
)

# Write the model
cfa.model <- write_lavaan(latent = latent)

fit.wlsmv <- cfa(model = cfa.model, 
               data = df2,
               estimator = "WLSMV",
               ordered = TRUE,
               rotation = "oblimin")

fit_metrics_robust <- c("chisq.scaled", "df", "pvalue.scaled", 
                         "cfi.robust", "tli.robust", "rmsea.robust", 
                        "rmsea.ci.lower.robust","rmsea.ci.upper.robust",
                        "srmr")

mfit1 <- 
  fitmeasures(fit.wlsmv, fit_metrics_robust) %>% 
  rbind() %>% 
  as.data.frame() %>% 
  mutate(across(where(is.numeric),~ round(.x, 3))) %>%
  rename(Chi2 = chisq.scaled,
         p = pvalue.scaled,
         CFI = cfi.robust,
         TLI = tli.robust,
         RMSEA = rmsea.robust,
         CI_low = rmsea.ci.lower.robust,
         CI_high = rmsea.ci.upper.robust,
         SRMR = srmr) %>% 
  add_column(Model = "WLSMV", .before = "Chi2")

kbl_rise(mfit1)

	Model	Chi2	df	p	CFI	TLI	RMSEA	CI_low	CI_high	SRMR
.	WLSMV	106.048	5	0	0.982	0.965	0.106	0.089	0.125	0.022

Code

modificationIndices(fit.wlsmv,
                    standardized = T) %>% 
  as.data.frame(row.names = NULL) %>% 
  #filter(mi > 30,
  #       op == "~~") %>% 
  arrange(desc(mi)) %>% 
  mutate(across(where(is.numeric),~ round(.x, 3))) %>%
  kbl_rise(fontsize = 14, tbl_width = 75)

lhs	op	rhs	mi	epc	sepc.lv	sepc.all	sepc.nox
aaq_1	~~	aaq_2	21.870	-0.090	-0.090	-0.213	-0.213
aaq_2	~~	aaq_7	15.519	0.080	0.080	0.258	0.258
aaq_6	~~	aaq_7	13.347	-0.069	-0.069	-0.229	-0.229
aaq_1	~~	aaq_7	7.687	0.057	0.057	0.177	0.177
aaq_5	~~	aaq_6	5.321	0.048	0.048	0.172	0.172
aaq_5	~~	aaq_7	4.966	-0.044	-0.044	-0.205	-0.205
aaq_1	~~	aaq_6	1.245	0.023	0.023	0.055	0.055
aaq_2	~~	aaq_6	0.490	0.014	0.014	0.035	0.035
aaq_1	~~	aaq_5	0.347	0.012	0.012	0.040	0.040
aaq_2	~~	aaq_5	0.097	-0.006	-0.006	-0.022	-0.022

13 Software used

Code

pkgs <- cite_packages(cite.tidyverse = TRUE,
                      output = "table",
                      bib.file = "grateful-refs.bib",
                      include.RStudio = TRUE,
                      out.dir = getwd())
formattable(pkgs, 
            table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 80%"')

Package	Version	Citation
base	4.3.3	R Core Team (2024)
car	3.1.2	Fox and Weisberg (2019)
eRm	1.0.6	Mair and Hatzinger (2007b); Mair and Hatzinger (2007a); Hatzinger and Rusch (2009); Rusch, Maier, and Hatzinger (2013); Koller, Maier, and Hatzinger (2015); Debelak and Koller (2019)
formattable	0.2.1	Ren and Russell (2021)
ggrepel	0.9.5	Slowikowski (2024)
glue	1.7.0	Hester and Bryan (2024)
janitor	2.2.0	Firke (2023)
kableExtra	1.4.0	Zhu (2024)
knitr	1.45	Xie (2014); Xie (2015); Xie (2023)
labelled	2.12.0	Larmarange (2023)
lavaan	0.6.17	Rosseel (2012)
lavaanExtra	0.2.0	Thériault (2023)
matrixStats	1.2.0	Bengtsson (2023)
mirt	1.41	Chalmers (2012)
patchwork	1.2.0	Pedersen (2024)
psych	2.4.3	William Revelle (2024)
psychotree	0.16.0	Trepte and Verbeet (2010); Strobl, Wickelmaier, and Zeileis (2011); Strobl, Kopf, and Zeileis (2015); Komboz, Zeileis, and Strobl (2018); Wickelmaier and Zeileis (2018)
reshape	0.8.9	Wickham (2007)
RISEkbmRasch	0.1.33.3	Johansson (2024)
rmarkdown	2.26	Xie, Allaire, and Grolemund (2018); Xie, Dervieux, and Riederer (2020); Allaire et al. (2024)
tidyverse	2.0.0	Wickham et al. (2019)

Code

sessionInfo()

R version 4.3.3 (2024-02-29)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.4.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Stockholm
tzcode source: internal

attached base packages:
 [1] parallel  grid      stats4    stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] furrr_0.3.1           future_1.33.2         lavaanExtra_0.2.0    
 [4] lavaan_0.6-17         doParallel_1.0.17     iterators_1.0.14     
 [7] foreach_1.5.2         labelled_2.12.0       haven_2.5.4          
[10] knitr_1.45            readxl_1.4.3          grateful_0.2.4       
[13] RISEkbmRasch_0.1.33.3 hexbin_1.28.3         catR_3.17            
[16] glue_1.7.0            ggrepel_0.9.5         patchwork_1.2.0      
[19] reshape_0.8.9         matrixStats_1.2.0     psychotree_0.16-0    
[22] psychotools_0.7-3     partykit_1.2-20       mvtnorm_1.2-4        
[25] libcoin_1.0-10        psych_2.4.3           mirt_1.41            
[28] lattice_0.22-5        eRm_1.0-6             lubridate_1.9.3      
[31] forcats_1.0.0         stringr_1.5.1         dplyr_1.1.4          
[34] purrr_1.0.2           readr_2.1.5           tidyr_1.3.1          
[37] tibble_3.2.1          ggplot2_3.5.0         tidyverse_2.0.0      
[40] kableExtra_1.4.0      formattable_0.2.1     car_3.1-2            
[43] carData_3.0-5        

loaded via a namespace (and not attached):
 [1] mnormt_2.1.1          pbapply_1.7-2         gridExtra_2.3        
 [4] permute_0.9-7         rlang_1.1.3           magrittr_2.0.3       
 [7] snakecase_0.11.1      compiler_4.3.3        mgcv_1.9-1           
[10] systemfonts_1.0.5     vctrs_0.6.5           quadprog_1.5-8       
[13] pkgconfig_2.0.3       fastmap_1.1.1         labeling_0.4.3       
[16] inum_1.0-5            pbivnorm_0.6.0        utf8_1.2.4           
[19] rmarkdown_2.26        tzdb_0.4.0            xfun_0.43            
[22] jsonlite_1.8.8        Deriv_4.1.3           cluster_2.1.6        
[25] R6_2.5.1              stringi_1.8.3         parallelly_1.36.0    
[28] rpart_4.1.23          cellranger_1.1.0      Rcpp_1.0.12          
[31] Matrix_1.6-5          splines_4.3.3         timechange_0.2.0     
[34] tidyselect_1.2.0      rstudioapi_0.16.0     abind_1.4-5          
[37] yaml_2.3.8            vegan_2.6-4           codetools_0.2-19     
[40] dcurver_0.9.2         listenv_0.9.0         plyr_1.8.9           
[43] withr_3.0.0           evaluate_0.23         survival_3.5-8       
[46] xml2_1.3.6            pillar_1.9.0          renv_1.0.3           
[49] generics_0.1.3        hms_1.1.3             munsell_0.5.0        
[52] scales_1.3.0          globals_0.16.2        janitor_2.2.0        
[55] tools_4.3.3           colorspace_2.1-0      nlme_3.1-164         
[58] Formula_1.2-5         cli_3.6.2             fansi_1.0.6          
[61] viridisLite_0.4.2     svglite_2.1.3         gtable_0.3.4         
[64] digest_0.6.34         GPArotation_2023.11-1 farver_2.1.1         
[67] htmlwidgets_1.6.4     htmltools_0.5.8.1     lifecycle_1.0.4      
[70] MASS_7.3-60.0.1

14 References

Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2024. rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.

Bengtsson, Henrik. 2023. matrixStats: Functions That Apply to Rows and Columns of Matrices (and to Vectors). https://CRAN.R-project.org/package=matrixStats.

Chalmers, R. Philip. 2012. “mirt: A Multidimensional Item Response Theory Package for the R Environment.” Journal of Statistical Software 48 (6): 1–29. https://doi.org/10.18637/jss.v048.i06.

Debelak, Rudolf, and Ingrid Koller. 2019. “Testing the Local Independence Assumption of the Rasch Model with Q3-Based Nonparametric Model Tests.” Applied Psychological Measurement 44. https://doi.org/10.1177/0146621619835501.

Firke, Sam. 2023. janitor: Simple Tools for Examining and Cleaning Dirty Data. https://CRAN.R-project.org/package=janitor.

Fox, John, and Sanford Weisberg. 2019. An R Companion to Applied Regression. Third. Thousand Oaks CA: Sage. https://socialsciences.mcmaster.ca/jfox/Books/Companion/.

Hatzinger, Reinhold, and Thomas Rusch. 2009. “IRT Models with Relaxed Assumptions in eRm: A Manual-Like Instruction.” Psychology Science Quarterly 51.

Hester, Jim, and Jennifer Bryan. 2024. glue: Interpreted String Literals. https://CRAN.R-project.org/package=glue.

Johansson, Magnus. 2024. RISEkbmRasch: Psychometric Analysis in r with Rasch Measurement Theory. https://github.com/pgmj/RISEkbmRasch.

Koller, Ingrid, Marco Maier, and Reinhold Hatzinger. 2015. “An Empirical Power Analysis of Quasi-Exact Tests for the Rasch Model: Measurement Invariance in Small Samples.” Methodology 11. https://doi.org/10.1027/1614-2241/a000090.

Komboz, Basil, Achim Zeileis, and Carolin Strobl. 2018. “Tree-Based Global Model Tests for Polytomous Rasch Models.” Educational and Psychological Measurement 78 (1): 128–66. https://doi.org/10.1177/0013164416664394.

Langer, Álvaro I., Fernando P. Ponce, Jorge L. Ordóñez-Carrasco, Reiner Fuentes-Ferrada, Scarlett Mac-Ginty, Jorge Gaete, and Daniel Núñez. 2024. “Psychometric Evidence of the Acceptance and Action Questionnaire-II (AAQ-II): An Item Response Theory Analysis in University Students from Chile.” BMC Psychology 12 (1): 111. https://doi.org/10.1186/s40359-024-01608-w.

Larmarange, Joseph. 2023. labelled: Manipulating Labelled Data. https://CRAN.R-project.org/package=labelled.

Mair, Patrick, and Reinhold Hatzinger. 2007a. “CML Based Estimation of Extended Rasch Models with the eRm Package in r.” Psychology Science 49. https://doi.org/10.18637/jss.v020.i09.

———. 2007b. “Extended Rasch Modeling: The eRm Package for the Application of IRT Models in r.” Journal of Statistical Software 20. https://doi.org/10.18637/jss.v020.i09.

Ong, Clarissa W., Benjamin G. Pierce, Douglas W. Woods, Michael P. Twohig, and Michael E. Levin. 2019. “The Acceptance and Action Questionnaire II: An Item Response Theory Analysis.” Journal of Psychopathology and Behavioral Assessment 41 (1): 123–34. https://doi.org/10.1007/s10862-018-9694-2.

Pedersen, Thomas Lin. 2024. patchwork: The Composer of Plots. https://CRAN.R-project.org/package=patchwork.

R Core Team. 2024. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Ren, Kun, and Kenton Russell. 2021. formattable: Create “Formattable” Data Structures. https://CRAN.R-project.org/package=formattable.

Rosseel, Yves. 2012. “lavaan: An R Package for Structural Equation Modeling.” Journal of Statistical Software 48 (2): 1–36. https://doi.org/10.18637/jss.v048.i02.

Rusch, Thomas, Marco Maier, and Reinhold Hatzinger. 2013. “Linear Logistic Models with Relaxed Assumptions in r.” In Algorithms from and for Nature and Life, edited by Berthold Lausen, Dirk van den Poel, and Alfred Ultsch. Studies in Classification, Data Analysis, and Knowledge Organization. New York: Springer. https://doi.org/10.1007/978-3-319-00035-0_34.

Slowikowski, Kamil. 2024. ggrepel: Automatically Position Non-Overlapping Text Labels with “ggplot2”. https://CRAN.R-project.org/package=ggrepel.

Strobl, Carolin, Julia Kopf, and Achim Zeileis. 2015. “Rasch Trees: A New Method for Detecting Differential Item Functioning in the Rasch Model.” Psychometrika 80 (2): 289–316. https://doi.org/10.1007/s11336-013-9388-3.

Strobl, Carolin, Florian Wickelmaier, and Achim Zeileis. 2011. “Accounting for Individual Differences in Bradley-Terry Models by Means of Recursive Partitioning.” Journal of Educational and Behavioral Statistics 36 (2): 135–53. https://doi.org/10.3102/1076998609359791.

Thériault, Rémi. 2023. “lavaanExtra: Convenience Functions for Lavaan.” Journal of Open Source Software 8 (90): 5701. https://doi.org/10.21105/joss.05701.

Trepte, Sabine, and Markus Verbeet, eds. 2010. Allgemeinbildung in Deutschland – Erkenntnisse Aus Dem SPIEGEL Studentenpisa-Test. Wiesbaden: VS Verlag.

Wickelmaier, Florian, and Achim Zeileis. 2018. “Using Recursive Partitioning to Account for Parameter Heterogeneity in Multinomial Processing Tree Models.” Behavior Research Methods 50 (3): 1217–33. https://doi.org/10.3758/s13428-017-0937-z.

Wickham, Hadley. 2007. “Reshaping Data with the Reshape Package.” Journal of Statistical Software 21 (12). https://www.jstatsoft.org/v21/i12/.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

William Revelle. 2024. psych: Procedures for Psychological, Psychometric, and Personality Research. Evanston, Illinois: Northwestern University. https://CRAN.R-project.org/package=psych.

Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2023. knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.

Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.

Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.

Zhu, Hao. 2024. kableExtra: Construct Complex Table with “kable” and Pipe Syntax. https://CRAN.R-project.org/package=kableExtra.

Reuse

CC BY 4.0