##########

#Analysis

##########

#Clear workspace (run if desired) 

rm(list = ls())

# #List of project directories

# dirs <- list(

#   data = "...",

#   analysis = "...")

# #Load required package

if (!require(psych))  { install.packages("psych") } ;  library(psych)

if (!require(lavaan)) { install.packages("lavaan") } ; library(lavaan)

if (!require(car)) { install.packages("car") } ; library(car)

if (!require(semTools)) { install.packages("semTools") } ; library(semTools)

#Load dataset

load(paste0(dirs$data, "GSE-3.rda"))

#########################################################################

#Split dataset between countries

GSE3_D <- subset(GSE3, subset=(COUN=="1"))

GSE3_UK <- subset(GSE3, subset=(COUN=="2"))

############################

#Step 1: Measurement models

############################

## Tau-congeneric measurement model

ASKU_MM <- 'ASKU =~ ASKU1 + ASKU2 + ASKU3'

# Joint

ASKU_MM.fit <- sem(ASKU_MM, data = GSE3, group = "COUN", estimator = "mlr", missing = "fiml", std.lv = TRUE)

summary(ASKU_MM.fit, standardized = TRUE, fit.measures = TRUE)

# Germany

ASKU_MM_DE.fit <- cfa(ASKU_MM, data = GSE3_D, estimator = "mlr", missing = "fiml", std.lv = TRUE)

summary(ASKU_MM_DE.fit, standardized = TRUE, fit.measures = TRUE)

# UK

ASKU_MM_UK.fit <- cfa(ASKU_MM, data = GSE3_UK, estimator = "mlr", missing = "fiml", std.lv = TRUE)

summary(ASKU_MM_UK.fit, standardized = TRUE, fit.measures = TRUE)

#########################################################################

## Essentially tau-equivalent model (equalized factor loadings within each country)

# Joint

ASKU_equ1 <- 'ASKU =~ c(a1,b1)*ASKU1 + c(a1,b1)*ASKU2 + c(a1,b1)*ASKU3'

ASKU_equ1.fit <- sem(ASKU_equ1, data = GSE3, group="COUN",estimator="MLR",std.lv=TRUE)

summary(ASKU_equ1.fit,standardized = TRUE, fit.measures = TRUE)

# Germany

ASKU_equ_DE <- 'ASKU =~ a1*ASKU1 + a1*ASKU2 + a1*ASKU3'

ASKU_equ_DE.fit <- cfa(ASKU_equ_DE, data = GSE3_D, estimator = "mlr", missing = "fiml", std.lv = TRUE)

summary(ASKU_equ_DE.fit, standardized = TRUE, fit.measures = TRUE)

# UK

ASKU_equ_UK <- 'ASKU =~ a1*ASKU1 + a1*ASKU2 + a1*ASKU3'

ASKU_equ_UK.fit <- cfa(ASKU_equ_UK, data = GSE3_UK, estimator = "mlr", missing = "fiml", std.lv = TRUE)

summary(ASKU_equ_UK.fit, standardized = TRUE, fit.measures = TRUE)

#########################################################################

## Tau-equivalent measurement model

ASKU_equ2 <- 'ASKU =~ a1*ASKU1 + a1*ASKU2 + a1*ASKU3

             # equal residual variance

              ASKU1~~b1*ASKU1

              ASKU2~~b1*ASKU2

              ASKU3~~b1*ASKU3'

# Germany

ASKU_equ2_DE.fit <- cfa(ASKU_equ2, data = GSE3_D, estimator = "mlr", missing = "fiml", std.lv = TRUE)

summary(ASKU_equ2_DE.fit, standardized = TRUE, fit.measures = TRUE)

# UK

ASKU_equ2_UK.fit <- cfa(ASKU_equ2, data = GSE3_UK, estimator = "mlr", missing = "fiml", std.lv = TRUE)

summary(ASKU_equ2_UK.fit, standardized = TRUE, fit.measures = TRUE)

#########################################################################

## Paralleltest measurement model (+ equalized intercepts)

ASKU_par <- 'ASKU =~ a1*ASKU1 + a1*ASKU2 + a1*ASKU3

             # equal residual variance

             ASKU1~~b1*ASKU1

             ASKU2~~b1*ASKU2

             ASKU3~~b1*ASKU3

             # equal intercepts

             ASKU1 ~ int2*1

             ASKU2 ~ int2*1

             ASKU3 ~ int2*1'

# Germany

ASKU_par_DE.fit <- cfa(ASKU_par, data = GSE3_D, estimator = "mlr", missing = "fiml", std.lv = TRUE)

summary(ASKU_par_DE.fit, standardized = TRUE, fit.measures = TRUE)

# UK

ASKU_par_UK.fit <- cfa(ASKU_par, data = GSE3_UK, estimator = "mlr", missing = "fiml", std.lv = TRUE)

summary(ASKU_par_UK.fit, standardized = TRUE, fit.measures = TRUE)

#####################################################

#Step 2: Scale values, selectivity, and correlations

#####################################################

# Germany

attach(GSE3_D)

describe((ASKU1+ASKU2+ASKU3)/3)

describe(ASKU1)

describe(ASKU2)

describe(ASKU3)

alpha(data.frame(ASKU1,ASKU2,ASKU3))

cor.test(ASKU1, ASKU2, use = "pairwise.complete.obs")

cor.test(ASKU1, ASKU3, use = "pairwise.complete.obs")

cor.test(ASKU2, ASKU3, use = "pairwise.complete.obs")

detach(GSE3_D)

# UK

attach(GSE3_UK)

describe((ASKU1+ASKU2+ASKU3)/3)

describe(ASKU1)

describe(ASKU2)

describe(ASKU3)

alpha(data.frame(ASKU1,ASKU2,ASKU3))

cor.test(ASKU1, ASKU2, use = "pairwise.complete.obs")

cor.test(ASKU1, ASKU3, use = "pairwise.complete.obs")

cor.test(ASKU2, ASKU3, use = "pairwise.complete.obs")

detach(GSE3_UK)

#####################

#Step 3: Reliability

#####################

## McDonald's omega (essentially tau-equivalent model)

ASKU_equ <- 'ASKU =~ a1*ASKU1 + a1*ASKU2 + a1*ASKU3'

# Germany

ASKU_equ_DE.fit <- cfa(ASKU_equ, data = GSE3_D, estimator = "mlr", missing = "fiml", std.lv = TRUE)

semTools::reliability(ASKU_equ_DE.fit)

# UK

ASKU_equ_UK.fit <- cfa(ASKU_equ, data = GSE3_UK, estimator = "mlr", missing = "fiml", std.lv = TRUE)

semTools::reliability(ASKU_equ_UK.fit)

#########################################################################

## Retest reliability

# Germany

attach(GSE3_D)

ASKU <- ASKU1+ASKU2+ASKU3

ASKUrt <- ASKU1rt+ASKU2rt+ASKU3rt

cor.test(ASKU, ASKUrt, use = "pairwise.complete.obs")

detach(GSE3_D)

# UK

attach(GSE3_UK)

ASKU <- ASKU1+ASKU2+ASKU3

ASKUrt <- ASKU1rt+ASKU2rt+ASKU3rt

cor.test(ASKU, ASKUrt, use = "pairwise.complete.obs")

detach(GSE3_UK)

############################

#Step 4: Criterion validity

############################

## Germany

attach(GSE3_D)

# BFI-2-XS

EXTR <- EXTR1R+EXTR2+EXTR3

AGRE <- AGRE1+AGRE2R+AGRE3

CONS <- CONS1R+CONS2R+CONS3

NEGA <- NEGA1+NEGA2+NEGA3R

OPEN <- OPEN1+OPEN2R+OPEN3

ASKU <- ASKU1+ASKU2+ASKU3

cor.test(ASKU, EXTR, use = "pairwise.complete.obs")

cor.test(ASKU, AGRE, use = "pairwise.complete.obs")

cor.test(ASKU, CONS, use = "pairwise.complete.obs")

cor.test(ASKU, NEGA, use = "pairwise.complete.obs")

cor.test(ASKU, OPEN, use = "pairwise.complete.obs")

# RSES

RSES <- RSES1+RSES2R+RSES3+RSES4+RSES5R+RSES6R+RSES7+RSES8R+RSES9R+RSES10

cor.test(ASKU, RSES, use = "pairwise.complete.obs")

# Health

cor.test(ASKU, HEAL, use = "pairwise.complete.obs")

# Employment status

EMPL.selfempl <- recode(EMPL, "2 = 2; 1 = 1; else = NA")

EMPL.unempl <- recode(EMPL, "3:4 = 2; 1:2 = 1; else = NA")

EMPL.retired <- recode(EMPL, "5 = 2; 8 = 2; 1:2 = 1; else = NA")

EMPL.student <- recode(EMPL, "6:7 = 2; 1:2 = 1; else = NA")

cor.test(ASKU, EMPL.unempl, use = "pairwise.complete.obs")

cor.test(ASKU, EMPL.selfempl, use = "pairwise.complete.obs")

cor.test(ASKU, EMPL.retired, use = "pairwise.complete.obs")

cor.test(ASKU, EMPL.student, use = "pairwise.complete.obs")

# Income

cor.test(ASKU, INCO, use = "pairwise.complete.obs")

# Education

cor.test(ASKU, SCHO, use = "pairwise.complete.obs")

# Age

cor.test(ASKU, AGE, use = "pairwise.complete.obs")

# Gender

cor.test(ASKU, SEX, use = "pairwise.complete.obs")

# I-8

URGE <- URGE1+URGE2

PREM <- PREM1+PREM2

PERS <- PERS1+PERS2

SENS <- SENS1+SENS2

cor.test(ASKU, URGE, use = "pairwise.complete.obs")

cor.test(ASKU, PREM, use = "pairwise.complete.obs")

cor.test(ASKU, PERS, use = "pairwise.complete.obs")

cor.test(ASKU, SENS, use = "pairwise.complete.obs")

# IE-4

ILOC <- ILOC1+ILOC2

ELOC <- ELOC1+ELOC2

cor.test(ASKU, ILOC, use = "pairwise.complete.obs")

cor.test(ASKU, ELOC, use = "pairwise.complete.obs")

# KUSIV3

KUSIV <- KUSI1+KUSI2R+KUSI3

cor.test(ASKU, KUSIV, use = "pairwise.complete.obs")

# L-1

cor.test(ASKU, LISA1, use = "pairwise.complete.obs")

# R-1

cor.test(ASKU, RISK1, use = "pairwise.complete.obs")

# PESS

IPEF <- IPEF1+IPEF2

EPEF <- EPEF1+EPEF2

cor.test(ASKU, IPEF, use = "pairwise.complete.obs")

cor.test(ASKU, EPEF, use = "pairwise.complete.obs")

# SOP2

SOP <- PESS1R+OPTI1

cor.test(ASKU, SOP, use = "pairwise.complete.obs")

# KSE-G

SDPQ <- SDPQ1+SDPQ2+SDPQ3

SDNQ <- SDNQ1+SDNQ2+SDNQ3

cor.test(ASKU, SDPQ, use = "pairwise.complete.obs")

cor.test(ASKU, SDNQ, use = "pairwise.complete.obs")

detach(GSE3_D)

#########################################################################

## UK

attach(GSE3_UK)

# BFI-2-XS

EXTR <- EXTR1R+EXTR2+EXTR3

AGRE <- AGRE1+AGRE2R+AGRE3

CONS <- CONS1R+CONS2R+CONS3

NEGA <- NEGA1+NEGA2+NEGA3R

OPEN <- OPEN1+OPEN2R+OPEN3

ASKU <- ASKU1+ASKU2+ASKU3

cor.test(ASKU, EXTR, use = "pairwise.complete.obs")

cor.test(ASKU, AGRE, use = "pairwise.complete.obs")

cor.test(ASKU, CONS, use = "pairwise.complete.obs")

cor.test(ASKU, NEGA, use = "pairwise.complete.obs")

cor.test(ASKU, OPEN, use = "pairwise.complete.obs")

# RSES

RSES <- RSES1+RSES2R+RSES3+RSES4+RSES5R+RSES6R+RSES7+RSES8R+RSES9R+RSES10

cor.test(ASKU, RSES, use = "pairwise.complete.obs")

# Health

cor.test(ASKU, HEAL, use = "pairwise.complete.obs")

# Employment status

EMPL.selfempl <- recode(EMPL, "2 = 2; 1 = 1; else = NA")

EMPL.unempl <- recode(EMPL, "3:4 = 2; 1:2 = 1; else = NA")

EMPL.retired <- recode(EMPL, "5 = 2; 8 = 2; 1:2 = 1; else = NA")

EMPL.student <- recode(EMPL, "6:7 = 2; 1:2 = 1; else = NA")

cor.test(ASKU, EMPL.unempl, use = "pairwise.complete.obs")

cor.test(ASKU, EMPL.selfempl, use = "pairwise.complete.obs")

cor.test(ASKU, EMPL.retired, use = "pairwise.complete.obs")

cor.test(ASKU, EMPL.student, use = "pairwise.complete.obs")

#Income

cor.test(ASKU, INCO, use = "pairwise.complete.obs")

#Education

cor.test(ASKU, SCHO, use = "pairwise.complete.obs")

# Age

cor.test(ASKU, AGE, use = "pairwise.complete.obs")

# Gender

cor.test(ASKU, SEX, use = "pairwise.complete.obs")

# I-8

URGE <- URGE1+URGE2

PREM <- PREM1+PREM2

PERS <- PERS1+PERS2

SENS <- SENS1+SENS2

cor.test(ASKU, URGE, use = "pairwise.complete.obs")

cor.test(ASKU, PREM, use = "pairwise.complete.obs")

cor.test(ASKU, PERS, use = "pairwise.complete.obs")

cor.test(ASKU, SENS, use = "pairwise.complete.obs")

# IE-4

ILOC <- ILOC1+ILOC2

ELOC <- ELOC1+ELOC2

cor.test(ASKU, ILOC, use = "pairwise.complete.obs")

cor.test(ASKU, ELOC, use = "pairwise.complete.obs")

# KUSIV3

KUSIV <- KUSI1+KUSI2R+KUSI3

cor.test(ASKU, KUSIV, use = "pairwise.complete.obs")

# L-1

cor.test(ASKU, LISA1, use = "pairwise.complete.obs")

# R-1

cor.test(ASKU, RISK1, use = "pairwise.complete.obs")

# PEKS

IPEF <- IPEF1+IPEF2

EPEF <- EPEF1+EPEF2

cor.test(ASKU, IPEF, use = "pairwise.complete.obs")

cor.test(ASKU, EPEF, use = "pairwise.complete.obs")

# SOP2

SOP <- PESS1R+OPTI1

cor.test(ASKU, SOP, use = "pairwise.complete.obs")

# KSE-G

SDPQ <- SDPQ1+SDPQ2+SDPQ3

SDNQ <- SDNQ1+SDNQ2+SDNQ3

cor.test(ASKU, SDPQ, use = "pairwise.complete.obs")

cor.test(ASKU, SDNQ, use = "pairwise.complete.obs")

detach(GSE3_UK)

################################

#Step 5: Measurement invariance

################################

## Measurement invariance based on essentially tau-equivalent model

# Configural invariance  = metric invariance

ASKU_metr <- 'ASKU = ~c(a1, b1)*ASKU1+c(a1,b1)*ASKU2+c(a1,b1)*ASKU3

              ASKU1 ~ c(1,1)*1

              ASKU2 ~ c(NA,NA)*1

              ASKU3 ~ c(NA,NA)*1

              ASKU ~ c(NA,NA)*1'

ASKU_metr.fit <- cfa(ASKU_metr, data = GSE3, group = "COUN",estimator="mlr", missing = "fiml")

summary(ASKU_metr.fit,standardized = T, fit.measures = T)

# Scalar invariance

ASKU_scalar <- 'ASKU=~c(a1,a1)*ASKU1+c(a1,a1)*ASKU2+c(a1,a1)*ASKU3

                ASKU1 ~ c(1,1)*1

                ASKU2 ~ c(c1,c1)*1

                ASKU3 ~ c(d1,d1)*1

                ASKU ~ c(NA,NA)*1'

ASKU_scalar.fit <- cfa(ASKU_scalar, data = GSE3, group = "COUN",estimator="mlr", missing = "fiml")

summary(ASKU_scalar.fit,standardized = T, fit.measures = T)

##########################

#Step 6: Reference values

##########################

#Quote 1: male, lower education, 18-29

#Quote 2: male, lower education, 30-49

#Quote 3: male, lower education, 50-69

#Quote 4: male, middle education, 18-29

#Quote 5: male, middle education, 30-49

#Quote 6: male, middle education, 50-69

#Quote 7: male, upper education, 18-29

#Quote 8: male, upper education, 30-49

#Quote 9: male, upper education, 50-69

#Quote 10: female, lower education, 18-29

#Quote 11: female, lower education, 30-49

#Quote 12: female, lower education, 50-69

#Quote 13: female, middle education, 18-29

#Quote 14: female, middle education, 30-49

#Quote 15: female, middle education, 50-69

#Quote 16: female, uppe# Strict invariance

#Quote 17: female, upper education, 30-49

#Quote 18: female, upper education, 50-69

# Germany

attach(GSE3_D)

ASKU <- (ASKU1+ASKU2+ASKU3)/3

tapply(ASKU, SEX, describe)

AGE1 <- subset(GSE3, AGE == 18 | AGE == 19 | AGE == 20 | AGE == 21 | AGE == 22 | AGE == 23 | AGE == 24 | AGE == 25 | AGE == 26

               | AGE == 27 | AGE == 28 | AGE == 29)

AGE2 <- subset(GSE3, AGE == 30 | AGE == 31 | AGE == 32 | AGE == 33 | AGE == 34 | AGE == 35 | AGE == 36 | AGE == 37 | AGE == 38

               | AGE == 39 | AGE == 40 | AGE == 41 | AGE == 42 | AGE == 43 | AGE == 44 | AGE == 45 | AGE == 46 | AGE == 47

               | AGE == 48 | AGE == 49)

AGE3 <- subset(GSE3, AGE == 50 | AGE == 51 | AGE == 52 | AGE == 53 | AGE == 54 | AGE == 55 | AGE == 56 | AGE == 57 | AGE == 58

               | AGE == 59 | AGE == 60 | AGE == 61 | AGE == 62 | AGE == 63 | AGE == 64 | AGE == 65 | AGE == 66 | AGE == 67

               | AGE == 68 | AGE == 69)

detach(GSE3_D)

attach(AGE1)

ASKU <- (ASKU1+ASKU2+ASKU3)/3

describe(ASKU)

detach(AGE1)

attach(AGE2)

ASKU <- (ASKU1+ASKU2+ASKU3)/3

describe(ASKU)

detach(AGE2)

attach(AGE3)

ASKU <- (ASKU1+ASKU2+ASKU3)/3

describe(ASKU)

detach(AGE3)

#########################################################################

# UK

attach(GSE3_UK)

ASKU <- (ASKU1+ASKU2+ASKU3)/3

tapply(ASKU, SEX, describe)

AGE1U <- subset(GSE3, AGE == 18 | AGE == 19 | AGE == 20 | AGE == 21 | AGE == 22 | AGE == 23 | AGE == 24 | AGE == 25 | AGE == 26

               | AGE == 27 | AGE == 28 | AGE == 29)

AGE2U <- subset(GSE3, AGE == 30 | AGE == 31 | AGE == 32 | AGE == 33 | AGE == 34 | AGE == 35 | AGE == 36 | AGE == 37 | AGE == 38

               | AGE == 39 | AGE == 40 | AGE == 41 | AGE == 42 | AGE == 43 | AGE == 44 | AGE == 45 | AGE == 46 | AGE == 47

               | AGE == 48 | AGE == 49)

AGE3U <- subset(GSE3, AGE == 50 | AGE == 51 | AGE == 52 | AGE == 53 | AGE == 54 | AGE == 55 | AGE == 56 | AGE == 57 | AGE == 58

               | AGE == 59 | AGE == 60 | AGE == 61 | AGE == 62 | AGE == 63 | AGE == 64 | AGE == 65 | AGE == 66 | AGE == 67

               | AGE == 68 | AGE == 69)

detach(GSE3_UK)

attach(AGE1U)

ASKU <- (ASKU1+ASKU2+ASKU3)/3

describe(ASKU)

detach(AGE1U)

attach(AGE2U)

ASKU <- (ASKU1+ASKU2+ASKU3)/3

describe(ASKU)

detach(AGE2U)

attach(AGE3U)

ASKU <- (ASKU1+ASKU2+ASKU3)/3

describe(ASKU)

detach(AGE3U)

################################

#Step 7: Descriptive statistics

################################

## Germany

attach(GSE3_D)

# Age

summary(AGE)

sd(AGE)

# Proportion women

table(SEX)

# Educational level

table(QUOT)

detach(GSE3_D)

#########################################################################

## UK

attach(GSE3_UK)

# Age

summary(AGE)

sd(AGE)

# Proportion women

table(SEX)

# Educational level

table(QUOT)

detach(GSE3_UK)