Full Dataset (8 + 4 Countries)
Multi-level Modeling

Display code

rm(list=ls())
options(digits = 2)

## 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, gt, lme4, car, lmerTest, 
               ggeffects, magrittr, broom, broom.mixed,
               backports, effects, interactions, plyr, sjPlot,
               modelsummary, patchwork)

## Import latest BCL (8 countries) dataset:
ds01 <- fread("~/Desktop/oxford/data/BCL/BCL02.csv")
ds01 <- as.data.table(ds01)

## Import latest Combined (4 countries) dataset:
ds02 <- fread("~/Desktop/oxford/data/cleanedds/combinedds.csv")
ds02 <- as.data.table(ds02)

## List of variables needed from both datasets:

list1 <- c("ID", "Country",  "Endorse_BCL", "Endorse_BBL",
           "IG_fusion", "IG_identification", "OG_bonds",
           "Empathic_concern", "Perspective_taking", 
           "Perceived_discrimination", "OG_hostility", 
           "OG_cooperation", "Fight_OG",
           "Age", "Female", "Married", "Wealth_level", 
           "Event_positive_affect", "Event_negative_affect", 
           "Event_episodic_recall", "Event_shared_perception", 
           "Event_reflection", "Event_transformative_indiv", 
           "Event_transformative_group")

## Variable: ID
ds01$counter <- 1:nrow(ds01)
ds01$country2 <- ifelse(ds01$Country == "Bangladesh", "BAN",
                 ifelse(ds01$Country == "Ghana", "GHA",
                 ifelse(ds01$Country == "Malawi", "MWI",
                 ifelse(ds01$Country == "Pakistan", "PAK01",
                 ifelse(ds01$Country == "Sierra Leone", "SLE",
                 ifelse(ds01$Country == "Tanzania", "TZA01",
                 ifelse(ds01$Country == "Uganda", "UGA01",
                 ifelse(ds01$Country == "USA", "USA", NA))))))))
ds01$ID <- paste0(ds01$country2,ds01$counter)

## Variable: Country
ds01$Country <- ds01$Country
ds02$Country <- ds02$Country

## Variable: Endorse_BCL
ds01$Endorse_BCL <- ds01$Endorse_BCL
ds02$Endorse_BCL <- ds02$BCL

## Variable: Endorse_BBL
ds01$Endorse_BBL <- ds01$Endorse_BBL
ds02$Endorse_BBL <- ds02$BBL

## Variable: IG_fusion
ds01$IG_fusion <- ds01$IG_Fusion
ds02$IG_fusion <- ds02$Ingroup_fusion

## Variable: IG_identification
ds01$IG_identification <- ds01$IG_Identification
ds02$IG_identification <- ds02$Ingroup_identification

## Variable: OG_bonds
ds01$OG_bonds <- ds01$OG_Bonds
ds02$OG_bonds <- ds02$Outgroup_bonds

## Variable: Empathic_concern
ds01$Empathic_concern <- ds01$Empathic_concern
ds02$Empathic_concern <- ds02$Empathic_concern

## Variable: Perspective_taking
ds01$Perspective_taking <- ds01$Perspective_taking
ds02$Perspective_taking <- ds02$Perspective_taking

## Variable: Perceived_discrimination: 
## Only available in ds02 (4 country combined dataset):
ds01$Perceived_discrimination <- NA
ds02$Perceived_discrimination <- ds02$Perceived_discrimination

## Variable: OG_hostility: 
## Only available in ds02 (4 country combined dataset):
ds01$OG_hostility <- NA
ds02$OG_hostility <- ds02$OG_hostility

## Variable: OG_cooperation: 
## Only available in ds02 (4 country combined dataset):
ds01$OG_cooperation <- NA
ds02$OG_cooperation <- ds02$OG_cooperation

## Variable: Fight_OG: 
## Only available in ds02 (4 country combined dataset):
ds01$Fight_OG <- NA
ds02$Fight_OG <- ds02$Fight_OG

## Variable: Age: 
ds01$Age <- as.numeric(ds01$Age)
ds02$Age <- as.numeric(ds02$Age)

## Variable: Female: 
ds01$Female <- ds01$Female
ds02$Female <- ifelse(ds02$gender == "Female", 1, 0)

## Variable: Married:
ds01$Married <- ds01$Married
ds02$Married <- ifelse(ds02$married == "Married", 1, 0)

## Variable: Wealth_level

## This variable needed a lot of work
## For ds01 (8 country ds), wealth_level is divided 1, 2, 3, 4 based on 
## quartiles of wealth_level responses which range from 0 to 100
## For ds02 (4 country ds), wealth level is divided into 1, 2, 3,4 based on
## SES responses (lower middle, middle, upper middle, upper)

ds01$Wealth_level <- ifelse(ds01$wealth_level %in% c(0:25), 1,
                     ifelse(ds01$wealth_level %in% c(26:50), 2,
                     ifelse(ds01$wealth_level %in% c(51:75), 3,
                     ifelse(ds01$wealth_level %in% c(76:100), 4, NA))))

ds01$Wealth_level <- factor(ds01$Wealth_level, 
                            levels = c("1", "2", "3", "4"))

ds01$Wealth_level <- as.numeric(ds01$Wealth_level)

ds02$Wealth_level <- ifelse(ds02$ses == "Lower middle", 1,
                     ifelse(ds02$ses == "Middle", 2,
                     ifelse(ds02$ses == "Upper middle", 3,
                     ifelse(ds02$ses == "Upper", 4, NA))))

ds02$Wealth_level <- factor(ds02$Wealth_level, 
                            levels = c("1", "2", "3", "4"))

ds02$Wealth_level <- as.numeric(ds02$Wealth_level)

## Variable: Event_positive affect
ds01$Event_positive_affect <- ds01$event_positive_affect
ds02$Event_positive_affect <- ds02$Event_positive_affect

## Variable: Event_negative affect
ds01$Event_negative_affect <- ds01$event_negative_affect
ds02$Event_negative_affect <- ds02$Event_negative_affect

## Variable: Event_episodic recall
ds01$Event_episodic_recall <- ds01$event_episodic_recall
ds02$Event_episodic_recall <- ds02$Event_episodic_recall

## Variable: Event_shared_perception
ds01$Event_shared_perception <- ds01$event_shared_perception
ds02$Event_shared_perception <- ds02$Event_shared_perception

## Variable: Event_reflection
ds01$Event_reflection <- ds01$event_event_reflection
ds02$Event_reflection <- ds02$Event_event_reflection

## Variable: Event_transformative_indiv
ds01$Event_transformative_indiv <- ds01$event_transformative_indiv
ds02$Event_transformative_indiv <- ds02$Event_transformative_indiv

## Variable: Event_transformative_group
ds01$Event_transformative_group <- ds01$event_transformative_group
ds02$Event_transformative_group <- ds02$Event_transformative_group

## Combine two datasets with needed variables only
ds01 <- ds01[, ..list1]
ds02 <- ds02[, ..list1]

ds <- rbind(ds01, ds02)

## ## Remove all objects in R except ds:
rm(list= ls()[!(ls() %in% c("ds"))])

df01 <- ds %>% drop_na(Endorse_BCL, Endorse_BBL, Country)

Section 1. Outcome: Endorsement of BCL

Note: Read Chapter 12 of Gelman & Hill (2016) for reference.

Unconditional means model

Also called “varying intercept model with no predictors” (Gelman and Hill, 2016, Chapter 12). Allows intercepts to randomly vary across countries:

Display code

## Varying intercept model with no predictors:
m00<- lmer(Endorse_BCL ~ 1 + (1 | Country), 
           data = ds)

summary(m00)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Endorse_BCL ~ 1 + (1 | Country)
   Data: ds

REML criterion at convergence: 12523

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.911 -0.565  0.256  0.715  1.686 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.196    0.443   
 Residual             2.442    1.563   
Number of obs: 3349, groups:  Country, 9

Fixed effects:
            Estimate Std. Error    df t value      Pr(>|t|)    
(Intercept)    5.035      0.153 8.071      33 0.00000000068 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random effects:

Variance for Intercept = 0.196. This is the variance of the means across level 1 categories (countries).

Display code

tab_model(m00)

	Endorse_BCL
Predictors	Estimates	CI	p
(Intercept)	5.04	4.74 – 5.33	<0.001
Random Effects
σ2	2.44
τ00 Country	0.20
ICC	0.07
N Country	9
Observations	3349
Marginal R2 / Conditional R2	0.000 / 0.074

We can see that ICC = 0.07. Lower ICC = low variance explained across groups. In this case, most of the variability is at individual-level (not group level). There is very little differing patterns between countries.

Random intercept model

Also called “varying intercept model with individual-level predictors” (Gelman and Hill, 2016, Chapter 12).

Display code

## Varying intercept models with individual-level predictors:

m01 <- lmer(Endorse_BCL~IG_fusion+IG_identification+OG_bonds+Empathic_concern+
              Perspective_taking+Age+Female+Married+Wealth_level+ 
              (1 | Country), 
            data = ds)

summary(m01)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
Endorse_BCL ~ IG_fusion + IG_identification + OG_bonds + Empathic_concern +  
    Perspective_taking + Age + Female + Married + Wealth_level +  
    (1 | Country)
   Data: ds

REML criterion at convergence: 8154

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.996 -0.538  0.123  0.654  3.798 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0326   0.181   
 Residual             1.4478   1.203   
Number of obs: 2523, groups:  Country, 8

Fixed effects:
                      Estimate  Std. Error          df t value
(Intercept)           1.109929    0.187006  250.177358    5.94
IG_fusion             0.109494    0.030774 2511.738528    3.56
IG_identification     0.367759    0.030388 2512.387934   12.10
OG_bonds              0.024270    0.014670 2509.199425    1.65
Empathic_concern      0.009422    0.027499 2438.639161    0.34
Perspective_taking    0.247633    0.029740 2443.966919    8.33
Age                   0.000364    0.002243 2512.793510    0.16
Female                0.001309    0.049590 2506.138158    0.03
Married               0.007210    0.055315 2509.822416    0.13
Wealth_level          0.077174    0.031997 2412.283088    2.41
                               Pr(>|t|)    
(Intercept)                0.0000000097 ***
IG_fusion                       0.00038 ***
IG_identification  < 0.0000000000000002 ***
OG_bonds                        0.09816 .  
Empathic_concern                0.73189    
Perspective_taking < 0.0000000000000002 ***
Age                             0.87119    
Female                          0.97895    
Married                         0.89630    
Wealth_level                    0.01594 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) IG_fsn IG_dnt OG_bnd Empth_ Prspc_ Age    Female Marrid
IG_fusion   -0.027                                                        
IG_idntfctn -0.063 -0.825                                                 
OG_bonds    -0.236 -0.075  0.117                                          
Empthc_cncr -0.338 -0.015 -0.076  0.028                                   
Prspctv_tkn -0.263 -0.124 -0.091 -0.085 -0.398                            
Age         -0.263 -0.026  0.028  0.033 -0.045 -0.078                     
Female      -0.093 -0.026  0.008  0.051 -0.053 -0.005  0.069              
Married      0.097 -0.009  0.014 -0.053 -0.066 -0.017 -0.422  0.007       
Wealth_levl -0.489  0.052 -0.052 -0.042  0.132  0.043  0.032  0.018  0.003

Display code

tab_model(m01)

	Endorse_BCL
Predictors	Estimates	CI	p
(Intercept)	1.11	0.74 – 1.48	<0.001
IG fusion	0.11	0.05 – 0.17	<0.001
IG identification	0.37	0.31 – 0.43	<0.001
OG bonds	0.02	-0.00 – 0.05	0.098
Empathic concern	0.01	-0.04 – 0.06	0.732
Perspective taking	0.25	0.19 – 0.31	<0.001
Age	0.00	-0.00 – 0.00	0.871
Female	0.00	-0.10 – 0.10	0.979
Married	0.01	-0.10 – 0.12	0.896
Wealth level	0.08	0.01 – 0.14	0.016
Random Effects
σ2	1.45
τ00 Country	0.03
ICC	0.02
N Country	8
Observations	2523
Marginal R2 / Conditional R2	0.363 / 0.377

Here, marginal R² is much higher compared to previous model. Adding individual-level predictors significantly increases explanatory power of the model. Again, evidence that most of the variation is at individual-level differences.

Display code

## Change class of all models so we can use stargazer():
class(m00) <- "lmerMod"
class(m01) <- "lmerMod"

## Tabulated results:
stargazer(m00, m01,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	Endorse_BCL
	(1)	(2)
IG_fusion		0.110***
		(0.031)
IG_identification		0.370***
		(0.030)
OG_bonds		0.024
		(0.015)
Empathic_concern		0.009
		(0.027)
Perspective_taking		0.250***
		(0.030)
Age		0.0004
		(0.002)
Female		0.001
		(0.050)
Married		0.007
		(0.055)
Wealth_level		0.077*
		(0.032)
Constant	5.000***	1.100***
	(0.150)	(0.190)
Observations	3,349	2,523
Log Likelihood	-6,262.000	-4,077.000
Akaike Inf. Crit.	12,529.000	8,178.000
Bayesian Inf. Crit.	12,547.000	8,248.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Histogram: Endorsement of BCL

Display code

summary(df01$Endorse_BCL)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0     4.3     5.7     5.1     6.3     7.0 

Display code

ggplot(data = df01, 
       aes(x = Endorse_BCL)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.10", 
                 xintercept = 5.10, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "Frequency",
       x = "Endorse_BCL score", 
       title = "Endorse_BCL")+
  theme_bw()

Display code

ggplot(data = df01, 
       aes(x = Endorse_BCL)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  labs(y = "Frequency",
       x = "Endorse_BCL score", 
       title = "Endorse_BCL")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

tbl01 <- aggregate(df01$Endorse_BCL, 
                    by=list(df01$Country),
                    FUN=mean)
tbl01$Country <- tbl01$Group.1
tbl01$Endorse_BCL <- tbl01$x
tbl01 <- tbl01[, 3:4]
tbl01

       Country Endorse_BCL
1   Bangladesh         4.3
2       Gambia         5.9
3        Ghana         4.9
4       Malawi         5.0
5     Pakistan         5.1
6 Sierra Leone         4.7
7     Tanzania         5.0
8       Uganda         5.6
9          USA         4.8

Display code

ggplot(data = df01, 
       aes(x = Endorse_BCL, 
           y = Country)) +
  geom_boxplot(color = "black",
               fill = "grey")+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.10", 
                 xintercept = 5.10, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "",
       x = "Endorse_BCL score", 
       title = "Endorse_BCL")+
  theme_bw()

Faceted plots: Endorse BCL vs Ingroup fusion

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = IG_fusion)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "IG_fusion score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Ingroup fusion")+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = IG_fusion, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(x = "IG_fusion score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Ingroup fusion")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = IG_fusion)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "IG_fusion score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Ingroup fusion")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Endorse BCL vs Ingroup identification

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = IG_identification))+
  xlim(1, 7)+
  ylim(1, 7)+
  geom_point() +
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "IG_identification score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Ingroup identification")+
#  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = IG_identification, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(x = "IG_identification score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Ingroup identification")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = IG_identification))+
  xlim(1, 7)+
  ylim(1, 7)+
  geom_point() +
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "IG_identification score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Ingroup identification")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Endorse BCL vs Outgroup bonds

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = OG_bonds)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "OG_bonds score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Outgroup bonds")+
#  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = OG_bonds, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(x = "OG_bonds score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Outgroup bonds")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = OG_bonds)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "OG_bonds score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Outgroup bonds")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Endorse BCL vs Empathetic concern

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = Empathic_concern)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "Empathic_concern score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Empathetic concern")+
#  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = Empathic_concern, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(x = "Empathic_concern score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Empathetic concern")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = Empathic_concern)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "Empathic_concern score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Empathetic concern")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Endorse BCL vs Perspective taking

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = Perspective_taking)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "Perspective_taking score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Perspective taking")+
#  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = Perspective_taking, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(x = "Perspective_taking score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Perspective taking")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BCL,
           x = Perspective_taking)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "Perspective_taking score",
       y = "Endorse_BCL score", 
       title = "Endorse_BCL vs Perspective taking")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Section 2. Outcome: Endorsement of BBL

Unconditional means model

Display code

## Varying intercept model with no predictors:

m02<- lmer(Endorse_BBL ~ 1 + (1 | Country), 
           data = ds)

summary(m02)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Endorse_BBL ~ 1 + (1 | Country)
   Data: ds

REML criterion at convergence: 13200

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.431 -0.596  0.185  0.763  1.792 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.203    0.451   
 Residual             2.990    1.729   
Number of obs: 3349, groups:  Country, 9

Fixed effects:
            Estimate Std. Error    df t value      Pr(>|t|)    
(Intercept)    4.762      0.156 8.419    30.5 0.00000000065 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Display code

tab_model(m02)

	Endorse_BBL
Predictors	Estimates	CI	p
(Intercept)	4.76	4.46 – 5.07	<0.001
Random Effects
σ2	2.99
τ00 Country	0.20
ICC	0.06
N Country	9
Observations	3349
Marginal R2 / Conditional R2	0.000 / 0.064

Random intercept model

Display code

## Varying intercept models with individual-level predictors:

m03 <- lmer(Endorse_BBL~IG_fusion+IG_identification+OG_bonds+Empathic_concern+
              Perspective_taking+Age+Female+Married+Wealth_level+ 
              (1 | Country), 
            data = ds)

summary(m03)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
Endorse_BBL ~ IG_fusion + IG_identification + OG_bonds + Empathic_concern +  
    Perspective_taking + Age + Female + Married + Wealth_level +  
    (1 | Country)
   Data: ds

REML criterion at convergence: 9325

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.671 -0.641  0.149  0.709  3.197 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.131    0.362   
 Residual             2.312    1.520   
Number of obs: 2520, groups:  Country, 8

Fixed effects:
                     Estimate Std. Error         df t value
(Intercept)           3.85805    0.25710   96.47321   15.01
IG_fusion             0.16349    0.03880 2506.98545    4.21
IG_identification     0.34222    0.03827 2507.83053    8.94
OG_bonds             -0.06075    0.01858 2509.71894   -3.27
Empathic_concern     -0.40816    0.03492 2501.53359  -11.69
Perspective_taking    0.06444    0.03773 2504.31918    1.71
Age                  -0.00697    0.00283 2508.03488   -2.46
Female               -0.07487    0.06277 2509.99094   -1.19
Married              -0.15589    0.06999 2509.62182   -2.23
Wealth_level          0.14831    0.04058 2502.90911    3.65
                               Pr(>|t|)    
(Intercept)        < 0.0000000000000002 ***
IG_fusion                      0.000026 ***
IG_identification  < 0.0000000000000002 ***
OG_bonds                        0.00109 ** 
Empathic_concern   < 0.0000000000000002 ***
Perspective_taking              0.08782 .  
Age                             0.01389 *  
Female                          0.23308    
Married                         0.02603 *  
Wealth_level                    0.00026 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) IG_fsn IG_dnt OG_bnd Empth_ Prspc_ Age    Female Marrid
IG_fusion   -0.021                                                        
IG_idntfctn -0.057 -0.825                                                 
OG_bonds    -0.211 -0.074  0.114                                          
Empthc_cncr -0.317 -0.016 -0.077  0.025                                   
Prspctv_tkn -0.249 -0.126 -0.089 -0.087 -0.389                            
Age         -0.241 -0.031  0.032  0.031 -0.045 -0.074                     
Female      -0.082 -0.030  0.014  0.047 -0.058 -0.002  0.066              
Married      0.083 -0.004  0.008 -0.049 -0.058 -0.018 -0.421  0.010       
Wealth_levl -0.449  0.051 -0.053 -0.040  0.134  0.044  0.026  0.012  0.009

Display code

tab_model(m03)

	Endorse_BBL
Predictors	Estimates	CI	p
(Intercept)	3.86	3.35 – 4.36	<0.001
IG fusion	0.16	0.09 – 0.24	<0.001
IG identification	0.34	0.27 – 0.42	<0.001
OG bonds	-0.06	-0.10 – -0.02	0.001
Empathic concern	-0.41	-0.48 – -0.34	<0.001
Perspective taking	0.06	-0.01 – 0.14	0.088
Age	-0.01	-0.01 – -0.00	0.014
Female	-0.07	-0.20 – 0.05	0.233
Married	-0.16	-0.29 – -0.02	0.026
Wealth level	0.15	0.07 – 0.23	<0.001
Random Effects
σ2	2.31
τ00 Country	0.13
ICC	0.05
N Country	8
Observations	2520
Marginal R2 / Conditional R2	0.206 / 0.249

Display code

## Change class of all models so we can use stargazer():
class(m02) <- "lmerMod"
class(m03) <- "lmerMod"

## Tabulated results:
stargazer(m02, m03,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	Endorse_BBL
	(1)	(2)
IG_fusion		0.160***
		(0.039)
IG_identification		0.340***
		(0.038)
OG_bonds		-0.061**
		(0.019)
Empathic_concern		-0.410***
		(0.035)
Perspective_taking		0.064
		(0.038)
Age		-0.007*
		(0.003)
Female		-0.075
		(0.063)
Married		-0.160*
		(0.070)
Wealth_level		0.150***
		(0.041)
Constant	4.800***	3.900***
	(0.160)	(0.260)
Observations	3,349	2,520
Log Likelihood	-6,600.000	-4,662.000
Akaike Inf. Crit.	13,206.000	9,349.000
Bayesian Inf. Crit.	13,224.000	9,419.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Histogram: Endorsement of BBL

Display code

summary(df01$Endorse_BBL)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0     3.3     5.0     4.7     6.0     7.0 

Display code

ggplot(data = df01, 
       aes(x = Endorse_BBL)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 4.70", 
                 xintercept = 4.70, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "Frequency",
       x = "Endorse_BBL score", 
       title = "Endorse_BBL")+
  theme_bw()

Display code

ggplot(data = df01, 
       aes(x = Endorse_BBL)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  labs(y = "Frequency",
       x = "Endorse_BBL score", 
       title = "Endorse_BBL")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

tbl02 <- aggregate(df01$Endorse_BBL, 
                    by=list(df01$Country),
                    FUN=mean)
tbl02$Country <- tbl02$Group.1
tbl02$Endorse_BBL <- tbl02$x
tbl02 <- tbl02[, 3:4]
tbl02

       Country Endorse_BBL
1   Bangladesh         4.4
2       Gambia         4.6
3        Ghana         5.3
4       Malawi         5.2
5     Pakistan         3.9
6 Sierra Leone         4.6
7     Tanzania         4.9
8       Uganda         5.0
9          USA         5.1

Display code

ggplot(data = df01, 
       aes(x = Endorse_BBL, 
           y = Country)) +
  geom_boxplot(color = "black",
               fill = "grey")+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 4.70", 
                 xintercept = 4.70, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "",
       x = "Endorse_BBL score", 
       title = "Endorse_BBL")+
  theme_bw()

Faceted plots: Endorse BBL vs Ingroup fusion

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = IG_fusion)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "IG_fusion score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Ingroup fusion")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = IG_fusion, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(x = "IG_fusion score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Ingroup fusion")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = IG_fusion)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "IG_fusion score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Ingroup fusion")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Endorse BBL vs Ingroup identification

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = IG_identification)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "IG_identification score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Ingroup identification")+
#  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = IG_identification, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(x = "IG_identification score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Ingroup identification")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = IG_identification)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "IG_identification score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Ingroup identification")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Endorse BBL vs Outgroup bonds

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = OG_bonds)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "OG_bonds score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Outgroup bonds")+
#  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = OG_bonds, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(x = "OG_bonds score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Outgroup bonds")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = OG_bonds)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "OG_bonds score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Outgroup bonds")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Endorse BBL vs Empathetic concern

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = Empathic_concern)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE)+
  ylim(1, 7)+
  labs(x = "Empathic_concern score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Empathetic concern")+
#  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = Empathic_concern, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(x = "Empathic concern",
       y = "Endorse_BBL score", 
       title = "Endorse BBL vs Empathic concern")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = Empathic_concern)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE)+
  ylim(1, 7)+
  labs(x = "Empathic_concern score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Empathetic concern")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Endorse BBL vs Perspective taking

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = Perspective_taking))+
  xlim(1, 7)+
  ylim(1, 7)+
  geom_point() +
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "Perspective_taking score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Perspective taking")+
#    facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = Perspective_taking, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(x = "Perspective_taking score",
       y = "Empathic_concern score", 
       title = "Empathic_concern vs Perspective taking")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = Endorse_BBL,
           x = Perspective_taking))+
  xlim(1, 7)+
  ylim(1, 7)+
  geom_point() +
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(x = "Perspective_taking score",
       y = "Endorse_BBL score", 
       title = "Endorse_BBL vs Perspective taking")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Section 3. Combined MLM results: BCL/BBL vs Group fusion/identification

Display code

## Tabulated results:
stargazer(m00, m01, m02, m03,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	Endorse_BCL		Endorse_BBL
	(1)	(2)	(3)	(4)
IG_fusion		0.110***		0.160***
		(0.031)		(0.039)
IG_identification		0.370***		0.340***
		(0.030)		(0.038)
OG_bonds		0.024		-0.061**
		(0.015)		(0.019)
Empathic_concern		0.009		-0.410***
		(0.027)		(0.035)
Perspective_taking		0.250***		0.064
		(0.030)		(0.038)
Age		0.0004		-0.007*
		(0.002)		(0.003)
Female		0.001		-0.075
		(0.050)		(0.063)
Married		0.007		-0.160*
		(0.055)		(0.070)
Wealth_level		0.077*		0.150***
		(0.032)		(0.041)
Constant	5.000***	1.100***	4.800***	3.900***
	(0.150)	(0.190)	(0.160)	(0.260)
Observations	3,349	2,523	3,349	2,520
Log Likelihood	-6,262.000	-4,077.000	-6,600.000	-4,662.000
Akaike Inf. Crit.	12,529.000	8,178.000	13,206.000	9,349.000
Bayesian Inf. Crit.	12,547.000	8,248.000	13,224.000	9,419.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Section 4. Outcome: Ingroup fusion

Unconditional means model

Display code

## Varying intercept model with no predictors:

m04 <- lmer(IG_fusion ~ 1 + (1 | Country), 
            data = ds)

summary(m04)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: IG_fusion ~ 1 + (1 | Country)
   Data: ds

REML criterion at convergence: 12972

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.903 -0.438  0.322  0.725  1.525 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.231    0.48    
 Residual             2.735    1.65    
Number of obs: 3367, groups:  Country, 9

Fixed effects:
            Estimate Std. Error    df t value     Pr(>|t|)    
(Intercept)    5.159      0.165 7.723    31.2 0.0000000021 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Display code

tab_model(m04)

	IG_fusion
Predictors	Estimates	CI	p
(Intercept)	5.16	4.83 – 5.48	<0.001
Random Effects
σ2	2.73
τ00 Country	0.23
ICC	0.08
N Country	9
Observations	3367
Marginal R2 / Conditional R2	0.000 / 0.078

Random intercept model

Display code

## Varying intercept models with individual-level predictors:

m05 <- lmer(IG_fusion~Event_positive_affect+Event_negative_affect+
              Event_episodic_recall+Event_shared_perception+Event_reflection+
              Event_transformative_indiv+Event_transformative_group+Age+Female+
              Married+Wealth_level+
              (1 | Country), data = ds)

summary(m05)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
IG_fusion ~ Event_positive_affect + Event_negative_affect + Event_episodic_recall +  
    Event_shared_perception + Event_reflection + Event_transformative_indiv +  
    Event_transformative_group + Age + Female + Married + Wealth_level +  
    (1 | Country)
   Data: ds

REML criterion at convergence: 9756

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.969 -0.559  0.057  0.606  4.924 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0916   0.303   
 Residual             1.2221   1.105   
Number of obs: 3181, groups:  Country, 9

Fixed effects:
                             Estimate Std. Error         df t value
(Intercept)                   0.98570    0.15300   34.32360    6.44
Event_positive_affect         0.04237    0.01118 3167.56853    3.79
Event_negative_affect         0.01405    0.01021 3167.43265    1.38
Event_episodic_recall         0.34054    0.01798 3168.97904   18.94
Event_shared_perception       0.16285    0.01717 3168.92299    9.49
Event_reflection              0.14297    0.01811 3168.98836    7.89
Event_transformative_indiv    0.13732    0.01789 3168.98451    7.68
Event_transformative_group    0.03710    0.01692 3165.67974    2.19
Age                           0.00462    0.00181 3143.10092    2.55
Female                        0.06058    0.04031 3166.06214    1.50
Married                       0.03887    0.04564 3168.76190    0.85
Wealth_level                 -0.14450    0.02487 3167.31166   -5.81
                                       Pr(>|t|)    
(Intercept)                   0.000000221148967 ***
Event_positive_affect                   0.00015 ***
Event_negative_affect                   0.16889    
Event_episodic_recall      < 0.0000000000000002 ***
Event_shared_perception    < 0.0000000000000002 ***
Event_reflection              0.000000000000004 ***
Event_transformative_indiv    0.000000000000021 ***
Event_transformative_group              0.02840 *  
Age                                     0.01089 *  
Female                                  0.13300    
Married                                 0.39443    
Wealth_level                  0.000000006902577 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
                 (Intr) Evnt_pst_ Evnt_n_ Evnt_psd_ Evnt_s_ Evnt_r
Evnt_pstv_f      -0.213                                           
Evnt_ngtv_f      -0.267  0.549                                    
Evnt_psdc_r      -0.067 -0.113    -0.060                          
Evnt_shrd_p      -0.065  0.035    -0.080  -0.332                  
Evnt_rflctn      -0.029 -0.019    -0.092  -0.262    -0.125        
Evnt_trnsfrmtv_n -0.051 -0.138    -0.054  -0.276     0.008  -0.318
Evnt_trnsfrmtv_g -0.026 -0.221    -0.023   0.035    -0.412  -0.141
Age              -0.341  0.036    -0.006  -0.059     0.013  -0.052
Female           -0.108  0.017     0.019  -0.025    -0.018   0.009
Married           0.071  0.006    -0.026  -0.031    -0.029   0.003
Wealth_levl      -0.353 -0.061     0.007   0.011    -0.001   0.001
                 Evnt_trnsfrmtv_n Evnt_trnsfrmtv_g Age    Female Marrid
Evnt_pstv_f                                                            
Evnt_ngtv_f                                                            
Evnt_psdc_r                                                            
Evnt_shrd_p                                                            
Evnt_rflctn                                                            
Evnt_trnsfrmtv_n                                                       
Evnt_trnsfrmtv_g -0.246                                                
Age               0.041            0.033                               
Female           -0.050            0.060            0.011              
Married           0.024           -0.020           -0.433  0.012       
Wealth_levl       0.028           -0.021            0.029  0.004  0.030

Display code

tab_model(m05)

	IG_fusion
Predictors	Estimates	CI	p
(Intercept)	0.99	0.69 – 1.29	<0.001
Event positive affect	0.04	0.02 – 0.06	<0.001
Event negative affect	0.01	-0.01 – 0.03	0.169
Event episodic recall	0.34	0.31 – 0.38	<0.001
Event shared perception	0.16	0.13 – 0.20	<0.001
Event reflection	0.14	0.11 – 0.18	<0.001
Event transformative indiv	0.14	0.10 – 0.17	<0.001
Event transformative group	0.04	0.00 – 0.07	0.028
Age	0.00	0.00 – 0.01	0.011
Female	0.06	-0.02 – 0.14	0.133
Married	0.04	-0.05 – 0.13	0.394
Wealth level	-0.14	-0.19 – -0.10	<0.001
Random Effects
σ2	1.22
τ00 Country	0.09
ICC	0.07
N Country	9
Observations	3181
Marginal R2 / Conditional R2	0.553 / 0.584

Display code

## Change class of all models so we can use stargazer():
class(m04) <- "lmerMod"
class(m05) <- "lmerMod"

## Tabulated results:
stargazer(m04, m05,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	IG_fusion
	(1)	(2)
Event_positive_affect		0.042***
		(0.011)
Event_negative_affect		0.014
		(0.010)
Event_episodic_recall		0.340***
		(0.018)
Event_shared_perception		0.160***
		(0.017)
Event_reflection		0.140***
		(0.018)
Event_transformative_indiv		0.140***
		(0.018)
Event_transformative_group		0.037*
		(0.017)
Age		0.005*
		(0.002)
Female		0.061
		(0.040)
Married		0.039
		(0.046)
Wealth_level		-0.140***
		(0.025)
Constant	5.200***	0.990***
	(0.170)	(0.150)
Observations	3,367	3,181
Log Likelihood	-6,486.000	-4,878.000
Akaike Inf. Crit.	12,978.000	9,784.000
Bayesian Inf. Crit.	12,996.000	9,869.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Histogram: Ingroup fusion

Display code

df01 <- ds %>% drop_na(IG_fusion, IG_identification)

summary(df01$IG_fusion)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0     4.7     6.0     5.4     6.7     7.0 

Display code

ggplot(data = df01, 
       aes(x = IG_fusion)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.40", 
                 xintercept = 5.40, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "Frequency",
       x = "IG_fusion score", 
       title = "IG_fusion")+
  theme_bw()

Display code

ggplot(data = df01, 
       aes(x = IG_fusion)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  labs(y = "Frequency",
       x = "IG_fusion score", 
       title = "IG_fusion")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

tbl02 <- aggregate(df01$IG_fusion, 
                    by=list(df01$Country),
                    FUN=mean)
tbl02$Country <- tbl02$Group.1
tbl02$IG_fusion <- tbl02$x
tbl02 <- tbl02[, 3:4]
tbl02

       Country IG_fusion
1   Bangladesh       4.6
2       Gambia       5.8
3        Ghana       4.9
4       Malawi       4.9
5     Pakistan       5.4
6 Sierra Leone       4.4
7     Tanzania       5.4
8       Uganda       5.8
9          USA       5.1

Display code

ggplot(data = df01, 
       aes(x = IG_fusion, 
           y = Country)) +
  geom_boxplot(color = "black",
               fill = "grey")+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.40", 
                 xintercept = 5.40, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "",
       x = "IG_fusion score", 
       title = "IG_fusion")+
  theme_bw()

Faceted plots: Ingroup fusion vs Positive affect about event

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_positive_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_positive_affect score", 
       title = "IG_fusion vs Event_positive_affect")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_positive_affect, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_fusion score",
       x = "Event_positive_affect score", 
       title = "IG_fusion vs Event_positive_affect")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_positive_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_positive_affect score", 
       title = "IG_fusion vs Event_positive_affect")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup fusion vs Negative affect about event

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_negative_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_negative_affect score", 
       title = "IG_fusion vs Event_negative_affect")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_negative_affect, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_fusion score",
       x = "Event_negative_affect score", 
       title = "IG_fusion vs Event_negative_affect")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_negative_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_negative_affect score", 
       title = "IG_fusion vs Event_negative_affect")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup fusion vs Episodic recall about event

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_episodic_recall)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_episodic_recall score", 
       title = "IG_fusion vs Event_episodic_recall")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_episodic_recall, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_fusion score",
       x = "Event_episodic_recall score", 
       title = "IG_fusion vs Event_episodic_recall")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_episodic_recall)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_episodic_recall score", 
       title = "IG_fusion vs Event_episodic_recall")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup fusion vs Shared perception about event

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_shared_perception)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_shared_perception score", 
       title = "IG_fusion vs Event_shared_perception")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_shared_perception, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_fusion score",
       x = "Event_shared_perception score", 
       title = "IG_fusion vs Event_shared_perception")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_shared_perception)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_shared_perception score", 
       title = "IG_fusion vs Event_shared_perception")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup fusion vs Reflection of event

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_reflection)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_reflection score", 
       title = "IG_fusion vs Event_reflection")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_reflection, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_fusion score",
       x = "Event_reflection score", 
       title = "IG_fusion vs Event_reflection")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_reflection)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_reflection score", 
       title = "IG_fusion vs Event_reflection")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup fusion vs Transformative event for individual

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_transformative_indiv)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_transformative_individual score", 
       title = "IG_fusion vs Event_transformative_individual")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_transformative_indiv, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_fusion score",
       x = "Event_transformative_individual score", 
       title = "IG_fusion vs Event_transformative_individual")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_transformative_indiv)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_transformative_individual score", 
       title = "IG_fusion vs Event_transformative_individual")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup fusion vs Transformative event for group

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_transformative_group)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_transformative_group score", 
       title = "IG_fusion vs Event_transformative_group")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_transformative_group, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_fusion score",
       x = "Event_transformative_group score", 
       title = "IG_fusion vs Event_transformative_group")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_fusion,
           x = Event_transformative_group)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_fusion score",
       x = "Event_transformative_group score", 
       title = "IG_fusion vs Event_transformative_group")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Section 5. Outcome: Ingroup identification

Unconditional means model

Display code

## Varying intercept model with no predictors:

m06 <- lmer(IG_identification ~ 1 + (1 | Country), 
            data = ds)

summary(m06)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: IG_identification ~ 1 + (1 | Country)
   Data: ds

REML criterion at convergence: 13102

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.913 -0.384  0.314  0.709  1.391 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.211    0.459   
 Residual             2.856    1.690   
Number of obs: 3363, groups:  Country, 9

Fixed effects:
            Estimate Std. Error    df t value      Pr(>|t|)    
(Intercept)    5.287      0.159 7.840    33.3 0.00000000099 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Display code

tab_model(m06)

	IG_identification
Predictors	Estimates	CI	p
(Intercept)	5.29	4.98 – 5.60	<0.001
Random Effects
σ2	2.86
τ00 Country	0.21
ICC	0.07
N Country	9
Observations	3363
Marginal R2 / Conditional R2	0.000 / 0.069

Random intercept model

Display code

## Varying intercept models with individual-level predictors:

m07 <- lmer(IG_identification~Event_positive_affect+Event_negative_affect+
              Event_episodic_recall+Event_shared_perception+Event_reflection+
              Event_transformative_indiv+Event_transformative_group+Age+Female+
              Married+Wealth_level+
              (1 | Country), data = ds)

summary(m07)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: IG_identification ~ Event_positive_affect + Event_negative_affect +  
    Event_episodic_recall + Event_shared_perception + Event_reflection +  
    Event_transformative_indiv + Event_transformative_group +  
    Age + Female + Married + Wealth_level + (1 | Country)
   Data: ds

REML criterion at convergence: 9726

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.489 -0.552  0.061  0.576  4.912 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0606   0.246   
 Residual             1.2165   1.103   
Number of obs: 3177, groups:  Country, 9

Fixed effects:
                              Estimate  Std. Error          df t value
(Intercept)                   1.121408    0.141025   52.228130    7.95
Event_positive_affect         0.021585    0.011150 3164.818208    1.94
Event_negative_affect        -0.013572    0.010198 3164.867264   -1.33
Event_episodic_recall         0.403666    0.017938 3164.282511   22.50
Event_shared_perception       0.209466    0.017170 3162.081281   12.20
Event_reflection              0.087462    0.018085 3163.878031    4.84
Event_transformative_indiv    0.170666    0.017816 3164.231134    9.58
Event_transformative_group   -0.006992    0.016918 3163.262464   -0.41
Age                           0.000873    0.001807 3105.288036    0.48
Female                        0.064235    0.040289 3163.756395    1.59
Married                      -0.014824    0.045538 3160.878991   -0.33
Wealth_level                 -0.105825    0.024801 3154.119175   -4.27
                                       Pr(>|t|)    
(Intercept)                       0.00000000015 ***
Event_positive_affect                     0.053 .  
Event_negative_affect                     0.183    
Event_episodic_recall      < 0.0000000000000002 ***
Event_shared_perception    < 0.0000000000000002 ***
Event_reflection                  0.00000138770 ***
Event_transformative_indiv < 0.0000000000000002 ***
Event_transformative_group                0.679    
Age                                       0.629    
Female                                    0.111    
Married                                   0.745    
Wealth_level                      0.00002039865 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
                 (Intr) Evnt_pst_ Evnt_n_ Evnt_psd_ Evnt_s_ Evnt_r
Evnt_pstv_f      -0.231                                           
Evnt_ngtv_f      -0.290  0.549                                    
Evnt_psdc_r      -0.073 -0.113    -0.062                          
Evnt_shrd_p      -0.069  0.037    -0.079  -0.334                  
Evnt_rflctn      -0.032 -0.015    -0.090  -0.264    -0.123        
Evnt_trnsfrmtv_n -0.054 -0.139    -0.051  -0.271     0.007  -0.322
Evnt_trnsfrmtv_g -0.029 -0.223    -0.025   0.035    -0.413  -0.141
Age              -0.371  0.037    -0.004  -0.060     0.010  -0.049
Female           -0.118  0.017     0.021  -0.028    -0.016   0.012
Married           0.075  0.008    -0.024  -0.030    -0.028  -0.002
Wealth_levl      -0.382 -0.064     0.005   0.012    -0.003   0.001
                 Evnt_trnsfrmtv_n Evnt_trnsfrmtv_g Age    Female Marrid
Evnt_pstv_f                                                            
Evnt_ngtv_f                                                            
Evnt_psdc_r                                                            
Evnt_shrd_p                                                            
Evnt_rflctn                                                            
Evnt_trnsfrmtv_n                                                       
Evnt_trnsfrmtv_g -0.246                                                
Age               0.038            0.037                               
Female           -0.050            0.059            0.012              
Married           0.027           -0.023           -0.431  0.014       
Wealth_levl       0.027           -0.019            0.032  0.004  0.028

Display code

tab_model(m07)

	IG_identification
Predictors	Estimates	CI	p
(Intercept)	1.12	0.84 – 1.40	<0.001
Event positive affect	0.02	-0.00 – 0.04	0.053
Event negative affect	-0.01	-0.03 – 0.01	0.183
Event episodic recall	0.40	0.37 – 0.44	<0.001
Event shared perception	0.21	0.18 – 0.24	<0.001
Event reflection	0.09	0.05 – 0.12	<0.001
Event transformative indiv	0.17	0.14 – 0.21	<0.001
Event transformative group	-0.01	-0.04 – 0.03	0.679
Age	0.00	-0.00 – 0.00	0.629
Female	0.06	-0.01 – 0.14	0.111
Married	-0.01	-0.10 – 0.07	0.745
Wealth level	-0.11	-0.15 – -0.06	<0.001
Random Effects
σ2	1.22
τ00 Country	0.06
ICC	0.05
N Country	9
Observations	3177
Marginal R2 / Conditional R2	0.578 / 0.598

Display code

## Change class of all models so we can use stargazer():
class(m06) <- "lmerMod"
class(m07) <- "lmerMod"

## Tabulated results:
stargazer(m06, m07,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	IG_identification
	(1)	(2)
Event_positive_affect		0.022
		(0.011)
Event_negative_affect		-0.014
		(0.010)
Event_episodic_recall		0.400***
		(0.018)
Event_shared_perception		0.210***
		(0.017)
Event_reflection		0.087***
		(0.018)
Event_transformative_indiv		0.170***
		(0.018)
Event_transformative_group		-0.007
		(0.017)
Age		0.001
		(0.002)
Female		0.064
		(0.040)
Married		-0.015
		(0.046)
Wealth_level		-0.110***
		(0.025)
Constant	5.300***	1.100***
	(0.160)	(0.140)
Observations	3,363	3,177
Log Likelihood	-6,551.000	-4,863.000
Akaike Inf. Crit.	13,108.000	9,754.000
Bayesian Inf. Crit.	13,126.000	9,839.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Histogram: Ingroup identification

Display code

df01 <- ds %>% drop_na(IG_identification)

summary(df01$IG_identification)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0     4.7     6.0     5.4     6.7     7.0 

Display code

ggplot(data = df01, 
       aes(x = IG_identification)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.40", 
                 xintercept = 5.40, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "Frequency",
       x = "IG_identification score", 
       title = "IG_identification")+
  theme_bw()

Display code

ggplot(data = df01, 
       aes(x = IG_identification)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  labs(y = "Frequency",
       x = "IG_identification score", 
       title = "IG_identification")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

tbl02 <- aggregate(df01$IG_identification, 
                    by=list(df01$Country),
                    FUN=mean)
tbl02$Country <- tbl02$Group.1
tbl02$IG_identification <- tbl02$x
tbl02 <- tbl02[, 3:4]
tbl02

       Country IG_identification
1   Bangladesh               4.6
2       Gambia               6.0
3        Ghana               5.1
4       Malawi               5.1
5     Pakistan               5.4
6 Sierra Leone               4.6
7     Tanzania               5.4
8       Uganda               5.9
9          USA               5.3

Display code

ggplot(data = df01, 
       aes(x = IG_identification, 
           y = Country)) +
  geom_boxplot(color = "black",
               fill = "grey")+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 5.40", 
                 xintercept = 5.40, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "",
       x = "IG_identification score", 
       title = "IG_identification")+
  theme_bw()

Faceted plots: Ingroup identification vs Positive affect about event

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_positive_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_positive_affect score", 
       title = "IG_identification vs Event_positive_affect")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_positive_affect, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_identification score",
       x = "Event_positive_affect score", 
       title = "IG_identification vs Event_positive_affect")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_positive_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_positive_affect score", 
       title = "IG_identification vs Event_positive_affect")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup identification vs Negative affect about event

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_negative_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_negative_affect score", 
       title = "IG_identification vs Event_negative_affect")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_negative_affect, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_identification score",
       x = "Event_negative_affect score", 
       title = "IG_identification vs Event_negative_affect")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_negative_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_negative_affect score", 
       title = "IG_identification vs Event_negative_affect")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup identification vs Episodic recall about event

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_episodic_recall)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_episodic_recall score", 
       title = "IG_identification vs Event_episodic_recall")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_episodic_recall, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_identification score",
       x = "Event_episodic_recall score", 
       title = "IG_identification vs Event_episodic_recall")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_episodic_recall)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_episodic_recall score", 
       title = "IG_identification vs Event_episodic_recall")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup identification vs Shared perception about event

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_shared_perception)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_shared_perception score", 
       title = "IG_identification vs Event_shared_perception")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_shared_perception, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_identification score",
       x = "Event_shared_perception score", 
       title = "IG_identification vs Event_shared_perception")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_shared_perception)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_shared_perception score", 
       title = "IG_identification vs Event_shared_perception")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup identification vs Reflection of event

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_reflection)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_reflection score", 
       title = "IG_identification vs Event_reflection")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_reflection, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_identification score",
       x = "Event_reflection score", 
       title = "IG_identification vs Event_reflection")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_reflection)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_reflection score", 
       title = "IG_identification vs Event_reflection")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup identification vs Transformative event for individual

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_transformative_indiv)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_transformative_individual score", 
       title = "IG_identification vs Event_transformative_individual")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_transformative_indiv, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_identification score",
       x = "Event_transformative_individual score", 
       title = "IG_identification vs Event_transformative_individual")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_transformative_indiv)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_transformative_individual score", 
       title = "IG_identification vs Event_transformative_individual")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Ingroup identification vs Transformative event for group

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_transformative_group)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_transformative_group score", 
       title = "IG_identification vs Event_transformative_group")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_transformative_group, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "IG_identification score",
       x = "Event_transformative_group score", 
       title = "IG_identification vs Event_transformative_group")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_identification,
           x = Event_transformative_group)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "IG_identification score",
       x = "Event_transformative_group score", 
       title = "IG_identification vs Event_transformative_group")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Section 6. Outcome: Outgroup bonds

Unconditional means model

Display code

## Varying intercept model with no predictors:

m08 <- lmer(OG_bonds ~ 1 + (1 | Country), 
            data = ds)

summary(m08)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: OG_bonds ~ 1 + (1 | Country)
   Data: ds

REML criterion at convergence: 10499

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.571 -0.846 -0.145  0.695  2.476 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.114    0.337   
 Residual             2.774    1.665   
Number of obs: 2716, groups:  Country, 8

Fixed effects:
            Estimate Std. Error    df t value    Pr(>|t|)    
(Intercept)    3.166      0.129 6.880    24.6 0.000000058 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Display code

tab_model(m08)

	OG_bonds
Predictors	Estimates	CI	p
(Intercept)	3.17	2.91 – 3.42	<0.001
Random Effects
σ2	2.77
τ00 Country	0.11
ICC	0.04
N Country	8
Observations	2716
Marginal R2 / Conditional R2	0.000 / 0.039

Random intercept model

Display code

## Varying intercept models with individual-level predictors:

m09 <- lmer(OG_bonds~Event_positive_affect+Event_negative_affect+
              Event_episodic_recall+Event_shared_perception+Event_reflection+
              Event_transformative_indiv+Event_transformative_group+Age+Female+
              Married+Wealth_level+
              (1 | Country), data = ds)

summary(m09)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
OG_bonds ~ Event_positive_affect + Event_negative_affect + Event_episodic_recall +  
    Event_shared_perception + Event_reflection + Event_transformative_indiv +  
    Event_transformative_group + Age + Female + Married + Wealth_level +  
    (1 | Country)
   Data: ds

REML criterion at convergence: 9828

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0848 -0.8085 -0.0979  0.6932  2.7804 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0738   0.272   
 Residual             2.6223   1.619   
Number of obs: 2567, groups:  Country, 8

Fixed effects:
                             Estimate Std. Error         df t value
(Intercept)                   2.33772    0.21643  117.44629   10.80
Event_positive_affect         0.08160    0.01813 2554.93844    4.50
Event_negative_affect         0.15376    0.01692 2552.74567    9.09
Event_episodic_recall        -0.19476    0.02959 2546.67932   -6.58
Event_shared_perception      -0.03127    0.02795 2541.72313   -1.12
Event_reflection              0.07076    0.02957 2536.19102    2.39
Event_transformative_indiv    0.04696    0.02946 2535.39473    1.59
Event_transformative_group    0.08547    0.02771 2553.57445    3.08
Age                          -0.00398    0.00299 2554.50741   -1.33
Female                       -0.16525    0.06610 2552.91745   -2.50
Married                       0.16519    0.07404 2551.79720    2.23
Wealth_level                  0.09725    0.04195 2474.88527    2.32
                                       Pr(>|t|)    
(Intercept)                < 0.0000000000000002 ***
Event_positive_affect            0.000007061918 ***
Event_negative_affect      < 0.0000000000000002 ***
Event_episodic_recall            0.000000000057 ***
Event_shared_perception                  0.2633    
Event_reflection                         0.0168 *  
Event_transformative_indiv               0.1110    
Event_transformative_group               0.0021 ** 
Age                                      0.1834    
Female                                   0.0125 *  
Married                                  0.0258 *  
Wealth_level                             0.0205 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
                 (Intr) Evnt_pst_ Evnt_n_ Evnt_psd_ Evnt_s_ Evnt_r
Evnt_pstv_f      -0.268                                           
Evnt_ngtv_f      -0.327  0.574                                    
Evnt_psdc_r      -0.089 -0.110    -0.082                          
Evnt_shrd_p      -0.079  0.040    -0.085  -0.332                  
Evnt_rflctn      -0.053  0.007    -0.066  -0.246    -0.134        
Evnt_trnsfrmtv_n -0.072 -0.138    -0.043  -0.282     0.035  -0.329
Evnt_trnsfrmtv_g -0.010 -0.218    -0.028   0.036    -0.413  -0.137
Age              -0.385  0.029    -0.008  -0.060     0.000  -0.053
Female           -0.141  0.028     0.027  -0.017    -0.032   0.001
Married           0.062  0.008    -0.021  -0.027    -0.023   0.004
Wealth_levl      -0.438 -0.051     0.026  -0.001    -0.005   0.011
                 Evnt_trnsfrmtv_n Evnt_trnsfrmtv_g Age    Female Marrid
Evnt_pstv_f                                                            
Evnt_ngtv_f                                                            
Evnt_psdc_r                                                            
Evnt_shrd_p                                                            
Evnt_rflctn                                                            
Evnt_trnsfrmtv_n                                                       
Evnt_trnsfrmtv_g -0.264                                                
Age               0.048            0.042                               
Female           -0.048            0.051            0.052              
Married           0.023           -0.024           -0.432  0.008       
Wealth_levl       0.049           -0.038            0.045  0.029  0.030

Display code

tab_model(m09)

	OG_bonds
Predictors	Estimates	CI	p
(Intercept)	2.34	1.91 – 2.76	<0.001
Event positive affect	0.08	0.05 – 0.12	<0.001
Event negative affect	0.15	0.12 – 0.19	<0.001
Event episodic recall	-0.19	-0.25 – -0.14	<0.001
Event shared perception	-0.03	-0.09 – 0.02	0.263
Event reflection	0.07	0.01 – 0.13	0.017
Event transformative indiv	0.05	-0.01 – 0.10	0.111
Event transformative group	0.09	0.03 – 0.14	0.002
Age	-0.00	-0.01 – 0.00	0.183
Female	-0.17	-0.29 – -0.04	0.012
Married	0.17	0.02 – 0.31	0.026
Wealth level	0.10	0.01 – 0.18	0.021
Random Effects
σ2	2.62
τ00 Country	0.07
ICC	0.03
N Country	8
Observations	2567
Marginal R2 / Conditional R2	0.059 / 0.085

Display code

## Change class of all models so we can use stargazer():
class(m08) <- "lmerMod"
class(m09) <- "lmerMod"

## Tabulated results:
stargazer(m08, m09,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	OG_bonds
	(1)	(2)
Event_positive_affect		0.082***
		(0.018)
Event_negative_affect		0.150***
		(0.017)
Event_episodic_recall		-0.200***
		(0.030)
Event_shared_perception		-0.031
		(0.028)
Event_reflection		0.071*
		(0.030)
Event_transformative_indiv		0.047
		(0.029)
Event_transformative_group		0.085**
		(0.028)
Age		-0.004
		(0.003)
Female		-0.170*
		(0.066)
Married		0.170*
		(0.074)
Wealth_level		0.097*
		(0.042)
Constant	3.200***	2.300***
	(0.130)	(0.220)
Observations	2,716	2,567
Log Likelihood	-5,249.000	-4,914.000
Akaike Inf. Crit.	10,505.000	9,856.000
Bayesian Inf. Crit.	10,523.000	9,938.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Histogram: Outgroup bonds

Display code

df02 <- ds %>% drop_na(OG_bonds)

summary(df02$OG_bonds)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0     2.0     3.0     3.3     4.5     7.0 

Display code

ggplot(data = df02, 
       aes(x = OG_bonds)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 3.30", 
                 xintercept = 3.30, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "Frequency",
       x = "OG_bonds score", 
       title = "OG_bonds")+
  theme_bw()

Display code

ggplot(data = df01, 
       aes(x = OG_bonds)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  labs(y = "Frequency",
       x = "OG_bonds score", 
       title = "OG_bonds")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

tbl02 <- aggregate(df02$OG_bonds, 
                    by=list(df02$Country),
                    FUN=mean)
tbl02$Country <- tbl02$Group.1
tbl02$OG_bonds <- tbl02$x
tbl02 <- tbl02[, 3:4]
tbl02

       Country OG_bonds
1   Bangladesh      3.7
2       Gambia      2.9
3        Ghana      2.9
4       Malawi      2.8
5     Pakistan      3.6
6 Sierra Leone      2.8
7     Tanzania      3.4
8       Uganda      3.0

Display code

ggplot(data = df01, 
       aes(x = OG_bonds, 
           y = Country)) +
  geom_boxplot(color = "black",
               fill = "grey")+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 3.30", 
                 xintercept = 3.30, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(y = "",
       x = "OG_bonds score", 
       title = "OG_bonds")+
  theme_bw()

Faceted plots: Outgroup bonds vs Positive affect about event

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_positive_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_positive_affect score", 
       title = "OG_bonds vs Event_positive_affect")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_positive_affect, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "OG_bonds score",
       x = "Event_positive_affect score", 
       title = "OG_bonds vs Event_positive_affect")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_positive_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_positive_affect score", 
       title = "OG_bonds vs Event_positive_affect")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Outgroup bonds vs Negative affect about event

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_negative_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_negative_affect score", 
       title = "OG_bonds vs Event_negative_affect")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_negative_affect, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "OG_bonds score",
       x = "Event_negative_affect score", 
       title = "OG_bonds vs Event_negative_affect")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_negative_affect)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_negative_affect score", 
       title = "OG_bonds vs Event_negative_affect")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Outgroup bonds vs Episodic recall about event

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_episodic_recall)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_episodic_recall score", 
       title = "OG_bonds vs Event_episodic_recall")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_episodic_recall, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "OG_bonds score",
       x = "Event_episodic_recall score", 
       title = "OG_bonds vs Event_episodic_recall")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_episodic_recall)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_episodic_recall score", 
       title = "OG_bonds vs Event_episodic_recall")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Outgroup bonds vs Shared perception about event

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_shared_perception)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_shared_perception score", 
       title = "OG_bonds vs Event_shared_perception")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_shared_perception, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "OG_bonds score",
       x = "Event_shared_perception score", 
       title = "OG_bonds vs Event_shared_perception")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_shared_perception)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_shared_perception score", 
       title = "OG_bonds vs Event_shared_perception")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Outgroup bonds vs Reflection of event

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_reflection)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_reflection score", 
       title = "OG_bonds vs Event_reflection")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_reflection, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "OG_bonds score",
       x = "Event_reflection score", 
       title = "OG_bonds vs Event_reflection")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_reflection)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_reflection score", 
       title = "OG_bonds vs Event_reflection")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Outgroup bonds vs Transformative event for individual

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_transformative_indiv)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_transformative_individual score", 
       title = "OG_bonds vs Event_transformative_individual")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_transformative_indiv, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "OG_bonds score",
       x = "Event_transformative_individual score", 
       title = "OG_bonds vs Event_transformative_individual")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_transformative_indiv)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_transformative_individual score", 
       title = "OG_bonds vs Event_transformative_individual")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Faceted plots: Outgroup bonds vs Transformative event for group

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_transformative_group)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_transformative_group score", 
       title = "OG_bonds vs Event_transformative_group")+
  # facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_transformative_group, 
           col = Country, 
           linetype = Country)) +
  geom_point()+
  xlim(1, 7)+
 ylim(1, 7)+
  stat_smooth(method="lm", se = F) +
  labs(y = "OG_bonds score",
       x = "Event_transformative_group score", 
       title = "OG_bonds vs Event_transformative_group")+
 scale_color_colorblind()+
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_bonds,
           x = Event_transformative_group)) +
  geom_point()+
  xlim(1, 7)+
  ylim(1, 7)+
  stat_smooth(method="lm", 
              fullrange=TRUE) +
  labs(y = "OG_bonds score",
       x = "Event_transformative_group score", 
       title = "OG_bonds vs Event_transformative_group")+
  facet_wrap( ~ Country, nrow = 3) +
  theme_bw()

Section 7. Combined MLM results: Group fusion/identification vs Imagistic items

Display code

## Tabulated results:
stargazer(m04, m05, m06, m07, m08, m09,
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	IG_fusion		IG_identification		OG_bonds
	(1)	(2)	(3)	(4)	(5)	(6)
Event_positive_affect		0.042***		0.022		0.082***
		(0.011)		(0.011)		(0.018)
Event_negative_affect		0.014		-0.014		0.150***
		(0.010)		(0.010)		(0.017)
Event_episodic_recall		0.340***		0.400***		-0.200***
		(0.018)		(0.018)		(0.030)
Event_shared_perception		0.160***		0.210***		-0.031
		(0.017)		(0.017)		(0.028)
Event_reflection		0.140***		0.087***		0.071*
		(0.018)		(0.018)		(0.030)
Event_transformative_indiv		0.140***		0.170***		0.047
		(0.018)		(0.018)		(0.029)
Event_transformative_group		0.037*		-0.007		0.085**
		(0.017)		(0.017)		(0.028)
Age		0.005*		0.001		-0.004
		(0.002)		(0.002)		(0.003)
Female		0.061		0.064		-0.170*
		(0.040)		(0.040)		(0.066)
Married		0.039		-0.015		0.170*
		(0.046)		(0.046)		(0.074)
Wealth_level		-0.140***		-0.110***		0.097*
		(0.025)		(0.025)		(0.042)
Constant	5.200***	0.990***	5.300***	1.100***	3.200***	2.300***
	(0.170)	(0.150)	(0.160)	(0.140)	(0.130)	(0.220)
Observations	3,367	3,181	3,363	3,177	2,716	2,567
Log Likelihood	-6,486.000	-4,878.000	-6,551.000	-4,863.000	-5,249.000	-4,914.000
Akaike Inf. Crit.	12,978.000	9,784.000	13,108.000	9,754.000	10,505.000	9,856.000
Bayesian Inf. Crit.	12,996.000	9,869.000	13,126.000	9,839.000	10,523.000	9,938.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Display code

## Write final dataset file:
# fwrite(ds, file = "~/Desktop/oxford/data/BCL/BCL02.csv")

Section 8. Coefficient plots

Display code

bg <- list(geom_point(aes(y = term, x = estimate),
                      alpha = .25, size = 15, color = 'red'))

Outcome: Endorsement of BCL

Display code

mp01 <- modelplot(m01,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Endorsement of BCL", 
       x = "Coefficient estimates with 95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

mp01

Display code

cfp01 <- modelplot(m01,
                   coef_omit = 'SD (Observations)',
                   coef_rename = TRUE, 
                   background = bg)+
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant"))+
  scale_color_manual(values = c("red", "blue"))+
  #  xlim(-3, 3)+
  labs(title = "Endorsement of BCL", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

cfp01

Outcome: Endorsement of BBL

Display code

mp03 <- modelplot(m03,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Endorsement of BBL", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)


mp03

Display code

cfp03 <- modelplot(m03,
                   coef_omit = 'SD (Observations)',
                   coef_rename = TRUE, 
                   background = bg)+
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant"))+
  scale_color_manual(values = c("red", "blue"))+
  #  xlim(-3, 3)+
  labs(title = "Endorsement of BBL", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

cfp03

Combined: Endorsement of BCL & BBL

Display code

mp01a <- mp01 + theme(legend.position = "none")

ol1 <- mp01a / mp03

ol1

Outcome: Ingroup fusion

Display code

mp05 <- modelplot(m05,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Ingroup fusion", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

mp05

Display code

cfp05 <- modelplot(m05,
                   coef_omit = 'SD (Observations)',
                   coef_rename = TRUE, 
                   background = bg)+
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant"))+
  scale_color_manual(values = c("red", "blue"))+
  #  xlim(-3, 3)+
  labs(title = "Ingroup fusion", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

cfp05

Outcome: Ingroup identification

Display code

mp07 <- modelplot(m07,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Ingroup identification", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)


mp07

Display code

cfp07 <- modelplot(m03,
                   coef_omit = 'SD (Observations)',
                   coef_rename = TRUE, 
                   background = bg)+
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant"))+
  scale_color_manual(values = c("red", "blue"))+
  #  xlim(-3, 3)+
  labs(title = "Ingroup identification", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

cfp07

Combined: Ingroup Fusion & Identification

Display code

mp05a <- mp05 + theme(legend.position = "none")
ol2 <- mp05a / mp07
ol2

Outcome: Outcome bonds

Display code

mp09 <- modelplot(m09,
                  coef_omit = 'SD (Observations)',
                  coef_rename = TRUE) +
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant")) +
  scale_color_manual(values = c("red", "blue", "black"))+
#  xlim(-3, 3)+
  labs(title = "Outgroup bonds", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

mp09

Display code

cfp09 <- modelplot(m09,
                   coef_omit = 'SD (Observations)',
                   coef_rename = TRUE, 
                   background = bg)+
  aes(color = ifelse(p.value < 0.05, "Significant", "Not significant"))+
  scale_color_manual(values = c("red", "blue"))+
  #  xlim(-3, 3)+
  labs(title = "Outgroup bonds", 
       x = "Coefficient estimates with95% CI")+
  geom_vline(xintercept = 0, 
             linetype = 2)

cfp09

To Do

• Some countries/datasets have more variables than others. This makes it difficult to estimate all explanatory variables of interest (eg - perceived discrimination, negative/positive contact/etc)

• No group-level predictors in the data. For eg: country-level variables that might explain the variation in DVs. Eg: country-level trust/religiosity scores, GDP, HDI, etc. Multi-level modeling is more informative if there are theoretically relevant group-level predictors.

• Other DV/IV combinations

• Other comments/suggestions?