options(digits = 2)
rm(list=ls())
## Install "pacman" package if not installed
# (remove the # symbol from the line below):
# install.packages("pacman")
## Load R packages:
pacman::p_load(data.table, tidyverse, haven, labelled, vtable,
psych, scales, weights, clipr, forcats,
stargazer, ggthemes, ggcharts, geomtextpath,
corrplot, tm, readxl)
## Import dataset:
ds <- read_excel("~/Desktop/oxford/data/gambia/GAMBIA Data Work.xlsx")ds$Age1 <- as.numeric(gsub("\\D", "", ds$Age))
ds$age <- ifelse(ds$Age1 %in% 18:90, ds$Age1, NA)
summary(ds$age) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
18 27 32 37 45 90 8 ds %>% drop_na(age)%>%
ggplot(aes(x = age))+
geom_histogram(color = "black",
fill = "gray",
bins = 50)+
geom_textvline(label = "Mean = 37.00",
xintercept = 37.00,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Age",
y = "Frequency",
title = "Age distribution: Gambia")+
theme_bw()
ds$gender <- ifelse(ds$Gender == "Male", "Male",
ifelse(ds$Gender == "Female", "Female", NA))
lp02 <- ds %>% drop_na(gender) %>%
lollipop_chart(x = gender,
line_color = "black",
point_color = "black")+
labs(y = "Frequency",
x = "",
title = "Gender distribution: Gambia")+
theme_bw()
lp02
bp01 <- ds %>% drop_na(gender, age) %>%
ggplot(aes(y = age,
x = gender))+
geom_boxplot(fill = "grey")+
labs(y = "Age",
x = "",
title = "Age distribution by gender: Gambia")+
coord_flip()+
theme_bw()
bp01
prov <- as.data.frame(table(ds$Province))
prov$Var1 <- as.character(prov$Var1)
prov$province <- ifelse(prov$Freq < 4, "Other", prov$Var1)
l1 <- as.list(prov$province)
ds2 <- ds[, c("Province")]
ds2$province <- ifelse(ds2$Province %in% l1, ds2$Province, "Other")
ds2$province <- ifelse(ds2$province == "", NA, ds2$province)
ds2$province <- as.character(ds2$province)
lp03 <- ds2 %>% drop_na(province) %>%
lollipop_chart(x = province,
line_color = "black",
point_color = "black")+
labs(y = "Frequency",
x = "",
title = "Province distribution: Gambia")+
theme_bw()
lp03
ds$hh_wealth <- haven::as_factor(ds$Household_level_of_wealth)
## Convert to sentence case:
ds$hh_wealth <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$hh_wealth, perl = TRUE)
ds$hh_wealth <- factor(ds$hh_wealth,
levels=c("Much lower", "Slightly lower", "Average",
"Slightly higher", "Much higher"))
table(ds$hh_wealth)
Much lower Slightly lower Average Slightly higher Much higher
21 40 129 39 4 lp04 <- ds %>% drop_na(hh_wealth) %>%
lollipop_chart(x = hh_wealth,
line_color = "black",
point_color = "black")+
labs(y = "Frequency",
x = "",
title = "Household wealth distribution: Gambia")+
theme_bw()
lp04
ds$ses1 <- haven::as_factor(ds$Socioeconomic_status)
## Convert to sentence case:
ds$ses1 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ses1, perl = TRUE)
ds$ses <- factor(ds$ses1,
levels=c("Lower middle", "Middle", "Upper middle", "Upper"))
table(ds$ses)
Lower middle Middle Upper middle Upper
42 139 40 9 lp05 <- ds %>% drop_na(ses) %>%
lollipop_chart(x = ses,
line_color = "black",
point_color = "black")+
labs(y = "Frequency",
x = "",
title = "Socioeconomic status: Gambia")+
theme_bw()
lp05
ds$income <- as.numeric(gsub("\\D", "", ds$Annual_income))
summary(ds$income) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
300 64250 150000 800599 500000 54000000 35 ## Remove outliers (greater than 5 million):
ds <- ds[ds$income < 5000000, ]
summary(ds$income) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
300 60000 134000 402137 500000 4000000 35 ds %>% drop_na(income)%>%
ggplot(aes(x = income))+
geom_histogram(color = "black",
fill = "gray",
bins = 80)+
geom_textvline(label = "Mean = 113K",
xintercept = 112746,
vjust = 1.3,
lwd = 1.05,
linetype = 2)+
scale_x_continuous(labels = scales::unit_format(unit = "K",
scale = 1e-3))+
labs(x = "Annual income",
y = "Frequency",
title = "Annual income distribution: Gambia")+
theme_bw()
## HH Wealth:
ds$hh_wealth <- haven::as_factor(ds$Household_level_of_wealth)
## Convert to sentence case:
ds$hh_wealth <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$hh_wealth, perl = TRUE)
ds$hh_wealth <- factor(ds$hh_wealth,
levels=c("Much lower", "Slightly lower", "Average",
"Slightly higher", "Much higher"))
## SES:
ds$ses1 <- haven::as_factor(ds$Socioeconomic_status)
## Convert to sentence case:
ds$ses1 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ses1, perl = TRUE)
ds$ses <- factor(ds$ses1,
levels=c("Lower middle", "Middle", "Upper middle", "Upper"))
ds$fses <- as.numeric(ds$ses)
ds$fhhw <- as.numeric(ds$hh_wealth)
ds3 <- cbind.data.frame(ds$fses, ds$fhhw, ds$income)
ds3 <- na.omit(ds3)
ds3$Annual_income <- ds3$`ds$income`
ds3$Household_wealth <- ds3$`ds$fhhw`
ds3$SES <- ds3$`ds$fses`
mtx <- cor(ds3[, c(4:6)])
corrplot(mtx, method = "number", number.cex = 1.0,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
ds$jobnature <- haven::as_factor(ds$Job_Nature)
ds$jobnature <- as.character(ds$jobnature)
ds$jobnature <- ifelse(ds$jobnature == "", NA, ds$jobnature)
lp06 <- ds %>% drop_na(jobnature) %>%
lollipop_chart(x = jobnature,
line_color = "black",
point_color = "black")+
labs(y = "Frequency",
x = "",
title = "Nature of employment: Gambia")+
theme_bw()
lp06
ds$r1 <- ds$Religious_Affiliation
ds$religion <- ifelse(ds$r1 == "Christian-Catholic", "Christian: Catholic",
ifelse(ds$r1 == "Christian-Other", "Christian: Other",
ifelse(ds$r1 == "Christian-Protestant", "Christian: Protestant",
ifelse(ds$r1 == "Muslim-Sunni", "Muslim: Sunni",
ifelse(ds$r1 == "Muslim", "Muslim: Other",
ifelse(ds$r1 == "Other-Muslim", "Muslim: Other", NA))))))
table(ds$religion)
Christian: Catholic Christian: Other Christian: Protestant
12 4 4
Muslim: Other Muslim: Sunni
3 173 lp07 <- ds %>% drop_na(religion) %>%
lollipop_chart(x = religion,
line_color = "black",
point_color = "black")+
labs(y = "Frequency",
x = "",
title = "Religious affiliation: Gambia")+
theme_bw()
lp07
ds$married <- ifelse(ds$Marital_status == "Married", "Married",
ifelse(ds$Marital_status == "Single", "Unmarried",
ifelse(ds$Marital_status == "In a relationship (not married)", "Unmarried",
"Other")))
table(ds$married)
Married Other Unmarried
129 2 63 bp01 <- ds %>% drop_na(married) %>%
ggplot(aes(x = fct_rev(fct_infreq(married))))+
geom_bar(fill = "black")+
labs(x = "",
y = "Frequency",
title = "Marital status")+
coord_flip()+
theme_bw()
# bp01
lp08 <- ds %>% drop_na(married) %>%
lollipop_chart(x = married,
line_color = "black",
point_color = "black")+
labs(y = "Frequency",
x = "",
title = "Marital status: Gambia")+
theme_bw()
lp08
ds$children <- as.numeric(gsub("\\D", "", ds$Number_of_Children))
summary(ds$children) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0 2 4 4 6 25 104 ds %>% drop_na(children)%>%
ggplot(aes(x = children))+
geom_bar(color = "black",
fill = "gray",
width = 0.75)+
geom_textvline(label = "Mean = 4.00",
xintercept = 4.00,
vjust = 1.3,
lwd = 1.05,
linetype = 2)+
labs(x = "Number of children",
y = "Frequency",
title = "Number of children: Gambia")+
theme_bw()
Note: I did my best to clean up the variable, but the categories need to be checked by an expert
eth <- as.data.frame(table(ds$Ethnicity))
eth$Var1 <- as.character(eth$Var1)
eth$ethnicity <- ifelse(eth$Freq < 2, "Other", eth$Var1)
l1 <- as.list(eth$ethnicity)
ds2 <- ds[, c("Ethnicity")]
ds2$ethnicity <- ifelse(ds2$Ethnicity %in% l1, ds2$Ethnicity, "Other")
ds2$ethnicity <- ifelse(ds2$ethnicity == "Serere", "Serer",
ifelse(ds2$ethnicity == "", NA, ds2$ethnicity))
table(ds2$ethnicity)
Balanta Fula Jola Karonika Mandinka Manjago Other Serer
3 31 45 3 74 5 48 4
Wollof
19 lp08 <- ds2 %>% drop_na(ethnicity) %>%
lollipop_chart(x = ethnicity,
line_color = "black",
point_color = "black")+
labs(y = "Frequency",
x = "",
title = "Ethnic distribution: Gambia")+
theme_bw()
lp08
Note: This needs a lot of work from someone familiar with Gambian education. The spreadsheet to edit can be found here:
Spreadsheet for cleaning: Gambia
In addition, education and ethnicity for Gambia has not been edited by anyone. The spreadsheet for Uganda is here:
Spreadsheet for cleaning: Uganda
ds$education <- as.character(haven::as_factor(ds$Level_of_education))
table(ds$education)
Advance Diploma
1
Advance Diploma in Business Studies
1
Advance Diploma in Management
3
Advance Diploma in Nursing
1
Advance Diploma in Tourism Management
1
B.S.C
1
Bachelor Degree
1
Bachelor in Business Administration
1
Bachelor's Degree
2
Bachelor's Degree in Political Science
1
Bachelors of Edu.
1
Bsc
1
BSc
2
BSC
4
BSc (Arabic Education)
1
BSc (Hons) Management
1
BSc in Accounting
1
BSc in Economic
1
BSc in Economica
1
BSc in Education
1
BSc in Human Resource
1
BSc in Management
1
BSc in Political Science
1
Certificate in Management Studies
1
Certificate Level
1
Certificate on Animal Health Product
1
Chartered Banker
1
College
8
College (HTC)
1
College Certificate
1
Degree in Islamic Studies
1
Diploma
6
Diploma 2 Business Administration
1
Diploma in Community Policing
1
Diploma in ICT
2
Diploma in Management
3
Diploma in Theology
1
Diploma in Travel and Tourism
1
Diplomat
1
Form 3 (Old education System)
2
Form 3 (Old Education System)
1
Form 4
1
Form 5 (Old education System)
2
Form 5 (Old Education System)
1
Gambia College (PTC)
1
GCE Level(WASSCE)
1
Grade 8 ( Upper Basic)
1
Graduate Diploma
1
High School
4
High School (WASSCE)
1
High School Certificate
1
High School Level
1
High School-Islamic Education
1
Higher Diploma
1
Higher Diploma in Education (HDC)
1
Higher Teacher Certificate (HTC)
7
Higher Teacher's Certificate (HTC)
5
Hotel School
1
HTC
2
HTC Primary
1
Informal Education (Madrasa)
1
Isachetovs Degree
1
Islamic School
1
Islamic Traditional School
1
Junior Sec School
1
Junior Secondary
1
LLB
1
Masters Degree
2
Masters in Economics
1
Memorized the Quran (Informal Education)
1
MSC
1
None
1
O Level
2
Ordinary Levels
1
Post Graduate
1
Primary School
3
Primary Teacher Certificate (PTC)
2
Primary Teacher's Certificate (PTC)
1
PTC
1
Quranic High School Graduate
1
Quranic School
1
Secondary level
1
Secondary School
6
Senior School Graduation
1
Senior Secondary School
9
Senior Secondary(WASSCE Certificate)
1
Tertiary
2
Tertiary Education
10
Tertiary level
1
Traditional Quranic Schools
1
Udergraduate
1
Undergraduate Degree in Communications
1
University
7
University Degree
6
University of the Gambia
1
WASSCE
15 ## Four ingroup fusion items:
# I have a deep emotional bond with the [ingroup].
ds$IGF01 <- as.numeric(ds$`Q11.1 Group Bonds`)
# I am strong because of the [ingroup].
ds$IGF02 <- as.numeric(ds$`Q11.2 Group Bonds`)
# I make the [ingroup] strong.
ds$IGF03 <- as.numeric(ds$`Q11.3 Group Bonds`)
# I am one with the [ingroup]
ds$IGF04 <- as.numeric(ds$`Q11.4 Group Bonds`)
## Four outgroup fusion items:
# I have a deep emotional bond with the [outgroup].
ds$OGF01 <- as.numeric(ds$`Q11.5 Group Bonds`)
# I am strong because of the [outgroup].
ds$OGF02 <- as.numeric(ds$`Q11.6 Group Bonds`)
# I make the [outgroup] strong.
ds$OGF03 <- as.numeric(ds$`Q11.7 Group Bonds`)
# I am one with the [outgroup].
ds$OGF04 <- as.numeric(ds$`Q11.8 Group Bonds`)
## Four ingroup identification items:
# I identify with the [ingroup].
ds$IGI01 <- as.numeric(ds$`Q11.9 Group Bonds`)
# I have a lot in common with the [ingroup].
ds$IGI02 <- as.numeric(ds$`Q11.10 Group Bonds`)
# I connect with the values of the [ingroup].
ds$IGI03 <- as.numeric(ds$`Q11.11 Group Bonds`)
# I feel a sense of belonging with the [ingroup].
ds$IGI04 <- as.numeric(ds$`Q11.12 Group Bonds`)
## Four outgroup identification items:
# I identify with the [outgroup].
ds$OGI01 <- as.numeric(ds$`Q11.13 Group Bonds`)
# I have a lot in common with the [outgroup].
ds$OGI02 <- as.numeric(ds$`Q11.14 Group Bonds`)
# I connect with the values of the [outgroup].
ds$OGI03 <- as.numeric(ds$`Q11.15 Group Bonds`)
# I feel a sense of belonging with the [outgroup].
ds$OGI04 <- as.numeric(ds$`Q11.16 Group Bonds`)## Bonds dataframe:
bonds <- cbind.data.frame(ds$IGF01, ds$IGF02, ds$IGF03, ds$IGF04,
ds$IGI01, ds$IGI02, ds$IGI03, ds$IGI04,
ds$OGF01, ds$OGF02, ds$OGF03, ds$OGF04,
ds$OGI01, ds$OGI02, ds$OGI03, ds$OGI04)
names(bonds) <- sub('ds\\$', '', names(bonds))
bonds <- na.omit(bonds)
mtx1 <- cor(bonds[, c(1:16)])Shorthand for item names:
IG = Ingroup, OG = Outgroup
F = Fusion, I = Identification
So for example:
IGF01 = “ingroup fusion: item 1”
OGI04 = “outgroup identification: item 4”, and so on
corrplot(mtx1, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
The correlation plot above shows the correlation coefficients between the sixteen different items (eight for ingroup fusion/identification and eight for outgroup fusion/identification).
KMO(r=cor(bonds))Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = cor(bonds))
Overall MSA = 0.87
MSA for each item =
IGF01 IGF02 IGF03 IGF04 IGI01 IGI02 IGI03 IGI04 OGF01 OGF02 OGF03 OGF04 OGI01
0.85 0.83 0.79 0.75 0.79 0.83 0.83 0.86 0.87 0.86 0.89 0.90 0.91
OGI02 OGI03 OGI04
0.92 0.88 0.90 cortest.bartlett(bonds)$chisq
[1] 2048
$p.value
[1] 0
$df
[1] 120parallel <- fa.parallel(bonds)
Parallel analysis suggests that the number of factors = 2 and the number of components = 2 Scree plot suggests a two factor structure. We will proceed with promax rotation, which assumes that the items are inter-correlated (that is, not independent from each other).
fit02 <- factanal(bonds, 2, rotation="promax")
fit02
Call:
factanal(x = bonds, factors = 2, rotation = "promax")
Uniquenesses:
IGF01 IGF02 IGF03 IGF04 IGI01 IGI02 IGI03 IGI04 OGF01 OGF02 OGF03 OGF04 OGI01
0.43 0.47 0.67 0.71 0.66 0.49 0.27 0.35 0.47 0.22 0.23 0.31 0.27
OGI02 OGI03 OGI04
0.38 0.33 0.26
Loadings:
Factor1 Factor2
IGF01 0.748
IGF02 0.726
IGF03 0.576
IGF04 0.541
IGI01 0.577
IGI02 0.715
IGI03 0.855
IGI04 0.794
OGF01 0.736
OGF02 0.877
OGF03 0.871
OGF04 0.832
OGI01 0.836 -0.106
OGI02 0.797
OGI03 0.821
OGI04 0.861
Factor1 Factor2
SS loadings 5.52 3.94
Proportion Var 0.34 0.25
Cumulative Var 0.34 0.59
Factor Correlations:
Factor1 Factor2
Factor1 1.000 -0.137
Factor2 -0.137 1.000
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 398 on 89 degrees of freedom.
The p-value is 0.00000000000000000000000000000000000000011 We want to split the ingroup/outgroup fusion/identification items into distinct factors. We can examine what the three and four factor solutions looks like.
fit03 <- factanal(bonds, 3, rotation="promax")
fit03
Call:
factanal(x = bonds, factors = 3, rotation = "promax")
Uniquenesses:
IGF01 IGF02 IGF03 IGF04 IGI01 IGI02 IGI03 IGI04 OGF01 OGF02 OGF03 OGF04 OGI01
0.40 0.43 0.65 0.68 0.65 0.50 0.28 0.33 0.47 0.15 0.11 0.30 0.28
OGI02 OGI03 OGI04
0.25 0.16 0.19
Loadings:
Factor1 Factor2 Factor3
IGF01 0.760
IGF02 0.741 -0.122
IGF03 0.583
IGF04 0.541 0.195
IGI01 0.578
IGI02 0.710
IGI03 0.850
IGI04 0.791 0.182
OGF01 0.726
OGF02 0.915 -0.202
OGF03 0.927 -0.265
OGF04 0.840
OGI01 0.826 -0.105
OGI02 0.754 0.373
OGI03 0.782 0.418
OGI04 0.826 0.290
Factor1 Factor2 Factor3
SS loadings 5.49 3.96 0.64
Proportion Var 0.34 0.25 0.04
Cumulative Var 0.34 0.59 0.63
Factor Correlations:
Factor1 Factor2 Factor3
Factor1 1.0000 -0.1415 -0.0869
Factor2 -0.1415 1.0000 0.0389
Factor3 -0.0869 0.0389 1.0000
Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 240 on 75 degrees of freedom.
The p-value is 0.00000000000000000042 fit04 <- factanal(bonds, 4, rotation="promax")
fit04
Call:
factanal(x = bonds, factors = 4, rotation = "promax")
Uniquenesses:
IGF01 IGF02 IGF03 IGF04 IGI01 IGI02 IGI03 IGI04 OGF01 OGF02 OGF03 OGF04 OGI01
0.43 0.46 0.48 0.25 0.65 0.49 0.19 0.31 0.45 0.15 0.11 0.30 0.27
OGI02 OGI03 OGI04
0.25 0.16 0.18
Loadings:
Factor1 Factor2 Factor3 Factor4
IGF01 0.643 0.101 -0.109
IGF02 0.590 0.126 -0.149
IGF03 0.167 0.533 -0.163
IGF04 0.907 0.160
IGI01 0.615
IGI02 0.749
IGI03 1.027 -0.165
IGI04 0.814 0.195
OGF01 0.709 0.195
OGF02 0.912 -0.200
OGF03 0.929 -0.137 -0.256
OGF04 0.838
OGI01 0.834 -0.130
OGI02 0.744 0.161 0.397
OGI03 0.778 0.448
OGI04 0.834 0.123 0.334
Factor1 Factor2 Factor3 Factor4
SS loadings 5.47 3.48 1.28 0.729
Proportion Var 0.34 0.22 0.08 0.046
Cumulative Var 0.34 0.56 0.64 0.685
Factor Correlations:
Factor1 Factor2 Factor3 Factor4
Factor1 1.0000 -0.135 -0.0839 -0.0213
Factor2 -0.1354 1.000 0.3182 -0.6302
Factor3 -0.0839 0.318 1.0000 -0.2679
Factor4 -0.0213 -0.630 -0.2679 1.0000
Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 180 on 62 degrees of freedom.
The p-value is 0.00000000000019 The fourth item has been removed from each sub-scale:
## Remove 4th item from all subscales and retry:
## New dataframe:
bonds <- cbind.data.frame(ds$IGF01, ds$IGF02, ds$IGF03,
ds$IGI01, ds$IGI02, ds$IGI03,
ds$OGF01, ds$OGF02, ds$OGF03,
ds$OGI01, ds$OGI02, ds$OGI03)
names(bonds) <- sub('ds\\$', '', names(bonds))
bonds <- na.omit(bonds)
mtx1 <- cor(bonds[, c(1:12)])
corrplot(mtx1, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
parallel <- fa.parallel(bonds)
Parallel analysis suggests that the number of factors = 2 and the number of components = 2 # Two factor model
fit04 <- factanal(bonds, 2, rotation="promax")
fit04
Call:
factanal(x = bonds, factors = 2, rotation = "promax")
Uniquenesses:
IGF01 IGF02 IGF03 IGI01 IGI02 IGI03 OGF01 OGF02 OGF03 OGI01 OGI02 OGI03
0.37 0.41 0.69 0.65 0.50 0.37 0.49 0.16 0.19 0.26 0.45 0.40
Loadings:
Factor1 Factor2
IGF01 0.787
IGF02 0.775
IGF03 0.565
IGI01 0.587
IGI02 0.704
IGI03 0.792
OGF01 0.718
OGF02 0.920
OGF03 0.899
OGI01 0.839
OGI02 0.746
OGI03 0.771
Factor1 Factor2
SS loadings 4.03 3.02
Proportion Var 0.34 0.25
Cumulative Var 0.34 0.59
Factor Correlations:
Factor1 Factor2
Factor1 1.000 -0.146
Factor2 -0.146 1.000
Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 207 on 43 degrees of freedom.
The p-value is 0.000000000000000000000023 # Three factor model:
fit05 <- factanal(bonds, 3, rotation="promax")
fit05
Call:
factanal(x = bonds, factors = 3, rotation = "promax")
Uniquenesses:
IGF01 IGF02 IGF03 IGI01 IGI02 IGI03 OGF01 OGF02 OGF03 OGI01 OGI02 OGI03
0.378 0.406 0.678 0.636 0.489 0.354 0.482 0.049 0.203 0.277 0.188 0.202
Loadings:
Factor1 Factor2 Factor3
IGF01 0.782
IGF02 0.767
IGF03 0.559
IGI01 0.586
IGI02 0.721
IGI03 0.808 0.121
OGF01 0.629 0.188
OGF02 1.023 -0.149
OGF03 0.907
OGI01 0.763 -0.106 0.141
OGI02 0.485 0.607
OGI03 0.542 0.532
Factor1 Factor2 Factor3
SS loadings 3.40 3.05 0.771
Proportion Var 0.28 0.25 0.064
Cumulative Var 0.28 0.54 0.602
Factor Correlations:
Factor1 Factor2 Factor3
Factor1 1.000 -0.124 0.390
Factor2 -0.124 1.000 -0.168
Factor3 0.390 -0.168 1.000
Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 105 on 33 degrees of freedom.
The p-value is 0.0000000022 # Four factor model:
fit06 <- factanal(bonds, 4, rotation="promax")
fit06
Call:
factanal(x = bonds, factors = 4, rotation = "promax")
Uniquenesses:
IGF01 IGF02 IGF03 IGI01 IGI02 IGI03 OGF01 OGF02 OGF03 OGI01 OGI02 OGI03
0.188 0.318 0.693 0.579 0.078 0.393 0.483 0.026 0.217 0.266 0.192 0.200
Loadings:
Factor1 Factor2 Factor3 Factor4
IGF01 1.025 -0.186
IGF02 0.882
IGF03 0.341 0.241
IGI01 0.571 -0.107
IGI02 -0.158 1.079
IGI03 0.474 0.372 0.116
OGF01 0.646 0.162
OGF02 1.041 -0.202
OGF03 0.904
OGI01 0.773 -0.167 0.123
OGI02 0.535 0.108 0.571
OGI03 0.588 0.501
Factor1 Factor2 Factor3 Factor4
SS loadings 3.6 2.20 1.77 0.704
Proportion Var 0.3 0.18 0.15 0.059
Cumulative Var 0.3 0.48 0.63 0.687
Factor Correlations:
Factor1 Factor2 Factor3 Factor4
Factor1 1.000 -0.127 0.106 0.365
Factor2 -0.127 1.000 -0.727 -0.177
Factor3 0.106 -0.727 1.000 0.193
Factor4 0.365 -0.177 0.193 1.000
Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 37 on 24 degrees of freedom.
The p-value is 0.047 FA basically suggests that ingroup / outgroup fusion/identification are not distinct factors. Only real split is between ingroup (identification + fusion) as one factor and outgroup (identification + fusion) as another factor. However, I’ll replicate the analysis as before - two constructs for ingroup fusion / identification and a single “bonds” construct for outgroup fusion+identification. Here is the reliability and inter-correlations of the three sub-scales:
igfusion <- cbind.data.frame(ds$IGF01, ds$IGF02, ds$IGF03)
names(igfusion) <- sub('ds\\$', '', names(igfusion))
igfusion <- na.omit(igfusion)
mtx2 <- cor(igfusion[, c(1:3)])
corrplot(mtx2, method = "number", number.cex = 0.7,
col=c("darkred", "red",
"darkgrey", "blue", "darkblue"))
## Reliability:
alpha02 <- psych::alpha(igfusion)
alpha02
Reliability analysis
Call: psych::alpha(x = igfusion)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.72 0.78 0.73 0.54 3.5 0.034 5.8 1 0.46
95% confidence boundaries
lower alpha upper
Feldt 0.64 0.72 0.78
Duhachek 0.65 0.72 0.79
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
IGF01 0.62 0.63 0.46 0.46 1.7 0.054 NA 0.46
IGF02 0.50 0.58 0.41 0.41 1.4 0.057 NA 0.41
IGF03 0.81 0.85 0.74 0.74 5.7 0.023 NA 0.74
Item statistics
n raw.r std.r r.cor r.drop mean sd
IGF01 189 0.79 0.86 0.80 0.65 6.4 0.83
IGF02 189 0.85 0.88 0.83 0.65 6.0 1.28
IGF03 189 0.82 0.75 0.51 0.47 5.0 1.61
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
IGF01 0.01 0.00 0.01 0.02 0.06 0.33 0.58 0
IGF02 0.02 0.02 0.01 0.07 0.14 0.29 0.46 0
IGF03 0.05 0.05 0.04 0.19 0.22 0.27 0.18 0igidentification <- cbind.data.frame(ds$IGI01, ds$IGI02, ds$IGI03)
names(igidentification) <- sub('ds\\$', '', names(igidentification))
igidentification <- na.omit(igidentification)
mtx3 <- cor(igidentification[, c(1:3)])
corrplot(mtx3, method = "number", number.cex = 0.7,
col=c("darkred", "red",
"darkgrey", "blue", "darkblue"))
## Reliability:
alpha03 <- psych::alpha(igidentification)
alpha03
Reliability analysis
Call: psych::alpha(x = igidentification)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.8 0.81 0.74 0.58 4.1 0.025 6 0.98 0.58
95% confidence boundaries
lower alpha upper
Feldt 0.74 0.8 0.84
Duhachek 0.75 0.8 0.84
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
IGI01 0.77 0.77 0.62 0.62 3.3 0.033 NA 0.62
IGI02 0.68 0.70 0.54 0.54 2.3 0.043 NA 0.54
IGI03 0.73 0.74 0.58 0.58 2.8 0.038 NA 0.58
Item statistics
n raw.r std.r r.cor r.drop mean sd
IGI01 195 0.86 0.83 0.69 0.62 5.8 1.4
IGI02 195 0.85 0.87 0.77 0.68 6.0 1.1
IGI03 195 0.83 0.85 0.73 0.64 6.1 1.0
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
IGI01 0.03 0.03 0.02 0.05 0.11 0.51 0.26 0
IGI02 0.01 0.02 0.01 0.03 0.09 0.52 0.33 0
IGI03 0.02 0.01 0.01 0.03 0.07 0.50 0.37 0ogbonds <- cbind.data.frame(ds$OGF01, ds$OGF02, ds$OGF03,
ds$OGI01, ds$OGI02, ds$OGI03)
names(ogbonds) <- sub('ds\\$', '', names(ogbonds))
ogbonds <- na.omit(ogbonds)
mtx3 <- cor(ogbonds[, c(1:6)])
corrplot(mtx3, method = "number", number.cex = 0.7,
col=c("darkred", "red",
"darkgrey", "blue", "darkblue"))
## Reliability:
alpha04 <- psych::alpha(ogbonds)
alpha04
Reliability analysis
Call: psych::alpha(x = ogbonds)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.93 0.93 0.93 0.69 13 0.0086 2.9 1.5 0.66
95% confidence boundaries
lower alpha upper
Feldt 0.91 0.93 0.94
Duhachek 0.91 0.93 0.94
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
OGF01 0.93 0.93 0.93 0.73 13.2 0.0085 0.0067 0.70
OGF02 0.91 0.91 0.90 0.67 9.9 0.0112 0.0066 0.65
OGF03 0.91 0.91 0.91 0.67 10.3 0.0109 0.0072 0.67
OGI01 0.91 0.91 0.92 0.68 10.5 0.0107 0.0082 0.64
OGI02 0.92 0.92 0.91 0.69 11.3 0.0101 0.0092 0.67
OGI03 0.91 0.91 0.91 0.68 10.7 0.0106 0.0099 0.64
Item statistics
n raw.r std.r r.cor r.drop mean sd
OGF01 184 0.79 0.78 0.71 0.69 3.3 2.0
OGF02 184 0.90 0.90 0.89 0.85 2.7 1.8
OGF03 184 0.88 0.88 0.87 0.83 2.5 1.7
OGI01 184 0.87 0.88 0.85 0.81 2.8 1.7
OGI02 184 0.84 0.84 0.81 0.77 3.1 1.8
OGI03 184 0.87 0.87 0.84 0.81 2.9 1.8
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
OGF01 0.21 0.28 0.07 0.12 0.11 0.11 0.09 0
OGF02 0.31 0.35 0.03 0.12 0.07 0.10 0.02 0
OGF03 0.35 0.33 0.04 0.12 0.08 0.07 0.02 0
OGI01 0.25 0.37 0.07 0.08 0.13 0.09 0.02 0
OGI02 0.23 0.27 0.07 0.10 0.24 0.08 0.01 0
OGI03 0.26 0.32 0.05 0.09 0.17 0.08 0.02 0## BCL and BBL items:
# BCL_01:
# Seek out opportunities to bridge social divisions with their opponents, enemies, opposition groups, or other outgroups.
# Variables: Q6.1, Q7.1
# BCL_02:
# Demonstrate willingness to compromise with their opponents, enemies, opposition groups, or other outgroups.
# Variables: Q6.2, Q7.2
# BCL_03:
# Try to understand and empathize with their opponents, enemies, opposition groups, or other outgroups.
# Variables: Q6.3, Q7.3
# BBL_01:
# Represent the interests of the communities and groups that they belong to even at the cost of other groups.
# Variables: Q6.4, Q7.4
# BBL_02:
# Focus on building stronger connections within the communities and groups they belong to rather than building stronger relationships with other groups across boundaries.
# Variables: Q6.5, Q7.5
# BBL_03:
# Try to gain benefits for the communities and groups they belong to even at the expense of other groups.
# Variables: Q6.6, Q7.6
ds$ENDBCL01 <- as.numeric(ds$`Q6.1 Leadership Quality`)
ds$ENDBCL02 <- as.numeric(ds$`Q6.2 Leadership Quality`)
ds$ENDBCL03 <- as.numeric(ds$`Q6.3 Leadership Quality`)
ds$ENDBBL01 <- as.numeric(ds$`Q6.4 Leadership Quality`)
ds$ENDBBL02 <- as.numeric(ds$`Q6.5 Leadership Quality`)
ds$ENDBBL03 <- as.numeric(ds$`Q6.6 Leadership Quality`)
ds$EXPBCL01 <- as.numeric(ds$`Q7.1 Leadership Experience`)
ds$EXPBCL02 <- as.numeric(ds$`Q7.2 Leadership Experience`)
ds$EXPBCL03 <- as.numeric(ds$`Q7.3 Leadership Experience`)
ds$EXPBBL01 <- as.numeric(ds$`Q7.4 Leadership Experience`)
ds$EXPBBL02 <- as.numeric(ds$`Q7.5 Leadership Experience`)
ds$EXPBBL03 <- as.numeric(ds$`Q7.6 Leadership Experience`)
leadership <- cbind.data.frame(ds$ENDBCL01, ds$ENDBCL02, ds$ENDBCL03,
ds$ENDBBL01, ds$ENDBBL02, ds$ENDBBL03,
ds$EXPBCL01, ds$EXPBCL02, ds$EXPBCL03,
ds$EXPBBL01, ds$EXPBBL02, ds$EXPBBL03)
names(leadership) <- sub('ds\\$', '', names(leadership))
leadership <- na.omit(leadership)
mtx1 <- cor(leadership[, c(1:12)])Shorthand for item names:
END = Endorse, EXP = Experience
BCL = BCL, BBL = BBL
So for example:
ENDBCL01 = “Endorsement of BCL, item 1”
EXPBBL02 = “Experience of BBL, item 2”, and so on
corrplot(mtx1, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
The correlation plot above shows the correlation coefficients between the twelve different items (six for endorsement of BCL/BBL, and six for experience of BCL/BBL).
## Kaiser-Meyer-Olkin (KMO) test of factorability
KMO(r=cor(leadership))Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = cor(leadership))
Overall MSA = 0.73
MSA for each item =
ENDBCL01 ENDBCL02 ENDBCL03 ENDBBL01 ENDBBL02 ENDBBL03 EXPBCL01 EXPBCL02
0.63 0.63 0.63 0.78 0.82 0.77 0.57 0.68
EXPBCL03 EXPBBL01 EXPBBL02 EXPBBL03
0.72 0.85 0.79 0.79 ## Bartlett's test of sphericity
cortest.bartlett(leadership)$chisq
[1] 880
$p.value
[1] 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000012
$df
[1] 66parallel <- fa.parallel(leadership)
Parallel analysis suggests that the number of factors = 4 and the number of components = 3 fit03 <- factanal(leadership, 3, rotation="promax")
fit03
Call:
factanal(x = leadership, factors = 3, rotation = "promax")
Uniquenesses:
ENDBCL01 ENDBCL02 ENDBCL03 ENDBBL01 ENDBBL02 ENDBBL03 EXPBCL01 EXPBCL02
0.670 0.104 0.512 0.622 0.659 0.523 0.005 0.613
EXPBCL03 EXPBBL01 EXPBBL02 EXPBBL03
0.817 0.488 0.428 0.267
Loadings:
Factor1 Factor2 Factor3
ENDBCL01 0.307 0.393
ENDBCL02 -0.174 1.003
ENDBCL03 0.742 -0.160
ENDBBL01 0.585 0.176 -0.194
ENDBBL02 0.575
ENDBBL03 0.708
EXPBCL01 -0.141 1.049
EXPBCL02 0.137 0.251 0.413
EXPBCL03 0.227 0.279
EXPBBL01 0.725
EXPBBL02 0.760 -0.194 0.111
EXPBBL03 0.877
Factor1 Factor2 Factor3
SS loadings 3.15 1.89 1.52
Proportion Var 0.26 0.16 0.13
Cumulative Var 0.26 0.42 0.55
Factor Correlations:
Factor1 Factor2 Factor3
Factor1 1.000 0.407 0.280
Factor2 0.407 1.000 0.288
Factor3 0.280 0.288 1.000
Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 143 on 33 degrees of freedom.
The p-value is 0.0000000000000012 fit04 <- factanal(leadership, 4, rotation="promax")
fit04
Call:
factanal(x = leadership, factors = 4, rotation = "promax")
Uniquenesses:
ENDBCL01 ENDBCL02 ENDBCL03 ENDBBL01 ENDBBL02 ENDBBL03 EXPBCL01 EXPBCL02
0.638 0.083 0.522 0.603 0.613 0.510 0.005 0.404
EXPBCL03 EXPBBL01 EXPBBL02 EXPBBL03
0.302 0.480 0.431 0.279
Loadings:
Factor1 Factor2 Factor3 Factor4
ENDBCL01 0.370 0.419 -0.164
ENDBCL02 -0.133 0.974
ENDBCL03 0.697 -0.139
ENDBBL01 0.603 0.195 -0.168
ENDBBL02 0.620 0.109 -0.159
ENDBBL03 0.709
EXPBCL01 1.015
EXPBCL02 0.101 0.323 0.547
EXPBCL03 0.828
EXPBBL01 0.676 0.170
EXPBBL02 0.743 -0.164 0.114
EXPBBL03 0.848
Factor1 Factor2 Factor3 Factor4
SS loadings 3.00 1.66 1.40 1.07
Proportion Var 0.25 0.14 0.12 0.09
Cumulative Var 0.25 0.39 0.51 0.59
Factor Correlations:
Factor1 Factor2 Factor3 Factor4
Factor1 1.000 0.324 -0.216 -0.293
Factor2 0.324 1.000 -0.204 -0.349
Factor3 -0.216 -0.204 1.000 0.276
Factor4 -0.293 -0.349 0.276 1.000
Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 81 on 24 degrees of freedom.
The p-value is 0.000000036 Again, FA structure does not really support a four factor Endorse/Experience BCL/BBL structure. I will still split the sub-scales into those four constructs.
end_bcl <- cbind.data.frame(ds$ENDBCL01, ds$ENDBCL02, ds$ENDBCL03)
end_bbl <- cbind.data.frame(ds$ENDBBL01, ds$ENDBBL02, ds$ENDBBL03)
exp_bcl <- cbind.data.frame(ds$EXPBCL01, ds$EXPBCL02, ds$EXPBCL03)
exp_bbl <- cbind.data.frame(ds$EXPBBL01, ds$EXPBBL02, ds$EXPBBL03)
names(end_bcl) <- sub('ds\\$', '', names(end_bcl))
names(end_bbl) <- sub('ds\\$', '', names(end_bbl))
names(exp_bcl) <- sub('ds\\$', '', names(exp_bcl))
names(exp_bbl) <- sub('ds\\$', '', names(exp_bbl))
end_bcl <- na.omit(end_bcl)
end_bbl <- na.omit(end_bbl)
exp_bcl <- na.omit(exp_bcl)
exp_bbl <- na.omit(exp_bbl)mtx1 <- cor(end_bcl[, c(1:3)])
corrplot(mtx1, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
alph01 <- psych::alpha(end_bcl)
alph01
Reliability analysis
Call: psych::alpha(x = end_bcl)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.68 0.7 0.65 0.43 2.3 0.04 5.9 1.1 0.38
95% confidence boundaries
lower alpha upper
Feldt 0.6 0.68 0.75
Duhachek 0.6 0.68 0.76
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
ENDBCL01 0.78 0.79 0.65 0.65 3.66 0.031 NA 0.65
ENDBCL02 0.41 0.42 0.27 0.27 0.74 0.081 NA 0.27
ENDBCL03 0.55 0.55 0.38 0.38 1.25 0.063 NA 0.38
Item statistics
n raw.r std.r r.cor r.drop mean sd
ENDBCL01 196 0.74 0.70 0.42 0.36 5.9 1.6
ENDBCL02 196 0.85 0.86 0.78 0.62 5.8 1.4
ENDBCL03 196 0.78 0.81 0.70 0.54 5.8 1.2
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
ENDBCL01 0.05 0.03 0.01 0.03 0.08 0.33 0.48 0
ENDBCL02 0.02 0.05 0.03 0.03 0.14 0.39 0.34 0
ENDBCL03 0.02 0.02 0.01 0.04 0.20 0.39 0.32 0mtx2 <- cor(end_bbl[, c(1:3)])
corrplot(mtx2, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
alph02 <- psych::alpha(end_bbl)
alph02
Reliability analysis
Call: psych::alpha(x = end_bbl)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.77 0.77 0.7 0.53 3.4 0.029 4.6 1.6 0.52
95% confidence boundaries
lower alpha upper
Feldt 0.71 0.77 0.82
Duhachek 0.71 0.77 0.83
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
ENDBBL01 0.68 0.68 0.52 0.52 2.2 0.045 NA 0.52
ENDBBL02 0.75 0.75 0.61 0.61 3.1 0.035 NA 0.61
ENDBBL03 0.63 0.63 0.46 0.46 1.7 0.053 NA 0.46
Item statistics
n raw.r std.r r.cor r.drop mean sd
ENDBBL01 193 0.83 0.83 0.70 0.61 4.7 2.0
ENDBBL02 193 0.79 0.80 0.62 0.55 4.8 1.9
ENDBBL03 193 0.86 0.86 0.75 0.66 4.2 2.0
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
ENDBBL01 0.10 0.10 0.05 0.11 0.22 0.21 0.20 0
ENDBBL02 0.09 0.10 0.05 0.09 0.22 0.26 0.18 0
ENDBBL03 0.13 0.12 0.08 0.16 0.18 0.20 0.12 0mtx3 <- cor(exp_bcl[, c(1:3)])
corrplot(mtx3, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
alph03 <- psych::alpha(exp_bcl)
alph03
Reliability analysis
Call: psych::alpha(x = exp_bcl)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.68 0.68 0.65 0.42 2.2 0.04 5.4 1.2 0.51
95% confidence boundaries
lower alpha upper
Feldt 0.59 0.68 0.75
Duhachek 0.60 0.68 0.76
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
EXPBCL01 0.71 0.71 0.55 0.55 2.49 0.041 NA 0.55
EXPBCL02 0.32 0.32 0.19 0.19 0.47 0.097 NA 0.19
EXPBCL03 0.67 0.67 0.51 0.51 2.06 0.047 NA 0.51
Item statistics
n raw.r std.r r.cor r.drop mean sd
EXPBCL01 195 0.74 0.72 0.51 0.40 5.5 1.6
EXPBCL02 195 0.88 0.88 0.81 0.69 5.3 1.6
EXPBCL03 195 0.73 0.74 0.56 0.42 5.3 1.5
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
EXPBCL01 0.05 0.03 0.04 0.11 0.19 0.23 0.36 0
EXPBCL02 0.03 0.03 0.06 0.15 0.18 0.29 0.26 0
EXPBCL03 0.03 0.03 0.05 0.19 0.19 0.23 0.28 0mtx4 <- cor(exp_bbl[, c(1:3)])
corrplot(mtx4, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
alph04 <- psych::alpha(exp_bbl)
alph04
Reliability analysis
Call: psych::alpha(x = exp_bbl)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.83 0.83 0.77 0.62 4.9 0.021 4.3 1.6 0.63
95% confidence boundaries
lower alpha upper
Feldt 0.79 0.83 0.87
Duhachek 0.79 0.83 0.87
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
EXPBBL01 0.81 0.81 0.69 0.69 4.4 0.027 NA 0.69
EXPBBL02 0.78 0.78 0.63 0.63 3.5 0.032 NA 0.63
EXPBBL03 0.70 0.70 0.54 0.54 2.4 0.042 NA 0.54
Item statistics
n raw.r std.r r.cor r.drop mean sd
EXPBBL01 196 0.84 0.84 0.70 0.64 4.5 1.9
EXPBBL02 196 0.85 0.86 0.75 0.68 4.5 1.8
EXPBBL03 196 0.90 0.89 0.83 0.75 3.9 2.0
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
EXPBBL01 0.09 0.07 0.14 0.21 0.12 0.16 0.20 0
EXPBBL02 0.09 0.08 0.11 0.18 0.21 0.19 0.15 0
EXPBBL03 0.18 0.09 0.13 0.22 0.12 0.15 0.11 0The four subscales have varying reliability scores.
## Proper (usable) variables in the model:
ds$Age1 <- as.numeric(gsub("\\D", "", ds$Age))
ds$Age <- as.numeric(ifelse(ds$Age1 == "6465", 64.5,
ifelse(ds$Age1 == "806", 80.5,
ifelse(ds$Age1 < 18, NA, ds$Age1))))
ds$Female <- ifelse(ds$gender == "Female", 1, 0)
ds$Married <- ifelse(ds$married == "Married", 1, 0)
ds$`SES-` <- ds$ses
ds$Endorse_BCL <- as.numeric((ds$ENDBCL01+ds$ENDBCL02+ds$ENDBCL03)/3)
ds$Endorse_BBL <- as.numeric((ds$ENDBBL01+ds$ENDBBL02+ds$ENDBBL03)/3)
ds$Experience_BCL <- as.numeric((ds$EXPBCL01+ds$EXPBCL02+ds$EXPBCL03)/3)
ds$Experience_BBL <- as.numeric((ds$EXPBBL01+ds$EXPBBL02+ds$EXPBBL03)/3)
ds$IG_Fusion <- as.numeric((ds$IGF01+ds$IGF02+ds$IGF03)/3)
ds$IG_Identification <- as.numeric((ds$IGI01+ds$IGI02+ds$IGI03)/3)
ds$OG_Bonds <- as.numeric((ds$OGF01+ds$OGF02+ds$OGF03+
ds$OGI01+ds$OGI02+ds$OGI03)/6)## Four regression models predicting endorsement and experience of BCL / BBL:
lm01 <- lm(Endorse_BCL~IG_Fusion+IG_Identification+OG_Bonds+Age+Female+Married+`SES-`,
data = ds)
lm02 <- lm(Experience_BCL~IG_Fusion+IG_Identification+OG_Bonds+Age+Female+Married+`SES-`,
data = ds)
lm03 <- lm(Endorse_BBL~IG_Fusion+IG_Identification+OG_Bonds+Age+Female+Married+`SES-`,
data = ds)
lm04 <- lm(Experience_BBL~IG_Fusion+IG_Identification+OG_Bonds+Age+Female+Married+`SES-`,
data = ds)## Tabulated results:
stargazer(lm01, lm02,
lm03, lm04,
type = "html",
star.cutoffs = c(0.05, 0.01, 0.001),
out = "table1.html")htmltools::includeHTML("table1.html")| Dependent variable: | ||||
| Endorse_BCL | Experience_BCL | Endorse_BBL | Experience_BBL | |
| (1) | (2) | (3) | (4) | |
| IG_Fusion | 0.200 | 0.410*** | 0.036 | 0.250 |
| (0.100) | (0.110) | (0.160) | (0.160) | |
| IG_Identification | 0.170 | 0.160 | 0.230 | 0.100 |
| (0.110) | (0.120) | (0.170) | (0.170) | |
| OG_Bonds | 0.110* | 0.066 | 0.180* | 0.140 |
| (0.056) | (0.060) | (0.090) | (0.086) | |
| Age | 0.005 | 0.005 | -0.006 | -0.008 |
| (0.007) | (0.008) | (0.011) | (0.011) | |
| Female | 0.043 | -0.130 | 0.100 | -0.001 |
| (0.190) | (0.200) | (0.290) | (0.290) | |
| Married | 0.120 | 0.170 | 0.160 | 0.350 |
| (0.200) | (0.210) | (0.310) | (0.310) | |
| `SES-`Middle | -0.055 | -0.650** | -0.750* | -0.640 |
| (0.230) | (0.240) | (0.360) | (0.350) | |
| `SES-`Upper middle | -0.220 | -0.370 | -0.490 | 0.0001 |
| (0.280) | (0.290) | (0.430) | (0.420) | |
| `SES-`Upper | 0.290 | -0.960 | -1.800* | -0.360 |
| (0.490) | (0.520) | (0.750) | (0.750) | |
| Constant | 3.200*** | 2.100** | 3.300** | 2.400* |
| (0.730) | (0.780) | (1.200) | (1.100) | |
| Observations | 166 | 164 | 164 | 166 |
| R2 | 0.130 | 0.250 | 0.086 | 0.097 |
| Adjusted R2 | 0.076 | 0.200 | 0.033 | 0.045 |
| Residual Std. Error | 1.000 (df = 156) | 1.100 (df = 154) | 1.600 (df = 154) | 1.600 (df = 156) |
| F Statistic | 2.500* (df = 9; 156) | 5.700*** (df = 9; 154) | 1.600 (df = 9; 154) | 1.900 (df = 9; 156) |
| Note: | *p<0.05; **p<0.01; ***p<0.001 | |||
## Empathic concern
ds$empathic_concern_01 <- as.numeric(ds$`Q18.1 Empathy`)
ds$empathic_concern_02a <- as.numeric(ds$`Q18.2 Empathy`)
ds$empathic_concern_02 <- as.numeric(8 - ds$empathic_concern_02a)
ds$empathic_concern_03 <- as.numeric(ds$`Q18.3 Empathy`)
ds$empathic_concern <- (ds$empathic_concern_01+ds$empathic_concern_02+
ds$empathic_concern_03)/3
summary(ds$empathic_concern) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
3 5 6 6 6 7 36 ds %>% drop_na(empathic_concern)%>%
ggplot(aes(x = empathic_concern))+
geom_histogram(color = "black",
fill = "gray",
bins = 30)+
geom_textvline(label = "Mean = 5.80",
xintercept = 5.80,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Empathic concern score",
y = "Frequency",
title = "Empathic concern: Gambia")+
theme_bw()
ds$perspective_taking_01 <- as.numeric(ds$`Q18.4 Empathy`)
ds$perspective_taking_02 <- as.numeric(ds$`Q18.5 Empathy`)
ds$perspective_taking_03 <- as.numeric(ds$`Q18.6 Empathy`)
ds$perspective_taking_04 <- as.numeric(ds$`Q18.7 Empathy`)
ds$perspective_taking <- (ds$perspective_taking_01+ds$perspective_taking_02+
ds$perspective_taking_03+ds$perspective_taking_04)/4
summary(ds$perspective_taking) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 5 6 6 6 7 40 ds %>% drop_na(perspective_taking)%>%
ggplot(aes(x = perspective_taking))+
geom_histogram(color = "black",
fill = "gray",
bins = 30)+
geom_textvline(label = "Mean = 5.80",
xintercept = 5.80,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Perspective taking score",
y = "Frequency",
title = "Perspective taking: Gambia")+
theme_bw()
ds$og_coop_01 <- as.numeric(ds$`Q12.1 Outgroup`)
ds$og_coop_02 <- as.numeric(ds$`Q12.2 Outgroup`)
ds$og_cooperation <- (ds$og_coop_01+ds$og_coop_02)/2
summary(ds$og_cooperation) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 3 4 4 4 7 36 ds %>% drop_na(og_cooperation)%>%
ggplot(aes(x = og_cooperation))+
geom_histogram(color = "black",
fill = "gray",
bins = 30)+
geom_textvline(label = "Mean = 3.60",
xintercept = 3.60,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Outgroup cooperation score",
y = "Frequency",
title = "Outgroup cooperation: Gambia")+
theme_bw()
ds$history_discrimination <- as.numeric(ds$`Q12.3 Outgroup`)
summary(ds$history_discrimination) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 4 5 5 6 7 36 ds %>% drop_na(history_discrimination)%>%
ggplot(aes(x = history_discrimination))+
geom_histogram(color = "black",
fill = "gray",
bins = 20)+
geom_textvline(label = "Mean = 5.00",
xintercept = 5.00,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Perceived history of discrimination score",
y = "Frequency",
title = "Perceived history of discrimination: Gambia")+
theme_bw()
ds$og_host_01 <- as.numeric(ds$`Q12.4 Outgroup`)
ds$og_host_02 <- as.numeric(ds$`Q12.5 Outgroup`)
ds$og_hostility <- (ds$og_host_01+ds$og_host_02)/2
summary(ds$og_hostility) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 2 3 3 4 7 36 ds %>% drop_na(og_hostility)%>%
ggplot(aes(x = og_hostility))+
geom_histogram(color = "black",
fill = "gray",
bins = 30)+
geom_textvline(label = "Mean = 3.30",
xintercept = 3.30,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Outgroup hostility score",
y = "Frequency",
title = "Outgroup hostility: Gambia")+
theme_bw()
ds$fight_outgroup <- as.numeric(ds$`Q12.6 Outgroup`)
summary(ds$fight_outgroup) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 3 4 4 5 7 36 ds %>% drop_na(fight_outgroup)%>%
ggplot(aes(x = fight_outgroup))+
geom_histogram(color = "black",
fill = "gray",
bins = 30)+
geom_textvline(label = "Mean = 3.80",
xintercept = 3.80,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Willingness to fight outgroup score",
y = "Frequency",
title = "Willingness to fight outgroup: Gambia")+
theme_bw()
ds3 <- cbind.data.frame(ds$empathic_concern, ds$perspective_taking,
ds$og_cooperation, ds$history_discrimination,
ds$og_hostility, ds$fight_outgroup)
ds3 <- na.omit(ds3)
names(ds3) <- sub('ds\\$', '', names(ds3))
mtx <- cor(ds3[, c(1:6)])
corrplot(mtx, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
This set of regression models have additional predictor variables compared to the last set. It adds empathic concern, perspective taking, and perceived history of discrimination as additional predictors.
## Four regression models predicting endorsement and experience of BCL / BBL:
lm01 <- lm(Endorse_BCL~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`,
data = ds)
lm02 <- lm(Experience_BCL~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`,
data = ds)
lm03 <- lm(Endorse_BBL~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`,
data = ds)
lm04 <- lm(Experience_BBL~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`,
data = ds)## Tabulated results:
stargazer(lm01, lm02,
lm03, lm04,
type = "html",
star.cutoffs = c(0.05, 0.01, 0.001),
out = "table1.html")htmltools::includeHTML("table1.html")| Dependent variable: | ||||
| Endorse_BCL | Experience_BCL | Endorse_BBL | Experience_BBL | |
| (1) | (2) | (3) | (4) | |
| IG_Fusion | 0.160 | 0.360** | 0.011 | 0.210 |
| (0.097) | (0.110) | (0.160) | (0.160) | |
| IG_Identification | 0.200 | 0.190 | 0.280 | 0.100 |
| (0.110) | (0.120) | (0.180) | (0.180) | |
| OG_Bonds | 0.140* | 0.077 | 0.250* | 0.150 |
| (0.061) | (0.070) | (0.110) | (0.100) | |
| empathic_concern | 0.010 | -0.073 | -0.160 | -0.100 |
| (0.100) | (0.120) | (0.170) | (0.170) | |
| perspective_taking | 0.220* | 0.270** | 0.190 | 0.210 |
| (0.089) | (0.100) | (0.150) | (0.150) | |
| history_discrimination | 0.021 | -0.088 | -0.210* | -0.170 |
| (0.055) | (0.063) | (0.097) | (0.092) | |
| Age | 0.001 | 0.003 | -0.006 | -0.007 |
| (0.007) | (0.008) | (0.012) | (0.011) | |
| Female | 0.094 | -0.067 | 0.120 | 0.045 |
| (0.180) | (0.200) | (0.300) | (0.300) | |
| Married | 0.092 | 0.130 | 0.099 | 0.270 |
| (0.190) | (0.220) | (0.320) | (0.320) | |
| `SES-`Middle | 0.087 | -0.630* | -0.860* | -0.730* |
| (0.210) | (0.240) | (0.360) | (0.360) | |
| `SES-`Upper middle | 0.030 | -0.230 | -0.410 | -0.012 |
| (0.260) | (0.300) | (0.450) | (0.430) | |
| `SES-`Upper | 0.620 | -0.660 | -1.700* | -0.190 |
| (0.460) | (0.520) | (0.780) | (0.770) | |
| Constant | 1.800* | 1.600 | 3.800* | 2.800 |
| (0.860) | (0.990) | (1.500) | (1.400) | |
| Observations | 162 | 160 | 160 | 162 |
| R2 | 0.220 | 0.290 | 0.120 | 0.120 |
| Adjusted R2 | 0.150 | 0.230 | 0.050 | 0.048 |
| Residual Std. Error | 0.950 (df = 149) | 1.100 (df = 147) | 1.600 (df = 147) | 1.600 (df = 149) |
| F Statistic | 3.500*** (df = 12; 149) | 5.000*** (df = 12; 147) | 1.700 (df = 12; 147) | 1.700 (df = 12; 149) |
| Note: | *p<0.05; **p<0.01; ***p<0.001 | |||
This set of regression models have different outcome variables which are: outgroup cooperation, outgroup hostility, and willingness to fight outgroup.
lm01 <- lm(og_cooperation~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`,
data = ds)
lm02 <- lm(og_hostility~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`,
data = ds)
lm03 <- lm(fight_outgroup~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`,
data = ds)## Tabulated results:
stargazer(lm01, lm02,
lm03, type = "html",
star.cutoffs = c(0.05, 0.01, 0.001),
out = "table1.html")htmltools::includeHTML("table1.html")| Dependent variable: | |||
| og_cooperation | og_hostility | fight_outgroup | |
| (1) | (2) | (3) | |
| IG_Fusion | -0.049 | -0.012 | 0.150 |
| (0.084) | (0.150) | (0.180) | |
| IG_Identification | 0.140 | 0.028 | -0.013 |
| (0.093) | (0.170) | (0.200) | |
| OG_Bonds | 0.064 | 0.031 | -0.026 |
| (0.053) | (0.093) | (0.110) | |
| empathic_concern | 0.054 | -0.130 | 0.100 |
| (0.087) | (0.150) | (0.180) | |
| perspective_taking | -0.058 | 0.190 | 0.060 |
| (0.077) | (0.140) | (0.160) | |
| history_discrimination | 0.120** | -0.140 | -0.098 |
| (0.047) | (0.084) | (0.100) | |
| Age | 0.005 | -0.004 | 0.013 |
| (0.006) | (0.010) | (0.013) | |
| Female | -0.360* | -0.200 | -0.120 |
| (0.150) | (0.270) | (0.330) | |
| Married | 0.096 | 0.640* | 0.450 |
| (0.160) | (0.290) | (0.350) | |
| `SES-`Middle | -0.250 | -0.390 | -0.370 |
| (0.180) | (0.320) | (0.390) | |
| `SES-`Upper middle | -0.120 | -0.003 | -0.240 |
| (0.230) | (0.400) | (0.470) | |
| `SES-`Upper | -0.620 | -0.550 | -2.000* |
| (0.400) | (0.700) | (0.840) | |
| Constant | 2.300** | 3.500** | 2.200 |
| (0.750) | (1.300) | (1.600) | |
| Observations | 161 | 161 | 161 |
| R2 | 0.170 | 0.086 | 0.073 |
| Adjusted R2 | 0.110 | 0.012 | -0.002 |
| Residual Std. Error (df = 148) | 0.820 | 1.400 | 1.700 |
| F Statistic (df = 12; 148) | 2.600** | 1.200 | 0.980 |
| Note: | *p<0.05; **p<0.01; ***p<0.001 | ||
ds$religious_freedom_exp_01a <- as.numeric(ds$`Q17.1 Life experience`)
ds$religious_freedom_exp_01 <- (8 - ds$religious_freedom_exp_01a)
ds$religious_freedom_exp_02a <- as.numeric(ds$`Q17.2 Life experience`)
ds$religious_freedom_exp_02 <- (8 - ds$religious_freedom_exp_02a)
ds$religious_freedom_exp_03a <- as.numeric(ds$`Q17.3 Life experience`)
ds$religious_freedom_exp_03 <- (8 - ds$religious_freedom_exp_03a)
ds$religious_freedom_exp_04a <- as.numeric(ds$`Q17.4 Life experience`)
ds$religious_freedom_exp_04 <- (8 - ds$religious_freedom_exp_04a)
ds$religious_freedom_exp_05a <- as.numeric(ds$`Q17.1 Life experience`)
ds$religious_freedom_exp_05 <- (8 - ds$religious_freedom_exp_05a)
ds$religious_freedom_exp_06a <- as.numeric(ds$`Q17.1 Life experience`)
ds$religious_freedom_exp_06 <- (8 - ds$religious_freedom_exp_06a)
ds$exp_religious_freedom <- (ds$religious_freedom_exp_01+ds$religious_freedom_exp_02+
ds$religious_freedom_exp_03+ds$religious_freedom_exp_04+
ds$religious_freedom_exp_05+ds$religious_freedom_exp_06)/6
summary(ds$exp_religious_freedom) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
2 6 7 6 7 7 38 ds %>% drop_na(exp_religious_freedom)%>%
ggplot(aes(x = exp_religious_freedom))+
geom_histogram(color = "black",
fill = "gray",
bins = 30)+
geom_textvline(label = "Mean = 6.20",
xintercept = 6.20,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Experience of religious freedom score",
y = "Frequency",
title = "Experience of religious freedom: Gambia")+
theme_bw()
ds$Religion <- ds$religion
ds2 <- ds %>%
drop_na(Religion, exp_religious_freedom)%>%
group_by(Religion) %>%
summarise(Exp_religious_freedom=mean(exp_religious_freedom),.groups = 'drop') %>%
as.data.frame()
ds2 Religion Exp_religious_freedom
1 Christian: Catholic 5.8
2 Christian: Other 5.8
3 Christian: Protestant 5.7
4 Muslim: Other 6.2
5 Muslim: Sunni 6.2ds %>% drop_na(religion, exp_religious_freedom)%>%
ggplot(aes(y = exp_religious_freedom,
x = religion))+
geom_boxplot()+
labs(x = "",
y = "Experience of religious freedom score",
title = "Experience of religious freedom: Gambia")+
coord_flip()+
theme_bw()
ds$pc01 <- ds$`Q15. Positve Contact`
## convert to sentence case:
ds$pc01 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$pc01, perl = TRUE)
ds$pc01 <- factor(ds$pc01,
levels = c("Never", "Very rarely", "Rarely",
"Sometimes", "Often", "Very often", "Always"))
table(ds$pc01)
Never Very rarely Rarely Sometimes Often Very often
1 10 20 75 38 27
Always
25 # never, very rarely, rarely
# sometimes, often, very often, always
ds$pc02 <- as.numeric(ds$pc01)
summary(ds$pc02) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 4 4 5 6 7 36 ds %>% drop_na(pc02)%>%
ggplot(aes(x = pc02))+
geom_histogram(color = "black",
fill = "gray",
bins = 30)+
geom_textvline(label = "Mean = 4.60",
xintercept = 4.60,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Positive contact score",
y = "Frequency",
title = "Positive contact with outgroup: Gambia")+
theme_bw()
ds$nc01 <- ds$`Q14 Negativ Outgroup`
## convert to sentence case:
ds$nc01 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$nc01, perl = TRUE)
ds$nc01 <- ifelse(ds$nc01 == "Ver rarely", "Very rarely", ds$nc01)
ds$nc01 <- factor(ds$nc01,
levels = c("Never", "Very rarely", "Rarely",
"Sometimes", "Often", "Very often", "Always"))
table(ds$nc01)
Never Very rarely Rarely Sometimes Often Very often
23 43 62 45 6 7
Always
10 # never, very rarely, rarely
# sometimes, often, very often, always
ds$nc02 <- as.numeric(ds$nc01)
summary(ds$nc02) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 2 3 3 4 7 36 ds %>% drop_na(pc02)%>%
ggplot(aes(x = pc02))+
geom_histogram(color = "black",
fill = "gray",
bins = 30)+
geom_textvline(label = "Mean = 3.20",
xintercept = 3.20,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Negative contact score",
y = "Frequency",
title = "Negative contact with outgroup: Gambia")+
theme_bw()
ds$intergroup_marriage_01 <- as.numeric(ds$`Q20.1 Marriage`)
ds$intergroup_marriage_02 <- as.numeric(ds$`Q20.2 Marriage`)
ds$intergroup_marriage_03 <- as.numeric(ds$`Q20.3 Marriage`)
ds$intergroup_marriage_04 <- as.numeric(ds$`Q20.4 Marriage`)
ig_marriage <- cbind.data.frame(ds$intergroup_marriage_01, ds$intergroup_marriage_02,
ds$intergroup_marriage_03, ds$intergroup_marriage_04)
names(ig_marriage) <- sub('ds\\$', '', names(ig_marriage))
ig_marriage <- na.omit(ig_marriage)
## Kaiser-Meyer-Olkin (KMO) test of factorability
KMO(r=cor(ig_marriage))Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = cor(ig_marriage))
Overall MSA = 0.82
MSA for each item =
intergroup_marriage_01 intergroup_marriage_02 intergroup_marriage_03
0.82 0.84 0.80
intergroup_marriage_04
0.82 ## Bartlett's test of sphericity
cortest.bartlett(ig_marriage)$chisq
[1] 514
$p.value
[1] 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000069
$df
[1] 6## Parallel test
parallel <- fa.parallel(ig_marriage)
Parallel analysis suggests that the number of factors = 2 and the number of components = 1 ## 1 factor solution:
fit01 <- factanal(ig_marriage, 1, rotation="promax")
fit01
Call:
factanal(x = ig_marriage, factors = 1, rotation = "promax")
Uniquenesses:
intergroup_marriage_01 intergroup_marriage_02 intergroup_marriage_03
0.29 0.35 0.23
intergroup_marriage_04
0.29
Loadings:
Factor1
intergroup_marriage_01 0.84
intergroup_marriage_02 0.81
intergroup_marriage_03 0.88
intergroup_marriage_04 0.84
Factor1
SS loadings 2.83
Proportion Var 0.71
Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 21 on 2 degrees of freedom.
The p-value is 0.000033 ## correlation between items
mtx1 <- cor(ig_marriage[, c(1:4)])
corrplot(mtx1, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
## Reliability analysis
psych::alpha(ig_marriage)
Reliability analysis
Call: psych::alpha(x = ig_marriage)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.9 0.91 0.89 0.71 9.7 0.011 5.9 1.1 0.7
95% confidence boundaries
lower alpha upper
Feldt 0.88 0.9 0.92
Duhachek 0.88 0.9 0.92
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se
intergroup_marriage_01 0.87 0.88 0.83 0.70 7.1 0.016
intergroup_marriage_02 0.88 0.89 0.85 0.72 7.9 0.015
intergroup_marriage_03 0.86 0.87 0.82 0.69 6.6 0.017
intergroup_marriage_04 0.88 0.88 0.84 0.72 7.5 0.014
var.r med.r
intergroup_marriage_01 0.0048 0.67
intergroup_marriage_02 0.0032 0.73
intergroup_marriage_03 0.0026 0.67
intergroup_marriage_04 0.0014 0.73
Item statistics
n raw.r std.r r.cor r.drop mean sd
intergroup_marriage_01 195 0.88 0.89 0.84 0.79 6.1 1.2
intergroup_marriage_02 195 0.86 0.87 0.81 0.76 6.0 1.2
intergroup_marriage_03 195 0.90 0.90 0.86 0.82 6.0 1.2
intergroup_marriage_04 195 0.89 0.88 0.82 0.78 5.5 1.5
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
intergroup_marriage_01 0.01 0.04 0.01 0.02 0.10 0.35 0.48 0
intergroup_marriage_02 0.01 0.02 0.02 0.03 0.14 0.39 0.40 0
intergroup_marriage_03 0.01 0.03 0.02 0.03 0.14 0.42 0.36 0
intergroup_marriage_04 0.04 0.03 0.02 0.11 0.18 0.31 0.30 0Based on the above output for factor analysis and reliability analysis, the “support for intergroup marriage” scale is highly reliable. Here is the visualization for the scale:
ds$endorse_intergroup_marriage <- (ds$intergroup_marriage_01+ds$intergroup_marriage_02+
ds$intergroup_marriage_03+ds$intergroup_marriage_04)/4
summary(ds$endorse_intergroup_marriage) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 6 6 6 7 7 37 ds %>% drop_na(endorse_intergroup_marriage)%>%
ggplot(aes(x = endorse_intergroup_marriage))+
geom_histogram(color = "black",
fill = "gray",
bins = 30)+
geom_textvline(label = "Mean = 5.90",
xintercept = 5.90,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Intergroup marriage endorsement score",
y = "Frequency",
title = "Support for intergroup marriage: Gambia")+
theme_bw()
## Feel outgroup: negative to positive:
ds$ogaf1 <- as.character(ds$`Q14.1 Feel Outgroup- Negative to Positive`)
## convert to sentence case:
ds$ogaf1 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf1, perl = TRUE)
ds$ogaf1 <- ifelse(ds$ogaf1 == "Very negative", "Very negative", ds$ogaf1)
ds$og_aff_01 <- factor(ds$ogaf1,
levels=c("Very negative", "Moderately negative",
"A little negative", "Neutral",
"A little positive", "Moderately positive",
"Very positive"))
## Feel Outgroup: Hostile to friendly:
ds$ogaf2 <- as.character(ds$`Q14.2 Feel Outgroup - Hostile to Friendly`)
## convert to sentence case:
ds$ogaf2 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf2, perl = TRUE)
ds$og_aff_02 <- factor(ds$ogaf2,
levels = c("Very hostile", "Moderately hostile",
"A little hostile", "Neutral",
"A little friendly", "Moderately friendly",
"Very friendly"))
table(ds$og_aff_02)
Very hostile Moderately hostile A little hostile Neutral
5 6 9 32
A little friendly Moderately friendly Very friendly
60 38 46 ## Feel Outgroup: Suspicious to trusting:
ds$ogaf3 <- as.character(ds$`Q14.3 Feel Outgroup - Suspicious to Trusting`)
## convert to sentence case:
ds$ogaf3 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf3, perl = TRUE)
# remove extra spaces:
ds$ogaf3 <- gsub("\\s+"," ",ds$ogaf3, perl = TRUE)
ds$ogaf3 <- ifelse(ds$ogaf3 == "Verytrusting", "Very trusting", ds$ogaf3)
ds$og_aff_03 <- factor(ds$ogaf3,
levels = c("Very suspicious", "Moderately suspicious",
"A little suspicious", "Neutral",
"A little trusting", "Moderately trusting",
"Very trusting"))
table(ds$og_aff_03)
Very suspicious Moderately suspicious A little suspicious
7 9 18
Neutral A little trusting Moderately trusting
51 51 25
Very trusting
33 ## Feel Outgroup: Contempt to respect:
ds$ogaf4 <- as.character(ds$`Q14.4 Feel Outgroup - Contempt to Respect`)
## convert to sentence case:
ds$ogaf4 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf4, perl = TRUE)
ds$og_aff_04 <- factor(ds$ogaf4,
levels = c("A lot of contempt", "Moderate contempt",
"A little contempt", "Neutral",
"A little respect", "Moderate respect",
"A lot of respect"))
table(ds$og_aff_04)
A lot of contempt Moderate contempt A little contempt Neutral
1 4 15 32
A little respect Moderate respect A lot of respect
53 34 56 ## Feel Outgroup: Concerned to Unconcerned:
ds$ogaf5 <- as.character(ds$`Q14.5 Feel Outgroup - Concerned to Unconcerned`)
## convert to sentence case:
ds$ogaf5 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf5, perl = TRUE)
# remove extra spaces:
ds$ogaf5 <- gsub("\\s+"," ",ds$ogaf5, perl = TRUE)
ds$ogaf5 <- ifelse(ds$ogaf5 == "Moderately unoncerned", "Moderately unconcerned", ds$ogaf5)
ds$og_aff_05 <- factor(ds$ogaf5,
levels = c("Very concerned", "Moderately concerned",
"A little concerned", "Neutral",
"A little unconcerned", "Moderately unconcerned",
"Very unconcerned"))
table(ds$og_aff_05)
Very concerned Moderately concerned A little concerned
24 17 24
Neutral A little unconcerned Moderately unconcerned
54 30 21
Very unconcerned
25 ## Feel Outgroup: Threatened to Relaxed:
ds$ogaf6 <- as.character(ds$`Q14.6 Feel Outgroup - Threatened to Relaxed`)
## convert to sentence case:
ds$ogaf6 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf6, perl = TRUE)
ds$og_aff_06 <- factor(ds$ogaf6,
levels = c("Very threatened", "Moderately threatened",
"A little threatened", "Neutral",
"A little relaxed", "Moderately relaxed",
"Very relaxed"))
table(ds$og_aff_06)
Very threatened Moderately threatened A little threatened
6 2 12
Neutral A little relaxed Moderately relaxed
60 39 23
Very relaxed
52 ds$NegativeToPositive <- as.numeric(ds$og_aff_01)
ds$HostileToFriendly <- as.numeric(ds$og_aff_02)
ds$SuspiciousToTrusting <- as.numeric(ds$og_aff_03)
ds$ContemptToRespect <- as.numeric(ds$og_aff_04)
ds$ConcernedToUnconcerned <- as.numeric(ds$og_aff_05)
ds$ThreatenedToRelaxed <- as.numeric(ds$og_aff_06)
og_affect <- cbind.data.frame(ds$NegativeToPositive, ds$HostileToFriendly,
ds$SuspiciousToTrusting, ds$ContemptToRespect,
ds$ConcernedToUnconcerned, ds$ThreatenedToRelaxed)
names(og_affect) <- sub('ds\\$', '', names(og_affect))
og_affect <- na.omit(og_affect)
mtx1 <- cor(og_affect[, c(1:6)])
corrplot(mtx1, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
## Kaiser-Meyer-Olkin (KMO) test of factorability
KMO(r=cor(og_affect))Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = cor(og_affect))
Overall MSA = 0.87
MSA for each item =
NegativeToPositive HostileToFriendly SuspiciousToTrusting
0.89 0.84 0.85
ContemptToRespect ConcernedToUnconcerned ThreatenedToRelaxed
0.89 0.87 0.88 ## Bartlett's test of sphericity
cortest.bartlett(og_affect)$chisq
[1] 545
$p.value
[1] 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000021
$df
[1] 15## Parallel test
parallel <- fa.parallel(og_affect)
Parallel analysis suggests that the number of factors = 1 and the number of components = 1 ## Reliability test
psych::alpha(og_affect)
Reliability analysis
Call: psych::alpha(x = og_affect)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.83 0.85 0.85 0.48 5.5 0.019 4.8 1.2 0.58
95% confidence boundaries
lower alpha upper
Feldt 0.8 0.83 0.87
Duhachek 0.8 0.83 0.87
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se
NegativeToPositive 0.81 0.82 0.82 0.48 4.6 0.023
HostileToFriendly 0.78 0.79 0.79 0.44 3.9 0.026
SuspiciousToTrusting 0.78 0.79 0.79 0.43 3.8 0.026
ContemptToRespect 0.80 0.81 0.81 0.46 4.3 0.024
ConcernedToUnconcerned 0.89 0.89 0.87 0.62 8.1 0.013
ThreatenedToRelaxed 0.79 0.80 0.81 0.45 4.1 0.025
var.r med.r
NegativeToPositive 0.0544 0.59
HostileToFriendly 0.0436 0.52
SuspiciousToTrusting 0.0434 0.52
ContemptToRespect 0.0495 0.56
ConcernedToUnconcerned 0.0052 0.62
ThreatenedToRelaxed 0.0572 0.54
Item statistics
n raw.r std.r r.cor r.drop mean sd
NegativeToPositive 193 0.76 0.76 0.69 0.62 4.6 1.8
HostileToFriendly 193 0.84 0.85 0.84 0.76 5.2 1.5
SuspiciousToTrusting 193 0.85 0.86 0.85 0.77 4.7 1.5
ContemptToRespect 193 0.77 0.79 0.74 0.67 5.3 1.4
ConcernedToUnconcerned 193 0.48 0.45 0.26 0.24 4.1 1.8
ThreatenedToRelaxed 193 0.81 0.81 0.77 0.71 5.1 1.5
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
NegativeToPositive 0.05 0.10 0.10 0.19 0.32 0.00 0.24 0
HostileToFriendly 0.03 0.03 0.05 0.17 0.30 0.19 0.24 0
SuspiciousToTrusting 0.04 0.05 0.09 0.26 0.26 0.13 0.17 0
ContemptToRespect 0.01 0.02 0.08 0.17 0.27 0.17 0.29 0
ConcernedToUnconcerned 0.12 0.09 0.12 0.27 0.16 0.11 0.12 0
ThreatenedToRelaxed 0.03 0.01 0.06 0.31 0.20 0.11 0.27 0Based on the output of factor analysis as well as reliability, it seems like one item (“concerned to unconcerned”) is problematic. Here is the scale without that item:
og_affect <- cbind.data.frame(ds$NegativeToPositive, ds$HostileToFriendly,
ds$SuspiciousToTrusting, ds$ContemptToRespect,
ds$ThreatenedToRelaxed)
names(og_affect) <- sub('ds\\$', '', names(og_affect))
og_affect <- na.omit(og_affect)
mtx1 <- cor(og_affect[, c(1:5)])
corrplot(mtx1, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
## Reliability analysis:
psych::alpha(og_affect)
Reliability analysis
Call: psych::alpha(x = og_affect)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.89 0.89 0.87 0.62 8.1 0.013 5 1.3 0.62
95% confidence boundaries
lower alpha upper
Feldt 0.86 0.89 0.91
Duhachek 0.86 0.89 0.91
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r
NegativeToPositive 0.88 0.88 0.86 0.65 7.6 0.014 0.0031
HostileToFriendly 0.85 0.85 0.82 0.59 5.8 0.018 0.0053
SuspiciousToTrusting 0.85 0.85 0.82 0.59 5.7 0.018 0.0041
ContemptToRespect 0.87 0.87 0.84 0.63 6.8 0.016 0.0047
ThreatenedToRelaxed 0.87 0.87 0.85 0.63 6.8 0.016 0.0070
med.r
NegativeToPositive 0.67
HostileToFriendly 0.58
SuspiciousToTrusting 0.59
ContemptToRespect 0.62
ThreatenedToRelaxed 0.65
Item statistics
n raw.r std.r r.cor r.drop mean sd
NegativeToPositive 193 0.80 0.78 0.70 0.66 4.6 1.8
HostileToFriendly 193 0.87 0.87 0.84 0.79 5.2 1.5
SuspiciousToTrusting 193 0.88 0.88 0.85 0.80 4.7 1.5
ContemptToRespect 193 0.80 0.82 0.75 0.70 5.3 1.4
ThreatenedToRelaxed 193 0.82 0.82 0.75 0.71 5.1 1.5
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
NegativeToPositive 0.05 0.10 0.10 0.19 0.32 0.00 0.24 0
HostileToFriendly 0.03 0.03 0.05 0.17 0.30 0.19 0.24 0
SuspiciousToTrusting 0.04 0.05 0.09 0.26 0.26 0.13 0.17 0
ContemptToRespect 0.01 0.02 0.08 0.17 0.27 0.17 0.29 0
ThreatenedToRelaxed 0.03 0.01 0.06 0.31 0.20 0.11 0.27 0Based on the above output for factor analysis and reliability analysis, the “outgroup affect” scale is highly reliable if we drop the problematic item (“concerned to unconcerned”). Here is the visualization for the scale after dropping the problematic item:
ds$og_affect <- (ds$NegativeToPositive+ds$HostileToFriendly+
ds$SuspiciousToTrusting+ds$ContemptToRespect+
ds$ThreatenedToRelaxed)/5
summary(ds$og_affect) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 4 5 5 6 7 39 ds %>% drop_na(og_affect)%>%
ggplot(aes(x = og_affect))+
geom_histogram(color = "black",
fill = "gray",
bins = 50)+
geom_textvline(label = "Mean = 5.00",
xintercept = 5.00,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Outgroup affect score",
y = "Frequency",
title = "Outgroup affect: Gambia")+
theme_bw()
This set of regression models have additional predictors: frequency of positive and negative contact, and interaction between perspective taking and perceived history of discrimination
ds$freq_positive_contact <- ds$pc02
ds$freq_negative_contact <- ds$nc02
ds$perspectiveXdiscrimination <- ds$perspective_taking*ds$history_discrimination
## Four regression models predicting endorsement and experience of BCL / BBL:
lm01 <- lm(Endorse_BCL~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`,
data = ds)
lm02 <- lm(Experience_BCL~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`,
data = ds)
lm03 <- lm(Endorse_BBL~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`,
data = ds)
lm04 <- lm(Experience_BBL~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`,
data = ds)## Tabulated results:
stargazer(lm01, lm02,
lm03, lm04,
type = "html",
star.cutoffs = c(0.05, 0.01, 0.001),
out = "table1.html")htmltools::includeHTML("table1.html")| Dependent variable: | ||||
| Endorse_BCL | Experience_BCL | Endorse_BBL | Experience_BBL | |
| (1) | (2) | (3) | (4) | |
| IG_Fusion | 0.170 | 0.370** | 0.0001 | 0.280 |
| (0.100) | (0.110) | (0.160) | (0.160) | |
| IG_Identification | 0.200 | 0.150 | 0.310 | 0.099 |
| (0.120) | (0.130) | (0.180) | (0.190) | |
| OG_Bonds | 0.130 | 0.039 | 0.240* | 0.100 |
| (0.067) | (0.075) | (0.110) | (0.110) | |
| freq_positive_contact | -0.019 | 0.010 | -0.250* | -0.018 |
| (0.069) | (0.077) | (0.110) | (0.110) | |
| freq_negative_contact | 0.061 | 0.160* | 0.300** | 0.210* |
| (0.058) | (0.065) | (0.092) | (0.094) | |
| empathic_concern | 0.038 | -0.027 | -0.160 | 0.016 |
| (0.110) | (0.120) | (0.170) | (0.180) | |
| perspective_taking | 0.097 | 0.300 | 0.120 | -0.160 |
| (0.230) | (0.250) | (0.370) | (0.360) | |
| history_discrimination | -0.100 | 0.002 | -0.190 | -0.540 |
| (0.260) | (0.290) | (0.410) | (0.420) | |
| perspectiveXdiscrimination | 0.024 | -0.014 | 0.012 | 0.072 |
| (0.045) | (0.050) | (0.073) | (0.073) | |
| Age | 0.001 | 0.005 | -0.011 | -0.007 |
| (0.008) | (0.008) | (0.012) | (0.012) | |
| Female | 0.066 | -0.120 | 0.019 | -0.054 |
| (0.180) | (0.200) | (0.290) | (0.290) | |
| Married | 0.069 | 0.018 | 0.049 | 0.230 |
| (0.200) | (0.230) | (0.320) | (0.320) | |
| `SES-`Middle | 0.076 | -0.670** | -0.920** | -0.780* |
| (0.220) | (0.250) | (0.350) | (0.360) | |
| `SES-`Upper middle | 0.001 | -0.320 | -0.400 | -0.150 |
| (0.270) | (0.300) | (0.440) | (0.440) | |
| `SES-`Upper | 0.600 | -0.550 | -1.100 | -0.300 |
| (0.500) | (0.560) | (0.800) | (0.810) | |
| Constant | 2.100 | 0.850 | 4.100* | 3.200 |
| (1.200) | (1.300) | (2.000) | (1.900) | |
| Observations | 160 | 158 | 158 | 160 |
| R2 | 0.230 | 0.330 | 0.230 | 0.170 |
| Adjusted R2 | 0.150 | 0.260 | 0.150 | 0.084 |
| Residual Std. Error | 0.950 (df = 144) | 1.100 (df = 142) | 1.500 (df = 142) | 1.500 (df = 144) |
| F Statistic | 2.900*** (df = 15; 144) | 4.600*** (df = 15; 142) | 2.800*** (df = 15; 142) | 2.000* (df = 15; 144) |
| Note: | *p<0.05; **p<0.01; ***p<0.001 | |||
Same predictors as section 14 above, but with different outcomes. This set of regression models has an additional outcome variable: affect towards outgroup.
## Four regression models:
lm01 <- lm(og_cooperation~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`,
data = ds)
lm02 <- lm(og_hostility~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`,
data = ds)
lm03 <- lm(fight_outgroup~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`,
data = ds)
lm04 <- lm(og_affect~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`,
data = ds)## Tabulated results:
stargazer(lm01, lm02,
lm03, lm04,
type = "html",
star.cutoffs = c(0.05, 0.01, 0.001),
out = "table1.html")htmltools::includeHTML("table1.html")| Dependent variable: | ||||
| og_cooperation | og_hostility | fight_outgroup | og_affect | |
| (1) | (2) | (3) | (4) | |
| IG_Fusion | -0.086 | -0.021 | 0.210 | 0.005 |
| (0.087) | (0.140) | (0.180) | (0.110) | |
| IG_Identification | 0.140 | 0.031 | -0.120 | 0.280* |
| (0.099) | (0.160) | (0.200) | (0.130) | |
| OG_Bonds | 0.086 | 0.017 | -0.160 | 0.140 |
| (0.058) | (0.095) | (0.120) | (0.076) | |
| freq_positive_contact | -0.059 | -0.240* | 0.080 | 0.240** |
| (0.060) | (0.098) | (0.120) | (0.080) | |
| freq_negative_contact | -0.017 | 0.300*** | 0.410*** | -0.051 |
| (0.050) | (0.082) | (0.100) | (0.065) | |
| empathic_concern | 0.003 | -0.140 | 0.220 | 0.050 |
| (0.094) | (0.150) | (0.190) | (0.120) | |
| perspective_taking | 0.270 | 0.220 | -0.210 | -0.210 |
| (0.190) | (0.320) | (0.390) | (0.250) | |
| history_discrimination | 0.540* | -0.010 | -0.330 | -0.096 |
| (0.220) | (0.370) | (0.450) | (0.290) | |
| perspectiveXdiscrimination | -0.072 | -0.008 | 0.046 | 0.050 |
| (0.039) | (0.064) | (0.079) | (0.051) | |
| Age | 0.004 | -0.013 | 0.012 | -0.001 |
| (0.006) | (0.011) | (0.013) | (0.008) | |
| Female | -0.340* | -0.290 | -0.260 | -0.023 |
| (0.160) | (0.260) | (0.320) | (0.200) | |
| Married | 0.140 | 0.640* | 0.210 | -0.028 |
| (0.170) | (0.280) | (0.350) | (0.230) | |
| `SES-`Middle | -0.330 | -0.430 | -0.250 | 0.015 |
| (0.190) | (0.310) | (0.390) | (0.260) | |
| `SES-`Upper middle | -0.130 | 0.032 | -0.300 | -0.380 |
| (0.230) | (0.380) | (0.470) | (0.310) | |
| `SES-`Upper | -0.410 | 0.130 | -1.700 | -0.820 |
| (0.430) | (0.710) | (0.870) | (0.570) | |
| Constant | 1.300 | 3.500* | 2.000 | 2.000 |
| (1.000) | (1.700) | (2.100) | (1.300) | |
| Observations | 159 | 159 | 159 | 156 |
| R2 | 0.200 | 0.210 | 0.170 | 0.320 |
| Adjusted R2 | 0.120 | 0.130 | 0.079 | 0.250 |
| Residual Std. Error | 0.820 (df = 143) | 1.300 (df = 143) | 1.700 (df = 143) | 1.100 (df = 140) |
| F Statistic | 2.400** (df = 15; 143) | 2.600** (df = 15; 143) | 1.900* (df = 15; 143) | 4.500*** (df = 15; 140) |
| Note: | *p<0.05; **p<0.01; ***p<0.001 | |||
ds$event_negative_affect <- as.numeric(ds$`Q9.1 Feelings about Event`)
table(ds$event_negative_affect)
1 2 3 4 5 6 7
14 18 7 7 31 51 68 ds$event_positive_affect_01 <- (8 - ds$event_negative_affect)
table(ds$event_positive_affect_01)
1 2 3 4 5 6 7
68 51 31 7 7 18 14 ds$event_positive_affect_02 <- as.numeric(ds$`Q9.2 Feelings about Event`)
table(ds$event_positive_affect_02)
1 2 3 4 5 6 7
59 42 11 10 20 38 14 ds$event_positive_affect <- ds$event_positive_affect_02
ds$event_episodic_recall_01 <- as.numeric(ds$`Q9.3 Feelings about Event`)
table(ds$event_episodic_recall_01)
1 2 3 4 5 6 7
1 3 3 22 41 79 46 ds$event_episodic_recall_02 <- as.numeric(ds$`Q9.4 Feelings about Event`)
table(ds$event_episodic_recall_02)
1 2 3 4 5 6 7
5 7 5 17 26 88 46 ds$event_shared_perception_01 <- as.numeric(ds$`Q9.5 Feelings about Event`)
table(ds$event_shared_perception_01)
1 2 3 4 5 6 7
2 5 2 19 26 77 65 ds$event_shared_perception_02 <- as.numeric(ds$`Q9.6 Feelings about Event`)
table(ds$event_shared_perception_02)
1 2 3 4 5 6 7
2 3 3 25 37 72 54 ds$event_event_reflection_01 <- as.numeric(ds$`Q9.9 Feelings about Event`)
table(ds$event_event_reflection_01)
1 2 3 4 5 6 7
8 41 23 36 28 43 15 ds$event_event_reflection_02 <- as.numeric(ds$`Q9.10 Feelings about Event`)
table(ds$event_event_reflection_02)
1 2 3 4 5 6 7
26 46 24 35 26 21 18 ds$event_transformative_indiv_01 <- as.numeric(ds$`Q9.7 Feelings about Event`)
ds$event_transformative_indiv_01 <- ifelse(ds$event_transformative_indiv_01 > 7, NA, ds$event_transformative_indiv_01)
table(ds$event_transformative_indiv_01)
1 2 3 4 5 6 7
27 17 9 17 22 57 44 ds$event_transformative_indiv_02 <- as.numeric(ds$`Q9.8 Feelings about Event`)
table(ds$event_transformative_indiv_02)
1 2 3 4 5 6 7
64 59 21 20 8 12 9 ds$event_transformative_group_01 <- as.numeric(ds$`Q9.11 Feelings about Event`)
table(ds$event_transformative_group_01)
1 2 3 4 5 6 7
11 5 3 12 32 80 53 ds$event_transformative_group_02 <- as.numeric(ds$`Q9.12 Feelings about Event`)
table(ds$event_transformative_group_02)
1 2 3 4 5 6 7
64 44 21 29 18 12 7 imagistic <- cbind.data.frame(ds$event_negative_affect, ds$event_positive_affect,
ds$event_episodic_recall_01, ds$event_episodic_recall_02,
ds$event_shared_perception_01, ds$event_shared_perception_02,
ds$event_event_reflection_01, ds$event_event_reflection_02,
ds$event_transformative_indiv_01, ds$event_transformative_indiv_02,
ds$event_transformative_group_01, ds$event_transformative_group_02)
imagistic <- na.omit(imagistic)
names(imagistic) <- sub('ds\\$event\\_', '', names(imagistic))
mtx <- cor(imagistic[, c(1:12)])
corrplot(mtx, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
## Kaiser-Meyer-Olkin (KMO) test of factorability
KMO(r=cor(imagistic))Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = cor(imagistic))
Overall MSA = 0.67
MSA for each item =
negative_affect positive_affect episodic_recall_01
0.61 0.56 0.73
episodic_recall_02 shared_perception_01 shared_perception_02
0.70 0.73 0.69
event_reflection_01 event_reflection_02 transformative_indiv_01
0.64 0.61 0.67
transformative_indiv_02 transformative_group_01 transformative_group_02
0.67 0.77 0.70 ## Bartlett's test of sphericity
cortest.bartlett(imagistic)$chisq
[1] 590
$p.value
[1] 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000003
$df
[1] 66## Parallel test
parallel <- fa.parallel(imagistic)
Parallel analysis suggests that the number of factors = 4 and the number of components = 3 Scree plot suggests there needs to be three factors, but theoretically we want six factors. Here is how the three factors break down:
fit03 <- factanal(imagistic, 3, rotation="promax")
fit03
Call:
factanal(x = imagistic, factors = 3, rotation = "promax")
Uniquenesses:
negative_affect positive_affect episodic_recall_01
0.42 0.52 0.63
episodic_recall_02 shared_perception_01 shared_perception_02
0.65 0.46 0.69
event_reflection_01 event_reflection_02 transformative_indiv_01
0.38 0.29 0.62
transformative_indiv_02 transformative_group_01 transformative_group_02
0.75 0.66 0.82
Loadings:
Factor1 Factor2 Factor3
negative_affect 0.294 0.754
positive_affect -0.190 -0.697
episodic_recall_01 0.602
episodic_recall_02 0.668 -0.203
shared_perception_01 0.711 0.390
shared_perception_02 0.533 0.274
event_reflection_01 0.815
event_reflection_02 -0.123 0.916 0.186
transformative_indiv_01 0.374 0.165 -0.305
transformative_indiv_02 0.400 -0.195
transformative_group_01 0.441 0.185 -0.101
transformative_group_02 0.368
Factor1 Factor2 Factor3
SS loadings 2.05 1.95 1.48
Proportion Var 0.17 0.16 0.12
Cumulative Var 0.17 0.33 0.46
Factor Correlations:
Factor1 Factor2 Factor3
Factor1 1.000 0.222 0.512
Factor2 0.222 1.000 0.170
Factor3 0.512 0.170 1.000
Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 102 on 33 degrees of freedom.
The p-value is 0.0000000056 The factor structure above makes sense but is a little incoherent. Here is what happens if we “force” six factors:
fit06 <- factanal(imagistic, 6, rotation="promax")
fit06
Call:
factanal(x = imagistic, factors = 6, rotation = "promax")
Uniquenesses:
negative_affect positive_affect episodic_recall_01
0.442 0.470 0.448
episodic_recall_02 shared_perception_01 shared_perception_02
0.542 0.456 0.466
event_reflection_01 event_reflection_02 transformative_indiv_01
0.005 0.370 0.613
transformative_indiv_02 transformative_group_01 transformative_group_02
0.584 0.005 0.422
Loadings:
Factor1 Factor2 Factor3 Factor4 Factor5 Factor6
negative_affect 0.106 -0.656 0.202
positive_affect 0.110 0.714
episodic_recall_01 0.797 -0.123
episodic_recall_02 -0.239 0.717
shared_perception_01 0.161 -0.198 0.532
shared_perception_02 0.135 0.721
event_reflection_01 1.105 -0.119 0.112
event_reflection_02 0.486 0.111 -0.295 0.270 -0.216
transformative_indiv_01 0.286 0.231 0.288 0.202
transformative_indiv_02 0.164 -0.112 0.575
transformative_group_01 1.035
transformative_group_02 -0.160 0.829
Factor1 Factor2 Factor3 Factor4 Factor5 Factor6
SS loadings 1.65 1.30 1.184 1.167 1.139 0.918
Proportion Var 0.14 0.11 0.099 0.097 0.095 0.076
Cumulative Var 0.14 0.24 0.344 0.441 0.536 0.613
Factor Correlations:
Factor1 Factor2 Factor3 Factor4 Factor5 Factor6
Factor1 1.0000 -0.4076 -0.0844 0.3952 -0.5587 -0.118
Factor2 -0.4076 1.0000 -0.0419 -0.5434 0.3297 0.281
Factor3 -0.0844 -0.0419 1.0000 0.0335 -0.0464 0.270
Factor4 0.3952 -0.5434 0.0335 1.0000 -0.3190 -0.488
Factor5 -0.5587 0.3297 -0.0464 -0.3190 1.0000 0.126
Factor6 -0.1183 0.2809 0.2702 -0.4879 0.1257 1.000
Test of the hypothesis that 6 factors are sufficient.
The chi square statistic is 24 on 9 degrees of freedom.
The p-value is 0.0049 FA solution still seems incoherent.
imagistic <- cbind.data.frame(ds$event_negative_affect, ds$event_positive_affect,
ds$event_episodic_recall_01, ds$event_episodic_recall_02,
ds$event_shared_perception_01, ds$event_shared_perception_02,
ds$event_event_reflection_01, ds$event_event_reflection_02,
ds$event_transformative_indiv_01, ds$event_transformative_indiv_02,
ds$event_transformative_group_01, ds$event_transformative_group_02)
## Reliability:
alph01 <- psych::alpha(imagistic)
alph01
Reliability analysis
Call: psych::alpha(x = imagistic)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.68 0.71 0.78 0.17 2.5 0.031 4.6 0.8 0.18
95% confidence boundaries
lower alpha upper
Feldt 0.62 0.68 0.74
Duhachek 0.62 0.68 0.74
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N
ds$event_negative_affect 0.71 0.73 0.77 0.19 2.6
ds$event_positive_affect 0.73 0.74 0.78 0.21 2.9
ds$event_episodic_recall_01 0.65 0.68 0.75 0.16 2.1
ds$event_episodic_recall_02 0.66 0.69 0.75 0.17 2.2
ds$event_shared_perception_01 0.66 0.69 0.75 0.17 2.2
ds$event_shared_perception_02 0.66 0.69 0.75 0.17 2.2
ds$event_event_reflection_01 0.63 0.67 0.73 0.16 2.0
ds$event_event_reflection_02 0.64 0.68 0.74 0.16 2.1
ds$event_transformative_indiv_01 0.64 0.68 0.74 0.16 2.1
ds$event_transformative_indiv_02 0.65 0.69 0.75 0.17 2.2
ds$event_transformative_group_01 0.64 0.68 0.75 0.16 2.1
ds$event_transformative_group_02 0.65 0.69 0.76 0.17 2.2
alpha se var.r med.r
ds$event_negative_affect 0.028 0.023 0.19
ds$event_positive_affect 0.026 0.019 0.20
ds$event_episodic_recall_01 0.034 0.032 0.17
ds$event_episodic_recall_02 0.033 0.032 0.18
ds$event_shared_perception_01 0.033 0.028 0.17
ds$event_shared_perception_02 0.033 0.031 0.17
ds$event_event_reflection_01 0.036 0.028 0.16
ds$event_event_reflection_02 0.035 0.029 0.17
ds$event_transformative_indiv_01 0.036 0.031 0.17
ds$event_transformative_indiv_02 0.034 0.031 0.18
ds$event_transformative_group_01 0.035 0.032 0.16
ds$event_transformative_group_02 0.034 0.032 0.18
Item statistics
n raw.r std.r r.cor r.drop mean sd
ds$event_negative_affect 196 0.24 0.27 0.172 0.036 5.3 1.9
ds$event_positive_affect 194 0.19 0.14 0.027 -0.035 3.3 2.2
ds$event_episodic_recall_01 195 0.53 0.58 0.528 0.436 5.7 1.1
ds$event_episodic_recall_02 194 0.48 0.52 0.454 0.354 5.6 1.4
ds$event_shared_perception_01 196 0.45 0.52 0.472 0.337 5.8 1.3
ds$event_shared_perception_02 196 0.46 0.52 0.466 0.347 5.7 1.2
ds$event_event_reflection_01 194 0.63 0.61 0.592 0.502 4.2 1.7
ds$event_event_reflection_02 196 0.60 0.58 0.555 0.449 3.6 1.9
ds$event_transformative_indiv_01 193 0.62 0.58 0.530 0.447 4.7 2.1
ds$event_transformative_indiv_02 193 0.54 0.51 0.448 0.402 2.6 1.7
ds$event_transformative_group_01 196 0.56 0.57 0.504 0.436 5.6 1.6
ds$event_transformative_group_02 195 0.51 0.49 0.415 0.358 2.8 1.8
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
ds$event_negative_affect 0.07 0.09 0.04 0.04 0.16 0.26 0.35 0.16
ds$event_positive_affect 0.30 0.22 0.06 0.05 0.10 0.20 0.07 0.16
ds$event_episodic_recall_01 0.01 0.02 0.02 0.11 0.21 0.41 0.24 0.16
ds$event_episodic_recall_02 0.03 0.04 0.03 0.09 0.13 0.45 0.24 0.16
ds$event_shared_perception_01 0.01 0.03 0.01 0.10 0.13 0.39 0.33 0.16
ds$event_shared_perception_02 0.01 0.02 0.02 0.13 0.19 0.37 0.28 0.16
ds$event_event_reflection_01 0.04 0.21 0.12 0.19 0.14 0.22 0.08 0.16
ds$event_event_reflection_02 0.13 0.23 0.12 0.18 0.13 0.11 0.09 0.16
ds$event_transformative_indiv_01 0.14 0.09 0.05 0.09 0.11 0.30 0.23 0.17
ds$event_transformative_indiv_02 0.33 0.31 0.11 0.10 0.04 0.06 0.05 0.17
ds$event_transformative_group_01 0.06 0.03 0.02 0.06 0.16 0.41 0.27 0.16
ds$event_transformative_group_02 0.33 0.23 0.11 0.15 0.09 0.06 0.04 0.16If we conceptualize “imagistic experience” as a single measure that is the sum of the six scales, here is the visualization of the measure:
ds$imagistic <- (((ds$event_positive_affect)+
(ds$event_negative_affect)+
((ds$event_episodic_recall_01+ds$event_episodic_recall_02)/2)+
((ds$event_shared_perception_01+ds$event_shared_perception_02)/2)+
((ds$event_event_reflection_01+ds$event_event_reflection_02)/2)+
((ds$event_transformative_indiv_01+ds$event_transformative_indiv_02)/2)+
((ds$event_transformative_group_01+ds$event_transformative_group_02)/2)))
summary(ds$imagistic) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
10 28 31 32 36 46 49 ds %>% drop_na(imagistic)%>%
ggplot(aes(x = imagistic))+
geom_histogram(color = "black",
fill = "gray",
bins = 50)+
geom_textvline(label = "Mean = 32.00",
xintercept = 32.00,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Imagistic experience",
y = "Frequency",
title = "Imagistic experience: Gambia")+
theme_bw()
Below are the individual Cronbach’s alpha for each of the seven subscales, followed by their individual visualizations, followed by the correlation between the seven sub-scales
## Positive affect:
negative_affect <- cbind.data.frame(ds$event_negative_affect)
negative_affect <- na.omit(negative_affect)
## Visualization:
ds$event_negative_affect <- (ds$event_negative_affect)
summary(ds$event_negative_affect) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 5 6 5 7 7 36 ds %>% drop_na(event_negative_affect)%>%
ggplot(aes(x = event_negative_affect))+
geom_histogram(color = "black",
fill = "gray",
bins = 50)+
geom_textvline(label = "Mean = 5.30",
xintercept = 5.30,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Negative affect about event",
y = "Frequency",
title = "Negative affect about event: Gambia")+
theme_bw()
## Positive affect:
positive_affect <- cbind.data.frame(ds$event_positive_affect)
positive_affect <- na.omit(positive_affect)
## Visualization:
ds$event_positive_affect <- (ds$event_positive_affect)
summary(ds$event_positive_affect) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 1 2 3 6 7 38 ds %>% drop_na(event_positive_affect)%>%
ggplot(aes(x = event_positive_affect))+
geom_histogram(color = "black",
fill = "gray",
bins = 50)+
geom_textvline(label = "Mean = 3.30",
xintercept = 3.30,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Positive affect about event",
y = "Frequency",
title = "Positive affect about event: Gambia")+
theme_bw()
## Episodic recall:
episodic_recall <- cbind.data.frame(ds$event_episodic_recall_01, ds$event_episodic_recall_02)
episodic_recall <- na.omit(episodic_recall)
names(episodic_recall) <- sub('ds\\$event\\_', '', names(episodic_recall))
mtx <- cor(episodic_recall[, c(1:2)])
## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
## Reliability:
alph01 <- psych::alpha(episodic_recall)
summary(alph01)
Reliability analysis
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.62 0.63 0.46 0.46 1.7 0.053 5.6 1.1 0.46## Visualization:
ds$event_episodic_recall <- ((ds$event_episodic_recall_01+ds$event_episodic_recall_02)/2)
summary(ds$event_episodic_recall) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
2 5 6 6 6 7 39 ds %>% drop_na(event_episodic_recall)%>%
ggplot(aes(x = event_episodic_recall))+
geom_histogram(color = "black",
fill = "gray",
bins = 50)+
geom_textvline(label = "Mean = 5.60",
xintercept = 5.60,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Episodic recall of event",
y = "Frequency",
title = "Episodic recall of event: Gambia")+
theme_bw()
## Shared perception:
shared_perception <- cbind.data.frame(ds$event_shared_perception_01, ds$event_shared_perception_02)
shared_perception <- na.omit(shared_perception)
names(shared_perception) <- sub('ds\\$event\\_', '', names(shared_perception))
mtx <- cor(shared_perception[, c(1:2)])
## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
## Reliability:
alph01 <- psych::alpha(shared_perception)
summary(alph01)
Reliability analysis
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.67 0.67 0.51 0.51 2.1 0.047 5.7 1.1 0.51## Visualization:
ds$event_shared_perception <- ((ds$event_shared_perception_01+ds$event_shared_perception_02)/2)
summary(ds$event_shared_perception) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
2 5 6 6 6 7 36 ds %>% drop_na(event_shared_perception)%>%
ggplot(aes(x = event_shared_perception))+
geom_histogram(color = "black",
fill = "gray",
bins = 50)+
geom_textvline(label = "Mean = 5.70",
xintercept = 5.70,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Shared perception of event",
y = "Frequency",
title = "Shared perception of event: Gambia")+
theme_bw()
## Reflection:
event_reflection <- cbind.data.frame(ds$event_event_reflection_01, ds$event_event_reflection_02)
event_reflection <- na.omit(event_reflection)
names(event_reflection) <- sub('ds\\$event\\_', '', names(event_reflection))
mtx <- cor(event_reflection[, c(1:2)])
## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
## Reliability:
alph01 <- psych::alpha(event_reflection)
summary(alph01)
Reliability analysis
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.77 0.77 0.62 0.62 3.3 0.033 3.9 1.6 0.62## Visualization:
ds$event_event_reflection <- ((ds$event_event_reflection_01+ds$event_event_reflection_02)/2)
summary(ds$event_event_reflection) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 2 4 4 5 7 38 ds %>% drop_na(event_event_reflection)%>%
ggplot(aes(x = event_event_reflection))+
geom_histogram(color = "black",
fill = "gray",
bins = 50)+
geom_textvline(label = "Mean = 3.90",
xintercept = 3.90,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Reflection of event",
y = "Frequency",
title = "Reflection of event: Gambia")+
theme_bw()
## Reflection:
event_transformative_indiv <- cbind.data.frame(ds$event_transformative_indiv_01, ds$event_transformative_indiv_02)
event_transformative_indiv <- na.omit(event_transformative_indiv)
names(event_transformative_indiv) <- sub('ds\\$event\\_', '', names(event_transformative_indiv))
mtx <- cor(event_transformative_indiv[, c(1:2)])
## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
## Reliability:
psych::alpha(event_transformative_indiv)
Reliability analysis
Call: psych::alpha(x = event_transformative_indiv)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.46 0.47 0.31 0.31 0.89 0.076 3.7 1.6 0.31
95% confidence boundaries
lower alpha upper
Feldt 0.29 0.46 0.60
Duhachek 0.31 0.46 0.61
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se
transformative_indiv_01 0.37 0.31 0.095 0.31 0.44 NA
transformative_indiv_02 0.26 0.31 0.095 0.31 0.44 NA
var.r med.r
transformative_indiv_01 0 0.31
transformative_indiv_02 0 0.31
Item statistics
n raw.r std.r r.cor r.drop mean sd
transformative_indiv_01 190 0.85 0.81 0.45 0.31 4.7 2.1
transformative_indiv_02 190 0.77 0.81 0.45 0.31 2.6 1.8
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
transformative_indiv_01 0.14 0.08 0.05 0.09 0.11 0.30 0.23 0
transformative_indiv_02 0.33 0.30 0.11 0.11 0.04 0.06 0.05 0## Visualization:
ds$event_transformative_indiv <- ((ds$event_transformative_indiv_01+ds$event_transformative_indiv_02)/2)
summary(ds$event_transformative_indiv) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 2 4 4 4 7 42 ds %>% drop_na(event_transformative_indiv)%>%
ggplot(aes(x = event_transformative_indiv))+
geom_histogram(color = "black",
fill = "gray",
bins = 50)+
geom_textvline(label = "Mean = 3.70",
xintercept = 3.70,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Transformative event score",
y = "Frequency",
title = "Transformative event for individual: Gambia")+
theme_bw()
event_transformative_group <- cbind.data.frame(ds$event_transformative_group_01, ds$event_transformative_group_02)
event_transformative_group <- na.omit(event_transformative_group)
names(event_transformative_group) <- sub('ds\\$event\\_', '', names(event_transformative_group))
mtx <- cor(event_transformative_group[, c(1:2)])
## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
## Reliability:
alph01 <- psych::alpha(event_transformative_group)
summary(alph01)
Reliability analysis
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.41 0.41 0.26 0.26 0.71 0.084 4.2 1.3 0.26## Visualization:
ds$event_transformative_group <- ((ds$event_transformative_group_01+ds$event_transformative_group_02)/2)
summary(ds$event_transformative_group) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 4 4 4 5 7 37 ds %>% drop_na(event_transformative_group)%>%
ggplot(aes(x = event_transformative_group))+
geom_histogram(color = "black",
fill = "gray",
bins = 50)+
geom_textvline(label = "Mean = 4.20",
xintercept = 4.20,
vjust = 1.1,
lwd = 1.05,
linetype = 2)+
labs(x = "Transformative event score",
y = "Frequency",
title = "Transformative event for group: Gambia")+
theme_bw()
imagistic_subscales <- cbind.data.frame(ds$event_positive_affect, ds$event_negative_affect,
ds$event_episodic_recall,
ds$event_shared_perception, ds$event_event_reflection,
ds$event_transformative_indiv, ds$event_transformative_group)
imagistic_subscales <- na.omit(imagistic_subscales)
names(imagistic_subscales) <- sub('ds\\$event\\_', '', names(imagistic_subscales))
mtx <- cor(imagistic_subscales[, c(1:6)])
## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
col=c("white", "darkred", "red",
"darkgrey", "blue", "darkblue"))
## Four regression models predicting IG/OG Fusion/Identification:
lm01 <- lm(IG_Fusion~event_positive_affect+event_negative_affect+event_episodic_recall+event_shared_perception+event_event_reflection+
event_transformative_indiv+event_transformative_group+Age+Female+Married+`SES-`,
data = ds)
lm02 <- lm(IG_Identification~event_positive_affect+event_negative_affect+event_episodic_recall+event_shared_perception+event_event_reflection+
event_transformative_indiv+event_transformative_group+Age+Female+Married+`SES-`,
data = ds)
lm03 <- lm(OG_Bonds~event_positive_affect+event_negative_affect+event_episodic_recall+event_shared_perception+event_event_reflection+
event_transformative_indiv+event_transformative_group+Age+Female+Married+`SES-`,
data = ds)## Tabulated results:
stargazer(lm01, lm02,
lm03,
type = "html",
star.cutoffs = c(0.05, 0.01, 0.001),
out = "table1.html")htmltools::includeHTML("table1.html")| Dependent variable: | |||
| IG_Fusion | IG_Identification | OG_Bonds | |
| (1) | (2) | (3) | |
| event_positive_affect | 0.027 | 0.110** | 0.025 |
| (0.039) | (0.037) | (0.064) | |
| event_negative_affect | 0.068 | 0.081 | 0.089 |
| (0.046) | (0.044) | (0.073) | |
| event_episodic_recall | 0.160* | 0.130 | -0.300* |
| (0.078) | (0.074) | (0.120) | |
| event_shared_perception | 0.069 | 0.150 | 0.037 |
| (0.082) | (0.080) | (0.130) | |
| event_event_reflection | 0.024 | -0.030 | 0.270** |
| (0.049) | (0.047) | (0.081) | |
| event_transformative_indiv | 0.081 | 0.110* | -0.120 |
| (0.055) | (0.053) | (0.091) | |
| event_transformative_group | 0.018 | -0.096 | 0.047 |
| (0.066) | (0.063) | (0.110) | |
| Age | -0.001 | -0.008 | 0.0004 |
| (0.006) | (0.006) | (0.010) | |
| Female | 0.055 | -0.110 | -0.560* |
| (0.160) | (0.150) | (0.260) | |
| Married | 0.001 | -0.100 | -0.032 |
| (0.170) | (0.170) | (0.280) | |
| `SES-`Middle | -0.250 | 0.021 | 0.079 |
| (0.200) | (0.190) | (0.320) | |
| `SES-`Upper middle | -0.170 | -0.079 | -0.610 |
| (0.240) | (0.230) | (0.390) | |
| `SES-`Upper | 0.710 | 0.870* | -0.230 |
| (0.430) | (0.410) | (0.690) | |
| Constant | 3.900*** | 4.100*** | 3.300*** |
| (0.520) | (0.500) | (0.840) | |
| Observations | 168 | 173 | 164 |
| R2 | 0.190 | 0.210 | 0.220 |
| Adjusted R2 | 0.120 | 0.150 | 0.150 |
| Residual Std. Error | 0.880 (df = 154) | 0.860 (df = 159) | 1.400 (df = 150) |
| F Statistic | 2.800** (df = 13; 154) | 3.300*** (df = 13; 159) | 3.300*** (df = 13; 150) |
| Note: | *p<0.05; **p<0.01; ***p<0.001 | ||
Social perception of religious freedom
Display code
Display code