CFA.Rmd

---
title: "HOMEWORK 3: CFA (version 2)"
author: "Maria Tikhonova"
output: 
  html_document:
    keep_md: no
    toc: yes
editor_options: 
  chunk_output_type: inline
---

# Load packages and data

```{r}
# setwd()
```

## Load packages

```{r, warning=FALSE, message=FALSE}
library(lavaan)
library(semPlot)
library(rio)
library(car)
library(dplyr)
```

## Import data

```{r}
pisa <- import("PISA_MR.sav")
```

Version 2: items ST161Q01HA, ST161Q02HA, ST161Q03HA, ST161Q06HA, ST161Q07HA, ST161Q08HA

```{r}
pisa <- pisa[c('ST161Q01HA', 'ST161Q02HA', 'ST161Q03HA', 'ST161Q06HA', 'ST161Q07HA', 'ST161Q08HA')]

pisa <- rename(pisa,
       "item_1" = 'ST161Q01HA',
       "item_2" = 'ST161Q02HA',
       "item_3" = 'ST161Q03HA',       
       "item_6" = 'ST161Q06HA',     
       "item_7" = 'ST161Q07HA',
       "item_8" = 'ST161Q08HA')
```

**Codebook**

Response options:  
Strongly disagree (1)  
Disagree (2)  
Agree (3)  
Strongly agree (4)  

item_1: I am a good reader  
item_2: I am able to understand difficult texts  
item_3: I read fluently  
item_6: I have always had difficulty with reading  
item_7: I have to read a text several times before completely understanding it.  
item_8: I find it difficult to answer questions about a text.  

Look at the data structure

```{r}
glimpse(pisa)
```

Delete missing values

```{r}
pisa <- na.omit(pisa)
```

Descriptive statistics

```{r}
cat("\nSummary\n")
summary(pisa)

cat("\nFrequency table\n")
table_freq <- sapply(pisa[, 1:6], table)
table_freq
table_freq <- as.data.frame(table_freq)

hist(pisa[,1:6])
```

Next, we will use function `cfa` from `lavaan` for model parameters estimation.

# MODEL SELECTION

## STEP 1. FIRST ATTEMPT

```{r}
model1 <- "F1 =~ item_1 + item_2 + item_3 + item_6 + item_7 + item_8"

fit_mod1 <- cfa(model = model1, data = pisa, estimator = "WLSMV", ordered = TRUE, std.lv = T)
summary(fit_mod1, fit=T, std=T)
```

```{r}
# factor scores only
head(standardizedSolution(fit_mod1, type = "std.all"), 6)
```

```{r}
# fit indexes
fitmeasures(fit_mod1, c("cfi.scaled", "tli.scaled", "rmsea.scaled", "srmr"))
```

```{r, warning = FALSE, message = FALSE}
# reliability
options(scipen = 999)
semTools::reliability(fit_mod1)
```

## STEP 2. RECODING

Recoding items

```{r}
pisa$item_6r <- car::recode(pisa$item_6, '1=4;2=3;3=2;4=1')
pisa$item_7r <- car::recode(pisa$item_7, '1=4;2=3;3=2;4=1')
pisa$item_8r <- car::recode(pisa$item_8, '1=4;2=3;3=2;4=1')
```

Creating the second model

```{r}
model2 <- "F1 =~ item_1 + item_2 + item_3 + item_6r + item_7r + item_8r"

fit_mod2 <- cfa(model = model2, data = pisa, estimator = "WLSMV", ordered = TRUE, std.lv = T)
summary(fit_mod2, fit=T, std=T)
```

```{r}
# factor scores only
head(standardizedSolution(fit_mod2, type = "std.all"), 6)
```

```{r}
# fit indexes
fitmeasures(fit_mod2, c("cfi.scaled", "tli.scaled", "rmsea.scaled", "srmr"))
```

```{r, warning = FALSE, message = FALSE}
# reliability
semTools::reliability(fit_mod2)
```

## STEP 3. DROP THE ITEM

```{r}
model3 <- "F1 =~ item_1 + item_2 + item_6r + item_7r + item_8r"

fit_mod3 <- cfa(model = model3, data = pisa, estimator = "WLSMV", ordered = TRUE, std.lv = T)
summary(fit_mod3, fit=T, std=T)
```

```{r}
# factor scores only
head(standardizedSolution(fit_mod3, type = "std.all"), 5)
```

```{r}
# fit indexes
fitmeasures(fit_mod3, c("cfi.scaled", "tli.scaled", "rmsea.scaled", "srmr"))
```

```{r, warning = FALSE, message = FALSE}
# reliability
semTools::reliability(fit_mod3)
```

## STEP 4. MODIFICATION INDEXES

```{r}
modindices(fit_mod3, minimum.value = 10, sort. = TRUE)
```

```{r}
model4 <- "F1 =~ item_1 + item_2 + item_6r + item_7r + item_8r

item_1 ~~ item_2"

fit_mod4 <- cfa(model = model4, data = pisa, estimator = "WLSMV", ordered = TRUE, std.lv = T)
summary(fit_mod4, fit=T, std=T)
```

```{r}
# factor scores only
head(standardizedSolution(fit_mod4, type = "std.all"), 6)
```

```{r}
# fit indexes
fitmeasures(fit_mod4, c("cfi.scaled", "tli.scaled", "rmsea.scaled", "srmr"))
```

```{r, warning = FALSE, message = FALSE}
# reliability
semTools::reliability(fit_mod4)
```

No further mi

```{r}
modindices(fit_mod4, minimum.value = 10, sort. = TRUE)
# no modindices
```

## STEP 5. DIMENSIONALITY CHANGE

```{r}
model5 <- "
F1 =~ item_1 + item_2
F2 =~ item_6r + item_7r + item_8r
"

fit_mod5 <- cfa(model = model5, data = pisa, estimator = "WLSMV", ordered = TRUE, std.lv = T)
summary(fit_mod5, fit=T, std=T)
```

```{r}
# factor scores only
head(standardizedSolution(fit_mod5, type = "std.all"), 5)
```

```{r}
# fit indexes
fitmeasures(fit_mod5, c("cfi.scaled", "tli.scaled", "rmsea.scaled", "srmr"))
```

```{r, warning = FALSE, message = FALSE}
# reliability
semTools::reliability(fit_mod5)
```

No further mi

```{r}
modindices(fit_mod5, minimum.value = 10, sort. = TRUE)
# no modindeces
```


--------------------------------------------------------------------------------

# Model comparison

fit_mod3: one-factor model, no MI; ("baseline")  
fit_mod4: one-factor model, MI;  
fit_mod5: two-factor model, no MI;

```{r}
anova34 <- anova(fit_mod3, fit_mod4, method = "satorra.2000")
anova34
```

```{r}
anova35 <- anova(fit_mod3, fit_mod5, method = "satorra.2000")
anova35
```

--------------------------------------------------------------------------------

# Plots

## Model 1

```{r}
lavaanPlot::lavaanPlot(model = fit_mod1, coefs = TRUE, covs = TRUE, stand = TRUE)
```

## Model 5

```{r}
lavaanPlot::lavaanPlot(model = fit_mod5, coefs = TRUE, covs = TRUE, stand = TRUE)
```

--------------------------------------------------------------------------------

# Tables

AI was used to create these tables.

## Fit indexes for all models

```{r}
# list of models
fits_all <- list(
  mod1 = fit_mod1,
  mod2 = fit_mod2,
  mod3 = fit_mod3,
  mod4 = fit_mod4,
  mod5 = fit_mod5
)

# extract values needed
get_table <- function(fit, name){
  data.frame(
    Model = name,
    Chisq_scaled = round(fitmeasures(fit, "chisq.scaled"), 3),
    Df_scaled    = fitmeasures(fit, "df.scaled"),
    CFI_scaled   = round(fitmeasures(fit, "cfi.scaled"), 3),
    TLI_scaled   = round(fitmeasures(fit, "tli.scaled"), 3),
    RMSEA_scaled = round(fitmeasures(fit, "rmsea.scaled"), 3),
    SRMR         = round(fitmeasures(fit, "srmr"), 3)
  )
}
table_fit <- do.call(rbind,
                  mapply(get_table,
                         fits_all,
                         names(fits_all),
                         SIMPLIFY = FALSE))
knitr::kable(table_fit)
```


## Loading table

```{r}
extract_loadings <- function(fit, model_name){
  pe <- parameterEstimates(fit, standardized = TRUE)
  
  loadings <- pe[pe$op == "=~", c("rhs", "std.all", "pvalue")]
  
  # for asterisks
  stars <- ifelse(loadings$pvalue < .001, "***",
           ifelse(loadings$pvalue < .01,  "**",
           ifelse(loadings$pvalue < .05,  "*", "")))
  
  # round
  values <- paste0(round(loadings$std.all, 3), stars)
  
  df <- data.frame(
    Item = loadings$rhs,
    value = values
  )
  
  colnames(df)[2] <- model_name
  
  df
}

# extract loadings
loads <- list(
load1 = extract_loadings(fit_mod1, "Mod1"),
load2 = extract_loadings(fit_mod2, "Mod2"),
load3 = extract_loadings(fit_mod3, "Mod3"),
load4 = extract_loadings(fit_mod4, "Mod4"),
load5 = extract_loadings(fit_mod5, "Mod5")
)

# combine by Item
table_load <- Reduce(function(x, y) merge(x, y, by = "Item", all = TRUE),
                         loads)

# round and NA
table_load[is.na(table_load)] <- "×"
knitr::kable(table_load)
```

## All Reliabilities in one output

```{r, warning = FALSE, message = FALSE}
round(semTools::reliability(fit_mod1, what = "omega3"), 3)
round(semTools::reliability(fit_mod2, what = "omega3"), 3)
round(semTools::reliability(fit_mod3, what = "omega3"), 3)
round(semTools::reliability(fit_mod4, what = "omega3"), 3)
round(semTools::reliability(fit_mod5, what = "omega3"), 3)
```