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(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)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.17 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)
)
describe(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 describe(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 EFA <- factanal(d, factors = 2, rotation = "promax")
print(EFA, digits=3, cutoff=.4, sort=F)
Call:
factanal(x = d, factors = 2, rotation = "promax")
Uniquenesses:
Hiring Discrimination Seen as Untrustworthy Seen as Dangerous
0.538 0.246 0.322
Devalued Looked Down On Seen as Less Intelligent
0.145 0.085 0.355
Loadings:
Factor1 Factor2
Hiring Discrimination 0.541
Seen as Untrustworthy 0.860
Seen as Dangerous 0.859
Devalued 0.864
Looked Down On 1.014
Seen as Less Intelligent 0.661
Factor1 Factor2
SS loadings 2.505 1.547
Proportion Var 0.418 0.258
Cumulative Var 0.418 0.675
Factor Correlations:
Factor1 Factor2
Factor1 1.000 -0.782
Factor2 -0.782 1.000
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 18.69 on 4 degrees of freedom.
The p-value is 0.000906 d <- subset(d, select = -c(`Seen as Dangerous`))
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 Devalued
0.531 0.518 0.129
Looked Down On Seen as Less Intelligent
0.119 0.360
Loadings:
Factor1
Hiring Discrimination 0.685
Seen as Untrustworthy 0.694
Devalued 0.933
Looked Down On 0.939
Seen as Less Intelligent 0.800
Factor1
SS loadings 3.343
Proportion Var 0.669
Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 23.77 on 5 degrees of freedom.
The p-value is 0.000241 describe(subset(df1, select=c(grep("Q9_",colnames(df1))))) vars n mean sd median trimmed mad min max range skew kurtosis se
Q9_1 1 420 2.19 1.11 2 2.06 1.48 1 5 4 0.83 -0.08 0.05
Q9_2 2 420 2.11 1.09 2 1.97 1.48 1 5 4 0.91 0.08 0.05
Q9_3 3 419 2.89 1.34 3 2.87 1.48 1 5 4 0.16 -1.24 0.07
Q9_4 4 419 1.70 1.03 1 1.48 0.00 1 5 4 1.64 2.11 0.05
Q9_5 5 420 2.11 1.18 2 1.95 1.48 1 5 4 0.95 -0.05 0.06
Q9_6 6 420 2.22 1.31 2 2.04 1.48 1 5 4 0.84 -0.50 0.06
Q9_7 7 420 2.90 1.31 3 2.88 1.48 1 5 4 0.25 -1.21 0.06
Q9_8 8 420 1.94 1.15 2 1.74 1.48 1 5 4 1.20 0.50 0.06
Q9_9 9 419 2.38 1.31 2 2.23 1.48 1 5 4 0.63 -0.73 0.06
Q9_10 10 419 2.42 1.30 2 2.28 1.48 1 5 4 0.59 -0.79 0.06
Q9_11 11 420 1.87 1.13 1 1.65 0.00 1 5 4 1.32 0.88 0.06
Q9_12 12 420 2.10 1.17 2 1.93 1.48 1 5 4 0.97 0.01 0.06
Q9_13 13 420 2.49 1.28 2 2.37 1.48 1 5 4 0.50 -0.87 0.06
Q9_14 14 420 2.55 1.36 2 2.43 1.48 1 5 4 0.50 -1.01 0.07
Q9_15 15 420 2.77 1.34 2 2.71 1.48 1 5 4 0.31 -1.14 0.07
Q9_16 16 420 2.80 1.37 3 2.74 1.48 1 5 4 0.24 -1.23 0.07d <- na.omit(subset(df1, select=c(grep("Q9_", colnames(df1)))))
labels <- c(
Q9_1 = "Clarity: Difficulty Making Sense of Feelings",
Q9_2 = "Clarity: Confused About Feelings",
Q9_3 = "Goals: Difficulty Working",
Q9_4 = "Impulse: Out of Control",
Q9_5 = "Strategies: Will Remain Upset",
Q9_6 = "Strategies: Will Feel Depressed",
Q9_7 = "Goals: Difficulty Focusing",
Q9_8 = "Impulse: Feels Out of Control",
Q9_9 = "Nonacceptance: Ashamed of Feelings",
Q9_10 = "Nonacceptance: Feels Weak",
Q9_11 = "Impulse: Difficulty Controlling Behavior",
Q9_12 = "Strategies: Nothing Will Help",
Q9_13 = "Nonacceptance: Irritated at Self",
Q9_14 = "Nonacceptance: Feels Bad About Self",
Q9_15 = "Goals: Can't Think of Anything Else",
Q9_16 = "Strategies: Emotions Overwhelming"
)
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 item 16, was crossloading on Goals subscale (should be Strategies).
d <- na.omit(subset(df1, select=c(grep("Q9_", colnames(df1)))))
# specify model: latent variables =~ observed indicators
model <- '
clarity =~ Q9_1 + Q9_2
goals =~ Q9_3 + Q9_7 + Q9_15
impulse =~ Q9_4 + Q9_8 + Q9_11
strategies =~ Q9_5 + Q9_6 + Q9_12
nonacceptance =~ Q9_9 + Q9_10 + Q9_13 + Q9_14
'
# 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.17 ended normally after 35 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 40
Number of observations 416
Model Test User Model:
Test statistic 273.199
Degrees of freedom 80
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 5704.653
Degrees of freedom 105
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.965
Tucker-Lewis Index (TLI) 0.955
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -7389.664
Loglikelihood unrestricted model (H1) -7253.064
Akaike (AIC) 14859.327
Bayesian (BIC) 15020.555
Sample-size adjusted Bayesian (SABIC) 14893.624
Root Mean Square Error of Approximation:
RMSEA 0.076
90 Percent confidence interval - lower 0.066
90 Percent confidence interval - upper 0.086
P-value H_0: RMSEA <= 0.050 0.000
P-value H_0: RMSEA >= 0.080 0.273
Standardized Root Mean Square Residual:
SRMR 0.032
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
clarity =~
Q9_1 0.991 0.045 22.237 0.000 0.991 0.897
Q9_2 1.024 0.043 23.991 0.000 1.024 0.943
goals =~
Q9_3 1.184 0.052 22.756 0.000 1.184 0.888
Q9_7 1.186 0.050 23.573 0.000 1.186 0.907
Q9_15 1.169 0.053 22.209 0.000 1.169 0.875
impulse =~
Q9_4 0.918 0.040 22.777 0.000 0.918 0.893
Q9_8 1.032 0.045 22.730 0.000 1.032 0.892
Q9_11 0.951 0.046 20.622 0.000 0.951 0.839
strategies =~
Q9_5 0.993 0.048 20.609 0.000 0.993 0.838
Q9_6 1.122 0.053 21.312 0.000 1.122 0.856
Q9_12 0.990 0.048 20.819 0.000 0.990 0.844
nonacceptance =~
Q9_9 1.130 0.052 21.774 0.000 1.130 0.864
Q9_10 1.107 0.052 21.421 0.000 1.107 0.855
Q9_13 1.096 0.051 21.506 0.000 1.096 0.857
Q9_14 1.207 0.053 22.736 0.000 1.207 0.887
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
clarity ~~
goals 0.616 0.035 17.657 0.000 0.616 0.616
impulse 0.590 0.037 16.096 0.000 0.590 0.590
strategies 0.700 0.031 22.719 0.000 0.700 0.700
nonacceptance 0.667 0.032 21.062 0.000 0.667 0.667
goals ~~
impulse 0.716 0.029 24.950 0.000 0.716 0.716
strategies 0.845 0.020 41.859 0.000 0.845 0.845
nonacceptance 0.721 0.028 25.877 0.000 0.721 0.721
impulse ~~
strategies 0.744 0.028 26.555 0.000 0.744 0.744
nonacceptance 0.640 0.034 19.059 0.000 0.640 0.640
strategies ~~
nonacceptance 0.800 0.023 34.134 0.000 0.800 0.800
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.Q9_1 0.239 0.034 6.976 0.000 0.239 0.196
.Q9_2 0.131 0.033 3.941 0.000 0.131 0.111
.Q9_3 0.377 0.037 10.317 0.000 0.377 0.212
.Q9_7 0.304 0.033 9.322 0.000 0.304 0.178
.Q9_15 0.420 0.039 10.852 0.000 0.420 0.235
.Q9_4 0.213 0.024 9.072 0.000 0.213 0.202
.Q9_8 0.273 0.030 9.135 0.000 0.273 0.204
.Q9_11 0.381 0.034 11.292 0.000 0.381 0.297
.Q9_5 0.417 0.037 11.296 0.000 0.417 0.297
.Q9_6 0.457 0.042 10.755 0.000 0.457 0.266
.Q9_12 0.396 0.036 11.146 0.000 0.396 0.288
.Q9_9 0.435 0.039 11.272 0.000 0.435 0.254
.Q9_10 0.451 0.039 11.513 0.000 0.451 0.269
.Q9_13 0.434 0.038 11.457 0.000 0.434 0.265
.Q9_14 0.396 0.038 10.462 0.000 0.396 0.214
clarity 1.000 1.000 1.000
goals 1.000 1.000 1.000
impulse 1.000 1.000 1.000
strategies 1.000 1.000 1.000
nonacceptance 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)
)
describe(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.692
Factor2 0.683 1.000 0.721
Factor3 0.692 0.721 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.000 -0.67 0.679
Factor2 -0.670 1.00 -0.700
Factor3 0.679 -0.70 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 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(m_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
Q40_4 1 87 2.79 1.04 3 2.86 1.48 1 4 3 -0.39 -1.04 0.11
Q42_4 2 99 2.49 0.97 3 2.49 1.48 1 4 3 -0.05 -1.01 0.10
Q44_4 3 94 2.69 0.89 3 2.74 1.48 1 4 3 -0.35 -0.61 0.09
Q46_4 4 90 2.46 0.93 3 2.44 1.48 1 4 3 -0.08 -0.90 0.10
Q50_4 5 99 2.91 0.89 3 3.00 0.00 1 4 3 -0.67 -0.18 0.09
Q54_4 6 96 2.58 0.99 3 2.60 1.48 1 4 3 -0.20 -1.02 0.10
Q58_4 7 96 2.76 0.88 3 2.82 1.48 1 4 3 -0.35 -0.57 0.09
Q62_4 8 88 2.11 0.99 2 2.03 1.48 1 4 3 0.34 -1.06 0.11
Q66_4 9 96 1.92 1.03 2 1.78 1.48 1 4 3 0.73 -0.77 0.11
Q68_4 10 93 2.16 1.08 2 2.08 1.48 1 4 3 0.41 -1.16 0.11
Q74_4 11 98 1.65 0.87 1 1.54 0.00 1 4 3 1.00 -0.25 0.09
Q78_4 12 92 2.20 0.95 2 2.15 1.48 1 4 3 0.14 -1.12 0.10
Q82_4 13 89 2.26 0.90 2 2.21 1.48 1 4 3 0.22 -0.77 0.10
Q86_4 14 92 2.33 0.94 2 2.28 1.48 1 4 3 0.03 -1.00 0.10
Q90_4 15 98 1.84 0.97 2 1.70 1.48 1 4 3 0.86 -0.41 0.10labels <- 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")
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()`).
Very imperfect due to structure of the data.
d2 <- d
set.seed(42) # for reproducibility
# mi <- mifa(
# data = d2,
# m = 30,
# ci = "boot",
# n_pc = 1:10,
# print = F
# )
# save(mi, file = "mi.RData")
load("mi.RData") # restores 'mi' object back into environment
# Use pooled covariance matrix from mifa, converted to correlation matrix
cov_pooled <- mi$cov_combined
cor_pooled <- cov2cor(cov_pooled)
ev <- eigen(cor_pooled) # get eigenvalues
ap <- parallel(subject = nrow(d2), var = ncol(cor_pooled), 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)
plotnScree(nS) # shows us the scree plot, look for the elbows
# EFA using pooled covariance matrix
EFA <- fa(
r = cov_pooled,
nfactors = 3, # update based on scree
n.obs = nrow(d2),
fm = "ml",
rotate = "promax"
)Loading required namespace: GPArotationprint(EFA$loadings, sort = TRUE, cutoff = .4)
Loadings:
ML1 ML2 ML3
Trauma2 0.601
Psychiatric Symptoms 1.012 -0.618
Introversion 0.694
Friendship 0.588
Narcissism 0.609
Trauma1 0.552
Toxic People 0.834
Manipulativeness 0.503
Love 0.528
Anxiety & Emotion Regulation 0.526
Depression 0.488
Intelligence 0.477
Attraction 0.474
Childhood 0.487
Family 0.425
ML1 ML2 ML3
SS loadings 3.004 2.097 1.658
Proportion Var 0.200 0.140 0.111
Cumulative Var 0.200 0.340 0.451# EFA using pooled covariance matrix
EFA <- fa(
r = cov_pooled,
nfactors = 2, # update based on scree
n.obs = nrow(d2),
fm = "ml",
rotate = "promax"
)
print(EFA$loadings, sort = TRUE, cutoff = .4)
Loadings:
ML1 ML2
Narcissism 0.603
Depression 0.556
Intelligence 0.648
Trauma1 0.851
Toxic People 0.738
Trauma2 0.539
Psychiatric Symptoms 0.812
Introversion 0.805
Attraction 0.568
Childhood 0.570
Friendship 0.580
Family 0.514
Love 0.486
Manipulativeness 0.468
Anxiety & Emotion Regulation 0.441
ML1 ML2
SS loadings 3.089 3.043
Proportion Var 0.206 0.203
Cumulative Var 0.206 0.409a_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")
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()`).
Very imperfect due to structure of the data.
d2 <- d
set.seed(42) # for reproducibility
# mi2 <- mifa(
# data = d2,
# m = 30,
# ci = "boot",
# n_pc = 1:10,
# print = F
# )
# save(mi2, file = "mi2.RData")
load("mi2.RData") # restores 'mi2' object back into environment
# Use pooled covariance matrix from mifa, converted to correlation matrix
cov_pooled <- mi2$cov_combined
cor_pooled <- cov2cor(cov_pooled)
ev <- eigen(cor_pooled) # get eigenvalues
ap <- parallel(subject = nrow(d2), var = ncol(cor_pooled), 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)
plotnScree(nS) # shows us the scree plot, look for the elbows
# EFA using pooled covariance matrix
EFA <- fa(
r = cov_pooled,
nfactors = 3, # update based on scree
n.obs = nrow(d2),
fm = "ml",
rotate = "promax"
)
print(EFA$loadings, sort = TRUE, cutoff = .4)
Loadings:
ML1 ML2 ML3
Trauma2 0.692
Introversion 0.681
Manipulativeness 0.734
Attraction 0.572
Childhood 0.788
Love 0.539
Toxic People 0.854
Psychiatric Symptoms 0.735
Narcissism 0.534
Depression 0.772
Intelligence 0.419 0.473
Trauma1 0.484
Friendship 0.463
Anxiety & Emotion Regulation 0.430
Family 0.493
ML1 ML2 ML3
SS loadings 3.294 1.905 1.509
Proportion Var 0.220 0.127 0.101
Cumulative Var 0.220 0.347 0.447# Identify target columns
target_cols <- grep(
paste(c(m_list, 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
Q40_4 6 87 2.79 1.04 3 2.86 1.48 1 4 3 -0.39 -1.04 0.11
Q42_4 7 99 2.49 0.97 3 2.49 1.48 1 4 3 -0.05 -1.01 0.10
Q44_4 8 94 2.69 0.89 3 2.74 1.48 1 4 3 -0.35 -0.61 0.09
Q46_4 9 90 2.46 0.93 3 2.44 1.48 1 4 3 -0.08 -0.90 0.10
Q50_4 10 99 2.91 0.89 3 3.00 0.00 1 4 3 -0.67 -0.18 0.09
Q52_4 11 88 2.40 0.99 2 2.38 1.48 1 4 3 0.07 -1.07 0.11
Q56_4 12 88 2.95 0.82 3 3.00 1.48 1 4 3 -0.42 -0.39 0.09
Q60_4 13 96 3.14 0.79 3 3.19 1.48 1 4 3 -0.49 -0.56 0.08
Q64_4 14 88 2.36 1.00 2 2.33 1.48 1 4 3 0.13 -1.07 0.11
Q70_4 15 91 2.23 0.92 2 2.18 1.48 1 4 3 0.21 -0.88 0.10
Q54_4 16 96 2.58 0.99 3 2.60 1.48 1 4 3 -0.20 -1.02 0.10
Q58_4 17 96 2.76 0.88 3 2.82 1.48 1 4 3 -0.35 -0.57 0.09
Q62_4 18 88 2.11 0.99 2 2.03 1.48 1 4 3 0.34 -1.06 0.11
Q66_4 19 96 1.92 1.03 2 1.78 1.48 1 4 3 0.73 -0.77 0.11
Q68_4 20 93 2.16 1.08 2 2.08 1.48 1 4 3 0.41 -1.16 0.11
Q72_4 21 86 2.31 1.05 2 2.27 1.48 1 4 3 0.07 -1.30 0.11
Q76_4 22 92 2.62 0.99 3 2.65 1.48 1 4 3 -0.20 -1.03 0.10
Q80_4 23 95 2.64 0.98 3 2.68 1.48 1 4 3 -0.19 -0.99 0.10
Q84_4 24 92 2.48 0.88 3 2.47 1.48 1 4 3 -0.12 -0.77 0.09
Q88_4 25 86 2.45 0.97 3 2.44 1.48 1 4 3 -0.06 -1.02 0.10
Q74_4 26 98 1.65 0.87 1 1.54 0.00 1 4 3 1.00 -0.25 0.09
Q78_4 27 92 2.20 0.95 2 2.15 1.48 1 4 3 0.14 -1.12 0.10
Q82_4 28 89 2.26 0.90 2 2.21 1.48 1 4 3 0.22 -0.77 0.10
Q86_4 29 92 2.33 0.94 2 2.28 1.48 1 4 3 0.03 -1.00 0.10
Q90_4 30 98 1.84 0.97 2 1.70 1.48 1 4 3 0.86 -0.41 0.10labels <- c(
Q40_4 = "m_Love",
Q42_4 = "m_Narcissism",
Q44_4 = "m_Depression",
Q46_4 = "m_Intelligence",
Q50_4 = "m_Trauma1",
Q54_4 = "m_Trauma2",
Q58_4 = "m_Toxic People",
Q62_4 = "m_Psychiatric Symptoms",
Q66_4 = "m_Introversion",
Q68_4 = "m_Manipulativeness",
Q74_4 = "m_Attraction",
Q78_4 = "m_Childhood",
Q82_4 = "m_Friendship",
Q86_4 = "m_Anxiety & Emotion Regulation",
Q90_4 = "m_Family",
Q23_4 = "a_Love",
Q34_4 = "a_Narcissism",
Q36_4 = "a_Depression",
Q38_4 = "a_Intelligence",
Q48_4 = "a_Trauma1",
Q52_4 = "a_Trauma2",
Q56_4 = "a_Toxic People",
Q60_4 = "a_Psychiatric Symptoms",
Q64_4 = "a_Introversion",
Q70_4 = "a_Manipulativeness",
Q72_4 = "a_Attraction",
Q76_4 = "a_Childhood",
Q80_4 = "a_Friendship",
Q84_4 = "a_Anxiety & Emotion Regulation",
Q88_4 = "a_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) <- ifelse(names(d) %in% names(labels),
labels[names(d)],
names(d))
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()Warning: Removed 2760 rows containing non-finite outside the scale range
(`stat_bin()`).
Very imperfect due to structure of the data.
d2 <- d
set.seed(42) # for reproducibility
# mi3 <- mifa(
# data = d2,
# m = 30,
# ci = "boot",
# n_pc = 1:10,
# print = F
# )
# save(mi3, file = "mi3.RData")
load("mi3.RData") # restores 'mi' object back into environment
# Use pooled covariance matrix from mifa, converted to correlation matrix
cov_pooled <- mi3$cov_combined
cor_pooled <- cov2cor(cov_pooled)
ev <- eigen(cor_pooled) # get eigenvalues
ap <- parallel(subject = nrow(d2), var = ncol(cor_pooled), 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)
plotnScree(nS) # shows us the scree plot, look for the elbows
# EFA using pooled covariance matrix
EFA <- fa(
r = cov_pooled,
nfactors = 3, # update based on scree
n.obs = nrow(d2),
fm = "ml",
rotate = "promax"
)
print(EFA$loadings, sort = TRUE, cutoff = .4)
Loadings:
ML1 ML2 ML3
a_Trauma2 0.558
m_Trauma2 0.520
m_Psychiatric Symptoms 0.727 -0.452
m_Introversion 0.700
m_Attraction 0.594
m_Friendship 0.530
m_Family 0.604
a_Love 0.574
a_Toxic People 0.600
a_Friendship 0.558
m_Narcissism 0.577
m_Intelligence 0.606
m_Toxic People 0.718
a_Narcissism
a_Depression 0.410
a_Intelligence 0.494
a_Trauma1 0.449
m_Love 0.407
m_Depression
m_Trauma1 0.479
a_Psychiatric Symptoms 0.494
a_Introversion 0.494
a_Manipulativeness 0.403
m_Manipulativeness 0.401 0.465
a_Attraction
a_Childhood 0.426
a_Anxiety & Emotion Regulation
a_Family
m_Childhood 0.490
m_Anxiety & Emotion Regulation 0.483
ML1 ML2 ML3
SS loadings 4.194 3.082 2.616
Proportion Var 0.140 0.103 0.087
Cumulative Var 0.140 0.243 0.330describe(subset(df2, select=c(grep("Q29_",colnames(df2))))) vars n mean sd median trimmed mad min max range skew kurtosis se
Q29_1 1 184 2.86 0.66 3 2.86 0.00 1 4 3 -0.51 0.69 0.05
Q29_2 2 184 2.95 0.67 3 2.96 0.00 1 4 3 -0.27 0.07 0.05
Q29_3 3 184 1.77 0.80 2 1.68 1.48 1 4 3 0.76 -0.14 0.06
Q29_4 4 184 2.02 0.95 2 1.93 1.48 1 4 3 0.53 -0.72 0.07
Q29_5 5 184 2.27 0.95 2 2.21 1.48 1 4 3 0.14 -0.98 0.07
Q29_6 6 184 2.28 0.96 2 2.23 1.48 1 4 3 0.15 -1.00 0.07
Q29_7 7 184 2.09 0.91 2 2.03 1.48 1 4 3 0.34 -0.86 0.07
Q29_8 8 184 1.87 0.88 2 1.78 1.48 1 4 3 0.63 -0.60 0.07d <- subset(df2, select=c(grep("Q29_", colnames(df2))))
labels <- c(
Q29_1 = "Ordinary People Unites",
Q29_2 = "Ordinary People Good",
Q29_3 = "Scientists Advantage",
Q29_4 = "Scientists Cahoots",
Q29_5 = "People Influence",
Q29_6 = "People Decisions",
Q29_7 = "Rely Experience",
Q29_8 = "Rely Common Sense"
)
# 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()
d <- subset(df2, select=c(grep("Q29_", colnames(df2))))
# Specify CFA model
model <- '
ppl =~ Q29_1 + Q29_2
eli =~ Q29_3 + Q29_4
dec =~ Q29_5 + Q29_6
tru =~ Q29_7 + Q29_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.17 ended normally after 29 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 22
Number of observations 184
Model Test User Model:
Test statistic 24.826
Degrees of freedom 14
P-value (Chi-square) 0.036
Model Test Baseline Model:
Test statistic 490.930
Degrees of freedom 28
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.977
Tucker-Lewis Index (TLI) 0.953
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -1594.557
Loglikelihood unrestricted model (H1) -1582.144
Akaike (AIC) 3233.114
Bayesian (BIC) 3303.843
Sample-size adjusted Bayesian (SABIC) 3234.163
Root Mean Square Error of Approximation:
RMSEA 0.065
90 Percent confidence interval - lower 0.016
90 Percent confidence interval - upper 0.106
P-value H_0: RMSEA <= 0.050 0.251
P-value H_0: RMSEA >= 0.080 0.301
Standardized Root Mean Square Residual:
SRMR 0.036
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
ppl =~
Q29_1 0.511 0.205 2.497 0.013 0.511 0.774
Q29_2 0.252 0.110 2.296 0.022 0.252 0.376
eli =~
Q29_3 0.640 0.059 10.910 0.000 0.640 0.803
Q29_4 0.704 0.070 10.125 0.000 0.704 0.746
dec =~
Q29_5 0.808 0.067 12.023 0.000 0.808 0.856
Q29_6 0.801 0.068 11.703 0.000 0.801 0.835
tru =~
Q29_7 0.747 0.063 11.824 0.000 0.747 0.820
Q29_8 0.692 0.061 11.268 0.000 0.692 0.785
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
ppl ~~
eli -0.067 0.109 -0.614 0.539 -0.067 -0.067
dec 0.152 0.111 1.363 0.173 0.152 0.152
tru 0.154 0.115 1.346 0.178 0.154 0.154
eli ~~
dec 0.532 0.074 7.213 0.000 0.532 0.532
tru 0.726 0.062 11.735 0.000 0.726 0.726
dec ~~
tru 0.628 0.064 9.811 0.000 0.628 0.628
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.Q29_1 0.175 0.206 0.848 0.396 0.175 0.401
.Q29_2 0.385 0.064 6.013 0.000 0.385 0.858
.Q29_3 0.225 0.049 4.633 0.000 0.225 0.355
.Q29_4 0.396 0.066 6.000 0.000 0.396 0.444
.Q29_5 0.237 0.066 3.574 0.000 0.237 0.267
.Q29_6 0.279 0.067 4.148 0.000 0.279 0.303
.Q29_7 0.272 0.055 4.946 0.000 0.272 0.328
.Q29_8 0.297 0.051 5.834 0.000 0.297 0.383
ppl 1.000 1.000 1.000
eli 1.000 1.000 1.000
dec 1.000 1.000 1.000
tru 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)
)
describe(subset(df2, select=c(grep("Q30_",colnames(df2))))) vars n mean sd median trimmed mad min max range skew kurtosis se
Q30_1 1 184 2.92 0.91 3 3.01 1.48 1 4 3 -0.53 -0.53 0.07
Q30_2 2 184 2.85 0.94 3 2.94 1.48 1 4 3 -0.60 -0.48 0.07
Q30_3 3 184 2.45 1.02 3 2.44 1.48 1 4 3 -0.02 -1.14 0.08
Q30_4 4 184 3.02 0.99 3 3.15 1.48 1 4 3 -0.72 -0.54 0.07
Q30_5 5 184 2.77 0.93 3 2.83 1.48 1 4 3 -0.33 -0.76 0.07
Q30_6 6 184 2.49 0.92 2 2.49 1.48 1 4 3 0.03 -0.84 0.07
Q30_7 7 184 2.65 0.90 3 2.68 1.48 1 4 3 -0.10 -0.79 0.07
Q30_8 8 184 2.77 0.80 3 2.79 0.00 1 4 3 -0.33 -0.29 0.06
Q30_9 9 184 2.82 0.82 3 2.84 1.48 1 4 3 -0.30 -0.45 0.06
Q30_10 10 184 2.55 0.95 3 2.56 1.48 1 4 3 0.01 -0.95 0.07
Q30_11 11 184 2.57 0.87 3 2.59 1.48 1 4 3 -0.02 -0.70 0.06
Q30_12 12 184 2.89 0.82 3 2.93 1.48 1 4 3 -0.38 -0.36 0.06
Q30_13 13 184 2.79 0.91 3 2.86 1.48 1 4 3 -0.38 -0.64 0.07
Q30_14 14 184 2.60 0.97 3 2.63 1.48 1 4 3 -0.15 -0.97 0.07
Q30_15 15 184 2.83 0.93 3 2.91 1.48 1 4 3 -0.44 -0.66 0.07d <- subset(df2, select=c(grep("Q30_", colnames(df2))))
labels <- c(
Q30_1 = "Breakdown Recognize",
Q30_2 = "Breakdown Values",
Q30_3 = "Breakdown Change",
Q30_4 = "Breakdown Morals",
Q30_5 = "Breakdown Trust",
Q30_6 = "Disintegration Morals",
Q30_7 = "Disintegration Selfish",
Q30_8 = "Disintegration Trust Rely",
Q30_9 = "Disintegration Trust",
Q30_10 = "Disintegration Care",
Q30_11 = "Disintegration Communal",
Q30_12 = "Disintegration End Justifies",
Q30_13 = "Disintegration Rules Fading",
Q30_14 = "Disintegration Good Evil",
Q30_15 = "Disintegration No Regard"
)
# 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()
Initial model had CFI/TLI in 80s. Dropped lower-loading items until fit indices improved. Dropped 1 and 3 from Breakdown, and 7, 9, 10, 11, and 12 from Disintegration
d <- subset(df2, select=c(grep("Q30_", colnames(df2))))
# Specify CFA model
model <- '
breakdown =~ Q30_2 + Q30_4 + Q30_5
disintegration =~ Q30_6 + Q30_8 + Q30_13 + Q30_14 + Q30_15
'
# 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.17 ended normally after 25 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 17
Number of observations 184
Model Test User Model:
Test statistic 46.082
Degrees of freedom 19
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 854.680
Degrees of freedom 28
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.967
Tucker-Lewis Index (TLI) 0.952
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -1559.581
Loglikelihood unrestricted model (H1) -1536.540
Akaike (AIC) 3153.162
Bayesian (BIC) 3207.816
Sample-size adjusted Bayesian (SABIC) 3153.973
Root Mean Square Error of Approximation:
RMSEA 0.088
90 Percent confidence interval - lower 0.056
90 Percent confidence interval - upper 0.121
P-value H_0: RMSEA <= 0.050 0.028
P-value H_0: RMSEA >= 0.080 0.685
Standardized Root Mean Square Residual:
SRMR 0.037
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
breakdown =~
Q30_2 0.693 0.063 11.064 0.000 0.693 0.737
Q30_4 0.779 0.064 12.243 0.000 0.779 0.793
Q30_5 0.702 0.061 11.460 0.000 0.702 0.756
disintegration =~
Q30_6 0.699 0.059 11.876 0.000 0.699 0.764
Q30_8 0.556 0.053 10.464 0.000 0.556 0.697
Q30_13 0.712 0.057 12.440 0.000 0.712 0.788
Q30_14 0.810 0.059 13.623 0.000 0.810 0.837
Q30_15 0.722 0.059 12.238 0.000 0.722 0.780
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
breakdown ~~
disintegration 0.926 0.028 32.597 0.000 0.926 0.926
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.Q30_2 0.405 0.051 7.963 0.000 0.405 0.457
.Q30_4 0.360 0.050 7.206 0.000 0.360 0.372
.Q30_5 0.371 0.048 7.749 0.000 0.371 0.429
.Q30_6 0.349 0.042 8.245 0.000 0.349 0.417
.Q30_8 0.327 0.038 8.687 0.000 0.327 0.514
.Q30_13 0.309 0.039 8.006 0.000 0.309 0.378
.Q30_14 0.279 0.038 7.316 0.000 0.279 0.299
.Q30_15 0.336 0.042 8.097 0.000 0.336 0.392
breakdown 1.000 1.000 1.000
disintegration 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)
)
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_6,Q39_8)), na.rm = TRUE)) %>%
mutate(reject = rowMeans(across(c(Q40_1,Q40_2,Q40_4,Q40_5,Q40_6)), na.rm = TRUE)) %>%
mutate(ders_clarity = rowMeans(across(c(Q9_1,Q9_2)), na.rm = TRUE)) %>%
mutate(ders_goals = rowMeans(across(c(Q9_3,Q9_7,Q9_15)), na.rm = TRUE)) %>%
mutate(ders_impulse = rowMeans(across(c(Q9_4,Q9_8,Q9_11)), na.rm = TRUE)) %>%
mutate(ders_strat = rowMeans(across(c(Q9_5,Q9_6,Q9_12)), na.rm = TRUE)) %>%
mutate(ders_nonacc = rowMeans(across(c(Q9_9,Q9_10,Q9_13,Q9_14)), 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)), 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))
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_6,Q7_8)), na.rm = TRUE)) %>%
mutate(reject = rowMeans(across(c(Q8_1,Q8_2,Q8_4,Q8_5,Q8_6)), na.rm = TRUE)) %>%
mutate(ders_clarity = rowMeans(across(c(Q9_1,Q9_2)), na.rm = TRUE)) %>%
mutate(ders_goals = rowMeans(across(c(Q9_3,Q9_7,Q9_15)), na.rm = TRUE)) %>%
mutate(ders_impulse = rowMeans(across(c(Q9_4,Q9_8,Q9_11)), na.rm = TRUE)) %>%
mutate(ders_strat = rowMeans(across(c(Q9_5,Q9_6,Q9_12)), na.rm = TRUE)) %>%
mutate(ders_nonacc = rowMeans(across(c(Q9_9,Q9_10,Q9_13,Q9_14)), 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)), 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(socsupp = rowMeans(across(c(Q31_1,Q31_2,Q31_5,Q31_10,Q31_3,Q31_4,Q31_8,Q31_11,Q31_6,Q31_7,Q31_9,Q31_12)), na.rm = TRUE)) %>%
mutate(socsupp_sig = rowMeans(across(c(Q31_1,Q31_2,Q31_5,Q31_10)), na.rm = TRUE)) %>%
mutate(socsupp_fam = rowMeans(across(c(Q31_3,Q31_4,Q31_8,Q31_11)), na.rm = TRUE)) %>%
mutate(socsupp_frd = rowMeans(across(c(Q31_6,Q31_7,Q31_9,Q31_12)), na.rm = TRUE)) %>%
mutate(scipop = rowMeans(across(c(Q29_1,Q29_2,Q29_3,Q29_4,Q29_5,Q29_6,Q29_7,Q29_8)), na.rm = TRUE)) %>%
mutate(scipop_ppl = rowMeans(across(c(Q29_1,Q29_2)), na.rm = TRUE)) %>%
mutate(scipop_eli = rowMeans(across(c(Q29_3,Q29_4)), na.rm = TRUE)) %>%
mutate(scipop_dec = rowMeans(across(c(Q29_5,Q29_6)), na.rm = TRUE)) %>%
mutate(scipop_tru = rowMeans(across(c(Q29_7,Q29_8)), na.rm = TRUE)) %>%
mutate(anomie = rowMeans(across(c(Q30_2,Q30_4,Q30_5,Q30_6,Q30_8,Q30_13,Q30_14,Q30_15)), na.rm = TRUE)) %>%
mutate(anomie_break = rowMeans(across(c(Q30_2,Q30_4,Q30_5)), na.rm = TRUE)) %>%
mutate(anomie_dis = rowMeans(across(c(Q30_6,Q30_8,Q30_13,Q30_14,Q30_15)), 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,
ders_clarity, ders_goals, ders_impulse, ders_strat, ders_nonacc,
idm_dp, idm_da, idm_ib, idm),
~ as.numeric(scale(.)),
.names = "{.col}_z"))
# Standardize composite variables in df2
df2 <- df2 %>%
mutate(across(c(stress_helpless, stress_efficacy, discrim, reject,
ders_clarity, ders_goals, ders_impulse, ders_strat, ders_nonacc,
idm_dp, idm_da, idm_ib, idm,
mis, acc, mis_cred, acc_cred, mis_dang, acc_dang,
socsupp, socsupp_sig, socsupp_fam, socsupp_frd,
scipop, scipop_ppl, scipop_eli, scipop_dec, scipop_tru,
anomie, anomie_break, anomie_dis),
~ 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 Stress, minority stress, and ER predict IDM (sample 1)
Discrimination, rejection sensitivity, and emotion dysregulation (strat) were significant. Students were significantly lower in IDM than non-students.
Experiences with discrimination, rejection sensitivity, and limited access to emotion regulation strategies are the most consistent predictors of inequality-driven mistrust. Although stress (self-efficacy) significantly predicted distrust of authority, it was not significant for the other subscales or for the overall measure, suggesting that IDM is a product of minority stress and not a product of general stress. Most of the emotion dysregulation subscales were not significant, although clarity significantly predicted distrust of authority and there were some borderline significant effects of difficulties engaging in goal-directed behavior on defensive processing. These findings suggest that minority stress contributes to inequality-driven mistrust by depleting access to emotion regulation strategies, undermining participants’ ability to manage the emotional weight of minority stress and to appraise social systems objectively.
reg4 <- lm(data=df1, idm_z ~ stress_helpless_z + stress_efficacy_z + discrim_z + reject_z + ders_clarity_z + ders_goals_z + ders_impulse_z + ders_strat_z + ders_nonacc_z + samp + parent_edu)
car::vif(reg4)stress_helpless_z stress_efficacy_z discrim_z reject_z
2.271912 1.728547 1.657134 1.483820
ders_clarity_z ders_goals_z ders_impulse_z ders_strat_z
2.108225 2.876729 2.123539 3.747133
ders_nonacc_z samp parent_edu
2.517875 1.593332 1.434601 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 + ders_clarity_z + ders_goals_z + ders_impulse_z +
ders_strat_z + ders_nonacc_z + samp + parent_edu, data = df1)
Residuals:
Min 1Q Median 3Q Max
-2.17549 -0.52898 -0.03084 0.42733 2.06355
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.333644 0.110190 3.028 0.00262 **
stress_helpless_z -0.039587 0.053264 -0.743 0.45778
stress_efficacy_z 0.001758 0.046477 0.038 0.96984
discrim_z 0.293925 0.045511 6.458 3.03e-10 ***
reject_z 0.211841 0.042932 4.934 1.18e-06 ***
ders_clarity_z 0.010149 0.051245 0.198 0.84311
ders_goals_z 0.011055 0.059854 0.185 0.85355
ders_impulse_z 0.021852 0.051400 0.425 0.67096
ders_strat_z 0.291718 0.068255 4.274 2.40e-05 ***
ders_nonacc_z 0.031592 0.056002 0.564 0.57299
sampSONA -0.490994 0.089448 -5.489 7.13e-08 ***
parent_edu -0.035158 0.038349 -0.917 0.35980
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7214 on 406 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.4937, Adjusted R-squared: 0.48
F-statistic: 35.99 on 11 and 406 DF, p-value: < 2.2e-16Effect of minority stress on IDM mediated by ER (sample 1)
Overall, the mediation model provided mixed support for the hypothesis. The direct effects are significant for both predictors, with rejection sensitivity showing a particularly strong direct effect, but the indirect effects are generally weak. The mediation is only partial, and the lower bounds for the rejection sensitivity mediation are approaching zero. This suggests that limited access to strategies partially explains the impact of discrimination and rejection sensitivity on inequality-driven mistrust, but that other factors explain the bulk of the effect.
# Define the model
model <- '
# Direct effects on Y (c-prime paths)
idm_z ~ c1*discrim_z + c2*reject_z + samp + parent_edu
# Effects of predictors on mediator (a paths)
ders_strat_z ~ a1*discrim_z + a2*reject_z + samp + parent_edu
# Effect of mediator on Y controlling for predictors (b path)
idm_z ~ b*ders_strat_z
# Covariance between predictors
discrim_z ~~ reject_z
# Indirect effects
indirect_discrim := a1*b
indirect_reject := a2*b
# Test whether indirect effects differ
indirect_diff := indirect_discrim - indirect_reject
# Total effects
total_discrim := c1 + (a1*b)
total_reject := c2 + (a2*b)
'
# Fit the model
# fit <- sem(model, data = df1, se = "bootstrap", bootstrap = 5000)
# saveRDS(fit, file = "mediation_fit1.rds")
fit1 <- readRDS("mediation_fit1.rds")
# Summarize results
summary(fit1, fit.measures = TRUE, ci = TRUE)lavaan 0.6-21 ended normally after 12 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 14
Used Total
Number of observations 418 420
Model Test User Model:
Test statistic 49.310
Degrees of freedom 4
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 504.959
Degrees of freedom 14
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.908
Tucker-Lewis Index (TLI) 0.677
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -2144.267
Loglikelihood unrestricted model (H1) NA
Akaike (AIC) 4316.533
Bayesian (BIC) 4373.030
Sample-size adjusted Bayesian (SABIC) 4328.604
Root Mean Square Error of Approximation:
RMSEA 0.165
90 Percent confidence interval - lower 0.125
90 Percent confidence interval - upper 0.207
P-value H_0: RMSEA <= 0.050 0.000
P-value H_0: RMSEA >= 0.080 1.000
Standardized Root Mean Square Residual:
SRMR 0.076
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
idm_z ~
discrim_z (c1) 0.296 0.045 6.570 0.000 0.210 0.387
reject_z (c2) 0.210 0.046 4.518 0.000 0.119 0.299
samp -0.496 0.085 -5.801 0.000 -0.665 -0.326
parent_ed -0.030 0.037 -0.819 0.413 -0.101 0.043
ders_strat_z ~
discrim_z (a1) 0.316 0.054 5.800 0.000 0.209 0.423
reject_z (a2) 0.135 0.061 2.221 0.026 0.011 0.254
samp -0.037 0.109 -0.334 0.738 -0.249 0.182
parent_ed 0.055 0.047 1.154 0.249 -0.038 0.146
idm_z ~
drs_strt_ (b) 0.317 0.042 7.586 0.000 0.237 0.402
Covariances:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
discrim_z ~~
reject_z 0.462 0.056 8.322 0.000 0.357 0.571
Variances:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
.idm_z 0.507 0.035 14.668 0.000 0.435 0.569
.ders_strat_z 0.843 0.060 13.974 0.000 0.719 0.957
discrim_z 0.996 0.061 16.392 0.000 0.877 1.114
reject_z 1.002 0.066 15.219 0.000 0.874 1.130
Defined Parameters:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
indirect_dscrm 0.100 0.024 4.231 0.000 0.059 0.151
indirect_rejct 0.043 0.019 2.208 0.027 0.004 0.081
indirect_diff 0.057 0.035 1.658 0.097 -0.004 0.131
total_discrim 0.397 0.044 9.090 0.000 0.312 0.483
total_reject 0.253 0.048 5.280 0.000 0.159 0.344lay <- get_layout(
"discrim_z", "", "",
"", "ders_strat_z", "idm_z",
"reject_z", "", "",
rows = 3
)
graph_sem(fit1, layout = lay)
Stress, minority stress, and ER predict IDM (sample 2, replication from sample 1)
High VIF score (> 5) for the ders_strat variable. Removing it from the models did not significantly change findings so it was retained.
Rejection sensitivity and stress (helplessness) were significant. In sample 1, discrimination and emotion dysregulation (strat) were also significant. Students are still lower in IDM than non-students.
Some of the central findings from H1 did not replicate. Discrimination and emotion dysregulation (strat) were no longer significant, while rejection sensitivity remained a strong predictor of IDM. The effect of stress (helplessness) was new. The lack of significant findings for discrimination and emotion dysregulation could be due to the smaller sample size.
reg8 <- lm(data=df2, idm_z ~ stress_helpless_z + stress_efficacy_z + discrim_z + reject_z + ders_clarity_z + ders_goals_z + ders_impulse_z + ders_strat_z + ders_nonacc_z + samp + parent_edu)
car::vif(reg8)stress_helpless_z stress_efficacy_z discrim_z reject_z
3.308621 1.915228 1.986450 1.757206
ders_clarity_z ders_goals_z ders_impulse_z ders_strat_z
2.965829 3.415925 3.592802 5.516620
ders_nonacc_z samp parent_edu
3.504264 1.747683 1.530543 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 + ders_clarity_z + ders_goals_z + ders_impulse_z +
ders_strat_z + ders_nonacc_z + samp + parent_edu, data = df2)
Residuals:
Min 1Q Median 3Q Max
-1.80172 -0.50084 -0.00658 0.40452 2.73346
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.26061 0.17880 1.458 0.14688
stress_helpless_z 0.20659 0.09927 2.081 0.03899 *
stress_efficacy_z -0.03946 0.07561 -0.522 0.60249
discrim_z 0.08019 0.07793 1.029 0.30496
reject_z 0.48230 0.07485 6.443 1.24e-09 ***
ders_clarity_z -0.08676 0.09578 -0.906 0.36633
ders_goals_z 0.05042 0.10176 0.495 0.62092
ders_impulse_z -0.10454 0.10399 -1.005 0.31625
ders_strat_z 0.04710 0.12963 0.363 0.71682
ders_nonacc_z 0.09379 0.10313 0.909 0.36444
sampSONA -0.47451 0.15967 -2.972 0.00341 **
parent_edu -0.03798 0.06438 -0.590 0.55601
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7351 on 164 degrees of freedom
(8 observations deleted due to missingness)
Multiple R-squared: 0.4943, Adjusted R-squared: 0.4604
F-statistic: 14.57 on 11 and 164 DF, p-value: < 2.2e-16Effect of minority stress on IDM mediated by ER (sample 2, replication from sample 1)
Overall, the mediation model provided mixed support for the hypothesis, but are broadly consistent from sample 1. The effect of rejection sensitivity on IDM and the effect of discrimination on limited access to strategies are supported. The indirect effects are weakly supported, with lower bounds approaching zero. There is a pattern of full vs partial mediation, although the effects are small. One thing to keep in mind is the small sample (n2 = 184, n1 = 420)
# Fit the model
# fit <- sem(model, data = df2, se = "bootstrap", bootstrap = 5000)
# saveRDS(fit, file = "mediation_fit2.rds")
fit2 <- readRDS("mediation_fit2.rds")
# Summarize results
summary(fit2, fit.measures = TRUE, ci = TRUE)lavaan 0.6-21 ended normally after 12 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 14
Used Total
Number of observations 176 184
Model Test User Model:
Test statistic 12.570
Degrees of freedom 4
P-value (Chi-square) 0.014
Model Test Baseline Model:
Test statistic 253.770
Degrees of freedom 14
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.964
Tucker-Lewis Index (TLI) 0.875
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -875.666
Loglikelihood unrestricted model (H1) NA
Akaike (AIC) 1779.333
Bayesian (BIC) 1823.719
Sample-size adjusted Bayesian (SABIC) 1779.385
Root Mean Square Error of Approximation:
RMSEA 0.110
90 Percent confidence interval - lower 0.045
90 Percent confidence interval - upper 0.182
P-value H_0: RMSEA <= 0.050 0.062
P-value H_0: RMSEA >= 0.080 0.809
Standardized Root Mean Square Residual:
SRMR 0.045
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
idm_z ~
discrim_z (c1) 0.125 0.075 1.665 0.096 -0.024 0.268
reject_z (c2) 0.473 0.068 6.935 0.000 0.337 0.605
samp -0.360 0.164 -2.202 0.028 -0.670 -0.032
parent_ed -0.035 0.064 -0.551 0.582 -0.163 0.088
ders_strat_z ~
discrim_z (a1) 0.391 0.082 4.756 0.000 0.233 0.552
reject_z (a2) 0.239 0.085 2.820 0.005 0.071 0.398
samp 0.256 0.153 1.669 0.095 -0.053 0.549
parent_ed -0.093 0.066 -1.415 0.157 -0.222 0.037
idm_z ~
drs_strt_ (b) 0.163 0.078 2.090 0.037 0.017 0.323
Covariances:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
discrim_z ~~
reject_z 0.547 0.104 5.279 0.000 0.351 0.757
Variances:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
.idm_z 0.540 0.064 8.429 0.000 0.403 0.654
.ders_strat_z 0.684 0.089 7.654 0.000 0.501 0.849
discrim_z 1.004 0.108 9.319 0.000 0.791 1.213
reject_z 0.963 0.107 8.961 0.000 0.750 1.175
Defined Parameters:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
indirect_dscrm 0.064 0.036 1.780 0.075 0.006 0.146
indirect_rejct 0.039 0.023 1.697 0.090 0.002 0.091
indirect_diff 0.025 0.031 0.800 0.424 -0.019 0.103
total_discrim 0.189 0.070 2.693 0.007 0.049 0.324
total_reject 0.512 0.066 7.784 0.000 0.384 0.638graph_sem(fit2, layout = lay)
# Pool the datasets and test the interaction
df_combined <- bind_rows(
df1 %>% mutate(sample = 1),
df2 %>% mutate(sample = 2)
)
df_combined$sample <- as.factor(df_combined$sample)
interaction_model <- lm(idm_z ~ stress_helpless_z + stress_efficacy_z + discrim_z + reject_z + ders_clarity_z + ders_goals_z + ders_impulse_z + ders_strat_z + ders_nonacc_z + samp + sample + parent_edu,
data = df_combined)
summary(interaction_model)
Call:
lm(formula = idm_z ~ stress_helpless_z + stress_efficacy_z +
discrim_z + reject_z + ders_clarity_z + ders_goals_z + ders_impulse_z +
ders_strat_z + ders_nonacc_z + samp + sample + parent_edu,
data = df_combined)
Residuals:
Min 1Q Median 3Q Max
-2.29272 -0.50726 -0.02999 0.47038 2.84169
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.2978513 0.0964063 3.090 0.002100 **
stress_helpless_z 0.0152548 0.0473009 0.323 0.747186
stress_efficacy_z -0.0192621 0.0396716 -0.486 0.627477
discrim_z 0.2398834 0.0395591 6.064 2.39e-09 ***
reject_z 0.2700060 0.0372801 7.243 1.40e-12 ***
ders_clarity_z 0.0009856 0.0451490 0.022 0.982591
ders_goals_z 0.0373833 0.0517371 0.723 0.470240
ders_impulse_z -0.0246665 0.0459618 -0.537 0.591699
ders_strat_z 0.2209344 0.0603925 3.658 0.000277 ***
ders_nonacc_z 0.0411414 0.0497049 0.828 0.408173
sampSONA -0.4807169 0.0783569 -6.135 1.58e-09 ***
sample2 -0.0760445 0.0668917 -1.137 0.256077
parent_edu -0.0249619 0.0331409 -0.753 0.451633
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7343 on 581 degrees of freedom
(10 observations deleted due to missingness)
Multiple R-squared: 0.4714, Adjusted R-squared: 0.4605
F-statistic: 43.17 on 12 and 581 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: 11 x 5
term b1 se1 b2 se2
<chr> <dbl> <dbl> <dbl> <dbl>
1 stress_helpless_z -0.0396 0.0533 0.207 0.0993
2 stress_efficacy_z 0.00176 0.0465 -0.0395 0.0756
3 discrim_z 0.294 0.0455 0.0802 0.0779
4 reject_z 0.212 0.0429 0.482 0.0749
5 ders_clarity_z 0.0101 0.0512 -0.0868 0.0958
6 ders_goals_z 0.0111 0.0599 0.0504 0.102
7 ders_impulse_z 0.0219 0.0514 -0.105 0.104
8 ders_strat_z 0.292 0.0683 0.0471 0.130
9 ders_nonacc_z 0.0316 0.0560 0.0938 0.103
10 sampSONA -0.491 0.0894 -0.475 0.160
11 parent_edu -0.0352 0.0383 -0.0380 0.0644append_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
0.3294 -0.6587 3.3413 -0.6587 15.3413
tau^2 (estimated amount of total heterogeneity): 0.0240 (SE = 0.0429)
tau (square root of estimated tau^2 value): 0.1548
I^2 (total heterogeneity / total variability): 79.06%
H^2 (total variability / sampling variability): 4.77
Test for Heterogeneity:
Q(df = 1) = 4.7746, p-val = 0.0289
Model Results:
estimate se zval pval ci.lb ci.ub
0.0692 0.1223 0.5664 0.5711 -0.1704 0.3089
---
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.7417 -3.4834 0.5166 -3.4834 12.5166
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0056)
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.2156, p-val = 0.6424
Model Results:
estimate se zval pval ci.lb ci.ub
-0.0095 0.0396 -0.2410 0.8095 -0.0871 0.0681
---
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.4707 -0.9413 3.0587 -0.9413 15.0587
tau^2 (estimated amount of total heterogeneity): 0.0188 (SE = 0.0323)
tau (square root of estimated tau^2 value): 0.1370
I^2 (total heterogeneity / total variability): 82.17%
H^2 (total variability / sampling variability): 5.61
Test for Heterogeneity:
Q(df = 1) = 5.6092, p-val = 0.0179
Model Results:
estimate se zval pval ci.lb ci.ub
0.1964 0.1065 1.8451 0.0650 -0.0122 0.4051 .
---
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.2353 -0.4706 3.5294 -0.4706 15.5294
tau^2 (estimated amount of total heterogeneity): 0.0329 (SE = 0.0517)
tau (square root of estimated tau^2 value): 0.1812
I^2 (total heterogeneity / total variability): 89.82%
H^2 (total variability / sampling variability): 9.82
Test for Heterogeneity:
Q(df = 1) = 9.8238, p-val = 0.0017
Model Results:
estimate se zval pval ci.lb ci.ub
0.3401 0.1350 2.5185 0.0118 0.0754 0.6048 *
---
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 == "ders_clarity_z"],
meta_inputs$b2[meta_inputs$term == "ders_clarity_z"]),
sei = c(meta_inputs$se1[meta_inputs$term == "ders_clarity_z"],
meta_inputs$se2[meta_inputs$term == "ders_clarity_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.2495 -2.4990 1.5010 -2.4990 13.5010
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0083)
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.7960, p-val = 0.3723
Model Results:
estimate se zval pval ci.lb ci.ub
-0.0114 0.0452 -0.2527 0.8005 -0.1000 0.0771
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs <- append_meta(meta_inputs, "ders_clarity_z", meta)
dat <- data.frame(
yi = c(meta_inputs$b1[meta_inputs$term == "ders_goals_z"],
meta_inputs$b2[meta_inputs$term == "ders_goals_z"]),
sei = c(meta_inputs$se1[meta_inputs$term == "ders_goals_z"],
meta_inputs$se2[meta_inputs$term == "ders_goals_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.5086 -3.0173 0.9827 -3.0173 12.9827
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0099)
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.1112, p-val = 0.7388
Model Results:
estimate se zval pval ci.lb ci.ub
0.0212 0.0516 0.4104 0.6815 -0.0799 0.1223
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs <- append_meta(meta_inputs, "ders_goals_z", meta)
dat <- data.frame(
yi = c(meta_inputs$b1[meta_inputs$term == "ders_impulse_z"],
meta_inputs$b2[meta_inputs$term == "ders_impulse_z"]),
sei = c(meta_inputs$se1[meta_inputs$term == "ders_impulse_z"],
meta_inputs$se2[meta_inputs$term == "ders_impulse_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.9960 -1.9920 2.0080 -1.9920 14.0080
tau^2 (estimated amount of total heterogeneity): 0.0013 (SE = 0.0113)
tau (square root of estimated tau^2 value): 0.0355
I^2 (total heterogeneity / total variability): 15.76%
H^2 (total variability / sampling variability): 1.19
Test for Heterogeneity:
Q(df = 1) = 1.1871, p-val = 0.2759
Model Results:
estimate se zval pval ci.lb ci.ub
-0.0090 0.0543 -0.1660 0.8682 -0.1154 0.0974
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs <- append_meta(meta_inputs, "ders_impulse_z", meta)
dat <- data.frame(
yi = c(meta_inputs$b1[meta_inputs$term == "ders_strat_z"],
meta_inputs$b2[meta_inputs$term == "ders_strat_z"]),
sei = c(meta_inputs$se1[meta_inputs$term == "ders_strat_z"],
meta_inputs$se2[meta_inputs$term == "ders_strat_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.3357 -0.6714 3.3286 -0.6714 15.3286
tau^2 (estimated amount of total heterogeneity): 0.0192 (SE = 0.0423)
tau (square root of estimated tau^2 value): 0.1385
I^2 (total heterogeneity / total variability): 64.13%
H^2 (total variability / sampling variability): 2.79
Test for Heterogeneity:
Q(df = 1) = 2.7881, p-val = 0.0950
Model Results:
estimate se zval pval ci.lb ci.ub
0.1942 0.1198 1.6218 0.1048 -0.0405 0.4290
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs <- append_meta(meta_inputs, "ders_strat_z", meta)
dat <- data.frame(
yi = c(meta_inputs$b1[meta_inputs$term == "ders_nonacc_z"],
meta_inputs$b2[meta_inputs$term == "ders_nonacc_z"]),
sei = c(meta_inputs$se1[meta_inputs$term == "ders_nonacc_z"],
meta_inputs$se2[meta_inputs$term == "ders_nonacc_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.4297 -2.8595 1.1405 -2.8595 13.1405
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0097)
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.2809, p-val = 0.5961
Model Results:
estimate se zval pval ci.lb ci.ub
0.0458 0.0492 0.9297 0.3525 -0.0507 0.1422
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs <- append_meta(meta_inputs, "ders_nonacc_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.1218 -2.2435 1.7565 -2.2435 13.7565
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0237)
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.0081, p-val = 0.9282
Model Results:
estimate se zval pval ci.lb ci.ub
-0.4871 0.0780 -6.2413 <.0001 -0.6400 -0.3341 ***
---
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
2.0180 -4.0361 -0.0361 -4.0361 11.9639
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0040)
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.0014, p-val = 0.9699
Model Results:
estimate se zval pval ci.lb ci.ub
-0.0359 0.0329 -1.0896 0.2759 -0.1005 0.0287
---
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) | t² | SE(t²) | I² | Q | Q p | ß (pooled) | SE (pooled) | p | 95% CI LL | 95% CI UL |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Stress (Helplessness) | -0.04 | 0.05 | 0.21 | 0.10 | 0.02 | 0.04 | 79.06 | 4.77 | 0.029 | 0.07 | 0.12 | 0.571 | -0.17 | 0.31 |
| Stress (Efficacy) | 0.00 | 0.05 | -0.04 | 0.08 | 0.00 | 0.01 | 0.00 | 0.22 | 0.642 | -0.01 | 0.04 | 0.810 | -0.09 | 0.07 |
| Discrimination | 0.29 | 0.05 | 0.08 | 0.08 | 0.02 | 0.03 | 82.17 | 5.61 | 0.018 | 0.20 | 0.11 | 0.065 | -0.01 | 0.41 |
| Rejection Sensitivity | 0.21 | 0.04 | 0.48 | 0.07 | 0.03 | 0.05 | 89.82 | 9.82 | 0.002 | 0.34 | 0.14 | 0.012 | 0.08 | 0.60 |
| ER: Clarity | 0.01 | 0.05 | -0.09 | 0.10 | 0.00 | 0.01 | 0.00 | 0.80 | 0.372 | -0.01 | 0.05 | 0.800 | -0.10 | 0.08 |
| ER: Goals | 0.01 | 0.06 | 0.05 | 0.10 | 0.00 | 0.01 | 0.00 | 0.11 | 0.739 | 0.02 | 0.05 | 0.682 | -0.08 | 0.12 |
| ER: Impulse | 0.02 | 0.05 | -0.10 | 0.10 | 0.00 | 0.01 | 15.76 | 1.19 | 0.276 | -0.01 | 0.05 | 0.868 | -0.12 | 0.10 |
| ER: Strategies | 0.29 | 0.07 | 0.05 | 0.13 | 0.02 | 0.04 | 64.13 | 2.79 | 0.095 | 0.19 | 0.12 | 0.105 | -0.04 | 0.43 |
| ER: Non-Acceptance | 0.03 | 0.06 | 0.09 | 0.10 | 0.00 | 0.01 | 0.00 | 0.28 | 0.596 | 0.05 | 0.05 | 0.353 | -0.05 | 0.14 |
| Student Status | -0.49 | 0.09 | -0.47 | 0.16 | 0.00 | 0.02 | 0.00 | 0.01 | 0.928 | -0.49 | 0.08 | 0.000 | -0.64 | -0.33 |
| Parental Education | -0.04 | 0.04 | -0.04 | 0.06 | 0.00 | 0.00 | 0.00 | 0.00 | 0.970 | -0.04 | 0.03 | 0.276 | -0.10 | 0.03 |
# Extract parameters from both lavaan models
s1_med <- parameterEstimates(fit1, boot.ci.type = "bca.simple") %>%
filter(label %in% c("a1", "a2", "b", "c1", "c2",
"indirect_discrim", "indirect_reject",
"total_discrim", "total_reject"))
s2_med <- parameterEstimates(fit2, boot.ci.type = "bca.simple") %>%
filter(label %in% c("a1", "a2", "b", "c1", "c2",
"indirect_discrim", "indirect_reject",
"total_discrim", "total_reject"))
# Join into single dataframe
meta_inputs_med <- left_join(
s1_med %>% select(label, est, se) %>% rename(b1 = est, se1 = se),
s2_med %>% select(label, est, se) %>% rename(b2 = est, se2 = se),
by = "label"
) %>%
mutate(vi1 = se1^2,
vi2 = se2^2) %>%
rename(term = label)
# Preview
print(meta_inputs_med) term b1 se1 b2 se2 vi1 vi2
1 c1 0.296 0.045 0.125 0.075 0.002 0.006
2 c2 0.210 0.046 0.473 0.068 0.002 0.005
3 a1 0.316 0.054 0.391 0.082 0.003 0.007
4 a2 0.135 0.061 0.239 0.085 0.004 0.007
5 b 0.317 0.042 0.163 0.078 0.002 0.006
6 indirect_discrim 0.100 0.024 0.064 0.036 0.001 0.001
7 indirect_reject 0.043 0.019 0.039 0.023 0.000 0.001
8 total_discrim 0.397 0.044 0.189 0.070 0.002 0.005
9 total_reject 0.253 0.048 0.512 0.066 0.002 0.004# --- a1 path (discrim -> ders_strat) ---
dat_a1 <- data.frame(
yi = c(meta_inputs_med$b1[meta_inputs_med$term == "a1"],
meta_inputs_med$b2[meta_inputs_med$term == "a1"]),
sei = c(meta_inputs_med$se1[meta_inputs_med$term == "a1"],
meta_inputs_med$se2[meta_inputs_med$term == "a1"])
)
dat_a1$vi <- dat_a1$sei^2
meta_a1 <- rma(yi, vi, data = dat_a1, method = "REML")
summary(meta_a1)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
1.4493 -2.8986 1.1014 -2.8986 13.1014
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0069)
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.5889, p-val = 0.4428
Model Results:
estimate se zval pval ci.lb ci.ub
0.3387 0.0454 7.4617 <.0001 0.2497 0.4277 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs_med <- append_meta(meta_inputs_med, "a1", meta_a1)
# --- a2 path (reject -> ders_strat) ---
dat_a2 <- data.frame(
yi = c(meta_inputs_med$b1[meta_inputs_med$term == "a2"],
meta_inputs_med$b2[meta_inputs_med$term == "a2"]),
sei = c(meta_inputs_med$se1[meta_inputs_med$term == "a2"],
meta_inputs_med$se2[meta_inputs_med$term == "a2"])
)
dat_a2$vi <- dat_a2$sei^2
meta_a2 <- rma(yi, vi, data = dat_a2, method = "REML")
summary(meta_a2)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
1.1928 -2.3855 1.6145 -2.3855 13.6145
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0077)
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.9910, p-val = 0.3195
Model Results:
estimate se zval pval ci.lb ci.ub
0.1704 0.0494 3.4484 0.0006 0.0735 0.2672 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs_med <- append_meta(meta_inputs_med, "a2", meta_a2)
# --- b path (ders_strat -> idm) ---
dat_b <- data.frame(
yi = c(meta_inputs_med$b1[meta_inputs_med$term == "b"],
meta_inputs_med$b2[meta_inputs_med$term == "b"]),
sei = c(meta_inputs_med$se1[meta_inputs_med$term == "b"],
meta_inputs_med$se2[meta_inputs_med$term == "b"])
)
dat_b$vi <- dat_b$sei^2
meta_b <- rma(yi, vi, data = dat_b, method = "REML")
summary(meta_b)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
0.7950 -1.5901 2.4099 -1.5901 14.4099
tau^2 (estimated amount of total heterogeneity): 0.0080 (SE = 0.0169)
tau (square root of estimated tau^2 value): 0.0896
I^2 (total heterogeneity / total variability): 67.26%
H^2 (total variability / sampling variability): 3.05
Test for Heterogeneity:
Q(df = 1) = 3.0539, p-val = 0.0805
Model Results:
estimate se zval pval ci.lb ci.ub
0.2541 0.0760 3.3435 0.0008 0.1051 0.4030 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs_med <- append_meta(meta_inputs_med, "b", meta_b)
# --- c1 path (discrim -> idm direct) ---
dat_c1 <- data.frame(
yi = c(meta_inputs_med$b1[meta_inputs_med$term == "c1"],
meta_inputs_med$b2[meta_inputs_med$term == "c1"]),
sei = c(meta_inputs_med$se1[meta_inputs_med$term == "c1"],
meta_inputs_med$se2[meta_inputs_med$term == "c1"])
)
dat_c1$vi <- dat_c1$sei^2
meta_c1 <- rma(yi, vi, data = dat_c1, method = "REML")
summary(meta_c1)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
0.6929 -1.3858 2.6142 -1.3858 14.6142
tau^2 (estimated amount of total heterogeneity): 0.0108 (SE = 0.0207)
tau (square root of estimated tau^2 value): 0.1039
I^2 (total heterogeneity / total variability): 73.73%
H^2 (total variability / sampling variability): 3.81
Test for Heterogeneity:
Q(df = 1) = 3.8060, p-val = 0.0511
Model Results:
estimate se zval pval ci.lb ci.ub
0.2214 0.0849 2.6079 0.0091 0.0550 0.3879 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs_med <- append_meta(meta_inputs_med, "c1", meta_c1)
# --- c2 path (reject -> idm direct) ---
dat_c2 <- data.frame(
yi = c(meta_inputs_med$b1[meta_inputs_med$term == "c2"],
meta_inputs_med$b2[meta_inputs_med$term == "c2"]),
sei = c(meta_inputs_med$se1[meta_inputs_med$term == "c2"],
meta_inputs_med$se2[meta_inputs_med$term == "c2"])
)
dat_c2$vi <- dat_c2$sei^2
meta_c2 <- rma(yi, vi, data = dat_c2, method = "REML")
summary(meta_c2)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
0.2618 -0.5237 3.4763 -0.5237 15.4763
tau^2 (estimated amount of total heterogeneity): 0.0313 (SE = 0.0490)
tau (square root of estimated tau^2 value): 0.1768
I^2 (total heterogeneity / total variability): 90.18%
H^2 (total variability / sampling variability): 10.18
Test for Heterogeneity:
Q(df = 1) = 10.1823, p-val = 0.0014
Model Results:
estimate se zval pval ci.lb ci.ub
0.3368 0.1316 2.5590 0.0105 0.0788 0.5947 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs_med <- append_meta(meta_inputs_med, "c2", meta_c2)
# --- Indirect effect: discrimination ---
dat_ind_discrim <- data.frame(
yi = c(meta_inputs_med$b1[meta_inputs_med$term == "indirect_discrim"],
meta_inputs_med$b2[meta_inputs_med$term == "indirect_discrim"]),
sei = c(meta_inputs_med$se1[meta_inputs_med$term == "indirect_discrim"],
meta_inputs_med$se2[meta_inputs_med$term == "indirect_discrim"])
)
dat_ind_discrim$vi <- dat_ind_discrim$sei^2
meta_ind_discrim <- rma(yi, vi, data = dat_ind_discrim, method = "REML")
summary(meta_ind_discrim)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
2.2154 -4.4307 -0.4307 -4.4307 11.5693
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0013)
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.7211, p-val = 0.3958
Model Results:
estimate se zval pval ci.lb ci.ub
0.0891 0.0198 4.5106 <.0001 0.0504 0.1278 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs_med <- append_meta(meta_inputs_med, "indirect_discrim", meta_ind_discrim)
# --- Indirect effect: rejection ---
dat_ind_reject <- data.frame(
yi = c(meta_inputs_med$b1[meta_inputs_med$term == "indirect_reject"],
meta_inputs_med$b2[meta_inputs_med$term == "indirect_reject"]),
sei = c(meta_inputs_med$se1[meta_inputs_med$term == "indirect_reject"],
meta_inputs_med$se2[meta_inputs_med$term == "indirect_reject"])
)
dat_ind_reject$vi <- dat_ind_reject$sei^2
meta_ind_reject <- rma(yi, vi, data = dat_ind_reject, method = "REML")
summary(meta_ind_reject)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
2.9241 -5.8482 -1.8482 -5.8482 10.1518
tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0006)
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.0174, p-val = 0.8950
Model Results:
estimate se zval pval ci.lb ci.ub
0.0412 0.0148 2.7813 0.0054 0.0122 0.0702 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs_med <- append_meta(meta_inputs_med, "indirect_reject", meta_ind_reject)
# --- Total effect: discrimination ---
dat_tot_discrim <- data.frame(
yi = c(meta_inputs_med$b1[meta_inputs_med$term == "total_discrim"],
meta_inputs_med$b2[meta_inputs_med$term == "total_discrim"]),
sei = c(meta_inputs_med$se1[meta_inputs_med$term == "total_discrim"],
meta_inputs_med$se2[meta_inputs_med$term == "total_discrim"])
)
dat_tot_discrim$vi <- dat_tot_discrim$sei^2
meta_tot_discrim <- rma(yi, vi, data = dat_tot_discrim, method = "REML")
summary(meta_tot_discrim)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
0.4998 -0.9996 3.0004 -0.9996 15.0004
tau^2 (estimated amount of total heterogeneity): 0.0181 (SE = 0.0305)
tau (square root of estimated tau^2 value): 0.1347
I^2 (total heterogeneity / total variability): 84.15%
H^2 (total variability / sampling variability): 6.31
Test for Heterogeneity:
Q(df = 1) = 6.3087, p-val = 0.0120
Model Results:
estimate se zval pval ci.lb ci.ub
0.3001 0.1035 2.8984 0.0038 0.0972 0.5030 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs_med <- append_meta(meta_inputs_med, "total_discrim", meta_tot_discrim)
# --- Total effect: rejection ---
dat_tot_reject <- data.frame(
yi = c(meta_inputs_med$b1[meta_inputs_med$term == "total_reject"],
meta_inputs_med$b2[meta_inputs_med$term == "total_reject"]),
sei = c(meta_inputs_med$se1[meta_inputs_med$term == "total_reject"],
meta_inputs_med$se2[meta_inputs_med$term == "total_reject"])
)
dat_tot_reject$vi <- dat_tot_reject$sei^2
meta_tot_reject <- rma(yi, vi, data = dat_tot_reject, method = "REML")
summary(meta_tot_reject)
Random-Effects Model (k = 2; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
0.2770 -0.5540 3.4460 -0.5540 15.4460
tau^2 (estimated amount of total heterogeneity): 0.0303 (SE = 0.0476)
tau (square root of estimated tau^2 value): 0.1742
I^2 (total heterogeneity / total variability): 90.16%
H^2 (total variability / sampling variability): 10.17
Test for Heterogeneity:
Q(df = 1) = 10.1672, p-val = 0.0014
Model Results:
estimate se zval pval ci.lb ci.ub
0.3785 0.1296 2.9193 0.0035 0.1244 0.6326 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1meta_inputs_med <- append_meta(meta_inputs_med, "total_reject", meta_tot_reject)term_labels <- c(
"c1" = "C1 (Discrimination > IDM)",
"c2" = "C2 (Rejection > IDM)",
"a1" = "A1 (Discrimination > ER: Strategies)",
"a2" = "A2 (Rejection > ER: Strategies)",
"b" = "B (ER: Strategies > IDM)",
"indirect_discrim" = "Indirect Effect: Discrimination > IDM",
"indirect_reject" = "Indirect Effect: Rejection > IDM",
"total_discrim" = "Total Effect: Discrimination",
"total_reject" = "Total Effect: Rejection"
)
# 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_med %>%
select(-vi1, -vi2) %>%
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)| Predictor | ß (S1) | SE (S1) | ß (S2) | SE (S2) | t² | SE(t²) | I² | Q | Q p | ß (pooled) | SE (pooled) | p | 95% CI LL | 95% CI UL |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| C1 (Discrimination > IDM) | 0.30 | 0.05 | 0.13 | 0.08 | 0.01 | 0.02 | 73.73 | 3.81 | 0.051 | 0.22 | 0.08 | 0.009 | 0.06 | 0.39 |
| C2 (Rejection > IDM) | 0.21 | 0.05 | 0.47 | 0.07 | 0.03 | 0.05 | 90.18 | 10.18 | 0.001 | 0.34 | 0.13 | 0.010 | 0.08 | 0.59 |
| A1 (Discrimination > ER: Strategies) | 0.32 | 0.05 | 0.39 | 0.08 | 0.00 | 0.01 | 0.00 | 0.59 | 0.443 | 0.34 | 0.05 | 0.000 | 0.25 | 0.43 |
| A2 (Rejection > ER: Strategies) | 0.14 | 0.06 | 0.24 | 0.08 | 0.00 | 0.01 | 0.00 | 0.99 | 0.320 | 0.17 | 0.05 | 0.001 | 0.07 | 0.27 |
| B (ER: Strategies > IDM) | 0.32 | 0.04 | 0.16 | 0.08 | 0.01 | 0.02 | 67.26 | 3.05 | 0.081 | 0.25 | 0.08 | 0.001 | 0.11 | 0.40 |
| Indirect Effect: Discrimination > IDM | 0.10 | 0.02 | 0.06 | 0.04 | 0.00 | 0.00 | 0.00 | 0.72 | 0.396 | 0.09 | 0.02 | 0.000 | 0.05 | 0.13 |
| Indirect Effect: Rejection > IDM | 0.04 | 0.02 | 0.04 | 0.02 | 0.00 | 0.00 | 0.00 | 0.02 | 0.895 | 0.04 | 0.01 | 0.005 | 0.01 | 0.07 |
| Total Effect: Discrimination | 0.40 | 0.04 | 0.19 | 0.07 | 0.02 | 0.03 | 84.15 | 6.31 | 0.012 | 0.30 | 0.10 | 0.004 | 0.10 | 0.50 |
| Total Effect: Rejection | 0.25 | 0.05 | 0.51 | 0.07 | 0.03 | 0.05 | 90.16 | 10.17 | 0.001 | 0.38 | 0.13 | 0.004 | 0.12 | 0.63 |
IDM predicts misinformation acceptance, accurate information acceptance, and increased susceptability to misinformation (sample 2)
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-05# # Plot residuals separately by group
# df2 %>%
# mutate(resid = residuals(reg1)) %>%
# ggplot(aes(x = resid, fill = samp)) +
# geom_density(alpha = 0.5) +
# theme_minimal()
#
# df2 %>%
# mutate(resid = residuals(reg1)) %>%
# group_by(samp) %>%
# summarise(
# shapiro_p = shapiro.test(resid)$p.value
# )reg2 <- 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-09Minority stress predicts misinformation acceptance, mediated by IDM (sample 2)
Condition (misinformation exposure) will predict less social support, higher science populism, higher anomie, controlling for student status, parent education status, and IDM (sample 2) Susceptibility to misinformation acceptance will be associated with lower perceived social support, higher science populism, and higher anomie, even when controlling for misinformation exposure, controlling for student status, parent education status, and IDM (sample 2)
Social Support
Histograms
Code
Code
CFA
Code
Code