EFA and CFA - Love languages (Study 3)
Libraries used
library(kableExtra) # for tables
library(psych)
library(ggplot2)
library(MVN) # for mvn()
library(parameters) # efa_to_cfa()
library(lavaan) # for cfa()
library(car) # for qqPlot()
library(semTools) # for AVE, compRelSEM(), discriminantValidity
library(apaTables) # for apa.cor.table()Data
love <- read.csv("love_cloud.csv")EFA and CFA datasets
Randomly split the original dataset into two, one for the EFA, and one for the CFA.
set.seed(243)
n <- floor(nrow(love)/2)
w_efa <- sample(1:nrow(love), n)
w_cfa <- setdiff(1:nrow(love), w_efa)
# leave only item variables (block 1)
w1 <- which(names(love) == "Q1")
w2 <- which(names(love) == "Q206")
vars <- names(love)[w1:w2]
aux_efa <- love[w_efa,vars]
aux_cfa <- love[w_cfa,vars]
f <- function(x){
1*(x == "Not at all") + 2*(x == "A little bit") + 3*(x == "Somewhat") + 4*(x == "Quite a bit") + 5*(x == "Very much")
}
b1 <- as.data.frame(f(as.matrix(aux_efa)))
b2 <- as.data.frame(f(as.matrix(aux_cfa)))Sample characteristics
Age
summary(love$age[w_efa])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18.00 33.00 41.50 44.04 55.00 85.00sd(love$age[w_efa], na.rm = T)## [1] 14.39224summary(love$age[w_cfa])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 19.00 32.00 41.00 43.55 53.25 87.00sd(love$age[w_cfa], na.rm = T)## [1] 15.04731Gender
summary(factor((love$gender[w_efa]))) ## Man Non-binary Woman
## 81 3 264summary(factor((love$gender[w_efa])))/length(w_efa)## Man Non-binary Woman
## 0.23275862 0.00862069 0.75862069summary(factor((love$gender[w_cfa]))) ## Man Non-binary Woman
## 87 3 258summary(factor((love$gender[w_cfa])))/length(w_cfa)## Man Non-binary Woman
## 0.25000000 0.00862069 0.74137931transgender
summary(factor((love$transgender[w_efa]))) ## No Yes
## 345 3summary(factor((love$transgender[w_efa])))/length(w_efa)## No Yes
## 0.99137931 0.00862069summary(factor((love$transgender[w_cfa]))) ## No
## 348summary(factor((love$transgender[w_cfa])))/length(w_cfa)## No
## 1Sexual orientation
summary(factor((love$orientation[w_efa])))## Asexual Bisexual Gay/Lesbian
## 9 25 9
## Heterosexual Pansexual Prefer not to answer
## 295 4 3
## Sexually Fluid Something else
## 2 1summary(factor((love$orientation[w_efa])))/length(w_efa)## Asexual Bisexual Gay/Lesbian
## 0.025862069 0.071839080 0.025862069
## Heterosexual Pansexual Prefer not to answer
## 0.847701149 0.011494253 0.008620690
## Sexually Fluid Something else
## 0.005747126 0.002873563summary(factor((love$orientation[w_cfa])))## Asexual Bisexual Gay/Lesbian
## 13 35 6
## Heterosexual Prefer not to answer Sexually Fluid
## 292 1 1summary(factor((love$orientation[w_cfa])))/length(w_cfa)## Asexual Bisexual Gay/Lesbian
## 0.037356322 0.100574713 0.017241379
## Heterosexual Prefer not to answer Sexually Fluid
## 0.839080460 0.002873563 0.002873563Race
summary(factor(c(love$race1[w_efa], love$race2[w_efa], love$race3[w_efa], love$race4[w_efa])))## African/Black Asian
## 35 15
## Hispanic/Latinx Native American
## 40 11
## Native Hawaiian/ Pacific Islander Prefer not to answer
## 2 1
## Something else White/ European
## 5 275
## NA's
## 1008summary(factor(c(love$race1[w_efa], love$race2[w_efa], love$race3[w_efa], love$race4[w_efa])))/length(w_efa)## African/Black Asian
## 0.100574713 0.043103448
## Hispanic/Latinx Native American
## 0.114942529 0.031609195
## Native Hawaiian/ Pacific Islander Prefer not to answer
## 0.005747126 0.002873563
## Something else White/ European
## 0.014367816 0.790229885
## NA's
## 2.896551724summary(factor(c(love$race1[w_cfa], love$race2[w_cfa], love$race3[w_cfa], love$race4[w_cfa])))## African/Black Asian
## 32 13
## Hispanic/Latinx Middle Eastern/North African
## 38 2
## Native American Native Hawaiian/ Pacific Islander
## 8 1
## Prefer not to answer Something else
## 1 3
## White/ European NA's
## 269 1025summary(factor(c(love$race1[w_cfa], love$race2[w_cfa], love$race3[w_cfa], love$race4[w_cfa])))/length(w_cfa)## African/Black Asian
## 0.091954023 0.037356322
## Hispanic/Latinx Middle Eastern/North African
## 0.109195402 0.005747126
## Native American Native Hawaiian/ Pacific Islander
## 0.022988506 0.002873563
## Prefer not to answer Something else
## 0.002873563 0.008620690
## White/ European NA's
## 0.772988506 2.945402299Rel length
mean(love$rel_length[w_efa], na.rm = T)## [1] 15.31063median(love$rel_length[w_efa], na.rm = T)## [1] 11.3sd(love$rel_length[w_efa], na.rm = T)## [1] 12.85766mean(love$rel_length[w_cfa], na.rm = T)## [1] 15.00403median(love$rel_length[w_cfa], na.rm = T)## [1] 10sd(love$rel_length[w_cfa], na.rm = T)## [1] 13.50026Rel type
summary(factor(love$rel_type[w_efa]))## Dating
## 8
## Engaged
## 36
## Married
## 242
## Partnered/committed, but not engaged or married
## 62summary(factor(love$rel_type[w_efa]))/length(w_efa)## Dating
## 0.02298851
## Engaged
## 0.10344828
## Married
## 0.69540230
## Partnered/committed, but not engaged or married
## 0.17816092summary(factor(love$rel_type[w_cfa]))## Dating
## 8
## Engaged
## 28
## Married
## 245
## Partnered/committed, but not engaged or married
## 67summary(factor(love$rel_type[w_cfa]))/length(w_cfa)## Dating
## 0.02298851
## Engaged
## 0.08045977
## Married
## 0.70402299
## Partnered/committed, but not engaged or married
## 0.19252874Living together
summary(factor(love$live[w_efa]))## No Yes
## 11 337summary(factor(love$live[w_efa]))/length(w_efa)## No Yes
## 0.0316092 0.9683908summary(factor(love$live[w_cfa]))## No Yes
## 8 340summary(factor(love$live[w_cfa]))/length(w_cfa)## No Yes
## 0.02298851 0.97701149EFA
Checking multivariate normality on the EFA dataset
Descriptive stats:
mvn1 <- mvn(b1)
mvn1_df <- data.frame(mvn1$Descriptives)
write.csv(mvn1_df, "mvn1_cloud.csv")
median(mvn1_df$Skew)## [1] -1.234714min(mvn1_df$Skew)## [1] -1.889561max(mvn1_df$Skew)## [1] -0.101065median(mvn1_df$Kurtosis)## [1] 0.9469358min(mvn1_df$Kurtosis)## [1] -1.347901max(mvn1_df$Kurtosis)## [1] 3.46337Number of potential outliers:
cutoff <- qchisq(1-0.001, ncol(b1))
mahal <- mahalanobis(b1, colMeans(b1), cov(b1))
s <- summary(mahal < cutoff); s## Mode FALSE TRUE
## logical 53 295Henze-Zirkler test for multivariate normality:
mvn1$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 3.489327 0 NOThe data is not multivariate normal (do principal axis factoring).
Correlations
cor1 <- cor(b1)
cor1_mat <- as.matrix(cor1)
diag(cor1_mat) <- NA
cor_min <- min(apply(cor1_mat, 2, min, na.rm = T)); cor_min## [1] 0.2115231cor_max <- max(apply(cor1_mat, 2, max, na.rm = T)); cor_max## [1] 0.7954706Correlations between two items varied from 0.2115231 to 0.7954706.
Factorability of the data (need KMO > 0.6)
X <- b1
# Kaiser-Meyer-Olkin (KMO) measure
KMO(r = cor1)## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor1)
## Overall MSA = 0.97
## MSA for each item =
## Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q11 Q12 Q13 Q14 Q15 Q16
## 0.97 0.98 0.98 0.97 0.98 0.98 0.98 0.98 0.98 0.98 0.98 0.98 0.98 0.98 0.97 0.97
## Q17 Q18 Q19 Q20 Q100 Q101 Q102 Q103 Q104 Q105 Q106 Q107 Q122 Q145 Q146 Q147
## 0.97 0.98 0.97 0.97 0.97 0.98 0.97 0.97 0.98 0.97 0.96 0.96 0.98 0.97 0.98 0.96
## Q148 Q149 Q150 Q151 Q152 Q161 Q162 Q163 Q164 Q165 Q166 Q167 Q168 Q169 Q179 Q180
## 0.96 0.98 0.97 0.97 0.97 0.98 0.97 0.97 0.97 0.97 0.97 0.96 0.97 0.96 0.96 0.97
## Q181 Q182 Q183 Q184 Q185 Q186 Q187 Q198 Q199 Q200 Q201 Q202 Q203 Q204 Q205 Q206
## 0.98 0.97 0.96 0.97 0.98 0.97 0.97 0.97 0.97 0.98 0.98 0.98 0.98 0.97 0.98 0.97Do Bartlett’s Test of Sphericity only if participants per item is low. Usually between 3:1 and 5:1. Small p-values (< 0.05) indicate that a factor analysis may be appropriate.
cortest.bartlett(cor1, n = nrow(X))## $chisq
## [1] 20004.16
##
## $p.value
## [1] 0
##
## $df
## [1] 2016A positive determinant means the factor analysis will probably run.
det(cor1)## [1] 1.798511e-27Number of factors to extract
#library(ggplot2)
fafit <- fa(cor1, nfactors = 20, n.obs = nrow(X), rotate = "oblimin", fm = "pa")## Loading required namespace: GPArotation## Warning in GPFoblq(L, Tmat = Tmat, normalize = normalize, eps = eps, maxit =
## maxit, : convergence not obtained in GPFoblq. 1000 iterations used.# n_factors <- length(fafit$e.values)
n_factors <- 10
scree <- data.frame(
Factor_n = as.factor(1:n_factors),
Eigenvalue = fafit$e.values[1:n_factors])
ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) +
geom_point() + geom_line() +
xlab("Number of factors") +
ylab("Initial eigenvalue") +
labs( title = "Scree Plot",
subtitle = "(Based on the unreduced correlation matrix)")
scree## Factor_n Eigenvalue
## 1 1 32.5230691
## 2 2 2.9010950
## 3 3 2.1184491
## 4 4 1.8335764
## 5 5 1.2721213
## 6 6 1.2508326
## 7 7 1.1160389
## 8 8 1.0281000
## 9 9 0.9734206
## 10 10 0.8653245Eigenvalues less than 1 reflect potentially unstable factors. No more than 8 factors.
Parallel analysis
parallel <- fa.parallel(cor1, n.obs = nrow(X), fa = "fa", fm = "pa")
## Parallel analysis suggests that the number of factors = 6 and the number of components = NA1 factor
Number of factors set to 1. After an initial elimination of items with loadings less than 0.32, we further eliminated items aiming to keep at least 2 items from each love language. This led to the following solution.
nfactors <- 1
remove <- c("Q1", "Q2", "Q3", "Q4", "Q6", "Q8", "Q9", "Q11", "Q12", "Q15", "Q16", "Q17", "Q18", "Q20", "Q100", "Q102", "Q103", "Q105", "Q106", "Q107", "Q122", "Q145", "Q146", "Q147", "Q148", "Q150", "Q151", "Q152", "Q161", "Q163", "Q166", "Q167", "Q169", "Q179", "Q180", "Q182", "Q183", "Q184", "Q185", "Q187", "Q198", "Q199", "Q202", "Q206")
#remove <- NA
cols <- setdiff(names(X), remove)
cor2 <- cor(X[,cols])
fa1 <- fa(r = cor2, n.obs = nrow(X[,cols]), nfactors = nfactors,
fm = "pa", max.iter = 100, rotate = "oblimin")
print(fa1$loadings, cutoff = 0.319)##
## Loadings:
## PA1
## Q5 0.729
## Q7 0.770
## Q10 0.823
## Q13 0.781
## Q14 0.733
## Q19 0.706
## Q101 0.766
## Q104 0.768
## Q149 0.748
## Q162 0.727
## Q164 0.796
## Q165 0.788
## Q168 0.771
## Q181 0.773
## Q186 0.807
## Q200 0.813
## Q201 0.758
## Q203 0.793
## Q204 0.740
## Q205 0.771
##
## PA1
## SS loadings 11.816
## Proportion Var 0.591fa1$RMSEA## RMSEA lower upper confidence
## 0.09037319 0.08326465 0.09794797 0.90000000fa1$TLI## [1] 0.9000017fa1$Vaccounted## PA1
## SS loadings 11.8160665
## Proportion Var 0.5908033# Factor
factors1 <- efa_to_cfa(fa1)
aux2 <- strsplit(factors1, " ")[[1]]
aux3 <- aux2[grep("Q", aux2)]
F1 <- aux3[1:length(aux3)]
# Alpha
alpha_F1 <- psych::alpha(X[,F1]); alpha_F1$total[1]## raw_alpha
## 0.96564322 factors
Number of factors set to 2. After an initial elimination of items with loadings less than 0.32 and/or with cross loadings, approximate simple structure was achieved with 41 items in the first factor and 9 items in the second one. To further reduce the number of items, we excluded items with lower factor loadings, lower communality, and higher complexity. This led to 15 items in F1 and 5 in F2. Of the 15 items in F1, 5 were words of affirmation, 4 were acts of service, 3 were quality time, and 3 were gifts. This led to the following 2-factor solution :
nfactors <- 2
remove <- c("Q1", "Q2", "Q4", "Q5", "Q6", "Q7", "Q8", "Q9", "Q11", "Q14", "Q16", "Q17", "Q18", "Q19", "Q20", "Q103", "Q104", "Q105", "Q106", "Q122", "Q146", "Q147", "Q148", "Q150", "Q151", "Q152", "Q161", "Q162", "Q164", "Q166", "Q167", "Q168", "Q181", "Q182", "Q183", "Q185", "Q186", "Q187", "Q198", "Q199", "Q202", "Q203", "Q205", "Q206")
#remove <- NA
cols <- setdiff(names(X), remove)
cor2 <- cor(X[,cols])
fa1 <- fa(r = cor2, n.obs = nrow(X[,cols]), nfactors = nfactors,
fm = "pa", max.iter = 100, rotate = "oblimin")
print(fa1$loadings, cutoff = 0.319)##
## Loadings:
## PA1 PA2
## Q3 0.820
## Q10 0.783
## Q12 0.770
## Q13 0.847
## Q15 0.844
## Q100 0.785
## Q101 0.717
## Q102 0.788
## Q107 0.618
## Q145 0.836
## Q149 0.778
## Q163 0.657
## Q165 0.627
## Q169 0.756
## Q179 0.792
## Q180 0.926
## Q184 0.690
## Q200 0.730
## Q201 0.726
## Q204 0.694
##
## PA1 PA2
## SS loadings 8.526 3.231
## Proportion Var 0.426 0.162
## Cumulative Var 0.426 0.588fa1$RMSEA## RMSEA lower upper confidence
## 0.07682237 0.06903262 0.08505370 0.90000000fa1$TLI## [1] 0.9215219fa1$Vaccounted## PA1 PA2
## SS loadings 8.6930667 3.3984906
## Proportion Var 0.4346533 0.1699245
## Cumulative Var 0.4346533 0.6045779
## Proportion Explained 0.7189369 0.2810631
## Cumulative Proportion 0.7189369 1.0000000# Factors
factors2 <- efa_to_cfa(fa1)
aux2 <- strsplit(factors2, " ")[[1]]
aux3 <- aux2[grep("Q", aux2)]
w <- grep("PA", aux3)
F1 <- aux3[1:w]
F2 <- aux3[(w+1):length(aux3)]
F1[length(F1)] <- substr(F1[length(F1)], 1, 4)
# Alphas
alpha_F1 <- psych::alpha(X[,F1]); alpha_F1$total[1]## raw_alpha
## 0.953899alpha_F2 <- psych::alpha(X[,F2]); alpha_F2$total[1]## raw_alpha
## 0.8986201# Correlation between factors
cor(apply(X[,F1], 1, mean), apply(X[,F2], 1, mean))## [1] 0.61959513 factors
When the number of factors was set to 3, the initial loadings suggested a 2-factor model. Specifically, there were no loadings above 0.32, only two cross-loadings below 0.4. However, working on the 4-factor solution led to three factors after an initial elimination of items with loadings less than 0.32 and/or with cross loadings. We then worked on a 3-factor solution starting with the 49 items from this 4-factor solution. We further decreased the number of items with lower loadings, cross-loadings, low communality, and conceptually inconsistent. Approximate simple structure was achieved with 11 items in the first factor, 11 items in the second one, and 10 items in the third one. To further reduce the number of items, we excluded items with lower factor loadings, lower communality, and higher complexity. This led to 7 items in F1 and 6 in F2 and 7 items in F3. Of the 7 items in F1, 5 were words of affirmation and 2 were quality time. Of the 7 items in F3, 5 were acts of service and 2 were gifts. This led to the following 3-factor solution:
nfactors <- 3
remove <- c("Q1", "Q5", "Q6", "Q107", "Q7", "Q9", "Q10", "Q11", "Q12", "Q13", "Q14", "Q15", "Q16", "Q17", "Q19", "Q20", "Q102", "Q103", "Q104", "Q106", "Q122", "Q146", "Q147", "Q148", "Q150", "Q151", "Q152", "Q161", "Q162", "Q165", "Q166", "Q167", "Q168", "Q181", "Q182", "Q183", "Q185", "Q187", "Q201", "Q202", "Q203", "Q204", "Q205", "Q206")
# remove <- NA
cols <- setdiff(names(X), remove)
cor2 <- cor(X[,cols])
fa1 <- fa(r = cor2, n.obs = nrow(X[,cols]), nfactors = nfactors,
fm = "pa", max.iter = 100, rotate = "oblimin")
print(fa1$loadings, cutoff = 0.319)##
## Loadings:
## PA1 PA2 PA3
## Q2 0.846
## Q3 0.627
## Q4 0.787
## Q8 0.724
## Q18 0.477
## Q100 0.770
## Q101 0.488
## Q105 0.565
## Q145 0.852
## Q149 0.650
## Q163 0.568
## Q164 0.473
## Q169 0.755
## Q179 0.767
## Q180 0.919
## Q184 0.623
## Q186 0.747
## Q198 0.788
## Q199 0.798
## Q200 0.948
##
## PA1 PA2 PA3
## SS loadings 3.828 3.787 3.340
## Proportion Var 0.191 0.189 0.167
## Cumulative Var 0.191 0.381 0.548fa1$RMSEA## RMSEA lower upper confidence
## 0.07101185 0.06254932 0.07990382 0.90000000fa1$TLI## [1] 0.9326798fa1$Vaccounted## PA1 PA2 PA3
## SS loadings 4.4887954 4.2424761 3.9301884
## Proportion Var 0.2244398 0.2121238 0.1965094
## Cumulative Var 0.2244398 0.4365636 0.6330730
## Proportion Explained 0.3545243 0.3350701 0.3104056
## Cumulative Proportion 0.3545243 0.6895944 1.0000000# Factors
factors3 <- efa_to_cfa(fa1)
aux2 <- strsplit(factors3, " ")[[1]]
aux3 <- aux2[grep("Q", aux2)]
w <- grep("PA", aux3)
F1 <- aux3[1:w[1]]
F2 <- aux3[(w[1]+1):w[2]]
F3 <- aux3[(w[2]+1):length(aux3)]
F1[length(F1)] <- substr(F1[length(F1)], 1, 4)
F2[length(F2)] <- substr(F2[length(F2)], 1, 4)
# Alphas
alpha_F1 <- psych::alpha(X[,F1]); alpha_F1$total[1]## raw_alpha
## 0.9226617alpha_F2 <- psych::alpha(X[,F2]); alpha_F2$total[1]## raw_alpha
## 0.9147622alpha_F3 <- psych::alpha(X[,F3]); alpha_F3$total[1]## raw_alpha
## 0.897103# Correlation between factors
F1m <- apply(X[,F1], 1, mean)
F2m <- apply(X[,F2], 1, mean)
F3m <- apply(X[,F3], 1, mean)
factors <- data.frame(F1m, F2m, F3m)
cor(factors)## F1m F2m F3m
## F1m 1.0000000 0.672156 0.7543078
## F2m 0.6721560 1.000000 0.5895530
## F3m 0.7543078 0.589553 1.00000004 factors
The initial work on a 4-factor solution suggested a 3-factor model. However, working on the 5-factor solution led to four factors after an initial elimination of items with loadings less than 0.32 and/or with cross loadings (the fifth factor had only 3 items and has highly correlated with one of the four factors). We then worked on a 4-factor solution starting with the 33 items from this 5-factor solution. Approximate simple structure was achieved with 10 items in the first factor, 9 items in the second one, 8 items in the third one, and 6 items in the fourth one. To further reduce the number of items, we excluded items with lower factor loadings, lower communality, and higher complexity. This led to 6 items in F1, 4 in F2, 4 items in F3, and 4 items in F4. Of the 6 items in F1, 4 were words of affirmation and 2 were quality time. The other factors mapped to the remaining love languages. This led to the following 4-factor solution:
nfactors <- 4
remove <- c("Q1", "Q5", "Q6", "Q7", "Q8", "Q9", "Q10", "Q11", "Q12", "Q13", "Q14", "Q15", "Q16", "Q101", "Q102", "Q103", "Q104", "Q105", "Q106", "Q122", "Q146", "Q147", "Q148", "Q150", "Q151", "Q152", "Q161", "Q162", "Q165", "Q166", "Q167", "Q168", "Q181", "Q183", "Q185", "Q187", "Q199", "Q201", "Q203", "Q204", "Q205", "Q206", "Q17", "Q3", "Q107", "Q182")
#remove <- NA
cols <- setdiff(names(X), remove)
cor2 <- cor(X[,cols])
fa1 <- fa(r = cor2, n.obs = nrow(X[,cols]), nfactors = nfactors,
fm = "pa", max.iter = 100, rotate = "oblimin")
print(fa1$loadings, cutoff = 0.319)##
## Loadings:
## PA1 PA2 PA4 PA3
## Q2 0.781
## Q4 0.685
## Q18 0.555
## Q19 0.750
## Q20 0.788
## Q100 0.707
## Q145 0.772
## Q149 0.617
## Q163 0.637
## Q164 0.506
## Q169 0.758
## Q179 0.735
## Q180 0.937
## Q184 0.580
## Q186 0.728
## Q198 0.619
## Q200 0.920
## Q202 0.593
##
## PA1 PA2 PA4 PA3
## SS loadings 2.904 2.722 2.146 1.934
## Proportion Var 0.161 0.151 0.119 0.107
## Cumulative Var 0.161 0.313 0.432 0.539fa1$RMSEA## RMSEA lower upper confidence
## 0.04837669 0.03628481 0.06044600 0.90000000fa1$TLI## [1] 0.971023fa1$Vaccounted## PA1 PA2 PA4 PA3
## SS loadings 3.5786543 3.1462183 2.7001419 2.5238021
## Proportion Var 0.1988141 0.1747899 0.1500079 0.1402112
## Cumulative Var 0.1988141 0.3736040 0.5236119 0.6638231
## Proportion Explained 0.2994986 0.2633079 0.2259757 0.2112177
## Cumulative Proportion 0.2994986 0.5628066 0.7887823 1.0000000# Factors
factors4_2 <- efa_to_cfa(fa1)
aux2 <- strsplit(factors4_2, " ")[[1]]
aux3 <- aux2[grep("Q", aux2)]
w <- grep("PA", aux3)
F1 <- aux3[1:w[1]]
F2 <- aux3[(w[1]+1):w[2]]
F3 <- aux3[(w[2]+1):w[3]]
F4 <- aux3[(w[3]+1):length(aux3)]
F1[length(F1)] <- substr(F1[length(F1)], 1, 4)
F2[length(F2)] <- substr(F2[length(F2)], 1, 4)
F3[length(F3)] <- substr(F3[length(F3)], 1, 4)
# Alphas
alpha_F1 <- psych::alpha(X[,F1]); alpha_F1$total[1]## raw_alpha
## 0.907973alpha_F2 <- psych::alpha(X[,F2]); alpha_F2$total[1]## raw_alpha
## 0.9073911alpha_F3 <- psych::alpha(X[,F3]); alpha_F3$total[1]## raw_alpha
## 0.8737413alpha_F4 <- psych::alpha(X[,F4]); alpha_F4$total[1]## raw_alpha
## 0.8686325# Correlation between factors
F1m <- apply(X[,F1], 1, mean)
F2m <- apply(X[,F2], 1, mean)
F3m <- apply(X[,F3], 1, mean)
F4m <- apply(X[,F4], 1, mean)
factors <- data.frame(F1m, F2m, F3m, F4m)
round(cor(factors),2)## F1m F2m F3m F4m
## F1m 1.00 0.64 0.70 0.69
## F2m 0.64 1.00 0.48 0.64
## F3m 0.70 0.48 1.00 0.68
## F4m 0.69 0.64 0.68 1.005 factors (25 items)
The initial work on a 5-factor solution suggested a 4-factor model. However, working on the 6-factor solution led to five factors after an initial elimination of items with loadings less than 0.32 and/or with cross loadings. To further reduce the number of items, we excluded items with lower factor loadings, lower communality, and higher complexity. This led to two solutions, one with 25 items (5 in each factor) and one with 20 items (4 in each factor).
nfactors <- 5
remove <- c("Q5", "Q6", "Q7", "Q8", "Q9", "Q10", "Q11", "Q12", "Q13", "Q101", "Q102", "Q103", "Q104", "Q105", "Q122", "Q146", "Q147", "Q148", "Q150", "Q151", "Q152", "Q161", "Q162", "Q163", "Q165", "Q166", "Q167", "Q181", "Q183", "Q201", "Q203", "Q204", "Q205", "Q206", "Q182", "Q185", "Q184", "Q168", "Q182", "Q185", "Q3") # solution 1 from 6 factors
#remove <- NA
cols <- setdiff(names(X), remove)
cor2 <- cor(X[,cols])
fa1 <- fa(r = cor2, n.obs = nrow(X[,cols]), nfactors = nfactors,
fm = "pa", max.iter = 100, rotate = "oblimin")
print(fa1$loadings, cutoff = 0.319)##
## Loadings:
## PA1 PA2 PA4 PA3 PA5
## Q1 0.407
## Q2 0.753
## Q4 0.685
## Q14 0.780
## Q15 0.636
## Q16 0.649
## Q17 0.554
## Q18 0.550
## Q19 0.711
## Q20 0.692
## Q100 0.665
## Q106 0.694
## Q107 0.586
## Q145 0.805
## Q149 0.573
## Q164 0.415
## Q169 0.778
## Q179 0.685
## Q180 0.847
## Q186 0.752
## Q187 0.692
## Q198 0.675
## Q199 0.676
## Q200 0.742
## Q202 0.620
##
## PA1 PA2 PA4 PA3 PA5
## SS loadings 2.843 2.919 2.251 2.139 2.166
## Proportion Var 0.114 0.117 0.090 0.086 0.087
## Cumulative Var 0.114 0.230 0.321 0.406 0.493fa1$RMSEA## RMSEA lower upper confidence
## 0.05696234 0.04931259 0.06491407 0.90000000fa1$TLI## [1] 0.9451039fa1$Vaccounted## PA1 PA2 PA4 PA3 PA5
## SS loadings 3.7237489 3.5215641 3.1251538 3.0237646 2.9108697
## Proportion Var 0.1489500 0.1408626 0.1250062 0.1209506 0.1164348
## Cumulative Var 0.1489500 0.2898125 0.4148187 0.5357693 0.6522040
## Proportion Explained 0.2283794 0.2159793 0.1916672 0.1854490 0.1785251
## Cumulative Proportion 0.2283794 0.4443587 0.6360259 0.8214749 1.0000000# Factors
factors5_1 <- efa_to_cfa(fa1)
aux2 <- strsplit(factors5_1, " ")[[1]]
aux3 <- aux2[grep("Q", aux2)]
w <- grep("PA", aux3)
F1 <- aux3[1:w[1]]
F2 <- aux3[(w[1]+1):w[2]]
F3 <- aux3[(w[2]+1):w[3]]
F4 <- aux3[(w[3]+1):w[4]]
F5 <- aux3[(w[4]+1):length(aux3)]
F1[length(F1)] <- substr(F1[length(F1)], 1, 4)
F2[length(F2)] <- substr(F2[length(F2)], 1, 4)
F3[length(F3)] <- substr(F3[length(F3)], 1, 4)
F4[length(F4)] <- substr(F4[length(F4)], 1, 4)
# Alphas
alpha_F1 <- psych::alpha(X[,F1]); alpha_F1$total[1]## raw_alpha
## 0.9240188alpha_F2 <- psych::alpha(X[,F2]); alpha_F2$total[1]## raw_alpha
## 0.8986201alpha_F3 <- psych::alpha(X[,F3]); alpha_F3$total[1]## raw_alpha
## 0.8920414alpha_F4 <- psych::alpha(X[,F4]); alpha_F4$total[1]## raw_alpha
## 0.8802723alpha_F5 <- psych::alpha(X[,F5]); alpha_F5$total[1]## raw_alpha
## 0.8718099# Correlation between factors
F1m <- apply(X[,F1], 1, mean)
F2m <- apply(X[,F2], 1, mean)
F3m <- apply(X[,F3], 1, mean)
F4m <- apply(X[,F4], 1, mean)
F5m <- apply(X[,F5], 1, mean)
factors <- data.frame(F1m, F2m, F3m, F4m, F5m)
round(cor(factors), 2)## F1m F2m F3m F4m F5m
## F1m 1.00 0.62 0.74 0.68 0.66
## F2m 0.62 1.00 0.64 0.64 0.50
## F3m 0.74 0.64 1.00 0.69 0.69
## F4m 0.68 0.64 0.69 1.00 0.71
## F5m 0.66 0.50 0.69 0.71 1.00nfactors <- 5
remove <- c("Q1", "Q3", "Q5", "Q6", "Q7", "Q8", "Q9", "Q10", "Q11", "Q12", "Q13", "Q101", "Q102", "Q103", "Q104", "Q105", "Q107", "Q122", "Q146", "Q148", "Q150", "Q151", "Q152", "Q161", "Q162", "Q163", "Q165", "Q166", "Q167", "Q181", "Q182", "Q183", "Q185", "Q199", "Q201", "Q203", "Q204", "Q205", "Q206", "Q168", "Q147", "Q17", "Q164", "Q187") # solution 2 from 6 factors
#remove <- NA
cols <- setdiff(names(X), remove)
cor2 <- cor(X[,cols])
fa1 <- fa(r = cor2, n.obs = nrow(X[,cols]), nfactors = nfactors,
fm = "pa", max.iter = 100, rotate = "oblimin")
print(fa1$loadings, cutoff = 0.319)##
## Loadings:
## PA2 PA5 PA1 PA4 PA3
## Q2 0.798
## Q4 0.708
## Q14 0.721
## Q15 0.613
## Q16 0.648
## Q18 0.471
## Q19 0.691
## Q20 0.699
## Q100 0.685
## Q106 0.753
## Q145 0.787
## Q149 0.590
## Q169 0.737
## Q179 0.739
## Q180 0.903
## Q184 0.551
## Q186 0.747
## Q198 0.578
## Q200 0.760
## Q202 0.552
##
## PA2 PA5 PA1 PA4 PA3
## SS loadings 2.640 2.163 1.953 1.989 1.615
## Proportion Var 0.132 0.108 0.098 0.099 0.081
## Cumulative Var 0.132 0.240 0.338 0.437 0.518fa1$RMSEA## RMSEA lower upper confidence
## 0.04765279 0.03631801 0.05897779 0.90000000fa1$TLI## [1] 0.9688624fa1$Vaccounted## PA2 PA5 PA1 PA4 PA3
## SS loadings 3.1128909 2.8084939 2.6650009 2.6299627 2.3070079
## Proportion Var 0.1556445 0.1404247 0.1332500 0.1314981 0.1153504
## Cumulative Var 0.1556445 0.2960692 0.4293193 0.5608174 0.6761678
## Proportion Explained 0.2301863 0.2076773 0.1970665 0.1944756 0.1705943
## Cumulative Proportion 0.2301863 0.4378635 0.6349301 0.8294057 1.0000000# Factors
factors5_2 <- efa_to_cfa(fa1)
aux2 <- strsplit(factors5_2, " ")[[1]]
aux3 <- aux2[grep("Q", aux2)]
w <- grep("PA", aux3)
F1 <- aux3[1:w[1]]
F2 <- aux3[(w[1]+1):w[2]]
F3 <- aux3[(w[2]+1):w[3]]
F4 <- aux3[(w[3]+1):w[4]]
F5 <- aux3[(w[4]+1):length(aux3)]
F1[length(F1)] <- substr(F1[length(F1)], 1, 4)
F2[length(F2)] <- substr(F2[length(F2)], 1, 4)
F3[length(F3)] <- substr(F3[length(F3)], 1, 4)
F4[length(F4)] <- substr(F4[length(F4)], 1, 4)
# Alphas
alpha_F1 <- psych::alpha(X[,F1]); alpha_F1$total[1]## raw_alpha
## 0.9073911alpha_F2 <- psych::alpha(X[,F2]); alpha_F2$total[1]## raw_alpha
## 0.8737413alpha_F3 <- psych::alpha(X[,F3]); alpha_F3$total[1]## raw_alpha
## 0.8870889alpha_F4 <- psych::alpha(X[,F4]); alpha_F4$total[1]## raw_alpha
## 0.8765027alpha_F5 <- psych::alpha(X[,F5]); alpha_F5$total[1]## raw_alpha
## 0.8686325# Correlation between factors
F1m <- apply(X[,F1], 1, mean)
F2m <- apply(X[,F2], 1, mean)
F3m <- apply(X[,F3], 1, mean)
F4m <- apply(X[,F4], 1, mean)
F5m <- apply(X[,F5], 1, mean)
factors <- data.frame(F1m, F2m, F3m, F4m, F5m)
round(cor(factors), 2)## F1m F2m F3m F4m F5m
## F1m 1.00 0.48 0.60 0.60 0.64
## F2m 0.48 1.00 0.68 0.65 0.68
## F3m 0.60 0.68 1.00 0.70 0.67
## F4m 0.60 0.65 0.70 1.00 0.66
## F5m 0.64 0.68 0.67 0.66 1.00CFA
The six models from the EFA will be tested with a CFA.
Checking multivariate normality on the CFA dataset
Descriptive stats:
mvn2 <- mvn(b2)
mvn2_df <- data.frame(mvn2$Descriptives)
write.csv(mvn2_df, "mvn2.csv")
median(mvn2_df$Skew)## [1] -1.198413min(mvn2_df$Skew)## [1] -2.086128max(mvn2_df$Skew)## [1] 0.07729608median(mvn2_df$Kurtosis)## [1] 0.665868min(mvn2_df$Kurtosis)## [1] -1.474451max(mvn2_df$Kurtosis)## [1] 4.025911Number of potential outliers:
cutoff <- qchisq(1-0.001, ncol(b2))
mahal <- mahalanobis(b2, colMeans(b2), cov(b2))
s <- summary(mahal < cutoff); s## Mode FALSE TRUE
## logical 51 297Henze-Zirkler test for multivariate normality:
mvn2$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 3.474456 0 NOThe data is not multivariate normal (use a robust method, “MLR”).
Create table to store measures of model fit
models <- data.frame(model = NA, chisq = NA, df = NA, p = NA, CFI = NA, TLI = NA, ecvi = NA, srmr = NA, RMSEA = NA, RMSEA_lower = NA, RMSEA_upper = NA)One-factor model
#factors1
model <- '
F1 =~ Q5 + Q7 + Q10 + Q13 + Q14 + Q19 + Q101 + Q104 + Q149 + Q162 + Q164 + Q165 + Q168 + Q181 + Q186 + Q200 + Q201 + Q203 + Q204 + Q205'
fit <- cfa(model, estimator = "MLR", data = b2)
fm <- fitMeasures(fit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "cfi.robust", "tli.robust", "ecvi", "srmr", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust"))
models[1,1] <- "1 factor"
models[1,2:ncol(models)] <- fm
summary(fit,
#fit.measures = TRUE,
standardized = TRUE
#,rsq = TRUE
)## lavaan 0.6-12 ended normally after 30 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 40
##
## Number of observations 348
##
## Model Test User Model:
## Standard Robust
## Test Statistic 710.799 400.159
## Degrees of freedom 170 170
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.776
## Yuan-Bentler correction (Mplus variant)
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~
## Q5 1.000 0.770 0.747
## Q7 1.107 0.064 17.182 0.000 0.852 0.729
## Q10 1.086 0.070 15.539 0.000 0.835 0.801
## Q13 1.028 0.082 12.505 0.000 0.791 0.769
## Q14 1.026 0.072 14.317 0.000 0.790 0.817
## Q19 1.190 0.096 12.335 0.000 0.916 0.740
## Q101 1.038 0.093 11.130 0.000 0.799 0.718
## Q104 1.133 0.080 14.127 0.000 0.872 0.830
## Q149 1.115 0.072 15.395 0.000 0.858 0.803
## Q162 1.301 0.103 12.641 0.000 1.001 0.784
## Q164 1.137 0.080 14.227 0.000 0.875 0.784
## Q165 1.027 0.073 13.995 0.000 0.790 0.809
## Q168 1.141 0.067 17.010 0.000 0.878 0.731
## Q181 1.102 0.064 17.126 0.000 0.848 0.780
## Q186 1.102 0.066 16.627 0.000 0.848 0.791
## Q200 1.134 0.071 16.048 0.000 0.873 0.820
## Q201 1.199 0.077 15.675 0.000 0.923 0.838
## Q203 1.247 0.086 14.567 0.000 0.960 0.861
## Q204 1.213 0.081 14.990 0.000 0.934 0.836
## Q205 1.124 0.082 13.761 0.000 0.865 0.798
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Q5 0.468 0.056 8.428 0.000 0.468 0.442
## .Q7 0.638 0.085 7.525 0.000 0.638 0.468
## .Q10 0.390 0.052 7.570 0.000 0.390 0.359
## .Q13 0.434 0.053 8.201 0.000 0.434 0.409
## .Q14 0.311 0.040 7.721 0.000 0.311 0.333
## .Q19 0.693 0.069 10.007 0.000 0.693 0.452
## .Q101 0.600 0.076 7.907 0.000 0.600 0.484
## .Q104 0.344 0.040 8.596 0.000 0.344 0.311
## .Q149 0.404 0.054 7.469 0.000 0.404 0.355
## .Q162 0.630 0.087 7.275 0.000 0.630 0.386
## .Q164 0.481 0.060 7.986 0.000 0.481 0.386
## .Q165 0.330 0.057 5.810 0.000 0.330 0.346
## .Q168 0.672 0.089 7.530 0.000 0.672 0.466
## .Q181 0.464 0.081 5.700 0.000 0.464 0.392
## .Q186 0.431 0.056 7.710 0.000 0.431 0.375
## .Q200 0.373 0.053 7.083 0.000 0.373 0.328
## .Q201 0.362 0.043 8.443 0.000 0.362 0.298
## .Q203 0.322 0.034 9.515 0.000 0.322 0.259
## .Q204 0.376 0.038 9.825 0.000 0.376 0.301
## .Q205 0.426 0.043 9.860 0.000 0.426 0.363
## F1 0.592 0.095 6.215 0.000 1.000 1.000fm## chisq.scaled df.scaled pvalue.scaled
## 400.159 170.000 0.000
## cfi.robust tli.robust ecvi
## 0.932 0.924 2.272
## srmr rmsea.robust rmsea.ci.lower.robust
## 0.039 0.083 0.073
## rmsea.ci.upper.robust
## 0.094res1 <- residuals(fit, type = "cor")$cov
res1[upper.tri(res1, diag = T)] <- NA
v1 <- as.vector(res1)
v2 <- v1[!is.na(v1)]
qqPlot(v2, id = F)
Two-factor model
#factors2
model <- '
F1 =~ Q3 + Q10 + Q12 + Q13 + Q15 + Q100 + Q101 + Q102 + Q149 + Q163 + Q165 + Q184 + Q200 + Q201 + Q204
F2 =~ Q107 + Q145 + Q169 + Q179 + Q180'
fit <- cfa(model, estimator = "MLR", data = b2)
fm <- fitMeasures(fit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "cfi.robust", "tli.robust", "ecvi", "srmr", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust"))
models[2,1] <- "2 factors"
models[2,2:ncol(models)] <- fm
summary(fit,
#fit.measures = TRUE,
standardized = TRUE,
rsq = TRUE
)## lavaan 0.6-12 ended normally after 39 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 41
##
## Number of observations 348
##
## Model Test User Model:
## Standard Robust
## Test Statistic 533.165 348.568
## Degrees of freedom 169 169
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.530
## Yuan-Bentler correction (Mplus variant)
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~
## Q3 1.000 0.685 0.784
## Q10 1.231 0.081 15.269 0.000 0.843 0.808
## Q12 1.200 0.058 20.552 0.000 0.822 0.843
## Q13 1.179 0.085 13.885 0.000 0.808 0.785
## Q15 1.157 0.073 15.820 0.000 0.792 0.802
## Q100 1.180 0.082 14.456 0.000 0.808 0.814
## Q101 1.129 0.084 13.483 0.000 0.773 0.695
## Q102 1.123 0.086 12.988 0.000 0.769 0.757
## Q149 1.311 0.084 15.673 0.000 0.898 0.840
## Q163 1.097 0.099 11.053 0.000 0.751 0.747
## Q165 1.167 0.073 15.920 0.000 0.799 0.818
## Q184 1.314 0.087 15.023 0.000 0.900 0.768
## Q200 1.227 0.096 12.816 0.000 0.840 0.789
## Q201 1.349 0.087 15.433 0.000 0.924 0.839
## Q204 1.333 0.092 14.487 0.000 0.913 0.817
## F2 =~
## Q107 1.000 0.810 0.652
## Q145 1.398 0.113 12.376 0.000 1.132 0.804
## Q169 1.349 0.103 13.120 0.000 1.092 0.858
## Q179 1.380 0.120 11.465 0.000 1.118 0.849
## Q180 1.441 0.114 12.646 0.000 1.167 0.869
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 ~~
## F2 0.387 0.065 5.993 0.000 0.698 0.698
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Q3 0.295 0.033 8.817 0.000 0.295 0.386
## .Q10 0.377 0.059 6.437 0.000 0.377 0.347
## .Q12 0.275 0.041 6.651 0.000 0.275 0.289
## .Q13 0.407 0.060 6.758 0.000 0.407 0.384
## .Q15 0.349 0.040 8.812 0.000 0.349 0.357
## .Q100 0.333 0.041 8.162 0.000 0.333 0.338
## .Q101 0.640 0.083 7.708 0.000 0.640 0.517
## .Q102 0.442 0.055 7.984 0.000 0.442 0.427
## .Q149 0.335 0.043 7.792 0.000 0.335 0.294
## .Q163 0.446 0.061 7.353 0.000 0.446 0.441
## .Q165 0.316 0.051 6.176 0.000 0.316 0.331
## .Q184 0.564 0.078 7.195 0.000 0.564 0.411
## .Q200 0.429 0.064 6.699 0.000 0.429 0.378
## .Q201 0.360 0.037 9.746 0.000 0.360 0.296
## .Q204 0.414 0.048 8.581 0.000 0.414 0.332
## .Q107 0.886 0.086 10.326 0.000 0.886 0.575
## .Q145 0.701 0.105 6.676 0.000 0.701 0.354
## .Q169 0.428 0.052 8.235 0.000 0.428 0.264
## .Q179 0.484 0.073 6.655 0.000 0.484 0.279
## .Q180 0.440 0.058 7.538 0.000 0.440 0.244
## F1 0.469 0.079 5.961 0.000 1.000 1.000
## F2 0.656 0.111 5.926 0.000 1.000 1.000
##
## R-Square:
## Estimate
## Q3 0.614
## Q10 0.653
## Q12 0.711
## Q13 0.616
## Q15 0.643
## Q100 0.662
## Q101 0.483
## Q102 0.573
## Q149 0.706
## Q163 0.559
## Q165 0.669
## Q184 0.589
## Q200 0.622
## Q201 0.704
## Q204 0.668
## Q107 0.425
## Q145 0.646
## Q169 0.736
## Q179 0.721
## Q180 0.756fm## chisq.scaled df.scaled pvalue.scaled
## 348.568 169.000 0.000
## cfi.robust tli.robust ecvi
## 0.951 0.945 1.768
## srmr rmsea.robust rmsea.ci.lower.robust
## 0.037 0.068 0.058
## rmsea.ci.upper.robust
## 0.079res1 <- residuals(fit, type = "cor")$cov
res1[upper.tri(res1, diag = T)] <- NA
v1 <- as.vector(res1)
v2 <- v1[!is.na(v1)]
qqPlot(v2, id = F)
Three-factor model
#factors3
model <- '
F1 =~ Q163 + Q164 + Q184 + Q186 + Q198 + Q199 + Q200
F2 =~ Q8 + Q105 + Q145 + Q169 + Q179 + Q180
F3 =~ Q2 + Q3 + Q4 + Q18 + Q100 + Q101 + Q149'
fit <- cfa(model, estimator = "MLR", data = b2)
fm <- fitMeasures(fit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "cfi.robust", "tli.robust", "ecvi", "srmr", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust"))
models[3,1] <- "3 factors"
models[3,2:ncol(models)] <- fm
summary(fit,
#fit.measures = TRUE,
standardized = TRUE
#,rsq = TRUE
)## lavaan 0.6-12 ended normally after 43 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 43
##
## Number of observations 348
##
## Model Test User Model:
## Standard Robust
## Test Statistic 540.491 352.037
## Degrees of freedom 167 167
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.535
## Yuan-Bentler correction (Mplus variant)
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~
## Q163 1.000 0.753 0.749
## Q164 1.163 0.076 15.354 0.000 0.876 0.784
## Q184 1.195 0.091 13.090 0.000 0.901 0.768
## Q186 1.160 0.084 13.806 0.000 0.874 0.815
## Q198 1.330 0.118 11.247 0.000 1.002 0.866
## Q199 1.323 0.111 11.975 0.000 0.997 0.874
## Q200 1.177 0.104 11.331 0.000 0.886 0.832
## F2 =~
## Q8 1.000 0.800 0.776
## Q105 0.935 0.068 13.778 0.000 0.747 0.694
## Q145 1.375 0.097 14.156 0.000 1.099 0.781
## Q169 1.369 0.084 16.379 0.000 1.095 0.860
## Q179 1.397 0.091 15.278 0.000 1.117 0.849
## Q180 1.452 0.091 15.969 0.000 1.161 0.865
## F3 =~
## Q2 1.000 0.840 0.744
## Q3 0.819 0.061 13.412 0.000 0.688 0.787
## Q4 0.939 0.058 16.155 0.000 0.788 0.747
## Q18 0.995 0.076 13.045 0.000 0.836 0.619
## Q100 0.989 0.068 14.544 0.000 0.830 0.836
## Q101 0.955 0.077 12.479 0.000 0.802 0.721
## Q149 1.098 0.070 15.710 0.000 0.922 0.863
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 ~~
## F2 0.434 0.073 5.936 0.000 0.720 0.720
## F3 0.569 0.076 7.525 0.000 0.899 0.899
## F2 ~~
## F3 0.465 0.067 6.928 0.000 0.692 0.692
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Q163 0.443 0.064 6.958 0.000 0.443 0.438
## .Q164 0.480 0.057 8.430 0.000 0.480 0.385
## .Q184 0.563 0.074 7.570 0.000 0.563 0.410
## .Q186 0.386 0.052 7.454 0.000 0.386 0.336
## .Q198 0.336 0.047 7.144 0.000 0.336 0.251
## .Q199 0.306 0.041 7.504 0.000 0.306 0.235
## .Q200 0.349 0.048 7.287 0.000 0.349 0.308
## .Q8 0.422 0.048 8.872 0.000 0.422 0.398
## .Q105 0.602 0.074 8.096 0.000 0.602 0.519
## .Q145 0.773 0.109 7.086 0.000 0.773 0.390
## .Q169 0.422 0.053 8.010 0.000 0.422 0.260
## .Q179 0.484 0.070 6.918 0.000 0.484 0.280
## .Q180 0.454 0.062 7.305 0.000 0.454 0.252
## .Q2 0.568 0.074 7.657 0.000 0.568 0.446
## .Q3 0.291 0.034 8.433 0.000 0.291 0.381
## .Q4 0.492 0.063 7.780 0.000 0.492 0.442
## .Q18 1.126 0.090 12.534 0.000 1.126 0.617
## .Q100 0.297 0.041 7.308 0.000 0.297 0.301
## .Q101 0.595 0.080 7.455 0.000 0.595 0.481
## .Q149 0.291 0.038 7.583 0.000 0.291 0.255
## F1 0.567 0.099 5.723 0.000 1.000 1.000
## F2 0.639 0.091 6.996 0.000 1.000 1.000
## F3 0.705 0.090 7.808 0.000 1.000 1.000fm## chisq.scaled df.scaled pvalue.scaled
## 352.037 167.000 0.000
## cfi.robust tli.robust ecvi
## 0.947 0.939 1.800
## srmr rmsea.robust rmsea.ci.lower.robust
## 0.049 0.070 0.060
## rmsea.ci.upper.robust
## 0.080res1 <- residuals(fit, type = "cor")$cov
res1[upper.tri(res1, diag = T)] <- NA
v1 <- as.vector(res1)
v2 <- v1[!is.na(v1)]
qqPlot(v2, id = F)
Four-factor model
#factors4_2
model <- '
F1 =~ Q163 + Q164 + Q184 + Q186 + Q198 + Q200
F2 =~ Q145 + Q169 + Q179 + Q180
F3 =~ Q2 + Q4 + Q100 + Q149
F4 =~ Q18 + Q19 + Q20 + Q202'
fit <- cfa(model, estimator = "MLR", data = b2)
fm <- fitMeasures(fit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "cfi.robust", "tli.robust", "ecvi", "srmr", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust"))
models[4,1] <- "4 factors"
models[4,2:ncol(models)] <- fm
summary(fit,
#fit.measures = TRUE,
standardized = TRUE
#,rsq = TRUE
)## lavaan 0.6-12 ended normally after 49 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 42
##
## Number of observations 348
##
## Model Test User Model:
## Standard Robust
## Test Statistic 370.068 252.106
## Degrees of freedom 129 129
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.468
## Yuan-Bentler correction (Mplus variant)
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~
## Q163 1.000 0.771 0.767
## Q164 1.155 0.070 16.430 0.000 0.891 0.798
## Q184 1.182 0.085 13.852 0.000 0.912 0.778
## Q186 1.144 0.079 14.494 0.000 0.882 0.823
## Q198 1.270 0.106 11.937 0.000 0.979 0.846
## Q200 1.126 0.096 11.786 0.000 0.869 0.816
## F2 =~
## Q145 1.000 1.128 0.801
## Q169 0.958 0.050 19.055 0.000 1.080 0.849
## Q179 1.007 0.053 19.082 0.000 1.136 0.863
## Q180 1.034 0.049 21.265 0.000 1.167 0.869
## F3 =~
## Q2 1.000 0.839 0.744
## Q4 0.957 0.059 16.201 0.000 0.803 0.761
## Q100 0.997 0.069 14.394 0.000 0.836 0.842
## Q149 1.113 0.068 16.324 0.000 0.933 0.874
## F4 =~
## Q18 1.000 0.986 0.730
## Q19 1.061 0.079 13.417 0.000 1.045 0.845
## Q20 1.112 0.073 15.156 0.000 1.096 0.847
## Q202 0.990 0.070 14.232 0.000 0.975 0.758
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 ~~
## F2 0.611 0.086 7.111 0.000 0.702 0.702
## F3 0.576 0.076 7.622 0.000 0.890 0.890
## F4 0.604 0.072 8.369 0.000 0.795 0.795
## F2 ~~
## F3 0.610 0.080 7.640 0.000 0.645 0.645
## F4 0.714 0.099 7.239 0.000 0.642 0.642
## F3 ~~
## F4 0.666 0.075 8.910 0.000 0.806 0.806
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Q163 0.416 0.060 6.900 0.000 0.416 0.411
## .Q164 0.453 0.052 8.779 0.000 0.453 0.364
## .Q184 0.543 0.070 7.756 0.000 0.543 0.395
## .Q186 0.372 0.050 7.450 0.000 0.372 0.323
## .Q198 0.381 0.049 7.759 0.000 0.381 0.284
## .Q200 0.380 0.054 7.050 0.000 0.380 0.335
## .Q145 0.709 0.108 6.590 0.000 0.709 0.358
## .Q169 0.453 0.056 8.151 0.000 0.453 0.280
## .Q179 0.442 0.071 6.254 0.000 0.442 0.255
## .Q180 0.440 0.064 6.835 0.000 0.440 0.244
## .Q2 0.569 0.077 7.432 0.000 0.569 0.447
## .Q4 0.468 0.062 7.547 0.000 0.468 0.421
## .Q100 0.287 0.040 7.129 0.000 0.287 0.291
## .Q149 0.270 0.036 7.483 0.000 0.270 0.236
## .Q18 0.853 0.097 8.796 0.000 0.853 0.467
## .Q19 0.439 0.050 8.720 0.000 0.439 0.287
## .Q20 0.473 0.066 7.182 0.000 0.473 0.282
## .Q202 0.704 0.086 8.140 0.000 0.704 0.425
## F1 0.595 0.099 6.032 0.000 1.000 1.000
## F2 1.273 0.133 9.557 0.000 1.000 1.000
## F3 0.704 0.091 7.752 0.000 1.000 1.000
## F4 0.972 0.119 8.179 0.000 1.000 1.000fm## chisq.scaled df.scaled pvalue.scaled
## 252.106 129.000 0.000
## cfi.robust tli.robust ecvi
## 0.961 0.954 1.305
## srmr rmsea.robust rmsea.ci.lower.robust
## 0.037 0.063 0.052
## rmsea.ci.upper.robust
## 0.075res1 <- residuals(fit, type = "cor")$cov
res1[upper.tri(res1, diag = T)] <- NA
v1 <- as.vector(res1)
v2 <- v1[!is.na(v1)]
qqPlot(v2, id = F)
Five-factor models
#factors5_1
model <- '
F1 =~ Q186 + Q187 + Q198 + Q199 + Q200
F2 =~ Q107 + Q145 + Q169 + Q179 + Q180
F3 =~ Q14 + Q15 + Q16 + Q106 + Q164
F4 =~ Q17 + Q18 + Q19 + Q20 + Q202
F5 =~ Q1 + Q2 + Q4 + Q100 + Q149'
fit <- cfa(model, estimator = "MLR", data = b2)
fm <- fitMeasures(fit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "cfi.robust", "tli.robust", "ecvi", "srmr", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust"))
models[5,1] <- "5 factors (1)"
models[5,2:ncol(models)] <- fm
summary(fit,
#fit.measures = TRUE,
standardized = TRUE
#,rsq = TRUE
)## lavaan 0.6-12 ended normally after 65 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 60
##
## Number of observations 348
##
## Model Test User Model:
## Standard Robust
## Test Statistic 639.384 440.391
## Degrees of freedom 265 265
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.452
## Yuan-Bentler correction (Mplus variant)
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 =~
## Q186 1.000 0.878 0.819
## Q187 1.012 0.044 22.851 0.000 0.889 0.849
## Q198 1.172 0.070 16.824 0.000 1.029 0.889
## Q199 1.151 0.071 16.253 0.000 1.011 0.887
## Q200 1.030 0.075 13.728 0.000 0.904 0.849
## F2 =~
## Q107 1.000 0.811 0.653
## Q145 1.392 0.113 12.346 0.000 1.129 0.802
## Q169 1.350 0.103 13.123 0.000 1.095 0.860
## Q179 1.377 0.120 11.451 0.000 1.117 0.849
## Q180 1.436 0.114 12.644 0.000 1.165 0.868
## F3 =~
## Q14 1.000 0.818 0.846
## Q15 1.003 0.066 15.143 0.000 0.820 0.829
## Q16 0.933 0.054 17.311 0.000 0.763 0.798
## Q106 0.902 0.057 15.888 0.000 0.737 0.746
## Q164 1.114 0.073 15.280 0.000 0.911 0.816
## F4 =~
## Q17 1.000 0.848 0.781
## Q18 1.165 0.093 12.461 0.000 0.988 0.731
## Q19 1.223 0.085 14.403 0.000 1.037 0.838
## Q20 1.272 0.097 13.085 0.000 1.079 0.833
## Q202 1.168 0.077 15.254 0.000 0.991 0.770
## F5 =~
## Q1 1.000 0.734 0.669
## Q2 1.175 0.084 13.932 0.000 0.862 0.764
## Q4 1.086 0.099 10.941 0.000 0.797 0.755
## Q100 1.143 0.094 12.175 0.000 0.839 0.845
## Q149 1.259 0.101 12.470 0.000 0.923 0.864
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## F1 ~~
## F2 0.486 0.079 6.127 0.000 0.683 0.683
## F3 0.629 0.093 6.742 0.000 0.875 0.875
## F4 0.606 0.084 7.210 0.000 0.814 0.814
## F5 0.538 0.071 7.620 0.000 0.835 0.835
## F2 ~~
## F3 0.455 0.069 6.632 0.000 0.686 0.686
## F4 0.445 0.073 6.069 0.000 0.647 0.647
## F5 0.393 0.061 6.455 0.000 0.661 0.661
## F3 ~~
## F4 0.551 0.080 6.863 0.000 0.795 0.795
## F5 0.504 0.064 7.920 0.000 0.840 0.840
## F4 ~~
## F5 0.512 0.063 8.153 0.000 0.822 0.822
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Q186 0.379 0.058 6.573 0.000 0.379 0.329
## .Q187 0.306 0.052 5.913 0.000 0.306 0.279
## .Q198 0.280 0.041 6.855 0.000 0.280 0.209
## .Q199 0.278 0.039 7.090 0.000 0.278 0.214
## .Q200 0.317 0.051 6.193 0.000 0.317 0.279
## .Q107 0.884 0.085 10.351 0.000 0.884 0.573
## .Q145 0.706 0.105 6.712 0.000 0.706 0.356
## .Q169 0.422 0.052 8.169 0.000 0.422 0.260
## .Q179 0.484 0.073 6.641 0.000 0.484 0.279
## .Q180 0.444 0.059 7.535 0.000 0.444 0.247
## .Q14 0.267 0.051 5.177 0.000 0.267 0.285
## .Q15 0.305 0.047 6.544 0.000 0.305 0.312
## .Q16 0.331 0.043 7.645 0.000 0.331 0.362
## .Q106 0.432 0.043 10.129 0.000 0.432 0.443
## .Q164 0.417 0.058 7.156 0.000 0.417 0.334
## .Q17 0.460 0.055 8.288 0.000 0.460 0.390
## .Q18 0.849 0.092 9.203 0.000 0.849 0.465
## .Q19 0.456 0.048 9.431 0.000 0.456 0.298
## .Q20 0.511 0.066 7.744 0.000 0.511 0.305
## .Q202 0.674 0.081 8.292 0.000 0.674 0.407
## .Q1 0.663 0.064 10.400 0.000 0.663 0.552
## .Q2 0.529 0.072 7.371 0.000 0.529 0.416
## .Q4 0.478 0.063 7.634 0.000 0.478 0.430
## .Q100 0.283 0.038 7.411 0.000 0.283 0.287
## .Q149 0.289 0.037 7.725 0.000 0.289 0.253
## F1 0.772 0.116 6.630 0.000 1.000 1.000
## F2 0.658 0.111 5.928 0.000 1.000 1.000
## F3 0.668 0.092 7.230 0.000 1.000 1.000
## F4 0.719 0.107 6.744 0.000 1.000 1.000
## F5 0.538 0.077 6.983 0.000 1.000 1.000fm## chisq.scaled df.scaled pvalue.scaled
## 440.391 265.000 0.000
## cfi.robust tli.robust ecvi
## 0.963 0.958 2.182
## srmr rmsea.robust rmsea.ci.lower.robust
## 0.038 0.053 0.044
## rmsea.ci.upper.robust
## 0.061res1 <- residuals(fit, type = "cor")$cov
res1[upper.tri(res1, diag = T)] <- NA
v1 <- as.vector(res1)
v2 <- v1[!is.na(v1)]
qqPlot(v2, id = F)
#factors5_2
model <-
'Touch =~ Q145 + Q169 + Q179 + Q180
Service =~ Q2 + Q4 + Q100 + Q149
Words =~ Q184 + Q186 + Q198 + Q200
Time =~ Q14 + Q15 + Q16 + Q106
Gifts =~ Q18 + Q19 + Q20 + Q202'
fit <- cfa(model, estimator = "MLR", data = b2)
fm <- fitMeasures(fit, fit.measures = c("chisq.scaled", "df.scaled", "pvalue.scaled", "cfi.robust", "tli.robust", "ecvi","srmr", "rmsea.robust", "rmsea.ci.lower.robust", "rmsea.ci.upper.robust"))
models[6,1] <- "5 factors (2)"
models[6,2:ncol(models)] <- fm
summary(fit,
#fit.measures = TRUE,
standardized = TRUE
,rsq = TRUE
)## lavaan 0.6-12 ended normally after 59 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 50
##
## Number of observations 348
##
## Model Test User Model:
## Standard Robust
## Test Statistic 368.736 256.292
## Degrees of freedom 160 160
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.439
## Yuan-Bentler correction (Mplus variant)
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Touch =~
## Q145 1.000 1.127 0.801
## Q169 0.961 0.051 18.929 0.000 1.083 0.851
## Q179 1.007 0.053 19.024 0.000 1.135 0.863
## Q180 1.034 0.049 21.238 0.000 1.165 0.868
## Service =~
## Q2 1.000 0.847 0.751
## Q4 0.944 0.058 16.305 0.000 0.799 0.758
## Q100 0.996 0.070 14.297 0.000 0.843 0.849
## Q149 1.091 0.067 16.336 0.000 0.924 0.865
## Words =~
## Q184 1.000 0.901 0.769
## Q186 0.975 0.060 16.171 0.000 0.878 0.819
## Q198 1.127 0.058 19.585 0.000 1.016 0.877
## Q200 0.997 0.064 15.479 0.000 0.898 0.843
## Time =~
## Q14 1.000 0.837 0.865
## Q15 0.969 0.068 14.230 0.000 0.811 0.820
## Q16 0.913 0.054 16.890 0.000 0.763 0.799
## Q106 0.887 0.056 15.836 0.000 0.742 0.751
## Gifts =~
## Q18 1.000 0.983 0.728
## Q19 1.067 0.080 13.327 0.000 1.049 0.848
## Q20 1.113 0.073 15.176 0.000 1.094 0.845
## Q202 0.992 0.069 14.382 0.000 0.975 0.758
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Touch ~~
## Service 0.615 0.080 7.679 0.000 0.645 0.645
## Words 0.694 0.095 7.312 0.000 0.684 0.684
## Time 0.630 0.081 7.778 0.000 0.668 0.668
## Gifts 0.711 0.098 7.225 0.000 0.642 0.642
## Service ~~
## Words 0.650 0.083 7.829 0.000 0.852 0.852
## Time 0.585 0.075 7.850 0.000 0.826 0.826
## Gifts 0.672 0.075 8.969 0.000 0.807 0.807
## Words ~~
## Time 0.663 0.087 7.656 0.000 0.880 0.880
## Gifts 0.709 0.081 8.771 0.000 0.801 0.801
## Time ~~
## Gifts 0.643 0.082 7.861 0.000 0.782 0.782
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Q145 0.711 0.108 6.590 0.000 0.711 0.359
## .Q169 0.447 0.056 8.043 0.000 0.447 0.276
## .Q179 0.444 0.071 6.261 0.000 0.444 0.256
## .Q180 0.443 0.064 6.871 0.000 0.443 0.246
## .Q2 0.555 0.075 7.406 0.000 0.555 0.436
## .Q4 0.474 0.063 7.489 0.000 0.474 0.426
## .Q100 0.275 0.040 6.919 0.000 0.275 0.279
## .Q149 0.287 0.038 7.579 0.000 0.287 0.252
## .Q184 0.562 0.075 7.487 0.000 0.562 0.409
## .Q186 0.379 0.057 6.606 0.000 0.379 0.330
## .Q198 0.308 0.042 7.328 0.000 0.308 0.230
## .Q200 0.328 0.048 6.844 0.000 0.328 0.289
## .Q14 0.235 0.048 4.931 0.000 0.235 0.251
## .Q15 0.320 0.049 6.474 0.000 0.320 0.328
## .Q16 0.330 0.044 7.550 0.000 0.330 0.361
## .Q106 0.425 0.044 9.780 0.000 0.425 0.436
## .Q18 0.858 0.096 8.899 0.000 0.858 0.470
## .Q19 0.431 0.050 8.612 0.000 0.431 0.281
## .Q20 0.478 0.066 7.208 0.000 0.478 0.286
## .Q202 0.704 0.088 8.007 0.000 0.704 0.425
## Touch 1.270 0.133 9.529 0.000 1.000 1.000
## Service 0.717 0.091 7.905 0.000 1.000 1.000
## Words 0.812 0.109 7.454 0.000 1.000 1.000
## Time 0.700 0.094 7.407 0.000 1.000 1.000
## Gifts 0.966 0.118 8.182 0.000 1.000 1.000
##
## R-Square:
## Estimate
## Q145 0.641
## Q169 0.724
## Q179 0.744
## Q180 0.754
## Q2 0.564
## Q4 0.574
## Q100 0.721
## Q149 0.748
## Q184 0.591
## Q186 0.670
## Q198 0.770
## Q200 0.711
## Q14 0.749
## Q15 0.672
## Q16 0.639
## Q106 0.564
## Q18 0.530
## Q19 0.719
## Q20 0.714
## Q202 0.575fm## chisq.scaled df.scaled pvalue.scaled
## 256.292 160.000 0.000
## cfi.robust tli.robust ecvi
## 0.973 0.968 1.347
## srmr rmsea.robust rmsea.ci.lower.robust
## 0.036 0.050 0.038
## rmsea.ci.upper.robust
## 0.061res1 <- residuals(fit, type = "cor")$cov
res1[upper.tri(res1, diag = T)] <- NA
v1 <- as.vector(res1)
v2 <- v1[!is.na(v1)]
qqPlot(v2, id = F)
Comparative measures of model fit
models$chisq <- round(models$chisq, 2)
models$CFI <- round(models$CFI, 4)
models$TLI <- round(models$TLI, 4)
models$ecvi <- round(models$ecvi, 2)
models$srmr <- round(models$srmr, 3)
models$RMSEA <- round(models$RMSEA, 4)
models$RMSEA_lower <- round(models$RMSEA_lower, 4)
models$RMSEA_upper <- round(models$RMSEA_upper, 4)
write.csv(models, "models_CFA_study3.csv", row.names = F)
models## model chisq df p CFI TLI ecvi srmr RMSEA
## 1 1 factor 400.16 170 0.000000e+00 0.9322 0.9242 2.27 0.039 0.0831
## 2 2 factors 348.57 169 1.498801e-14 0.9512 0.9452 1.77 0.037 0.0683
## 3 3 factors 352.04 167 2.775558e-15 0.9467 0.9394 1.80 0.049 0.0699
## 4 4 factors 252.11 129 5.485111e-10 0.9613 0.9541 1.30 0.037 0.0634
## 5 5 factors (1) 440.39 265 7.077383e-11 0.9628 0.9579 2.18 0.038 0.0525
## 6 5 factors (2) 256.29 160 2.029545e-06 0.9734 0.9684 1.35 0.036 0.0499
## RMSEA_lower RMSEA_upper
## 1 0.0726 0.0937
## 2 0.0581 0.0785
## 3 0.0597 0.0801
## 4 0.0517 0.0750
## 5 0.0438 0.0611
## 6 0.0382 0.0610Models with 4 and 5 factors have a better CFI and TLI compared to models with 1 to 3 factors. RMSEA values are also better for these models. The smaller model with 5 factors has the best CFI, TLI, and RMSEA, while having a relatively small ECVI, which indicates that this model is likely to best predict new data.
Selected five-factor model - alphas and correlations
Touch_items <- c("Q145", "Q169", "Q179", "Q180")
Service_items <- c("Q2", "Q4", "Q100", "Q149")
Words_items <- c("Q184", "Q186", "Q198", "Q200")
Time_items <- c("Q14", "Q15", "Q16", "Q106")
Gifts_items <- c("Q18", "Q19", "Q20", "Q202")
Touch <- apply(b2[,Touch_items], 1, mean)
Service <- apply(b2[,Service_items], 1, mean)
Words <- apply(b2[,Words_items], 1, mean)
Time <- apply(b2[,Time_items], 1, mean)
Gifts <- apply(b2[,Gifts_items], 1, mean)
factors <- data.frame(Touch, Service, Words, Time, Gifts)
round(cor(factors),2) ## Touch Service Words Time Gifts
## Touch 1.00 0.57 0.61 0.60 0.59
## Service 0.57 1.00 0.76 0.72 0.71
## Words 0.61 0.76 1.00 0.78 0.70
## Time 0.60 0.72 0.78 1.00 0.68
## Gifts 0.59 0.71 0.70 0.68 1.00apa.cor.table(factors, "cor_cfa.doc")##
##
## Means, standard deviations, and correlations with confidence intervals
##
##
## Variable M SD 1 2 3 4
## 1. Touch 3.78 1.19
##
## 2. Service 4.15 0.91 .57**
## [.50, .64]
##
## 3. Words 4.22 0.97 .61** .76**
## [.54, .68] [.71, .80]
##
## 4. Time 4.29 0.84 .60** .72** .78**
## [.53, .66] [.66, .77] [.73, .82]
##
## 5. Gifts 3.73 1.10 .59** .71** .70** .68**
## [.52, .66] [.65, .76] [.64, .75] [.62, .73]
##
##
## Note. M and SD are used to represent mean and standard deviation, respectively.
## Values in square brackets indicate the 95% confidence interval.
## The confidence interval is a plausible range of population correlations
## that could have caused the sample correlation (Cumming, 2014).
## * indicates p < .05. ** indicates p < .01.
## alpha_Touch <- psych::alpha(b2[,Touch_items]); alpha_Touch$total[1]## raw_alpha
## 0.9087784alpha_Service <- psych::alpha(b2[,Service_items]); alpha_Service$total[1]## raw_alpha
## 0.8819057alpha_Words <- psych::alpha(b2[,Words_items]); alpha_Words$total[1]## raw_alpha
## 0.89382alpha_Time <- psych::alpha(b2[,Time_items]); alpha_Time$total[1]## raw_alpha
## 0.8844603alpha_Gifts <- psych::alpha(b2[,Gifts_items]); alpha_Gifts$total[1]## raw_alpha
## 0.8713219Five-factor model - AVE and composite reliability
AVE(fit)## Touch Service Words Time Gifts
## 0.713 0.647 0.684 0.655 0.630compRelSEM(fit)## Touch Service Words Time Gifts
## 0.908 0.874 0.899 0.880 0.872Five-factor model - are factors distinct enough?
discriminantValidity(fit, merge = TRUE, level = 0.95)## Some of the latent variable variances are estimated instead of fixed to 1. The model is re-estimated by scaling the latent variables by fixing their variances and freeing all factor loadings.## lhs op rhs est ci.lower ci.upper Df AIC BIC
## 1 Touch ~~ Service 0.6445072 0.5474660 0.7415483 164 17082.98 17260.18
## 2 Touch ~~ Words 0.6836805 0.5929424 0.7744186 164 17061.31 17238.51
## 3 Touch ~~ Time 0.6684051 0.5786492 0.7581611 164 17053.74 17230.94
## 4 Touch ~~ Gifts 0.6416311 0.5374729 0.7457893 164 17033.99 17211.19
## 5 Service ~~ Words 0.8523675 0.7823925 0.9223426 164 16744.09 16921.29
## 6 Service ~~ Time 0.8258140 0.7560666 0.8955613 164 16766.35 16943.55
## 7 Service ~~ Gifts 0.8073330 0.7300541 0.8846120 164 16774.21 16951.41
## 8 Words ~~ Time 0.8797581 0.8144989 0.9450172 164 16711.16 16888.36
## 9 Words ~~ Gifts 0.8007022 0.7371807 0.8642237 164 16797.66 16974.86
## 10 Time ~~ Gifts 0.7819998 0.7048329 0.8591666 164 16805.10 16982.31
## Chisq Chisq diff Df diff Pr(>Chisq)
## 1 826.5000 204.50508 4 4.038182e-43
## 2 804.8240 160.51139 4 1.135657e-33
## 3 797.2521 193.49720 4 9.391317e-41
## 4 777.5070 178.92163 4 1.270932e-37
## 5 487.6053 50.64783 4 2.644314e-10
## 6 509.8669 59.00080 4 4.703758e-12
## 7 517.7298 56.00995 4 1.995570e-11
## 8 454.6724 34.36693 4 6.266037e-07
## 9 541.1754 85.40784 4 1.242924e-17
## 10 548.6204 66.20917 4 1.431113e-13