library(psych)
library(kableExtra)
library(ggplot2)
library(tidyr)
library(dplyr)
library(nFactors)
# remotes::install_version("lavaan", version = "0.6.17")
library(lavaan)
library(semPlot)
library(semTools)
library(tidySEM)
library(stringr)
library(corrplot)
library(mifa) # MI covariance matrix for EFA
library(mice) # underlying imputation engine
library(sjPlot)
library(metafor)
library(broom)
# library(DiagrammeR)
library(lme4)
library(lmerTest)s1 <- read.csv(file="data/s1 - export.csv", header=T)
s2 <- read.csv(file="data/s2 - export.csv", header=T)
# stress - Q6
# everyday discrimination - Q7 & Q39
# expectation of rejection - Q8 & Q40
# ders - Q9
# idm - Q18
# condition
# samp
# misinformation acceptance
# accurate information acceptance
# social support - Q31
# scipop - Q29
# anomie - Q30
# demographics
df1 <- subset(s1, select=c(id, samp, grep("Q6_",colnames(s1)), grep("Q39_",colnames(s1)), grep("Q40_",colnames(s1)), grep("Q9_",colnames(s1)), grep("Q18_",colnames(s1))))
df2 <- subset(s2, select=c(id, condition, samp, grep("Q6_",colnames(s2)), grep("Q7_",colnames(s2)), grep("Q8_",colnames(s2)), grep("Q9_",colnames(s2)), grep("Q18_",colnames(s2)), grep("Q31_",colnames(s2)), grep("Q29_",colnames(s2)), grep("Q30_",colnames(s2))))
df1_attn <- subset(s1, select=c(id, Q39_7, Q18_10))
df2_attn <- subset(s2, select=c(id, Q18_12))
df_attn <- bind_rows(df1_attn, df2_attn)
df1 <- subset(df1, select=-c(Q18_10, Q39_7))
df2 <- subset(df2, select=-c(Q18_12))
a_list <- c("23","34","36","38","48","52","56","60","64","70","72","76","80","84","88")
m_list <- c("40","42","44","46","50","54","58","62","66","68","74","78","82","86","90")
# Identify target columns
target_cols <- grep(
paste(c(a_list, m_list), collapse = "|"),
names(s2),
value = TRUE
)
df2 <- cbind(df2, s2[, intersect(target_cols, colnames(s2))])
df1$study <- "s1"
df2$study <- "s2"Already cleaned in imported files.
df_attn$attn_violations <- rowSums(
cbind(
df_attn$Q39_7 != 3,
df_attn$Q18_10 != 2,
df_attn$Q18_12 != 4
),
na.rm = TRUE
)
rm(df_attn, df1_attn, df2_attn)describe(subset(df1, select=c(grep("Q6_",colnames(df1))))) vars n mean sd median trimmed mad min max range skew kurtosis se
Q6_1 1 420 3.00 0.97 3 3.01 1.48 1 5 4 -0.07 -0.26 0.05
Q6_2 2 420 3.12 1.17 3 3.15 1.48 1 5 4 -0.07 -0.74 0.06
Q6_3 3 420 3.77 1.08 4 3.88 1.48 1 5 4 -0.67 -0.20 0.05
Q6_4 4 420 3.36 1.03 3 3.38 1.48 1 5 4 -0.29 -0.34 0.05
Q6_5 5 420 3.03 0.96 3 3.02 1.48 1 5 4 0.06 -0.24 0.05
Q6_6 6 420 2.73 1.09 3 2.71 1.48 1 5 4 0.18 -0.57 0.05
Q6_7 7 420 3.24 0.97 3 3.24 1.48 1 5 4 -0.20 -0.31 0.05
Q6_8 8 420 3.18 1.01 3 3.16 1.48 1 5 4 -0.09 -0.49 0.05
Q6_9 9 420 3.00 1.07 3 3.00 1.48 1 5 4 0.00 -0.63 0.05
Q6_10 10 420 2.82 1.22 3 2.78 1.48 1 5 4 0.16 -0.87 0.06d <- subset(df1, select=c(grep("Q6_", colnames(df1))))
labels <- c(
Q6_1 = "Helplessness: Unexpected Upset",
Q6_2 = "Helplessness: Loss of Control",
Q6_3 = "Helplessness: Stressed",
Q6_4 = "Self-Efficacy: Confident",
Q6_5 = "Self-Efficacy: Going Your Way",
Q6_6 = "Helplessness: Unable to Cope",
Q6_7 = "Self-Efficacy: Control Irritations",
Q6_8 = "Self-Efficacy: On Top of Things",
Q6_9 = "Helplessness: Uncontrollable Anger",
Q6_10 = "Helplessness: Difficulties Piling Up"
)
# shorten labels to max 40 characters (you can adjust the number)
# labels <- str_trunc(labels, width = 60, side = "right")
# replace column names
names(d) <- labels
d_long <- stack(d)
names(d_long) <- c("value", "variable")
ggplot(d_long, aes(x = value)) +
geom_histogram(binwidth = 1, fill = "steelblue", color = "white") +
facet_wrap(~ variable, labeller = label_wrap_gen(width = 30)) +
theme_minimal()
Dropped 1 from Helpless to get CFI/TLI above .95
d <- subset(df1, select=c(grep("Q6_", colnames(df1))))
# specify model: latent variables =~ observed indicators
model <- '
helpless =~ Q6_2 + Q6_3 + Q6_6 + Q6_9 + Q6_10
efficacy =~ Q6_4 + Q6_5 + Q6_7 + Q6_8
'
# fit CFA
fit <- cfa(model, data = d, std.lv = TRUE)
# summary with fit indices and standardized loadings
summary(fit, fit.measures = TRUE, standardized = TRUE)lavaan 0.6-21 ended normally after 18 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 19
Number of observations 420
Model Test User Model:
Test statistic 75.807
Degrees of freedom 26
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 1893.172
Degrees of freedom 36
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.973
Tucker-Lewis Index (TLI) 0.963
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -4686.789
Loglikelihood unrestricted model (H1) -4648.885
Akaike (AIC) 9411.577
Bayesian (BIC) 9488.342
Sample-size adjusted Bayesian (SABIC) 9428.049
Root Mean Square Error of Approximation:
RMSEA 0.068
90 Percent confidence interval - lower 0.050
90 Percent confidence interval - upper 0.085
P-value H_0: RMSEA <= 0.050 0.049
P-value H_0: RMSEA >= 0.080 0.131
Standardized Root Mean Square Residual:
SRMR 0.035
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
helpless =~
Q6_2 0.908 0.050 18.210 0.000 0.908 0.776
Q6_3 0.773 0.048 16.203 0.000 0.773 0.714
Q6_6 0.884 0.046 19.319 0.000 0.884 0.808
Q6_9 0.709 0.048 14.640 0.000 0.709 0.661
Q6_10 1.043 0.049 21.192 0.000 1.043 0.859
efficacy =~
Q6_4 0.823 0.044 18.566 0.000 0.823 0.802
Q6_5 0.766 0.041 18.480 0.000 0.766 0.799
Q6_7 0.572 0.046 12.415 0.000 0.572 0.590
Q6_8 0.772 0.044 17.407 0.000 0.772 0.765
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
helpless ~~
efficacy -0.687 0.034 -20.337 0.000 -0.687 -0.687
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.Q6_2 0.544 0.046 11.888 0.000 0.544 0.398
.Q6_3 0.575 0.045 12.722 0.000 0.575 0.490
.Q6_6 0.415 0.037 11.234 0.000 0.415 0.347
.Q6_9 0.647 0.049 13.175 0.000 0.647 0.563
.Q6_10 0.387 0.040 9.641 0.000 0.387 0.262
.Q6_4 0.376 0.037 10.109 0.000 0.376 0.357
.Q6_5 0.332 0.033 10.189 0.000 0.332 0.361
.Q6_7 0.614 0.046 13.233 0.000 0.614 0.652
.Q6_8 0.421 0.038 11.066 0.000 0.421 0.414
helpless 1.000 1.000 1.000
efficacy 1.000 1.000 1.000# modindices(fit, sort. = T)
# basic path diagram
semPaths(
fit,
what = "std", # show standardized loadings
layout = "tree", # neat hierarchical layout
rotation = 2,
style = "lisrel", # clean style
residuals = FALSE, # hide residual arrows (optional)
intercepts = FALSE,
sizeMan = 6, # size of observed variables
sizeLat = 8, # size of latent variables
edge.label.cex = 0.9, # path label size
label.cex = 1.1, # node label size
nCharNodes = 0, # keep full variable names
mar = c(6, 6, 6, 6)
)
# Omega (default)
compRelSEM(fit)$helpless
Composite `helpless` is composed of observed variables:
Q6_2, Q6_3, Q6_6, Q6_9, Q6_10
True-score variance is represented by common factor(s):
helpless
Total variance of composite `helpless` determined from the unrestricted model.
The proportion attributable to "true" scores is its model-based estimate of reliability ("omega"):
[1] 0.874
$efficacy
Composite `efficacy` is composed of observed variables:
Q6_4, Q6_5, Q6_7, Q6_8
True-score variance is represented by common factor(s):
efficacy
Total variance of composite `efficacy` determined from the unrestricted model.
The proportion attributable to "true" scores is its model-based estimate of reliability ("omega"):
[1] 0.831# Alpha
compRelSEM(fit, tau.eq = TRUE)$helpless
Composite `helpless` is composed of observed variables:
Q6_2, Q6_3, Q6_6, Q6_9, Q6_10
Total variance of composite `helpless` determined from the unrestricted model.
Coefficient alpha would be:
[1] 0.877
$efficacy
Composite `efficacy` is composed of observed variables:
Q6_4, Q6_5, Q6_7, Q6_8
Total variance of composite `efficacy` determined from the unrestricted model.
Coefficient alpha would be:
[1] 0.827describe(subset(df1, select=c(grep("Q39_",colnames(df1))))) vars n mean sd median trimmed mad min max range skew kurtosis se
Q39_1 1 420 2.85 1.30 3 2.78 1.48 1 6 5 0.28 -0.47 0.06
Q39_2 2 420 2.82 1.27 3 2.75 1.48 1 6 5 0.35 -0.33 0.06
Q39_3 3 420 1.99 1.04 2 1.85 1.48 1 6 5 0.97 0.66 0.05
Q39_4 4 420 2.77 1.40 3 2.67 1.48 1 6 5 0.38 -0.74 0.07
Q39_5 5 420 2.00 1.22 2 1.79 1.48 1 6 5 1.18 0.76 0.06
Q39_6 6 419 2.15 1.17 2 2.00 1.48 1 6 5 0.88 0.20 0.06
Q39_8 7 420 3.10 1.29 3 3.10 1.48 1 6 5 0.05 -0.69 0.06
Q39_9 8 420 2.25 1.24 2 2.10 1.48 1 6 5 0.90 0.24 0.06
Q39_10 9 420 1.76 0.92 2 1.63 1.48 1 6 5 1.35 2.37 0.05d <- na.omit(subset(df1, select=c(grep("Q39_", colnames(df1)))))
labels <- c(
Q39_1 = "Less Courtesy",
Q39_2 = "Less Respect",
Q39_3 = "Poorer Service",
Q39_4 = "Seen as Not Smart",
Q39_5 = "Seen as Threatening",
Q39_6 = "Seen as Dishonest",
Q39_8 = "Inferior",
Q39_9 = "Insulted",
Q39_10 = "Threatened/Harassed"
)
names(d) <- labels
d_long <- stack(d)
names(d_long) <- c("value", "variable")
ggplot(d_long, aes(x = value)) +
geom_histogram(binwidth = 1, fill = "steelblue", color = "white") +
facet_wrap(~ variable, labeller = label_wrap_gen(width = 30)) +
theme_minimal()
ev <- eigen(cor(d)) # get eigenvalues
ap <- parallel(subject=nrow(d),var=ncol(d),rep=100,cent=.05) # run the parallel analysis, gives us another perspective on how many factors should be used in the model
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) # creates the scree plot
plotnScree(nS) # shows us the scree plot, look for the elbows
EFA <- factanal(d, factors = 1, rotation = "promax")
print(EFA, digits=3, cutoff=.4, sort=F)
Call:
factanal(x = d, factors = 1, rotation = "promax")
Uniquenesses:
Less Courtesy Less Respect Poorer Service Seen as Not Smart
0.175 0.151 0.684 0.542
Seen as Threatening Seen as Dishonest Inferior Insulted
0.804 0.658 0.504 0.678
Threatened/Harassed
0.716
Loadings:
Factor1
Less Courtesy 0.908
Less Respect 0.921
Poorer Service 0.562
Seen as Not Smart 0.677
Seen as Threatening 0.443
Seen as Dishonest 0.585
Inferior 0.705
Insulted 0.568
Threatened/Harassed 0.533
Factor1
SS loadings 4.089
Proportion Var 0.454
Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 373.2 on 27 degrees of freedom.
The p-value is 1.39e-62 d <- subset(d, select = -c(`Seen as Threatening`, `Threatened/Harassed`, `Insulted`, `Poorer Service`))
ev <- eigen(cor(d)) # get eigenvalues
ap <- parallel(subject=nrow(d),var=ncol(d),rep=100,cent=.05) # run the parallel analysis, gives us another perspective on how many factors should be used in the model
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) # creates the scree plot
plotnScree(nS) # shows us the scree plot, look for the elbows
EFA <- factanal(d, factors = 1, rotation = "promax")
print(EFA, digits=3, cutoff=.4, sort=F)
Call:
factanal(x = d, factors = 1, rotation = "promax")
Uniquenesses:
Less Courtesy Less Respect Seen as Not Smart Seen as Dishonest
0.132 0.111 0.595 0.712
Inferior
0.545
Loadings:
Factor1
Less Courtesy 0.931
Less Respect 0.943
Seen as Not Smart 0.637
Seen as Dishonest 0.536
Inferior 0.674
Factor1
SS loadings 2.904
Proportion Var 0.581
Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 125.02 on 5 degrees of freedom.
The p-value is 2.71e-25 RMSEA is high, likely due to low degrees of freedom.
End up with four items – treated less courtesy, treated less respect, seen as not smart, and inferior.
Kenny, D. A., Kaniskan, B., & McCoach, D. B. (2015). The performance of RMSEA in models with small degrees of freedom. Sociological methods & research, 44(3), 486-507.
# labels <- c(
# Q39_1 = "Less Courtesy",
# Q39_2 = "Less Respect",
# Q39_4 = "Seen as Not Smart",
# Q39_8 = "Inferior",
# )
d <- na.omit(subset(df2, select=c(grep("Q7_", colnames(df2)))))
# specify model: latent variables =~ observed indicators
model <- '
discrim =~ Q7_1 + Q7_2 + Q7_4 + Q7_8
'
# fit CFA
fit <- cfa(model, data = d, std.lv = TRUE)
# summary with fit indices and standardized loadings
summary(fit, fit.measures = TRUE, standardized = TRUE)lavaan 0.6-21 ended normally after 20 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 8
Number of observations 184
Model Test User Model:
Test statistic 10.197
Degrees of freedom 2
P-value (Chi-square) 0.006
Model Test Baseline Model:
Test statistic 473.394
Degrees of freedom 6
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.982
Tucker-Lewis Index (TLI) 0.947
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -1030.948
Loglikelihood unrestricted model (H1) -1025.850
Akaike (AIC) 2077.896
Bayesian (BIC) 2103.615
Sample-size adjusted Bayesian (SABIC) 2078.277
Root Mean Square Error of Approximation:
RMSEA 0.149
90 Percent confidence interval - lower 0.068
90 Percent confidence interval - upper 0.245
P-value H_0: RMSEA <= 0.050 0.026
P-value H_0: RMSEA >= 0.080 0.924
Standardized Root Mean Square Residual:
SRMR 0.046
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
discrim =~
Q7_1 1.256 0.079 15.979 0.000 1.256 0.924
Q7_2 1.371 0.079 17.451 0.000 1.371 0.974
Q7_4 0.911 0.092 9.924 0.000 0.911 0.659
Q7_8 0.655 0.086 7.603 0.000 0.655 0.529
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.Q7_1 0.269 0.059 4.543 0.000 0.269 0.146
.Q7_2 0.102 0.063 1.621 0.105 0.102 0.051
.Q7_4 1.081 0.117 9.234 0.000 1.081 0.565
.Q7_8 1.103 0.117 9.428 0.000 1.103 0.720
discrim 1.000 1.000 1.000# modindices(fit, sort. = T)
# basic path diagram
semPaths(
fit,
what = "std", # show standardized loadings
layout = "tree", # neat hierarchical layout
rotation = 2,
style = "lisrel", # clean style
residuals = FALSE, # hide residual arrows (optional)
intercepts = FALSE,
sizeMan = 6, # size of observed variables
sizeLat = 8, # size of latent variables
edge.label.cex = 0.9, # path label size
label.cex = 1.1, # node label size
nCharNodes = 0, # keep full variable names
mar = c(6, 6, 6, 6)
)
# Omega (default)
compRelSEM(fit)$discrim
Composite `discrim` is composed of observed variables:
Q7_1, Q7_2, Q7_4, Q7_8
True-score variance is represented by common factor(s):
discrim
Total variance of composite `discrim` determined from the unrestricted model.
The proportion attributable to "true" scores is its model-based estimate of reliability ("omega"):
[1] 0.856# Alpha
compRelSEM(fit, tau.eq = TRUE)$discrim
Composite `discrim` is composed of observed variables:
Q7_1, Q7_2, Q7_4, Q7_8
Total variance of composite `discrim` determined from the unrestricted model.
Coefficient alpha would be:
[1] 0.861describe(subset(df1, select=c(grep("Q40_",colnames(df1))))) vars n mean sd median trimmed mad min max range skew kurtosis se
Q40_1 1 420 1.69 0.89 1 1.55 0.00 1 4 3 1.09 0.19 0.04
Q40_2 2 420 1.45 0.71 1 1.30 0.00 1 4 3 1.45 1.22 0.03
Q40_3 3 420 1.38 0.71 1 1.20 0.00 1 4 3 1.85 2.62 0.03
Q40_4 4 420 1.80 0.92 2 1.68 1.48 1 4 3 0.80 -0.51 0.04
Q40_5 5 420 1.82 0.95 2 1.70 1.48 1 4 3 0.78 -0.61 0.05
Q40_6 6 420 1.85 0.93 2 1.76 1.48 1 4 3 0.65 -0.80 0.05d <- na.omit(subset(df1, select=c(grep("Q40_", colnames(df1)))))
labels <- c(
Q40_1 = "Hiring Discrimination",
Q40_2 = "Seen as Untrustworthy",
Q40_3 = "Seen as Dangerous",
Q40_4 = "Devalued",
Q40_5 = "Looked Down On",
Q40_6 = "Seen as Less Intelligent"
)
names(d) <- labels
d_long <- stack(d)
names(d_long) <- c("value", "variable")
ggplot(d_long, aes(x = value)) +
geom_histogram(binwidth = 1, fill = "steelblue", color = "white") +
facet_wrap(~ variable, labeller = label_wrap_gen(width = 30)) +
theme_minimal()
ev <- eigen(cor(d)) # get eigenvalues
ap <- parallel(subject=nrow(d),var=ncol(d),rep=100,cent=.05) # run the parallel analysis, gives us another perspective on how many factors should be used in the model
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) # creates the scree plot
plotnScree(nS) # shows us the scree plot, look for the elbows
EFA <- factanal(d, factors = 1, rotation = "promax")
print(EFA, digits=3, cutoff=.4, sort=F)
Call:
factanal(x = d, factors = 1, rotation = "promax")
Uniquenesses:
Hiring Discrimination Seen as Untrustworthy Seen as Dangerous
0.529 0.488 0.556
Devalued Looked Down On Seen as Less Intelligent
0.125 0.138 0.357
Loadings:
Factor1
Hiring Discrimination 0.686
Seen as Untrustworthy 0.716
Seen as Dangerous 0.666
Devalued 0.935
Looked Down On 0.928
Seen as Less Intelligent 0.802
Factor1
SS loadings 3.807
Proportion Var 0.635
Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 140.3 on 9 degrees of freedom.
The p-value is 8.97e-26 Dropped ‘person like you is dangerous and unpredictable’
d <- na.omit(subset(df2, select=c(grep("Q8_", colnames(df2)))))
# specify model: latent variables =~ observed indicators
model <- '
discrim =~ Q8_1 + Q8_2 + Q8_4 + Q8_5 + Q8_6
'
# fit CFA
fit <- cfa(model, data = d, std.lv = TRUE)
# summary with fit indices and standardized loadings
summary(fit, fit.measures = TRUE, standardized = TRUE)lavaan 0.6-21 ended normally after 26 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 10
Number of observations 184
Model Test User Model:
Test statistic 8.919
Degrees of freedom 5
P-value (Chi-square) 0.112
Model Test Baseline Model:
Test statistic 681.994
Degrees of freedom 10
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.994
Tucker-Lewis Index (TLI) 0.988
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -893.823
Loglikelihood unrestricted model (H1) -889.363
Akaike (AIC) 1807.645
Bayesian (BIC) 1839.795
Sample-size adjusted Bayesian (SABIC) 1808.122
Root Mean Square Error of Approximation:
RMSEA 0.065
90 Percent confidence interval - lower 0.000
90 Percent confidence interval - upper 0.134
P-value H_0: RMSEA <= 0.050 0.297
P-value H_0: RMSEA >= 0.080 0.423
Standardized Root Mean Square Residual:
SRMR 0.025
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
discrim =~
Q8_1 0.660 0.061 10.749 0.000 0.660 0.702
Q8_2 0.524 0.050 10.398 0.000 0.524 0.685
Q8_4 0.908 0.054 16.821 0.000 0.908 0.942
Q8_5 0.904 0.055 16.435 0.000 0.904 0.929
Q8_6 0.811 0.060 13.478 0.000 0.811 0.822
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.Q8_1 0.448 0.049 9.067 0.000 0.448 0.507
.Q8_2 0.311 0.034 9.116 0.000 0.311 0.531
.Q8_4 0.105 0.021 4.944 0.000 0.105 0.113
.Q8_5 0.129 0.023 5.686 0.000 0.129 0.137
.Q8_6 0.315 0.037 8.439 0.000 0.315 0.324
discrim 1.000 1.000 1.000# modindices(fit, sort. = T)
# basic path diagram
semPaths(
fit,
what = "std", # show standardized loadings
layout = "tree", # neat hierarchical layout
rotation = 2,
style = "lisrel", # clean style
residuals = FALSE, # hide residual arrows (optional)
intercepts = FALSE,
sizeMan = 6, # size of observed variables
sizeLat = 8, # size of latent variables
edge.label.cex = 0.9, # path label size
label.cex = 1.1, # node label size
nCharNodes = 0, # keep full variable names
mar = c(6, 6, 6, 6)
)
# Omega (default)
compRelSEM(fit)$discrim
Composite `discrim` is composed of observed variables:
Q8_1, Q8_2, Q8_4, Q8_5, Q8_6
True-score variance is represented by common factor(s):
discrim
Total variance of composite `discrim` determined from the unrestricted model.
The proportion attributable to "true" scores is its model-based estimate of reliability ("omega"):
[1] 0.914# Alpha
compRelSEM(fit, tau.eq = TRUE)$discrim
Composite `discrim` is composed of observed variables:
Q8_1, Q8_2, Q8_4, Q8_5, Q8_6
Total variance of composite `discrim` determined from the unrestricted model.
Coefficient alpha would be:
[1] 0.909describe(subset(df1, select=c(grep("Q18_",colnames(df1))))) vars n mean sd median trimmed mad min max range skew kurtosis se
Q18_1 1 420 1.96 0.86 2 1.90 1.48 1 4 3 0.48 -0.66 0.04
Q18_2 2 419 2.28 0.98 2 2.23 1.48 1 4 3 0.07 -1.11 0.05
Q18_3 3 420 2.46 1.01 3 2.45 1.48 1 4 3 -0.10 -1.11 0.05
Q18_4 4 420 2.49 1.02 3 2.49 1.48 1 4 3 -0.12 -1.12 0.05
Q18_5 5 420 2.37 0.99 2 2.33 1.48 1 4 3 0.06 -1.08 0.05
Q18_6 6 420 2.53 1.04 3 2.54 1.48 1 4 3 -0.17 -1.16 0.05
Q18_7 7 420 1.82 0.78 2 1.74 1.48 1 4 3 0.72 0.12 0.04
Q18_8 8 420 1.66 0.79 1 1.54 0.00 1 4 3 1.00 0.27 0.04
Q18_9 9 420 1.93 0.95 2 1.83 1.48 1 4 3 0.60 -0.79 0.05
Q18_11 10 420 2.30 1.02 2 2.24 1.48 1 4 3 0.16 -1.14 0.05
Q18_12 11 420 1.86 0.92 2 1.76 1.48 1 4 3 0.64 -0.76 0.05
Q18_13 12 420 2.79 0.87 3 2.84 1.48 1 4 3 -0.33 -0.56 0.04
Q18_14 13 420 2.76 0.87 3 2.82 1.48 1 4 3 -0.45 -0.39 0.04
Q18_15 14 420 2.24 0.91 2 2.19 1.48 1 4 3 0.21 -0.82 0.04
Q18_16 15 420 2.11 0.79 2 2.10 1.48 1 4 3 0.24 -0.50 0.04d <- na.omit(subset(df1, select=c(grep("Q18_", colnames(df1)))))
labels <- c(
Q18_1 = "Unfair Treatment Expected",
Q18_2 = "Information Untruthful",
Q18_3 = "Distrust from Experience",
Q18_4 = "Insincere Intentions",
Q18_5 = "Can't Trust Others",
Q18_6 = "Will Be Taken Advantage Of",
Q18_7 = "No One Would Help",
Q18_8 = "Others Out to Get Me",
Q18_9 = "Distrust Authority",
Q18_11 = "Officials Untrustworthy",
Q18_12 = "Treated Unjustly",
Q18_13 = "Prefer Self-Research",
Q18_14 = "Faith Leads to Hurt",
Q18_15 = "Question Why Told Things",
Q18_16 = "Ignore Others' Advice"
)
# shorten labels to max 40 characters (you can adjust the number)
# labels <- str_trunc(labels, width = 60, side = "right")
# replace column names
names(d) <- labels
d_long <- stack(d)
names(d_long) <- c("value", "variable")
ggplot(d_long, aes(x = value)) +
geom_histogram(binwidth = 1, fill = "steelblue", color = "white") +
facet_wrap(~ variable, labeller = label_wrap_gen(width = 30)) +
theme_minimal()
ev <- eigen(cor(d)) # get eigenvalues
ap <- parallel(subject=nrow(d),var=ncol(d),rep=100,cent=.05) # run the parallel analysis, gives us another perspective on how many factors should be used in the model
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) # creates the scree plot
plotnScree(nS) # shows us the scree plot, look for the elbows
EFA <- factanal(d, factors = 3, rotation = "promax")
print(EFA, digits=3, cutoff=.4, sort=T)
Call:
factanal(x = d, factors = 3, rotation = "promax")
Uniquenesses:
Unfair Treatment Expected Information Untruthful
0.546 0.648
Distrust from Experience Insincere Intentions
0.231 0.155
Can't Trust Others Will Be Taken Advantage Of
0.342 0.497
No One Would Help Others Out to Get Me
0.684 0.535
Distrust Authority Officials Untrustworthy
0.342 0.450
Treated Unjustly Prefer Self-Research
0.391 0.757
Faith Leads to Hurt Question Why Told Things
0.435 0.388
Ignore Others' Advice
0.692
Loadings:
Factor1 Factor2 Factor3
Information Untruthful 0.622
Others Out to Get Me 0.697
Distrust Authority 0.932
Officials Untrustworthy 0.821
Treated Unjustly 0.780
Distrust from Experience 0.926
Insincere Intentions 0.942
Can't Trust Others 0.632
Faith Leads to Hurt 0.714
Question Why Told Things 0.884
Ignore Others' Advice 0.571
Unfair Treatment Expected 0.481
Will Be Taken Advantage Of 0.433
No One Would Help
Prefer Self-Research 0.444
Factor1 Factor2 Factor3
SS loadings 3.557 2.454 1.898
Proportion Var 0.237 0.164 0.127
Cumulative Var 0.237 0.401 0.527
Factor Correlations:
Factor1 Factor2 Factor3
Factor1 1.000 0.683 0.721
Factor2 0.683 1.000 0.692
Factor3 0.721 0.692 1.000
Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 121.21 on 63 degrees of freedom.
The p-value is 1.48e-05 d <- subset(d, select = -c(`No One Would Help`))
ev <- eigen(cor(d)) # get eigenvalues
ap <- parallel(subject=nrow(d),var=ncol(d),rep=100,cent=.05) # run the parallel analysis, gives us another perspective on how many factors should be used in the model
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) # creates the scree plot
plotnScree(nS) # shows us the scree plot, look for the elbows
EFA <- factanal(d, factors = 3, rotation = "promax")
print(EFA, digits=3, cutoff=.4, sort=T)
Call:
factanal(x = d, factors = 3, rotation = "promax")
Uniquenesses:
Unfair Treatment Expected Information Untruthful
0.548 0.643
Distrust from Experience Insincere Intentions
0.230 0.155
Can't Trust Others Will Be Taken Advantage Of
0.342 0.498
Others Out to Get Me Distrust Authority
0.556 0.340
Officials Untrustworthy Treated Unjustly
0.432 0.393
Prefer Self-Research Faith Leads to Hurt
0.749 0.423
Question Why Told Things Ignore Others' Advice
0.413 0.688
Loadings:
Factor1 Factor2 Factor3
Information Untruthful 0.614
Others Out to Get Me 0.657
Distrust Authority 0.915
Officials Untrustworthy 0.823
Treated Unjustly 0.763
Distrust from Experience 0.919
Insincere Intentions 0.931
Can't Trust Others 0.628
Faith Leads to Hurt 0.722
Question Why Told Things 0.836
Ignore Others' Advice 0.567
Unfair Treatment Expected 0.470
Will Be Taken Advantage Of 0.431
Prefer Self-Research 0.451
Factor1 Factor2 Factor3
SS loadings 3.267 2.411 1.805
Proportion Var 0.233 0.172 0.129
Cumulative Var 0.233 0.406 0.535
Factor Correlations:
Factor1 Factor2 Factor3
Factor1 1.00 0.670 0.700
Factor2 0.67 1.000 0.679
Factor3 0.70 0.679 1.000
Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 85.57 on 52 degrees of freedom.
The p-value is 0.00231 # labels <- c(
# Q18_1 = "Unfair Treatment Expected",
# Q18_2 = "Information Untruthful",
# Q18_3 = "Distrust from Experience",
# Q18_4 = "Insincere Intentions",
# Q18_5 = "Can't Trust Others",
# Q18_6 = "Will Be Taken Advantage Of",
# Q18_8 = "Others Out to Get Me",
# Q18_9 = "Distrust Authority",
# Q18_10 = "Officials Untrustworthy",
# Q18_11 = "Treated Unjustly",
# Q18_13 = "Prefer Self-Research",
# Q18_14 = "Faith Leads to Hurt",
# Q18_15 = "Question Why Told Things",
# Q18_16 = "Ignore Others' Advice"
# )
d <- na.omit(subset(df2, select=c(grep("Q18_", colnames(df2)))))
# specify model: latent variables =~ observed indicators
model <- '
da =~ Q18_8 + Q18_9 + Q18_10 + Q18_11
ib =~ Q18_3 + Q18_4 + Q18_5 + Q18_6
dp =~ Q18_14 + Q18_15 + Q18_16 + Q18_13
'
# fit CFA
fit <- cfa(model, data = d, std.lv = TRUE)
# summary with fit indices and standardized loadings
summary(fit, fit.measures = TRUE, standardized = TRUE)lavaan 0.6-21 ended normally after 26 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 27
Number of observations 184
Model Test User Model:
Test statistic 132.881
Degrees of freedom 51
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 1278.294
Degrees of freedom 66
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.932
Tucker-Lewis Index (TLI) 0.913
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -2492.399
Loglikelihood unrestricted model (H1) -2425.958
Akaike (AIC) 5038.798
Bayesian (BIC) 5125.601
Sample-size adjusted Bayesian (SABIC) 5040.085
Root Mean Square Error of Approximation:
RMSEA 0.093
90 Percent confidence interval - lower 0.074
90 Percent confidence interval - upper 0.113
P-value H_0: RMSEA <= 0.050 0.000
P-value H_0: RMSEA >= 0.080 0.878
Standardized Root Mean Square Residual:
SRMR 0.059
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
da =~
Q18_8 0.594 0.060 9.900 0.000 0.594 0.680
Q18_9 0.810 0.061 13.295 0.000 0.810 0.841
Q18_10 0.784 0.072 10.872 0.000 0.784 0.729
Q18_11 0.802 0.068 11.742 0.000 0.802 0.771
ib =~
Q18_3 0.902 0.059 15.357 0.000 0.902 0.897
Q18_4 0.933 0.061 15.203 0.000 0.933 0.892
Q18_5 0.795 0.058 13.777 0.000 0.795 0.839
Q18_6 0.771 0.068 11.372 0.000 0.771 0.736
dp =~
Q18_14 0.764 0.066 11.627 0.000 0.764 0.792
Q18_15 0.555 0.064 8.606 0.000 0.555 0.625
Q18_16 0.577 0.063 9.165 0.000 0.577 0.657
Q18_13 0.485 0.070 6.890 0.000 0.485 0.518
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
da ~~
ib 0.711 0.047 15.232 0.000 0.711 0.711
dp 0.676 0.059 11.516 0.000 0.676 0.676
ib ~~
dp 0.755 0.047 15.932 0.000 0.755 0.755
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.Q18_8 0.410 0.049 8.368 0.000 0.410 0.538
.Q18_9 0.272 0.045 6.110 0.000 0.272 0.293
.Q18_10 0.542 0.068 7.960 0.000 0.542 0.468
.Q18_11 0.439 0.059 7.465 0.000 0.439 0.406
.Q18_3 0.197 0.031 6.339 0.000 0.197 0.195
.Q18_4 0.224 0.034 6.524 0.000 0.224 0.205
.Q18_5 0.267 0.034 7.753 0.000 0.267 0.297
.Q18_6 0.503 0.058 8.692 0.000 0.503 0.458
.Q18_14 0.346 0.057 6.031 0.000 0.346 0.372
.Q18_15 0.481 0.058 8.273 0.000 0.481 0.610
.Q18_16 0.438 0.055 8.016 0.000 0.438 0.568
.Q18_13 0.640 0.072 8.845 0.000 0.640 0.731
da 1.000 1.000 1.000
ib 1.000 1.000 1.000
dp 1.000 1.000 1.000# modindices(fit, sort. = T)
# basic path diagram
semPaths(
fit,
what = "std", # show standardized loadings
layout = "tree", # neat hierarchical layout
rotation = 2,
style = "lisrel", # clean style
residuals = FALSE, # hide residual arrows (optional)
intercepts = FALSE,
sizeMan = 6, # size of observed variables
sizeLat = 8, # size of latent variables
edge.label.cex = 0.9, # path label size
label.cex = 1.1, # node label size
nCharNodes = 0, # keep full variable names
mar = c(6, 6, 6, 6)
)
# model_1factor <- '
# oneidm =~ Q18_8 + Q18_3 + Q18_4 + Q18_5 + Q18_6 + Q18_14
# '
#
# fit_1factor <- cfa(model_1factor, data = df2, estimator = "ML")
# summary(fit_1factor, fit.measures = TRUE, standardized = TRUE)
model2 <- '
da =~ Q18_8 + Q18_9 + Q18_10 + Q18_11
ib =~ Q18_3 + Q18_4 + Q18_5 + Q18_6
'
# fit CFA
fit2 <- cfa(model2, data = d, std.lv = TRUE)
# summary with fit indices and standardized loadings
summary(fit2, fit.measures = TRUE, standardized = TRUE)lavaan 0.6-21 ended normally after 23 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 17
Number of observations 184
Model Test User Model:
Test statistic 62.767
Degrees of freedom 19
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 941.262
Degrees of freedom 28
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.952
Tucker-Lewis Index (TLI) 0.929
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -1646.279
Loglikelihood unrestricted model (H1) -1614.895
Akaike (AIC) 3326.557
Bayesian (BIC) 3381.211
Sample-size adjusted Bayesian (SABIC) 3327.368
Root Mean Square Error of Approximation:
RMSEA 0.112
90 Percent confidence interval - lower 0.082
90 Percent confidence interval - upper 0.143
P-value H_0: RMSEA <= 0.050 0.001
P-value H_0: RMSEA >= 0.080 0.959
Standardized Root Mean Square Residual:
SRMR 0.060
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
da =~
Q18_8 0.599 0.060 9.993 0.000 0.599 0.686
Q18_9 0.800 0.061 13.034 0.000 0.800 0.831
Q18_10 0.769 0.073 10.565 0.000 0.769 0.715
Q18_11 0.818 0.068 12.042 0.000 0.818 0.786
ib =~
Q18_3 0.902 0.059 15.309 0.000 0.902 0.896
Q18_4 0.940 0.061 15.380 0.000 0.940 0.899
Q18_5 0.794 0.058 13.735 0.000 0.794 0.837
Q18_6 0.761 0.068 11.148 0.000 0.761 0.726
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
da ~~
ib 0.713 0.046 15.353 0.000 0.713 0.713
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.Q18_8 0.404 0.049 8.282 0.000 0.404 0.530
.Q18_9 0.286 0.046 6.242 0.000 0.286 0.309
.Q18_10 0.565 0.070 8.041 0.000 0.565 0.489
.Q18_11 0.413 0.058 7.149 0.000 0.413 0.381
.Q18_3 0.199 0.032 6.233 0.000 0.199 0.196
.Q18_4 0.210 0.034 6.140 0.000 0.210 0.192
.Q18_5 0.268 0.035 7.712 0.000 0.268 0.299
.Q18_6 0.519 0.059 8.721 0.000 0.519 0.473
da 1.000 1.000 1.000
ib 1.000 1.000 1.000# Omega (default)
compRelSEM(fit)$da
Composite `da` is composed of observed variables:
Q18_8, Q18_9, Q18_10, Q18_11
True-score variance is represented by common factor(s):
da
Total variance of composite `da` determined from the unrestricted model.
The proportion attributable to "true" scores is its model-based estimate of reliability ("omega"):
[1] 0.851
$ib
Composite `ib` is composed of observed variables:
Q18_3, Q18_4, Q18_5, Q18_6
True-score variance is represented by common factor(s):
ib
Total variance of composite `ib` determined from the unrestricted model.
The proportion attributable to "true" scores is its model-based estimate of reliability ("omega"):
[1] 0.912
$dp
Composite `dp` is composed of observed variables:
Q18_14, Q18_15, Q18_16, Q18_13
True-score variance is represented by common factor(s):
dp
Total variance of composite `dp` determined from the unrestricted model.
The proportion attributable to "true" scores is its model-based estimate of reliability ("omega"):
[1] 0.747# Alpha
compRelSEM(fit, tau.eq = TRUE)$da
Composite `da` is composed of observed variables:
Q18_8, Q18_9, Q18_10, Q18_11
Total variance of composite `da` determined from the unrestricted model.
Coefficient alpha would be:
[1] 0.835
$ib
Composite `ib` is composed of observed variables:
Q18_3, Q18_4, Q18_5, Q18_6
Total variance of composite `ib` determined from the unrestricted model.
Coefficient alpha would be:
[1] 0.902
$dp
Composite `dp` is composed of observed variables:
Q18_14, Q18_15, Q18_16, Q18_13
Total variance of composite `dp` determined from the unrestricted model.
Coefficient alpha would be:
[1] 0.742m_list <- c("40","42","44","46","50","54","58","62","66","68","74","78","82","86","90")
# Identify target columns
target_cols <- grep(
paste(c(m_list), collapse = "|"),
names(df2),
value = TRUE
)
d <- subset(s2, select = target_cols)
d <- d[, grepl("_4$", colnames(d))]
labels <- c(
Q40_4 = "Love",
Q42_4 = "Narcissism",
Q44_4 = "Depression",
Q46_4 = "Intelligence",
Q50_4 = "Trauma1",
Q54_4 = "Trauma2",
Q58_4 = "Toxic People",
Q62_4 = "Psychiatric Symptoms",
Q66_4 = "Introversion",
Q68_4 = "Manipulativeness",
Q74_4 = "Attraction",
Q78_4 = "Childhood",
Q82_4 = "Friendship",
Q86_4 = "Anxiety & Emotion Regulation",
Q90_4 = "Family"
)
# shorten labels to max 40 characters (you can adjust the number)
# labels <- str_trunc(labels, width = 60, side = "right")
# replace column names
names(d) <- labels
d_long <- stack(d)
names(d_long) <- c("value", "variable")
desc <- data.frame(describe(d))
# Clean up and round
desc_table <- desc %>%
select(n, mean, sd, median, min, max, skew, kurtosis, se) %>%
round(2)
# Create table
desc_table %>%
kable(
format = "html",
caption = "Descriptive Statistics",
col.names = c("N", "Mean", "SD", "Median", "Min", "Max", "Skew", "Kurtosis", "SE")
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "left"
) %>%
column_spec(1, bold = TRUE)| N | Mean | SD | Median | Min | Max | Skew | Kurtosis | SE | |
|---|---|---|---|---|---|---|---|---|---|
| Love | 87 | 2.79 | 1.04 | 3 | 1 | 4 | -0.39 | -1.04 | 0.11 |
| Narcissism | 99 | 2.49 | 0.97 | 3 | 1 | 4 | -0.05 | -1.01 | 0.10 |
| Depression | 94 | 2.69 | 0.89 | 3 | 1 | 4 | -0.35 | -0.61 | 0.09 |
| Intelligence | 90 | 2.46 | 0.93 | 3 | 1 | 4 | -0.08 | -0.90 | 0.10 |
| Trauma1 | 99 | 2.91 | 0.89 | 3 | 1 | 4 | -0.67 | -0.18 | 0.09 |
| Trauma2 | 96 | 2.58 | 0.99 | 3 | 1 | 4 | -0.20 | -1.02 | 0.10 |
| Toxic People | 96 | 2.76 | 0.88 | 3 | 1 | 4 | -0.35 | -0.57 | 0.09 |
| Psychiatric Symptoms | 88 | 2.11 | 0.99 | 2 | 1 | 4 | 0.34 | -1.06 | 0.11 |
| Introversion | 96 | 1.92 | 1.03 | 2 | 1 | 4 | 0.73 | -0.77 | 0.11 |
| Manipulativeness | 93 | 2.16 | 1.08 | 2 | 1 | 4 | 0.41 | -1.16 | 0.11 |
| Attraction | 98 | 1.65 | 0.87 | 1 | 1 | 4 | 1.00 | -0.25 | 0.09 |
| Childhood | 92 | 2.20 | 0.95 | 2 | 1 | 4 | 0.14 | -1.12 | 0.10 |
| Friendship | 89 | 2.26 | 0.90 | 2 | 1 | 4 | 0.22 | -0.77 | 0.10 |
| Anxiety & Emotion Regulation | 92 | 2.33 | 0.94 | 2 | 1 | 4 | 0.03 | -1.00 | 0.10 |
| Family | 98 | 1.84 | 0.97 | 2 | 1 | 4 | 0.86 | -0.41 | 0.10 |
ggplot(d_long, aes(x = value)) +
geom_histogram(binwidth = 1, fill = "steelblue", color = "white") +
facet_wrap(~ variable, labeller = label_wrap_gen(width = 30)) +
theme_minimal()Warning: Removed 1353 rows containing non-finite outside the scale range
(`stat_bin()`).
a_list <- c("23","34","36","38","48","52","56","60","64","70","72","76","80","84","88")
# Identify target columns
target_cols <- grep(
paste(c(a_list), collapse = "|"),
names(df2),
value = TRUE
)
d <- subset(s2, select = target_cols)
d <- d[, grepl("_4$", colnames(d))]
describe(d) vars n mean sd median trimmed mad min max range skew kurtosis se
Q23_4 1 97 2.97 0.77 3 2.99 1.48 1 4 3 -0.22 -0.67 0.08
Q34_4 2 85 2.84 0.90 3 2.91 1.48 1 4 3 -0.46 -0.53 0.10
Q36_4 3 90 2.78 0.92 3 2.85 1.48 1 4 3 -0.41 -0.66 0.10
Q38_4 4 94 2.70 0.97 3 2.75 1.48 1 4 3 -0.22 -0.97 0.10
Q48_4 5 85 3.01 0.87 3 3.07 1.48 1 4 3 -0.46 -0.65 0.09
Q52_4 6 88 2.40 0.99 2 2.38 1.48 1 4 3 0.07 -1.07 0.11
Q56_4 7 88 2.95 0.82 3 3.00 1.48 1 4 3 -0.42 -0.39 0.09
Q60_4 8 96 3.14 0.79 3 3.19 1.48 1 4 3 -0.49 -0.56 0.08
Q64_4 9 88 2.36 1.00 2 2.33 1.48 1 4 3 0.13 -1.07 0.11
Q70_4 10 91 2.23 0.92 2 2.18 1.48 1 4 3 0.21 -0.88 0.10
Q72_4 11 86 2.31 1.05 2 2.27 1.48 1 4 3 0.07 -1.30 0.11
Q76_4 12 92 2.62 0.99 3 2.65 1.48 1 4 3 -0.20 -1.03 0.10
Q80_4 13 95 2.64 0.98 3 2.68 1.48 1 4 3 -0.19 -0.99 0.10
Q84_4 14 92 2.48 0.88 3 2.47 1.48 1 4 3 -0.12 -0.77 0.09
Q88_4 15 86 2.45 0.97 3 2.44 1.48 1 4 3 -0.06 -1.02 0.10labels <- c(
Q23_4 = "Love",
Q34_4 = "Narcissism",
Q36_4 = "Depression",
Q38_4 = "Intelligence",
Q48_4 = "Trauma1",
Q52_4 = "Trauma2",
Q56_4 = "Toxic People",
Q60_4 = "Psychiatric Symptoms",
Q64_4 = "Introversion",
Q70_4 = "Manipulativeness",
Q72_4 = "Attraction",
Q76_4 = "Childhood",
Q80_4 = "Friendship",
Q84_4 = "Anxiety & Emotion Regulation",
Q88_4 = "Family"
)
# shorten labels to max 40 characters (you can adjust the number)
# labels <- str_trunc(labels, width = 60, side = "right")
# replace column names
names(d) <- labels
d_long <- stack(d)
names(d_long) <- c("value", "variable")
desc <- data.frame(describe(d))
# Clean up and round
desc_table <- desc %>%
select(n, mean, sd, median, min, max, skew, kurtosis, se) %>%
round(2)
# Create table
desc_table %>%
kable(
format = "html",
caption = "Descriptive Statistics",
col.names = c("N", "Mean", "SD", "Median", "Min", "Max", "Skew", "Kurtosis", "SE")
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "left"
) %>%
column_spec(1, bold = TRUE)| N | Mean | SD | Median | Min | Max | Skew | Kurtosis | SE | |
|---|---|---|---|---|---|---|---|---|---|
| Love | 97 | 2.97 | 0.77 | 3 | 1 | 4 | -0.22 | -0.67 | 0.08 |
| Narcissism | 85 | 2.84 | 0.90 | 3 | 1 | 4 | -0.46 | -0.53 | 0.10 |
| Depression | 90 | 2.78 | 0.92 | 3 | 1 | 4 | -0.41 | -0.66 | 0.10 |
| Intelligence | 94 | 2.70 | 0.97 | 3 | 1 | 4 | -0.22 | -0.97 | 0.10 |
| Trauma1 | 85 | 3.01 | 0.87 | 3 | 1 | 4 | -0.46 | -0.65 | 0.09 |
| Trauma2 | 88 | 2.40 | 0.99 | 2 | 1 | 4 | 0.07 | -1.07 | 0.11 |
| Toxic People | 88 | 2.95 | 0.82 | 3 | 1 | 4 | -0.42 | -0.39 | 0.09 |
| Psychiatric Symptoms | 96 | 3.14 | 0.79 | 3 | 1 | 4 | -0.49 | -0.56 | 0.08 |
| Introversion | 88 | 2.36 | 1.00 | 2 | 1 | 4 | 0.13 | -1.07 | 0.11 |
| Manipulativeness | 91 | 2.23 | 0.92 | 2 | 1 | 4 | 0.21 | -0.88 | 0.10 |
| Attraction | 86 | 2.31 | 1.05 | 2 | 1 | 4 | 0.07 | -1.30 | 0.11 |
| Childhood | 92 | 2.62 | 0.99 | 3 | 1 | 4 | -0.20 | -1.03 | 0.10 |
| Friendship | 95 | 2.64 | 0.98 | 3 | 1 | 4 | -0.19 | -0.99 | 0.10 |
| Anxiety & Emotion Regulation | 92 | 2.48 | 0.88 | 3 | 1 | 4 | -0.12 | -0.77 | 0.09 |
| Family | 86 | 2.45 | 0.97 | 3 | 1 | 4 | -0.06 | -1.02 | 0.10 |
ggplot(d_long, aes(x = value)) +
geom_histogram(binwidth = 1, fill = "steelblue", color = "white") +
facet_wrap(~ variable, labeller = label_wrap_gen(width = 30)) +
theme_minimal()Warning: Removed 1407 rows containing non-finite outside the scale range
(`stat_bin()`).
mis_ICC <- ICC(subset(df2, select=c(Q40_4,Q42_4,Q44_4,Q46_4,Q50_4,Q54_4,Q58_4,Q62_4,Q66_4,Q68_4,Q74_4,Q78_4,Q82_4,Q86_4,Q90_4)), missing = TRUE)
mis_ICCCall: ICC(x = subset(df2, select = c(Q40_4, Q42_4, Q44_4, Q46_4, Q50_4,
Q54_4, Q58_4, Q62_4, Q66_4, Q68_4, Q74_4, Q78_4, Q82_4, Q86_4,
Q90_4)), missing = TRUE)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound
Single_raters_absolute ICC1 0.30 7.5 183 2576 1.5e-138 0.26
Single_random_raters ICC2 0.31 9.3 183 2562 1.4e-176 0.25
Single_fixed_raters ICC3 0.36 9.3 183 2562 1.4e-176 0.31
Average_raters_absolute ICC1k 0.87 7.5 183 2576 1.5e-138 0.84
Average_random_raters ICC2k 0.87 9.3 183 2562 1.4e-176 0.84
Average_fixed_raters ICC3k 0.89 9.3 183 2562 1.4e-176 0.87
upper bound
Single_raters_absolute 0.36
Single_random_raters 0.37
Single_fixed_raters 0.41
Average_raters_absolute 0.89
Average_random_raters 0.90
Average_fixed_raters 0.91
Number of subjects = 184 Number of Judges = 15
See the help file for a discussion of the other 4 McGraw and Wong estimates,acc_ICC <- ICC(subset(df2, select=c(Q23_4,Q34_4,Q36_4,Q38_4,Q48_4,Q52_4,Q56_4,Q60_4,Q64_4,Q70_4,Q72_4,Q76_4,Q80_4,Q84_4,Q88_4)), missing = TRUE)
acc_ICCCall: ICC(x = subset(df2, select = c(Q23_4, Q34_4, Q36_4, Q38_4, Q48_4,
Q52_4, Q56_4, Q60_4, Q64_4, Q70_4, Q72_4, Q76_4, Q80_4, Q84_4,
Q88_4)), missing = TRUE)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound
Single_raters_absolute ICC1 0.32 8.2 183 2576 3.4e-153 0.28
Single_random_raters ICC2 0.33 9.3 183 2562 2.1e-176 0.28
Single_fixed_raters ICC3 0.36 9.3 183 2562 2.1e-176 0.30
Average_raters_absolute ICC1k 0.88 8.2 183 2576 3.4e-153 0.85
Average_random_raters ICC2k 0.88 9.3 183 2562 2.1e-176 0.85
Average_fixed_raters ICC3k 0.89 9.3 183 2562 2.1e-176 0.87
upper bound
Single_raters_absolute 0.38
Single_random_raters 0.39
Single_fixed_raters 0.41
Average_raters_absolute 0.90
Average_random_raters 0.90
Average_fixed_raters 0.91
Number of subjects = 184 Number of Judges = 15
See the help file for a discussion of the other 4 McGraw and Wong estimates,mis_cred_ICC <- ICC(subset(df2, select=c(Q40_4,Q44_4,Q46_4,Q50_4,Q54_4,Q68_4,Q86_4)), missing = TRUE)
mis_cred_ICCCall: ICC(x = subset(df2, select = c(Q40_4, Q44_4, Q46_4, Q50_4, Q54_4,
Q68_4, Q86_4)), missing = TRUE)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.37 5.1 183 1104 2.1e-65 0.31 0.43
Single_random_raters ICC2 0.37 5.6 183 1098 2.9e-75 0.31 0.44
Single_fixed_raters ICC3 0.40 5.6 183 1098 2.9e-75 0.34 0.47
Average_raters_absolute ICC1k 0.80 5.1 183 1104 2.1e-65 0.76 0.84
Average_random_raters ICC2k 0.81 5.6 183 1098 2.9e-75 0.76 0.85
Average_fixed_raters ICC3k 0.82 5.6 183 1098 2.9e-75 0.78 0.86
Number of subjects = 184 Number of Judges = 7
See the help file for a discussion of the other 4 McGraw and Wong estimates,acc_cred_ICC <- ICC(subset(df2, select=c(Q23_4,Q34_4,Q36_4,Q38_4,Q48_4,Q70_4,Q84_4)), missing = TRUE)
acc_cred_ICCCall: ICC(x = subset(df2, select = c(Q23_4, Q34_4, Q36_4, Q38_4, Q48_4,
Q70_4, Q84_4)), missing = TRUE)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.32 4.4 183 1104 6.9e-53 0.26 0.39
Single_random_raters ICC2 0.33 4.9 183 1098 6.1e-63 0.27 0.40
Single_fixed_raters ICC3 0.36 4.9 183 1098 6.1e-63 0.30 0.43
Average_raters_absolute ICC1k 0.77 4.4 183 1104 6.9e-53 0.72 0.82
Average_random_raters ICC2k 0.78 4.9 183 1098 6.1e-63 0.72 0.82
Average_fixed_raters ICC3k 0.80 4.9 183 1098 6.1e-63 0.75 0.84
Number of subjects = 184 Number of Judges = 7
See the help file for a discussion of the other 4 McGraw and Wong estimates,mis_dang_ICC <- ICC(subset(df2, select=c(Q42_4,Q46_4,Q50_4,Q54_4,Q58_4,Q62_4,Q78_4,Q90_4)), missing = TRUE)
mis_dang_ICCCall: ICC(x = subset(df2, select = c(Q42_4, Q46_4, Q50_4, Q54_4, Q58_4,
Q62_4, Q78_4, Q90_4)), missing = TRUE)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.27 4.0 183 1288 7.8e-49 0.22 0.34
Single_random_raters ICC2 0.28 4.9 183 1281 6.5e-65 0.22 0.36
Single_fixed_raters ICC3 0.33 4.9 183 1281 6.5e-65 0.27 0.39
Average_raters_absolute ICC1k 0.75 4.0 183 1288 7.8e-49 0.69 0.80
Average_random_raters ICC2k 0.76 4.9 183 1281 6.5e-65 0.69 0.82
Average_fixed_raters ICC3k 0.79 4.9 183 1281 6.5e-65 0.75 0.84
Number of subjects = 184 Number of Judges = 8
See the help file for a discussion of the other 4 McGraw and Wong estimates,acc_dang_ICC <- ICC(subset(df2, select=c(Q38_4,Q48_4,Q52_4,Q56_4,Q76_4,Q88_4)), missing = TRUE)
acc_cred_ICCCall: ICC(x = subset(df2, select = c(Q23_4, Q34_4, Q36_4, Q38_4, Q48_4,
Q70_4, Q84_4)), missing = TRUE)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.32 4.4 183 1104 6.9e-53 0.26 0.39
Single_random_raters ICC2 0.33 4.9 183 1098 6.1e-63 0.27 0.40
Single_fixed_raters ICC3 0.36 4.9 183 1098 6.1e-63 0.30 0.43
Average_raters_absolute ICC1k 0.77 4.4 183 1104 6.9e-53 0.72 0.82
Average_random_raters ICC2k 0.78 4.9 183 1098 6.1e-63 0.72 0.82
Average_fixed_raters ICC3k 0.80 4.9 183 1098 6.1e-63 0.75 0.84
Number of subjects = 184 Number of Judges = 7
See the help file for a discussion of the other 4 McGraw and Wong estimates,df1 <- df1 %>%
mutate(stress_helpless = rowMeans(across(c(Q6_2,Q6_3,Q6_6,Q6_9,Q6_10)), na.rm = TRUE)) %>%
mutate(stress_efficacy = rowMeans(across(c(Q6_4,Q6_4,Q6_7,Q6_8)), na.rm = TRUE)) %>%
mutate(discrim = rowMeans(across(c(Q39_1,Q39_2,Q39_4,Q39_8)), na.rm = TRUE)) %>%
mutate(reject = rowMeans(across(c(Q40_1,Q40_2,Q40_4,Q40_5,Q40_6)), na.rm = TRUE)) %>%
mutate(idm_dp = rowMeans(across(c(Q18_13,Q18_14,Q18_15,Q18_16)), na.rm = TRUE)) %>%
mutate(idm_da = rowMeans(across(c(Q18_2,Q18_9,Q18_11,Q18_12,Q18_1,Q18_8)), na.rm = TRUE)) %>%
mutate(idm_ib = rowMeans(across(c(Q18_3,Q18_4,Q18_5,Q18_6)), na.rm = TRUE)) %>%
mutate(idm = rowMeans(across(c(idm_dp,idm_da,idm_ib)), na.rm = TRUE)) %>%
mutate(idm2 = rowMeans(across(c(Q18_3,Q18_4,Q18_5,Q18_6,Q18_8,Q18_14)), na.rm = TRUE))
df2 <- df2 %>%
mutate(stress_helpless = rowMeans(across(c(Q6_2,Q6_3,Q6_6,Q6_9,Q6_10)), na.rm = TRUE)) %>%
mutate(stress_efficacy = rowMeans(across(c(Q6_4,Q6_4,Q6_7,Q6_8)), na.rm = TRUE)) %>%
mutate(discrim = rowMeans(across(c(Q7_1,Q7_2,Q7_4,Q7_8)), na.rm = TRUE)) %>%
mutate(reject = rowMeans(across(c(Q8_1,Q8_2,Q8_4,Q8_5,Q8_6)), na.rm = TRUE)) %>%
mutate(idm_dp = rowMeans(across(c(Q18_13,Q18_14,Q18_15,Q18_16)), na.rm = TRUE)) %>%
mutate(idm_da = rowMeans(across(c(Q18_2,Q18_9,Q18_10,Q18_11,Q18_1,Q18_2)), na.rm = TRUE)) %>%
mutate(idm_ib = rowMeans(across(c(Q18_3,Q18_4,Q18_5,Q18_6)), na.rm = TRUE)) %>%
mutate(idm = rowMeans(across(c(idm_dp,idm_da,idm_ib)), na.rm = TRUE)) %>%
mutate(idm2 = rowMeans(across(c(Q18_3,Q18_4,Q18_5,Q18_6,Q18_8,Q18_14)), na.rm = TRUE)) %>%
mutate(mis = rowMeans(across(c(Q40_4,Q42_4,Q44_4,Q46_4,Q50_4,Q54_4,Q58_4,Q62_4,Q66_4,Q68_4,Q74_4,Q78_4,Q82_4,Q86_4,Q90_4)), na.rm = TRUE)) %>%
mutate(acc = rowMeans(across(c(Q23_4,Q34_4,Q36_4,Q38_4,Q48_4,Q52_4,Q56_4,Q60_4,Q64_4,Q70_4,Q72_4,Q76_4,Q80_4,Q84_4,Q88_4)), na.rm = TRUE)) %>%
mutate(mis_cred = rowMeans(across(c(Q40_4,Q44_4,Q46_4,Q50_4,Q54_4,Q68_4,Q86_4)), na.rm = TRUE)) %>%
mutate(acc_cred = rowMeans(across(c(Q23_4,Q34_4,Q36_4,Q38_4,Q48_4,Q84_4)), na.rm = TRUE)) %>%
mutate(mis_dang = rowMeans(across(c(Q42_4,Q46_4,Q50_4,Q54_4,Q58_4,Q62_4,Q78_4,Q90_4)), na.rm = TRUE)) %>%
mutate(acc_dang = rowMeans(across(c(Q38_4,Q48_4,Q52_4,Q56_4,Q76_4,Q88_4)), na.rm = TRUE))
# Standardize composite variables in df1
df1 <- df1 %>%
mutate(across(c(stress_helpless, stress_efficacy, discrim, reject,
idm_dp, idm_da, idm_ib, idm, idm2),
~ as.numeric(scale(.)),
.names = "{.col}_z"))
# Standardize composite variables in df2
df2 <- df2 %>%
mutate(across(c(stress_helpless, stress_efficacy, discrim, reject,
idm_dp, idm_da, idm_ib, idm, idm2,
mis, acc, mis_cred, acc_cred, mis_dang, acc_dang),
~ as.numeric(scale(.)),
.names = "{.col}_z"))
s1 <- s1 %>%
mutate(
p1_edu = ifelse(Q36 == 6, NA, Q36),
p2_edu = ifelse(Q37 == 6, NA, Q37),
parent_edu = rowMeans(across(c(p1_edu, p2_edu)), na.rm = TRUE)
)
df1 <- df1 %>%
mutate(parent_edu = s1$parent_edu)
s2 <- s2 %>%
mutate(
p1_edu = ifelse(Q7 == 6, NA, Q7),
p2_edu = ifelse(Q8 == 6, NA, Q8),
parent_edu = rowMeans(across(c(p1_edu, p2_edu)), na.rm = TRUE)
)
df2 <- df2 %>%
mutate(parent_edu = s2$parent_edu)corrout1 <- corr.test(subset(df1, select=c(60:72)))
corrplot(
corrout1$r,
p.mat = corrout1$p, # add p-values
sig.level = 0.05, # hide correlations above this p-value
insig = "pch", # or "pch" to mark nonsignificant ones
method = "color",
type = "upper",
tl.col = "black",
tl.srt = 45,
addCoef.col = "black",
number.cex = 0.8)
corrout2 <- corr.test(subset(df2, select=c(216:242)))
corrplot(
corrout2$r,
p.mat = corrout2$p, # add p-values
sig.level = 0.05, # hide correlations above this p-value
insig = "pch", # or "pch" to mark nonsignificant ones
method = "color",
type = "upper",
tl.col = "black",
tl.srt = 45,
addCoef.col = "black",
number.cex = 0.8)
t.test(parent_edu ~ samp, data = df1)
Welch Two Sample t-test
data: parent_edu by samp
t = -12.26, df = 368.52, p-value < 2.2e-16
alternative hypothesis: true difference in means between group Prolific and group SONA is not equal to 0
95 percent confidence interval:
-1.3390379 -0.9688569
sample estimates:
mean in group Prolific mean in group SONA
2.546053 3.700000 t.test(parent_edu ~ samp, data = df2)
Welch Two Sample t-test
data: parent_edu by samp
t = -8.0798, df = 83.629, p-value = 4.312e-12
alternative hypothesis: true difference in means between group Prolific and group SONA is not equal to 0
95 percent confidence interval:
-1.5966197 -0.9658899
sample estimates:
mean in group Prolific mean in group SONA
2.577236 3.858491 H1 Stress, experiences with discrimination, and rejection sensitivity will predict inequality-driven mistrust
reg4 <- lm(data=df1, idm_z ~ stress_helpless_z + stress_efficacy_z + discrim_z + reject_z + samp + parent_edu)
car::vif(reg4)stress_helpless_z stress_efficacy_z discrim_z reject_z
1.745736 1.553836 1.497610 1.434414
samp parent_edu
1.518453 1.386463 plot_model(reg4, type="diag")[[1]]
[[2]]
[[3]]
[[4]]
plot(reg4, 5)
summary(reg4)
Call:
lm(formula = idm_z ~ stress_helpless_z + stress_efficacy_z +
discrim_z + reject_z + samp + parent_edu, data = df1)
Residuals:
Min 1Q Median 3Q Max
-2.18235 -0.48041 -0.06168 0.51653 2.56221
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.239343 0.115039 2.081 0.0381 *
stress_helpless_z 0.116206 0.049446 2.350 0.0192 *
stress_efficacy_z -0.092525 0.046666 -1.983 0.0481 *
discrim_z 0.316745 0.045863 6.906 1.9e-11 ***
reject_z 0.243337 0.044702 5.444 9.0e-08 ***
sampSONA -0.496822 0.092474 -5.373 1.3e-07 ***
parent_edu -0.003748 0.039925 -0.094 0.9253
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.764 on 411 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.4254, Adjusted R-squared: 0.417
F-statistic: 50.72 on 6 and 411 DF, p-value: < 2.2e-16H2 Hypothesis 1 will be replicated with a new sample and in a mini meta-analysis
reg8 <- lm(data=df2, idm_z ~ stress_helpless_z + stress_efficacy_z + discrim_z + reject_z + samp + parent_edu)
car::vif(reg8)stress_helpless_z stress_efficacy_z discrim_z reject_z
2.390220 1.831291 1.763868 1.517440
samp parent_edu
1.629979 1.444221 plot_model(reg8, type="diag")[[1]]
[[2]]
[[3]]
[[4]]
plot(reg8, 5)
summary(reg8)
Call:
lm(formula = idm_z ~ stress_helpless_z + stress_efficacy_z +
discrim_z + reject_z + samp + parent_edu, data = df2)
Residuals:
Min 1Q Median 3Q Max
-1.83626 -0.48588 0.00287 0.41192 2.81892
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.25310 0.17337 1.460 0.14617
stress_helpless_z 0.22915 0.08404 2.727 0.00707 **
stress_efficacy_z -0.05298 0.07364 -0.719 0.47290
discrim_z 0.05549 0.07355 0.754 0.45164
reject_z 0.47862 0.06928 6.909 9.53e-11 ***
sampSONA -0.46278 0.15358 -3.013 0.00298 **
parent_edu -0.03630 0.06229 -0.583 0.56079
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7321 on 169 degrees of freedom
(8 observations deleted due to missingness)
Multiple R-squared: 0.483, Adjusted R-squared: 0.4647
F-statistic: 26.32 on 6 and 169 DF, p-value: < 2.2e-16# Extract coefficients and SEs from both models
reg_s1 <- tidy(reg4) %>%
filter(term %in% c("discrim_z", "reject_z", "ders_strat_z", "stress_helpless_z", "stress_efficacy_z", "ders_clarity_z", "ders_goals_z", "ders_impulse_z", "ders_nonacc_z", "sampSONA", "parent_edu")) %>%
select(term, estimate, std.error) %>%
rename(b1 = estimate, se1 = std.error)
reg_s2 <- tidy(reg8) %>%
filter(term %in% c("discrim_z", "reject_z", "ders_strat_z", "stress_helpless_z", "stress_efficacy_z", "ders_clarity_z", "ders_goals_z", "ders_impulse_z", "ders_nonacc_z", "sampSONA", "parent_edu")) %>%
select(term, estimate, std.error) %>%
rename(b2 = estimate, se2 = std.error)
# Join them together
meta_inputs <- left_join(reg_s1, reg_s2, by = "term")
# Preview
print(meta_inputs)# A tibble: 6 × 5
term b1 se1 b2 se2
<chr> <dbl> <dbl> <dbl> <dbl>
1 stress_helpless_z 0.116 0.0494 0.229 0.0840
2 stress_efficacy_z -0.0925 0.0467 -0.0530 0.0736
3 discrim_z 0.317 0.0459 0.0555 0.0735
4 reject_z 0.243 0.0447 0.479 0.0693
5 sampSONA -0.497 0.0925 -0.463 0.154
6 parent_edu -0.00375 0.0399 -0.0363 0.0623append_meta <- function(inputs, term_name, meta_obj) {
inputs[inputs$term == term_name, "tau2"] <- meta_obj$tau2
inputs[inputs$term == term_name, "se.tau2"] <- meta_obj$se.tau2
inputs[inputs$term == term_name, "I2"] <- meta_obj$I2
inputs[inputs$term == term_name, "QE"] <- meta_obj$QE
inputs[inputs$term == term_name, "QEp"] <- meta_obj$QEp
inputs[inputs$term == term_name, "beta"] <- as.numeric(meta_obj$beta)
inputs[inputs$term == term_name, "SE"] <- as.numeric(meta_obj$se)
inputs[inputs$term == term_name, "pval"] <- as.numeric(meta_obj$pval)
inputs[inputs$term == term_name, "ci.lb"] <- as.numeric(meta_obj$ci.lb)
inputs[inputs$term == term_name, "ci.ub"] <- as.numeric(meta_obj$ci.ub)
return(inputs)
}
dat <- data.frame(
yi = c(meta_inputs$b1[meta_inputs$term == "stress_helpless_z"],
meta_inputs$b2[meta_inputs$term == "stress_helpless_z"]),
sei = c(meta_inputs$se1[meta_inputs$term == "stress_helpless_z"],
meta_inputs$se2[meta_inputs$term == "stress_helpless_z"])
)
dat$vi <- dat$sei^2
meta <- rma(yi, vi, data = dat, method = "REML")
summary(meta)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
1.1085 -2.2170 1.7830 -2.2170 13.7830
tau^2 (estimated amount of total heterogeneity): 0.0016 (SE = 0.0090)
tau (square root of estimated tau^2 value): 0.0403
I^2 (total heterogeneity / total variability): 25.47%
H^2 (total variability / sampling variability): 1.34
Test for Heterogeneity:
Q(df = 1) = 1.3417, p-val = 0.2467
Model Results:
estimate se zval pval ci.lb ci.ub
0.1522 0.0526 2.8919 0.0038 0.0491 0.2554 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs <- append_meta(meta_inputs, "stress_helpless_z", meta)
dat <- data.frame(
yi = c(meta_inputs$b1[meta_inputs$term == "stress_efficacy_z"],
meta_inputs$b2[meta_inputs$term == "stress_efficacy_z"]),
sei = c(meta_inputs$se1[meta_inputs$term == "stress_efficacy_z"],
meta_inputs$se2[meta_inputs$term == "stress_efficacy_z"])
)
dat$vi <- dat$sei^2
meta <- rma(yi, vi, data = dat, method = "REML")
summary(meta)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
1.7645 -3.5290 0.4710 -3.5290 12.4710
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0054)
tau (square root of estimated tau^2 value): 0
I^2 (total heterogeneity / total variability): 0.00%
H^2 (total variability / sampling variability): 1.00
Test for Heterogeneity:
Q(df = 1) = 0.2058, p-val = 0.6501
Model Results:
estimate se zval pval ci.lb ci.ub
-0.0812 0.0394 -2.0598 0.0394 -0.1585 -0.0039 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs <- append_meta(meta_inputs, "stress_efficacy_z", meta)
dat <- data.frame(
yi = c(meta_inputs$b1[meta_inputs$term == "discrim_z"],
meta_inputs$b2[meta_inputs$term == "discrim_z"]),
sei = c(meta_inputs$se1[meta_inputs$term == "discrim_z"],
meta_inputs$se2[meta_inputs$term == "discrim_z"])
)
dat$vi <- dat$sei^2
meta <- rma(yi, vi, data = dat, method = "REML")
summary(meta)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
0.2699 -0.5398 3.4602 -0.5398 15.4602
tau^2 (estimated amount of total heterogeneity): 0.0304 (SE = 0.0483)
tau (square root of estimated tau^2 value): 0.1743
I^2 (total heterogeneity / total variability): 88.99%
H^2 (total variability / sampling variability): 9.09
Test for Heterogeneity:
Q(df = 1) = 9.0851, p-val = 0.0026
Model Results:
estimate se zval pval ci.lb ci.ub
0.1924 0.1305 1.4749 0.1402 -0.0633 0.4482
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs <- append_meta(meta_inputs, "discrim_z", meta)
dat <- data.frame(
yi = c(meta_inputs$b1[meta_inputs$term == "reject_z"],
meta_inputs$b2[meta_inputs$term == "reject_z"]),
sei = c(meta_inputs$se1[meta_inputs$term == "reject_z"],
meta_inputs$se2[meta_inputs$term == "reject_z"])
)
dat$vi <- dat$sei^2
meta <- rma(yi, vi, data = dat, method = "REML")
summary(meta)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
0.3746 -0.7492 3.2508 -0.7492 15.2508
tau^2 (estimated amount of total heterogeneity): 0.0243 (SE = 0.0391)
tau (square root of estimated tau^2 value): 0.1558
I^2 (total heterogeneity / total variability): 87.72%
H^2 (total variability / sampling variability): 8.14
Test for Heterogeneity:
Q(df = 1) = 8.1432, p-val = 0.0043
Model Results:
estimate se zval pval ci.lb ci.ub
0.3550 0.1175 3.0218 0.0025 0.1247 0.5853 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs <- append_meta(meta_inputs, "reject_z", meta)
dat <- data.frame(
yi = c(meta_inputs$b1[meta_inputs$term == "sampSONA"],
meta_inputs$b2[meta_inputs$term == "sampSONA"]),
sei = c(meta_inputs$se1[meta_inputs$term == "sampSONA"],
meta_inputs$se2[meta_inputs$term == "sampSONA"])
)
dat$vi <- dat$sei^2
meta <- rma(yi, vi, data = dat, method = "REML")
summary(meta)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
1.1285 -2.2569 1.7431 -2.2569 13.7431
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0227)
tau (square root of estimated tau^2 value): 0
I^2 (total heterogeneity / total variability): 0.00%
H^2 (total variability / sampling variability): 1.00
Test for Heterogeneity:
Q(df = 1) = 0.0361, p-val = 0.8494
Model Results:
estimate se zval pval ci.lb ci.ub
-0.4878 0.0792 -6.1569 <.0001 -0.6430 -0.3325 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs <- append_meta(meta_inputs, "sampSONA", meta)
dat <- data.frame(
yi = c(meta_inputs$b1[meta_inputs$term == "parent_edu"],
meta_inputs$b2[meta_inputs$term == "parent_edu"]),
sei = c(meta_inputs$se1[meta_inputs$term == "parent_edu"],
meta_inputs$se2[meta_inputs$term == "parent_edu"])
)
dat$vi <- dat$sei^2
meta <- rma(yi, vi, data = dat, method = "REML")
summary(meta)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
1.9347 -3.8694 0.1306 -3.8694 12.1306
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0039)
tau (square root of estimated tau^2 value): 0
I^2 (total heterogeneity / total variability): 0.00%
H^2 (total variability / sampling variability): 1.00
Test for Heterogeneity:
Q(df = 1) = 0.1936, p-val = 0.6599
Model Results:
estimate se zval pval ci.lb ci.ub
-0.0132 0.0336 -0.3935 0.6939 -0.0791 0.0527
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs <- append_meta(meta_inputs, "parent_edu", meta)term_labels <- c(
"stress_helpless_z" = "Stress (Helplessness)",
"stress_efficacy_z" = "Stress (Efficacy)",
"discrim_z" = "Discrimination",
"reject_z" = "Rejection Sensitivity",
"ders_clarity_z" = "ER: Clarity",
"ders_goals_z" = "ER: Goals",
"ders_impulse_z" = "ER: Impulse",
"ders_strat_z" = "ER: Strategies",
"ders_nonacc_z" = "ER: Non-Acceptance",
"sampSONA" = "Student Status",
"parent_edu" = "Parental Education"
)
sig_rows <- which(meta_inputs$pval < .05)
trend_rows <- which(meta_inputs$pval >= .05 & meta_inputs$pval < .11)
replication_rows <- c(1) # <-- specify row numbers here
meta_inputs %>%
mutate(term = term_labels[term]) %>%
mutate(across(c(b1, se1, b2, se2, tau2, se.tau2, I2, QE, beta, SE, ci.lb, ci.ub),
~ round(., 2))) %>%
mutate(across(c(QEp, pval), ~ round(., 3))) %>%
mutate(across(where(is.numeric), ~ round(., 3))) %>%
rename(
"Predictor" = term,
"β (S1)" = b1,
"SE (S1)" = se1,
"β (S2)" = b2,
"SE (S2)" = se2,
"τ²" = tau2,
"SE(τ²)" = se.tau2,
"I²" = I2,
"Q" = QE,
"Q p" = QEp,
"β (pooled)" = beta,
"SE (pooled)" = SE,
"p" = pval,
"95% CI LL" = ci.lb,
"95% CI UL" = ci.ub
) %>%
kable(format = "html", align = "c") %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = T,
font_size = 10
) %>%
add_header_above(c(
" " = 5,
"Heterogeneity" = 5,
"Pooled Effect" = 5
)) %>%
column_spec(1, bold = TRUE) %>%
column_spec(11, italic = TRUE) %>%
row_spec(0, bold = TRUE) %>%
row_spec(sig_rows, bold = TRUE, background = "#d4edda") %>% # green for significant
row_spec(trend_rows, bold = FALSE, background = "#fff3cd") %>%
row_spec(replication_rows, bold = FALSE, background = "#f8d7da")| Predictor | β (S1) | SE (S1) | β (S2) | SE (S2) | τ² | SE(τ²) | I² | Q | Q p | β (pooled) | SE (pooled) | p | 95% CI LL | 95% CI UL |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Stress (Helplessness) | 0.12 | 0.05 | 0.23 | 0.08 | 0.00 | 0.01 | 25.47 | 1.34 | 0.247 | 0.15 | 0.05 | 0.004 | 0.05 | 0.26 |
| Stress (Efficacy) | -0.09 | 0.05 | -0.05 | 0.07 | 0.00 | 0.01 | 0.00 | 0.21 | 0.650 | -0.08 | 0.04 | 0.039 | -0.16 | 0.00 |
| Discrimination | 0.32 | 0.05 | 0.06 | 0.07 | 0.03 | 0.05 | 88.99 | 9.09 | 0.003 | 0.19 | 0.13 | 0.140 | -0.06 | 0.45 |
| Rejection Sensitivity | 0.24 | 0.04 | 0.48 | 0.07 | 0.02 | 0.04 | 87.72 | 8.14 | 0.004 | 0.36 | 0.12 | 0.003 | 0.12 | 0.59 |
| Student Status | -0.50 | 0.09 | -0.46 | 0.15 | 0.00 | 0.02 | 0.00 | 0.04 | 0.849 | -0.49 | 0.08 | 0.000 | -0.64 | -0.33 |
| Parental Education | 0.00 | 0.04 | -0.04 | 0.06 | 0.00 | 0.00 | 0.00 | 0.19 | 0.660 | -0.01 | 0.03 | 0.694 | -0.08 | 0.05 |
H3 Inequality-driven mistrust will predict misinformation acceptance, accurate information acceptance, and increased susceptibility to misinformation
IDM predicts misinformation acceptance (b = .14, p = .049) but not accurate information acceptance (b = .04, p = .549). Students are lower in misinformation (b = -.65, p < .001) and accurate information (b = -.71, p < .001) acceptance. IDM also predicts increased susceptability to misinformation acceptance (using residualized change approach; b = .11, p = .045). Students have lower susceptability to misinformation than non-students, but the difference is only borderline significant (b = -.23, p = .080).
IDM predicts misinformation acceptance and increased susceptability to misinformation, but not accurate information acceptance.
reg1 <- lm(data = df2, mis_z ~ idm_z + samp)
car::vif(reg1) idm_z samp
1.043019 1.043019 plot_model(reg1, type="diag")[[1]]
[[2]]`geom_smooth()` using formula = 'y ~ x'
[[3]]
[[4]]`geom_smooth()` using formula = 'y ~ x'
plot(reg1, 5)
summary(reg1)
Call:
lm(formula = mis_z ~ idm_z + samp, data = df2)
Residuals:
Min 1Q Median 3Q Max
-2.0573 -0.6790 -0.1536 0.6758 2.3577
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.20620 0.08413 2.451 0.0152 *
idm_z 0.13992 0.07081 1.976 0.0497 *
sampSONA -0.65416 0.15201 -4.304 2.75e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.938 on 181 degrees of freedom
Multiple R-squared: 0.1298, Adjusted R-squared: 0.1202
F-statistic: 13.5 on 2 and 181 DF, p-value: 3.441e-06reg1 <- lm(data = df2, mis_cred_z ~ idm_z + samp)
car::vif(reg1) idm_z samp
1.034972 1.034972 plot_model(reg1, type="diag")[[1]]
[[2]]`geom_smooth()` using formula = 'y ~ x'
[[3]]
[[4]]`geom_smooth()` using formula = 'y ~ x'
plot(reg1, 5)
summary(reg1)
Call:
lm(formula = mis_cred_z ~ idm_z + samp, data = df2)
Residuals:
Min 1Q Median 3Q Max
-2.28894 -0.66733 0.01516 0.80209 2.21841
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.20748 0.08632 2.404 0.0173 *
idm_z 0.13004 0.07191 1.808 0.0723 .
sampSONA -0.66705 0.15554 -4.289 2.99e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9389 on 172 degrees of freedom
(9 observations deleted due to missingness)
Multiple R-squared: 0.1285, Adjusted R-squared: 0.1184
F-statistic: 12.69 on 2 and 172 DF, p-value: 7.265e-06reg1 <- lm(data = df2, mis_dang_z ~ idm_z + samp)
car::vif(reg1) idm_z samp
1.047961 1.047961 plot_model(reg1, type="diag")[[1]]
[[2]]`geom_smooth()` using formula = 'y ~ x'
[[3]]
[[4]]`geom_smooth()` using formula = 'y ~ x'
plot(reg1, 5)
summary(reg1)
Call:
lm(formula = mis_dang_z ~ idm_z + samp, data = df2)
Residuals:
Min 1Q Median 3Q Max
-2.25962 -0.62586 -0.05038 0.52366 2.10238
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.17786 0.08530 2.085 0.038479 *
idm_z 0.15994 0.07201 2.221 0.027597 *
sampSONA -0.56356 0.15501 -3.636 0.000362 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9474 on 179 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.1123, Adjusted R-squared: 0.1024
F-statistic: 11.32 on 2 and 179 DF, p-value: 2.346e-05reg2 <- lm(data = df2, acc_z ~ idm_z + samp)
car::vif(reg2) idm_z samp
1.043019 1.043019 plot_model(reg2, type="diag")[[1]]
[[2]]`geom_smooth()` using formula = 'y ~ x'
[[3]]
[[4]]`geom_smooth()` using formula = 'y ~ x'
plot(reg2, 5)
summary(reg2)
Call:
lm(formula = acc_z ~ idm_z + samp, data = df2)
Residuals:
Min 1Q Median 3Q Max
-2.7085 -0.6091 0.1030 0.6058 1.9538
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.22461 0.08471 2.652 0.00872 **
idm_z 0.04272 0.07130 0.599 0.54982
sampSONA -0.71257 0.15305 -4.656 6.22e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9444 on 181 degrees of freedom
Multiple R-squared: 0.1178, Adjusted R-squared: 0.108
F-statistic: 12.08 on 2 and 181 DF, p-value: 1.187e-05reg2 <- lm(data = df2, acc_cred_z ~ idm_z + samp)
car::vif(reg2) idm_z samp
1.043019 1.043019 plot_model(reg2, type="diag")[[1]]
[[2]]`geom_smooth()` using formula = 'y ~ x'
[[3]]
[[4]]`geom_smooth()` using formula = 'y ~ x'
plot(reg2, 5)
summary(reg2)
Call:
lm(formula = acc_cred_z ~ idm_z + samp, data = df2)
Residuals:
Min 1Q Median 3Q Max
-2.2940 -0.6899 0.2544 0.5986 1.9294
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.09294 0.08916 1.042 0.2986
idm_z 0.03951 0.07505 0.526 0.5992
sampSONA -0.29484 0.16109 -1.830 0.0689 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9941 on 181 degrees of freedom
Multiple R-squared: 0.02263, Adjusted R-squared: 0.01183
F-statistic: 2.096 on 2 and 181 DF, p-value: 0.126reg2 <- lm(data = df2, acc_dang_z ~ idm_z + samp)
car::vif(reg2) idm_z samp
1.0322 1.0322 plot_model(reg2, type="diag")[[1]]
[[2]]`geom_smooth()` using formula = 'y ~ x'
[[3]]
[[4]]`geom_smooth()` using formula = 'y ~ x'
plot(reg2, 5)
summary(reg2)
Call:
lm(formula = acc_dang_z ~ idm_z + samp, data = df2)
Residuals:
Min 1Q Median 3Q Max
-2.2269 -0.6101 0.0968 0.7245 2.1597
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.209512 0.090508 2.315 0.0219 *
idm_z 0.001824 0.075484 0.024 0.9808
sampSONA -0.692930 0.165699 -4.182 4.76e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9535 on 159 degrees of freedom
(22 observations deleted due to missingness)
Multiple R-squared: 0.1021, Adjusted R-squared: 0.09085
F-statistic: 9.044 on 2 and 159 DF, p-value: 0.0001905Uses residualized change approach (Castro-Schilo & Grimm, 2018).
reg3 <- lm(data = df2, mis_z ~ acc_z + idm_z + samp)
car::vif(reg3) acc_z idm_z samp
1.133511 1.045087 1.167931 plot_model(reg3, type="diag")[[1]]
[[2]]`geom_smooth()` using formula = 'y ~ x'
[[3]]
[[4]]`geom_smooth()` using formula = 'y ~ x'
plot(reg3, 5)
summary(reg3)
Call:
lm(formula = mis_z ~ acc_z + idm_z + samp, data = df2)
Residuals:
Min 1Q Median 3Q Max
-2.06177 -0.52635 -0.03887 0.50343 1.87752
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.07146 0.06853 1.043 0.2984
acc_z 0.59987 0.05900 10.167 <2e-16 ***
idm_z 0.11430 0.05665 2.018 0.0451 *
sampSONA -0.22672 0.12855 -1.764 0.0795 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7497 on 180 degrees of freedom
Multiple R-squared: 0.4472, Adjusted R-squared: 0.438
F-statistic: 48.55 on 3 and 180 DF, p-value: < 2.2e-16reg3 <- lm(data = df2, mis_cred_z ~ acc_cred_z + idm_z + samp)
car::vif(reg3)acc_cred_z idm_z samp
1.023751 1.037512 1.053649 plot_model(reg3, type="diag")[[1]]
[[2]]`geom_smooth()` using formula = 'y ~ x'
[[3]]
[[4]]`geom_smooth()` using formula = 'y ~ x'
plot(reg3, 5)
summary(reg3)
Call:
lm(formula = mis_cred_z ~ acc_cred_z + idm_z + samp, data = df2)
Residuals:
Min 1Q Median 3Q Max
-2.41975 -0.51261 0.00346 0.63414 1.76366
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.16315 0.07287 2.239 0.026458 *
acc_cred_z 0.50734 0.05997 8.460 1.15e-14 ***
idm_z 0.10466 0.06063 1.726 0.086124 .
sampSONA -0.51820 0.13215 -3.921 0.000127 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7906 on 171 degrees of freedom
(9 observations deleted due to missingness)
Multiple R-squared: 0.3857, Adjusted R-squared: 0.3749
F-statistic: 35.78 on 3 and 171 DF, p-value: < 2.2e-16reg3 <- lm(data = df2, mis_dang_z ~ acc_dang_z + idm_z + samp)
car::vif(reg3)acc_dang_z idm_z samp
1.109669 1.036897 1.146297 plot_model(reg3, type="diag")[[1]]
[[2]]`geom_smooth()` using formula = 'y ~ x'
[[3]]
[[4]]`geom_smooth()` using formula = 'y ~ x'
plot(reg3, 5)
summary(reg3)
Call:
lm(formula = mis_dang_z ~ acc_dang_z + idm_z + samp, data = df2)
Residuals:
Min 1Q Median 3Q Max
-2.56046 -0.55236 0.04106 0.62231 1.94738
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.09610 0.08672 1.108 0.2695
acc_dang_z 0.42090 0.07460 5.642 7.7e-08 ***
idm_z 0.14944 0.07137 2.094 0.0379 *
sampSONA -0.24012 0.16541 -1.452 0.1486
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8955 on 156 degrees of freedom
(24 observations deleted due to missingness)
Multiple R-squared: 0.246, Adjusted R-squared: 0.2315
F-statistic: 16.97 on 3 and 156 DF, p-value: 1.376e-09H4 Experiences with discrimination and rejection sensitivity predict misinformation acceptance but not accurate information acceptance, mediated by inequality-driven mistrust
Overall, model is consistent with the hypothesized pathway but does not confirm it.
model <- '
# Direct effects on mis_z (c-prime paths)
mis_z ~ c1*discrim_z + c2*reject_z + samp
# Direct effects on acc_z (c-prime paths)
acc_z ~ c3*discrim_z + c4*reject_z + samp
# Effects of predictors on mediator (a paths)
idm_z ~ a1*discrim_z + a2*reject_z + samp
# Effect of mediator on outcomes (b paths)
mis_z ~ b1*idm_z
acc_z ~ b2*idm_z
# Residual covariance between outcomes
mis_z ~~ acc_z
# Contrast: does IDM differentially predict mis_z vs acc_z?
b_diff := b1 - b2
# Indirect effects on mis_z
indirect_discrim_mis := a1*b1
indirect_reject_mis := a2*b1
# Indirect effects on acc_z
indirect_discrim_acc := a1*b2
indirect_reject_acc := a2*b2
# Total effects on mis_z
total_discrim_mis := c1 + (a1*b1)
total_reject_mis := c2 + (a2*b1)
# Total effects on acc_z
total_discrim_acc := c3 + (a1*b2)
total_reject_acc := c4 + (a2*b2)
'
# Fit the model
# fit <- sem(model, data = df2, se = "bootstrap", bootstrap = 5000)
# saveRDS(fit, file = "mediation_fit10.rds")
fit6 <- readRDS("mediation_fit10.rds")
# Summarize results
summary(fit6, fit.measures = TRUE, ci = TRUE)lavaan 0.6.17 ended normally after 12 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 15
Number of observations 184
Model Test User Model:
Test statistic 0.000
Degrees of freedom 0
Model Test Baseline Model:
Test statistic 230.726
Degrees of freedom 12
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 1.000
Tucker-Lewis Index (TLI) 1.000
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -666.387
Loglikelihood unrestricted model (H1) NA
Akaike (AIC) 1362.774
Bayesian (BIC) 1410.998
Sample-size adjusted Bayesian (SABIC) 1363.489
Root Mean Square Error of Approximation:
RMSEA 0.000
90 Percent confidence interval - lower 0.000
90 Percent confidence interval - upper 0.000
P-value H_0: RMSEA <= 0.050 NA
P-value H_0: RMSEA >= 0.080 NA
Standardized Root Mean Square Residual:
SRMR 0.000
Parameter Estimates:
Standard errors Bootstrap
Number of requested bootstrap draws 5000
Number of successful bootstrap draws 5000
Regressions:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
mis_z ~
discrim_z (c1) 0.070 0.085 0.816 0.414 -0.099 0.241
reject_z (c2) -0.085 0.093 -0.914 0.361 -0.259 0.108
samp -0.683 0.157 -4.341 0.000 -0.985 -0.375
acc_z ~
discrim_z (c3) -0.053 0.089 -0.592 0.554 -0.226 0.120
reject_z (c4) -0.052 0.098 -0.528 0.598 -0.241 0.147
samp -0.686 0.159 -4.317 0.000 -1.000 -0.379
idm_z ~
discrim_z (a1) 0.164 0.072 2.282 0.022 0.018 0.306
reject_z (a2) 0.498 0.067 7.422 0.000 0.364 0.628
samp -0.371 0.131 -2.842 0.004 -0.627 -0.107
mis_z ~
idm_z (b1) 0.160 0.092 1.739 0.082 -0.022 0.340
acc_z ~
idm_z (b2) 0.098 0.092 1.075 0.283 -0.081 0.273
Covariances:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
.mis_z ~~
.acc_z 0.526 0.068 7.785 0.000 0.387 0.645
Variances:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
.mis_z 0.861 0.077 11.227 0.000 0.690 0.996
.acc_z 0.872 0.092 9.432 0.000 0.672 1.036
.idm_z 0.593 0.065 9.165 0.000 0.462 0.714
Defined Parameters:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
b_diff 0.061 0.097 0.634 0.526 -0.126 0.254
indrct_dscrm_m 0.026 0.020 1.288 0.198 -0.005 0.074
indrct_rjct_ms 0.079 0.047 1.705 0.088 -0.011 0.170
indrct_dscrm_c 0.016 0.017 0.952 0.341 -0.016 0.053
indrct_rjct_cc 0.049 0.046 1.061 0.289 -0.041 0.138
total_dscrm_ms 0.096 0.085 1.126 0.260 -0.071 0.267
total_rejct_ms -0.005 0.086 -0.062 0.950 -0.171 0.169
total_dscrm_cc -0.036 0.088 -0.414 0.679 -0.211 0.134
total_rejct_cc -0.003 0.095 -0.031 0.975 -0.186 0.184# lay3 <- get_layout(
# "discrim_z", "", "samp",
# "", "idm_z", "mis_z",
# "reject_z", "", "",
# rows = 3
# )
#
# graph_sem(fit6, layout = lay3)
#
# lay3 <- get_layout(
# "discrim_z", "", "samp",
# "", "idm_z", "acc_z",
# "reject_z", "", "",
# rows = 3
# )
#
# graph_sem(fit6, layout = lay3)
lay3 <- get_layout(
"discrim_z", "", "mis_z", "",
"", "idm_z", "", "samp",
"reject_z", "", "acc_z", "",
rows = 3
)
# Prepare the graph object using your existing fit and layout
graph_data <- prepare_graph(fit6, layout = lay3)
# Filter out edges with p >= .10
edges(graph_data) <- edges(graph_data)[edges(graph_data)$pval < .10, ]
edges(graph_data)$label[edges(graph_data)$from == "idm_z" & edges(graph_data)$to == "mis_z"] <- "0.16†"
edges(graph_data) <- edges(graph_data)[edges(graph_data)$curvature == 0 | is.na(edges(graph_data)$curvature) | edges(graph_data)$from != "mis_z" | edges(graph_data)$to != "acc_z", ]
# Inspect the current node labels
nodes(graph_data) name shape label x y node_xmin node_xmax node_ymin node_ymax show
1 acc_z rect acc_z 6 2 5.4 6.6 1.6 2.4 TRUE
2 discrim_z rect discrim_z 2 6 1.4 2.6 5.6 6.4 TRUE
3 idm_z rect idm_z 4 4 3.4 4.6 3.6 4.4 TRUE
4 mis_z rect mis_z 6 6 5.4 6.6 5.6 6.4 TRUE
5 reject_z rect reject_z 2 2 1.4 2.6 1.6 2.4 TRUE
6 samp rect samp 8 4 7.4 8.6 3.6 4.4 TRUE# Update labels
nodes(graph_data)$label[nodes(graph_data)$name == "discrim_z"] <- "Discrimination"
nodes(graph_data)$label[nodes(graph_data)$name == "reject_z"] <- "Rejection\nSensitivity"
nodes(graph_data)$label[nodes(graph_data)$name == "idm_z"] <- "Inequality-Driven\nMistrust"
nodes(graph_data)$label[nodes(graph_data)$name == "mis_z"] <- "Misinfo.\nAcceptance"
nodes(graph_data)$label[nodes(graph_data)$name == "acc_z"] <- "Accurate Info.\nAcceptance"
nodes(graph_data)$label[nodes(graph_data)$name == "samp"] <- "Student Status"
# Plot the filtered graph
plot(graph_data)Some edges involve nodes with argument 'show = FALSE'. These were dropped.Warning in (function (mapping = NULL, data = NULL, stat = "identity", position
= "identity", : Ignoring unknown parameters: `linewidth`Warning in (function (mapping = NULL, data = NULL, stat = "identity", position
= "identity", : Ignoring unknown parameters: `linewidth`Warning in (function (mapping = NULL, data = NULL, stat = "identity", position
= "identity", : Ignoring unknown parameters: `linewidth`