# Development and validation of the Attitudes towards Italian Mafias Scale (AIMS), Schepisi et al. #

# STUDY 1 #

# LOAD LIBRARIES

library(tidyverse)
library(ltm)
library(psych)
library(lavaan)
library(semPlot)
library(semTools)
library(mirt)
library(lavaanPlot)
library(dplyr)

# IMPORT THE DATASET; ITEMS ARE ALREADY REVERSED, EXCEPT FOR E3
df<-read.csv("data_matrix_study1.csv", sep=";", header=T)
str(df)

# CHANGE INTEGER VARIABLES INTO NUMERIC
df <- df %>% 
  mutate_at(c(2,7:82), as.numeric)

# CHANGE VARIABLES INTO FACTOR
df$SUB_ID <- factor(df$SUB_ID)
df$GENDER <- factor(df$GENDER)
df$BIRTH_REGION <- factor(df$BIRTH_REGION)
df$EDUCATION <- factor(df$EDUCATION)
df$INCOME <- factor(df$INCOME)

# CREATE DATAFRAME WITH ONLY AIMS ITEMS
df= subset(df,select=c(7:82))

# REVERSE ITEM E3
df[, c('E3')] <- 8 - df[, c('E3')]

# ITEMS SELECTION: STEP 1; EXCLUDE PROBLEMATIC VARIABLES
df.1<- subset( df, select = -c(E3,E10,B6,C31,C32,C33)) 


# ITEMS SELECTION: STEP 2; EXCLUDE VARIABLES WITH EXTREME RESPONSES AND RETAIN THE VARIABLES BELOW
df.2<-subset(df.1,select=c(E5, E7, E9, E11, B3, B4, B7, B12, B19, B20, B21, C1, C6, C7, C8, C10, C13, C15, C18, C23, C24, C25, C27, C29, C30, C34, C35, C37))

# TEST SAMPLE ADEQUACY 
cortest.bartlett(df.2)
KMO(df.2)
det(cor(df.2, use = "pairwise.complete.obs"))

# EXPLORATORY FACTOR ANALYSIS (EFA) #
outparallel2 <- fa.parallel(df.2, fa ="fa", quant = .99, fm = "pa")
Matrix2 <- cor(df.2, use = "pairwise.complete.obs")
scree(Matrix2, pc = F)
fa# source("http://personality-project.org/r/vss.r")

# EFA WITH 3 FACTORS #
outfa.3.2 <- fa(df.2, nfactors = 3, rotate = "oblimin", fm = "pa")
print(outfa.3.2, cut=.1, digits = 2, sort = T)

Matrix2 <- cor(df, use = "pairwise.complete.obs")# correlation matrix for EFA

main<-omega(Matrix2,3, fm = "ml", n.obs = 292)# EFA con ML
main<-omega(Matrix2,3, fm = "pa", n.obs = 292)#EFA con PA
rt<-main$schmid$sl #loadings, communalities, uniquenesses after Schmid-Leiman orthogonalization
eigen < - colSums(rt^2) #eigenvalues

# STUDY 2 #

#load libraries

library(tidyverse)
library(ltm)
library(psych)
library(lavaan)
library(semPlot)
library(semTools)
library(mirt)
library(lme4)
library(car)
library(semTable)
library(lavaanPlot)
library(effsize)

# IMPORT THE DATASET
df_original<-read.csv("data_matrix_study2.csv", sep=";", header=T)


# CHANGE VARIABLES INTO FACTOR
df_original$Gender_birth <- factor(df_original$Gender_birth)
df_original$Region_birth <- factor(df_original$Region_birth)
df_original$Donation <- factor(df_original$Donation)
df_original$Education <- factor(df_original$Education)
df_original$Income <- factor(df_original$Income)
df_original$sub <- factor(df_original$sub)
df_original$DRMR_REST<- factor(df_original$DRMR_REST)
df_original$CCS_REST<- factor(df_original$CCS_REST)

str(df_original)
# CHANGE INTEGER VARIABLES INTO NUMERIC
df_original <- df_original %>% 
  mutate_at(c(2,8,10:78), as.numeric)

df=df_original

str(df)

# RECODE CSS ITEMS SUCH THAT THEY GO IN THE SAME DIRECTION AS THE AIMS
# which is higher scores mean positive attitude towards Mafia and  crime
df[, c('E5', 'E7', 'E9','B3', 'B7', 'B12', 'B19', 'B21', 'C1', 'C15','C30','L1','L2','L6','L8','L9','T3','T5','T6','T8','P1','P2','P4','P6','TLV3','TLV6','ICO2','ICO6')] <- 8 - df[, c('E5', 'E7', 'E9','B3', 'B7', 'B12', 'B19', 'B21', 'C1', 'C15','C30','L1','L2','L6','L8','L9','T3','T5','T6','T8','P1','P2','P4','P6','TLV3','TLV6','ICO2','ICO6')]

# CREATE DATAFRAME WITH ONLY CSS AND AIMS ITEMS
df_scales = subset(df,select=c(10:78))

# CREATE DATAFRAME WITH ONLY CSS ITEMS
dfCSS= subset(df_scales,select = c(29:69))

# CREATE DATAFRAME WITH ONLY AIMS ITEMS
dfAIMS = subset(df_scales,select = c(1:28))

# EXCLUDE ITEMS THAT IN THE EFA IN STUDY 1 LOADED LESS THAN .3 INTO AT LEAST ONE OF THE 3 FACTORS:
dfAIMS2<-subset(dfAIMS, select=-c(C7,C18,C24,C27,C37))


# CONFIRMATORY FACTOR ANALYSIS (CFA) WITH BIFACTOR MODEL WITH THE RESULTING 23 ITEMS OF THE AIMS
# BIFACTOR MODEL WITH ONE GENERAL FACTOR AND THREE SPECIFIC SUB-FACTORS AS FOUND IN EFA IN STUDY 1
bifactor.model <- '
MafiaAttitude=~ B3 + B4 + B7+ B12 + B19 + B20 + B21 + E5 + E7 + E9 + E11 + C1 + C6 + C8 + C10 + C13 + C15 + C23 + C25 + C29 + C30 + C34 + C35 
Beha  =~  B3 + B4 + B7 + B12 + B19 + B20 + B21
EmoCog =~ E5 + E7 + E9 + E11 + C1 + C6 + C10 + C15 + C30
Cog =~   C8 + C13 + C23 + C25 + C29 + C34 + C35  
'

fit.bifactAIMS <- cfa(bifactor.model, data = dfAIMS2, std.lv=TRUE,
                      ordered = names(dfAIMS2), orthogonal = T)# DWLS as estimator for ordinal variables

fitMeasures(fit.bifactAIMS, fit.measures = c("chisq", "df","cfi", "tli", "srmr","rmsea"), output = "matrix")

summary(fit.bifactAIMS, standardized=TRUE)

# TEST PERFOMANCE OF BIFACTOR MODEL AGAINST ALTERNATIVE MODELS

# TRIDIMENSIONAL MODEL  
tri.model <- '
Beha  =~  B3 + B4 + B7 + B12 + B19 + B20 + B21
EmoCog =~ E5 + E7 + E9 + E11 + C1 + C6 + C10 + C15 + C30
Cog =~    C8 + C13 + C23 + C25 + C29 + C34 + C35 
'
fit.triAIMS <- cfa(tri.model, data = dfAIMS2, std.lv=TRUE,
                   ordered = names(dfAIMS2))

fitMeasures(fit.triAIMS, fit.measures = c("chisq", "df","cfi", "tli", "srmr","rmsea"), output = "matrix")

# UNIDIMENSIONAL MODEL 
uni.model <- '
Uni =~ B3 + B4 + B7+ B12 + B19 + B20 + B21 + E5 + E7 + E9 + E11 + C1 + C6 + C8 + C10 + C13 + C15 + C23 + C25 + C29 + C30 + C34 + C35
'
fit.uniAIMS <- cfa(uni.model, data = dfAIMS2, std.lv=TRUE,
                   ordered = names(dfAIMS2))

fitMeasures(fit.uniAIMS, fit.measures = c("chisq", "df","cfi", "tli", "srmr","rmsea"), output = "matrix")

# SECOND-ORDER MODEL
so.model <- '
Beha  =~  B3 + B4 + B7 + B12 + B19 + B20 + B21
EmoCog =~ E5 + E7 + E9 + E11 + C1 + C6 + C10 + C15 + C30
Cog =~    C8 + C13 + C23 + C25 + C29 + C34 + C35 

Gfactor =~ Beha + EmoCog + Cog
'
fit.soAIMS <- cfa(so.model, data = dfAIMS2, std.lv=TRUE,
                  ordered = names(dfAIMS2))

fitMeasures(fit.soAIMS, fit.measures = c("chisq", "df","cfi", "tli", "srmr","rmsea"), output = "matrix")


# MODEL COMPARISON BIFACTOR VS. UNIDIMENSIONAL VS. TRIDIMENSIONAL
anova(fit.uniAIMS, fit.bifactAIMS, fit.triAIMS)


# MODEL COMPARISON BIFACTOR VS. UNIDIMENSIONAL VS. SECOND-ORDER
anova(fit.uniAIMS, fit.bifactAIMS, fit.soAIMS)

# CHECK FOR BIFACTOR MODEL'S INTERNAL CONSISTENCY - MCDONAND'S OMEGA
omegaFromSem(fit.bifactAIMS, flip= FALSE)


# RE-FIT BIFACTOR MODEL WITHOUT ITEMS THAT LOAD LESS THAN .3 ONTO THE GENERAL FACTOR
# EXCLUDE C8, C23, C25, E9
bifactor.model3 <- '
MafiaAttitude=~ B3 + B4 + B7+ B12 + B19 + B20 + B21 + E5 + E7 +  E11 + C1 + C6 + C10 + C13 + C15 +  C29 + C30 + C34 + C35 
Beha  =~  B3 + B4 + B7 + B12 + B19 + B20 + B21
EmoCog =~ E5 + E7 + E11 + C1 + C6 + C10 + C15 + C30
Cog =~ C13 + C29 + C34 + C35  
'

fit.bifactAIMS_3 <- cfa(bifactor.model3, data = dfAIMS2, std.lv=TRUE,
                        ordered = names(dfAIMS2), orthogonal = T)# DWLS as estimator for ordinal variables

fitMeasures(fit.bifactAIMS_3, fit.measures = c("chisq", "df","cfi", "tli", "srmr","rmsea"), output = "matrix")

omegaFromSem(fit.bifactAIMS_3, flip= FALSE)
summary(fit.bifactAIMS_3)


# RE-FIT BIFACTOR MODEL WITHOUT ITEMS THAT LOAD LESS THAN .3 ONTO THE GENERAL FACTOR
# EXCLUDE C10
bifactor.model4 <- '
MafiaAttitude=~ B3 + B4 + B7+ B12 + B19 + B20 + B21 + E5 + E7 + E11 + C1 + C6 + C13 + C15 + C29 + C30 + C34 + C35 
Beha  =~  B3 + B4 + B7 + B12 + B19 + B20 + B21
EmoCog =~ E5 + E7 + E11 + C1 + C6 + C15 + C30
Cog =~ C13 + C29 + C34 + C35  
'

fit.bifactAIMS_4 <- cfa(bifactor.model4, data = dfAIMS2, std.lv=TRUE,
                        ordered = names(dfAIMS2), orthogonal = T)# DWLS as estimator for ordinal variables

fitMeasures(fit.bifactAIMS_4, fit.measures = c("chisq", "df","cfi", "tli", "srmr","rmsea"), output = "matrix")

summary(fit.bifactAIMS_4)

# CHECK FOR BIFACTOR MODEL'S INTERNAL CONSISTENCY
omegaFromSem(fit.bifactAIMS_4,flip=FALSE)

# TEST PERFOMANCE OF THE NEW 18-ITEM BIFACTOR MODEL AGAINST ALTERNATIVE MODELS

# TRIDIMENSIONAL MODEL  
tri.model <- '
Beha  =~  B3 + B4 + B7 + B12 + B19 + B20 + B21
EmoCog =~ E5 + E7 + E11 + C1 + C6 + C15 + C30
Cog =~   C13 + C29 + C34 + C35 
'
fit.triAIMS <- cfa(tri.model, data = dfAIMS2, std.lv=TRUE,
                   ordered = names(dfAIMS2))

fitMeasures(fit.triAIMS, fit.measures = c("chisq", "df","cfi", "tli", "srmr","rmsea"), output = "matrix")

# UNIDIMENSIONAL MODEL 
uni.model <- '
Uni =~ B3 + B4 + B7+ B12 + B19 + B20 + B21 + E5 + E7 + E11 + C1 + C6 + C13 + C15 + C29 + C30 + C34 + C35 
'
fit.uniAIMS <- cfa(uni.model, data = dfAIMS2, std.lv=TRUE,
                   ordered = names(dfAIMS2))

fitMeasures(fit.uniAIMS, fit.measures = c("chisq", "df","cfi", "tli", "srmr","rmsea"), output = "matrix")

# SECOND-ORDER MODEL
so.model <- '
Beha  =~  B3 + B4 + B7 + B12 + B19 + B20 + B21
EmoCog =~ E5 + E7 + E11 + C1 + C6 + C15 + C30
Cog =~   C13 + C29 + C34 + C35 

Gfactor =~ Beha + EmoCog + Cog
'
fit.soAIMS <- cfa(so.model, data = dfAIMS2, std.lv=TRUE,
                   ordered = names(dfAIMS2))

fitMeasures(fit.soAIMS, fit.measures = c("chisq", "df","cfi", "tli", "srmr","rmsea"), output = "matrix")


# MODEL COMPARISON BIFACTOR VS. UNIDIMENSIONAL VS. TRIDIMENSIONAL
anova(fit.uniAIMS, fit.bifactAIMS_4, fit.triAIMS)


# MODEL COMPARISON BIFACTOR VS. UNIDIMENSIONAL VS. SECOND-ORDER
anova(fit.uniAIMS, fit.bifactAIMS_4, fit.soAIMS)

# CONSTRUCT-DISCRIMINANT VALIDITY AGAINST CRIME SENTIMENT SCALE (CSS; GENDRAU ET AL., 1997)
# FIT AND COMPARE A DISCRIMINANT MODEL AGAINST CONVERGENT AND UNIDIMENSIONAL MODELS

# DISCRIMINANT MODEL
discr.model <- '
MafiaAttitude =~ B3 + B4 + B7+ B12 + B19 + B20 + B21 + E5 + E7 + E11 + C1 + C6 + C13 + C15 + C29 + C30 + C34 + C35 
Beha  =~  B3 + B4 + B7 + B12 + B19 + B20 + B21
EmoCog =~ E5 + E7 + E11 + C1 + C6 + C15 + C30
Cog =~   C13 + C29 + C34 + C35  

CSS =~ L1 + L2 +L3 + L4 + L5 + L6 + L7 + L8 + L9 + L10 + T1 + T2 + T3 + T4 + T5 + T6 + T7 + T8 + P1 + P2 + P3 + P4 + P5 + P6 + P7 + TLV1 + TLV2  + TLV3 + TLV4 + TLV5 + TLV6 + TLV7 + TLV8 + TLV9 + TLV10  + ICO1 + ICO2 + ICO3 + ICO4 + ICO5 + ICO6 
L =~ L1 + L2 +L3 + L4 + L5 + L6 + L7 + L8 + L9 + L10
T =~ T1 + T2 + T3 + T4 + T5 + T6 + T7 + T8
P =~ P1 + P2 + P3 + P4 + P5 + P6 + P7
TLV =~ TLV1 + TLV2  + TLV3 + TLV4 + TLV5 + TLV6 + TLV7 + TLV8 + TLV9 + TLV10
ICO =~ ICO1 + ICO2 + ICO3 + ICO4 + ICO5 + ICO6 

CSS ~~ MafiaAttitude
'
fit.discr <- cfa(discr.model, data = df, std.lv=TRUE, orthogonal = T, estimator='MLR')

fitMeasures(fit.discr, fit.measures = c("cfi", "tli", "srmr","rmsea"), output = "matrix")

summary(fit.discr2)

# UNIDIMENSIONAL MODEL
big.uni.model <- '
UNI =~ B3 + B4 + B7+ B12 + B19 + B20 + B21 + E5 + E7 + E11 + C1 + C6 + C13 + C15 + C29 + C30 + C34 + C35  + L1 + L2 +L3 + L4 + L5 + L6 + L7 + L8 + L9 + L10 + T1 + T2 + T3 + T4 + T5 + T6 + T7 + T8 + P1 + P2 + P3 + P4 + P5 + P6 + P7 + TLV1 + TLV2  + TLV3 + TLV4 + TLV5 + TLV6 + TLV7 + TLV8 + TLV9 + TLV10  + ICO1 + ICO2 + ICO3 + ICO4 + ICO5 + ICO6 
'
fit.big.uni <- cfa(big.uni.model, data = df, std.lv=TRUE,estimator='MLR')

fitMeasures(fit.big.uni, fit.measures = c("cfi", "tli", "srmr","rmsea"), output = "matrix") 

summary(fit.big.uni2)

# CONVERGENT MODEL
conv.model <- '
CONV =~ B3 + B4 + B7+ B12 + B19 + B20 + B21 + E5 + E7 + E11 + C1 + C6 + C13 + C15 + C29 + C30 + C34 + C35  + L1 + L2 +L3 + L4 + L5 + L6 + L7 + L8 + L9 + L10 + T1 + T2 + T3 + T4 + T5 + T6 + T7 + T8 + P1 + P2 + P3 + P4 + P5 + P6 + P7 + TLV1 + TLV2  + TLV3 + TLV4 + TLV5 + TLV6 + TLV7 + TLV8 + TLV9 + TLV10  + ICO1 + ICO2 + ICO3 + ICO4 + ICO5 + ICO6 
L =~ L1 + L2 +L3 + L4 + L5 + L6 + L7 + L8 + L9 + L10
Beha  =~  B3 + B4 + B7 + B12 + B19 + B20 + B21
EmoCog =~ E5 + E7 + E11 + C1 + C6 + C15 + C30
Cog =~   C13 + C29 + C34 + C35  
T =~ T1 + T2 + T3 + T4 + T5 + T6 + T7 + T8
P =~ P1 + P2 + P3 + P4 + P5 + P6 + P7
TLV =~ TLV1 + TLV2  + TLV3 + TLV4 + TLV5 + TLV6 + TLV7 + TLV8 + TLV9 + TLV10
ICO =~ ICO1 + ICO2 + ICO3 + ICO4 + ICO5 + ICO6 
'
fit.conv <- cfa(conv.model, data = df, std.lv=TRUE,estimator='MLR',orthogonal=T) 

fitMeasures(fit.conv, fit.measures = c("cfi", "tli", "srmr","rmsea"), output = "matrix") 

summary(fit.conv)

# MODEL COMPARISON WITH DISCRIMINANT, CONVERGENT AND UNIDIMENSIONAL MODELS
anova(fit.conv, fit.discr, fit.big.uni) 

# PLOT DISCRIMINANT, CONVERGENT AND UNIDIMENSIONAL MODELS
omegaFromSem(fit.discr,flip=TRUE)

omegaFromSem(fit.conv, flip=TRUE)

omegaFromSem(fit.big.uni,flip=TRUE)

# CONSTRUCT-CONVERGENT VALIDITY 
# CORRELATIONS BETWEEN AIMS AND CSS SCORES #

dfAIMS3<-subset(dfAIMS2, select=-c(E9,C8,C10,C23,C25))# CREATE DATAFRAME WITH THE FINAL 18-ITEM AIMS BY EXCLUDING OLD ITEMS

dfAIMS3$AIMS_mean=rowMeans(dfAIMS3)# CREATE MEAN AIMS SCORE

dfCSS$CSS_mean=rowMeans(dfCSS)# CREATE MEAN CSS SCORE

df_AIMS_CSS_scores=cbind(dfAIMS3$AIMS_mean,dfCSS$CSS_mean)

df_AIMS_CSS_means= as.data.frame(df_AIMS_CSS_scores)

colnames(df_AIMS_CSS_means)[1] <- "AIMS_mean"
colnames(df_AIMS_CSS_means)[2] <- "CSS_mean" 

corr_AIMS_CSS= cor.test(df_AIMS_CSS_means$AIMS_mean,df_AIMS_CSS_means$CSS_mean, method= "spearman",exact=FALSE)


# PREDICTIVE VALIDITY OF AIMS SCORES ON DONATION BEHAVIOR #
df= cbind(df,df_AIMS_CSS_means)

model.AIMS_Donation = glm(Donation ~ AIMS_mean, data = df, family=binomial)
Anova(model.AIMS_Donation,type="3",test="F")
summary(model.AIMS_Donation)
sd(df$AIMS_mean)
mean(df$AIMS_mean)

# PROBABILITY OF DONATING
exp(-1.4979+(2.534351+0.7912799)*0.5776)/(1+exp(-1.4979+(2.534351+0.7912799)*0.5776)) # +1sd

exp(-1.4979+(2.534351-0.7912799)*0.5776)/(1+exp(-1.4979+(2.534351-0.7912799)*0.5776)) # -1sd

# PREDICTIVE VALIDITY OF FACTOR SCORES 
fscoresAIMS <- lavaan::predict(fit.bifactAIMS_4)

df$fscoresAIMS <- fscoresAIMS[, 1]
df$fscoresBeha <- fscoresAIMS[, 2]
df$fscoresEmoCog <- fscoresAIMS[, 3]
df$fscoresCog <- fscoresAIMS[, 4]


model.AIMS_Donation.fsc = glm(Donation ~ fscoresAIMS, data = df, family=binomial)
Anova(model.AIMS_Donation.fsc,type="3")
summary(model.AIMS_Donation.fsc) 


# MEASUREMENT INVARIANCE - COMPARE DRMR (Apulia, Basilicata, Calabria, Campania, Sicily) VS. REST OF ITALY
dfAIMS.region=df# CREATE NEW DATAFRAME 
dfAIMS.region<-subset(df, select=c(DRMR_REST,B3, B4, B7,B12,B19,B20,B21,E5,E7,E11,C1,C6,C13,C15,C29,C30,C34,C35))

# CHECK RESPONSE FREQUENCY FOR EACH ITEM OF THE AIMS IN BOTH GROUPS TO AVOID EMPTY CELLS
table(dfAIMS.region$E5, dfAIMS.region$DRMR_REST)#  missing values at group 1 response cat 6 
table(dfAIMS.region$E7, dfAIMS.region$DRMR_REST)# OK
table(dfAIMS.region$E11, dfAIMS.region$DRMR_REST)# OK

table(dfAIMS.region$B3, dfAIMS.region$DRMR_REST)# OK
table(dfAIMS.region$B4, dfAIMS.region$DRMR_REST)# OK
table(dfAIMS.region$B7, dfAIMS.region$DRMR_REST)# OK
table(dfAIMS.region$B12, dfAIMS.region$DRMR_REST)# OK
table(dfAIMS.region$B19, dfAIMS.region$DRMR_REST)# OK
table(dfAIMS.region$B20, dfAIMS.region$DRMR_REST)# OK
table(dfAIMS.region$B21, dfAIMS.region$DRMR_REST)#  OK

table(dfAIMS.region$C1, dfAIMS.region$DRMR_REST)# missing values at group 2 response cat 7
table(dfAIMS.region$C6, dfAIMS.region$DRMR_REST)# OK
table(dfAIMS.region$C13, dfAIMS.region$DRMR_REST)# OK
table(dfAIMS.region$C15, dfAIMS.region$DRMR_REST)# missing values at group 1 response cat 6
table(dfAIMS.region$C29, dfAIMS.region$DRMR_REST)# OK
table(dfAIMS.region$C30, dfAIMS.region$DRMR_REST)# OK 
table(dfAIMS.region$C34, dfAIMS.region$DRMR_REST)# OK
table(dfAIMS.region$C35, dfAIMS.region$DRMR_REST)# OK


# COLLAPSE RESPONSE FROM 7 TO 6, MERGING CATEGORIES 6 AND 7 INTO 6 TO AVOID EMPTY CELLS
dfAIMS.region[,2:19]<- lapply(dfAIMS.region[ ,2:19] , 
                              FUN = function(x) recode(x, "1=1; 2=2; 3=3; 4=4;5=5; 6=6; 7=6"))

# FIT 18-ITEM BIFACTOR MODEL FOR MEASUREMENT INVARIANCE 
bifactor.model.region <- '
MafiaAttitude =~ B3 + B4 + B7+ B12 + B19 + B20 + B21 + E5 + E7 + E11 + C1 + C6 + C13 + C15 + C29 + C30 + C34 + C35 
Beha  =~  B3 + B4 + B7 + B12 + B19 + B20 + B21
EmoCog =~ E5 + E7 + E11 + C1 + C6  + C15 + C30
Cog =~   C13 + C29 + C34 + C35  
'

fit.bifactAIMS_Region <- cfa(bifactor.model.region, data = dfAIMS.region, std.lv=TRUE,orthogonal = T, 
                             ordered = names(dfAIMS.region[2:19]))

fitMeasures(fit.bifactAIMS_Region, fit.measures = c("cfi", "tli", "srmr","rmsea"), output = "matrix")

summary(fit.bifactAIMS_Region)


# CONFIGURAL INVARIANCE- BASELINE MODEL
model.config <-  measEq.syntax(configural.model = bifactor.model,
                               data = dfAIMS.region,
                               ordered = names(dfAIMS.region[2:19]),orthogonal=T,
                               parameterization = "delta", ID.fac = "std.lv",
                               ID.cat = "Wu.Estabrook.2016",
                               group = "DRMR_REST",group.equal = "configural")

model.baseline <- as.character(model.config)

fit.model.baseline <- cfa(model.baseline,data = dfAIMS.region,ordered = names(dfAIMS.region[2:19]),group = "DRMR_REST") 

fitmeasures(fit.model.baseline, fit.measures = c("chisq", "df", "tli", "cfi", "rmsea", "srmr"),
            output  = "matrix")

summary(fit.model.baseline)

# METRIC INVARIANCE - THRESHOLDS MODEL
model.thresholds <-   measEq.syntax(configural.model = bifactor.model,
                                    data = dfAIMS.region,
                                    ordered = names(dfAIMS.region[2:19]), orthogonal=T, 
                                    parameterization = "delta", ID.fac = "std.lv",
                                    ID.cat = "Wu.Estabrook.2016",
                                    group = "DRMR_REST",group.equal = c("thresholds"))

thresholds <- as.character(model.thresholds)

fit_thr <- cfa(thresholds, data = dfAIMS.region, ordered = names(dfAIMS.region[2:19]), group = "DRMR_REST")

fitmeasures(fit_thr, fit.measures = c("chisq", "df", "tli", "cfi", "rmsea", "srmr"),
            output ="matrix")

# MODEL COMPARISON BETWEEN BASELINE AND THRESHOLDS MODELS
lavTestLRT(fit.model.baseline, fit_thr)

summary(fit_thr)

# SCALAR INVARIANCE - THRESHOLDS AND LOADINGS MODEL
load_thre <-  measEq.syntax(configural.model = bifactor.model,
                            data = dfAIMS.region,
                            ordered = names(dfAIMS.region[2:19]),orthogonal=T, parameterization = "delta", ID.fac = "std.lv",
                            ID.cat = "Wu.Estabrook.2016",
                            group = "DRMR_REST",group.equal = c("thresholds", "loadings"))

load_thre <- as.character(load_thre)

fit_loth <- cfa(load_thre, ordered =  names(dfAIMS.region[2:19]), data = dfAIMS.region,group = "DRMR_REST")

fitmeasures(fit_loth, fit.measures = c("chisq", "df", "tli", "cfi", "rmsea", "srmr"),
            output ="matrix")

lavTestScore(fit_loth)

# MODEL COMPARISON BETWEEN THRESHOLDS AND THRESHOLDS/LOADINGS MODELS
lavTestLRT(fit_loth, fit_thr)


# MEAN SCORES COMPARISON BETWEEN DRMR VS. REST OF ITALY
df$AIMS_mean_collapsed=rowMeans(dfAIMS.region[,2:19])# mean calculated using scores from collapsed response categories
t.test(df$AIMS_mean_collapsed~df$DRMR_REST)
ES<-cohen.d(df$AIMS_mean_collapsed~df$DRMR_REST)# effect size

df$AIMS_mean_NOTcollapsed=rowMeans(dfAIMS3[,1:18])# mean calculated using scores from NOT collapsed response categories
t.test(df$AIMS_mean_NOTcollapsed~df$DRMR_REST)
ES<-cohen.d(df$AIMS_mean_NOTcollapsed~df$DRMR_REST)# effect size

t.test(df$fscoresAIMS~df$DRMR_REST)# mean calculated using the factor score of the AIMS general factor (created from not collapsed categories)
ES<-cohen.d(df$fscoresAIMS~df$DRMR_REST)# effect size

# CHECK FOR DIFFERENCE IN AGE BETWEEN DRMR AND REST OF ITALY
t.test(df_original$Age~df_original$DRMR_REST)

# CHECK FOR DIFFERENCE IN SOCIO-ECONOMIC STATUS BETWEEN DRMR AND REST OF ITALY
t.test(df_original$SES~df_original$DRMR_REST)

# CHECK FOR DIFFERENCE IN EDUCATION BETWEEN DRMR AND REST OF ITALY
chisq.test(df_original$Education,df_original$DRMR_REST)

# CHECK FOR DIFFERENCE IN INCOME BETWEEN DRMR AND REST OF ITALY
chisq.test(df_original$Income,df_original$DRMR_REST)# TEST INCOME DIFFERENCES

#CALCULATE GENERAL FACTOR SCORE ON BIFACTOR MODEL WITH COLLAPSED RESPONSE CATEGORIES
fscoresAIMS.region <- lavaan::predict(fit.bifactAIMS_Region)

df$fscoresAIMS.collapsed <- fscoresAIMS.region[, 1]

# Export dataframe for analysis in JASP #
write.table(df, file = "dfAIMS.xls",quote=F, sep=";", dec=".", na="", row.names=T, col.names=T)