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