BCL Cross-country Data
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)

## Import latest BCL dataset:
ds <- fread("~/Desktop/oxford/data/BCL/BCL01.csv")
ds <- as.data.table(ds)

ds$Empathic_concern <- ds$empathic_concern
ds$Perspective_taking <- ds$perspective_taking

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

Section 1. Outcome: Endorsement of BCL

Note: Read pg. 259 of Gelman & Hill (2016) for reference. Additional instructions are also available at this online guide.

Unconditional means model

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

Display code

## Varying intercept model with no predictors:
m00<- lmer(formula = 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: 7270

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.298 -0.714  0.174  0.817  1.485 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0587   0.242   
 Residual             3.0734   1.753   
Number of obs: 1832, groups:  Country, 8

Fixed effects:
            Estimate Std. Error     df t value       Pr(>|t|)    
(Intercept)   4.7333     0.0978 7.6887    48.4 0.000000000078 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random intercept models

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(formula = 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: 3933

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.693 -0.708  0.077  0.700  3.795 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0133   0.115   
 Residual             1.6374   1.280   
Number of obs: 1166, groups:  Country, 7

Fixed effects:
                     Estimate Std. Error         df t value
(Intercept)           0.15584    0.37906 1050.65432    0.41
IG_Fusion             0.04592    0.04879 1155.92858    0.94
IG_Identification     0.43698    0.04744 1155.99989    9.21
OG_Bonds             -0.01722    0.02055 1145.16377   -0.84
Empathic_concern      0.13935    0.07401 1155.99156    1.88
Perspective_taking    0.39042    0.05907 1152.78834    6.61
Age                  -0.00240    0.00525 1155.59856   -0.46
Female                0.09286    0.07928 1115.92948    1.17
Married               0.02951    0.09894 1128.49603    0.30
Wealth_level          0.00233    0.00190 1110.62211    1.23
                               Pr(>|t|)    
(Intercept)                        0.68    
IG_Fusion                          0.35    
IG_Identification  < 0.0000000000000002 ***
OG_Bonds                           0.40    
Empathic_concern                   0.06 .  
Perspective_taking       0.000000000059 ***
Age                                0.65    
Female                             0.24    
Married                            0.77    
Wealth_level                       0.22    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) IG_Fsn IG_Idn OG_Bnd Empth_ Prspc_ Age    Female Marrid
IG_Fusion   -0.012                                                        
IG_Idntfctn  0.079 -0.877                                                 
OG_Bonds    -0.134 -0.151  0.182                                          
Empthc_cncr -0.664  0.006 -0.089 -0.014                                   
Prspctv_tkn -0.380 -0.087 -0.156 -0.045 -0.153                            
Age         -0.395  0.023 -0.007 -0.023  0.034 -0.004                     
Female      -0.115 -0.009 -0.006  0.073  0.020  0.004  0.031              
Married      0.157 -0.039  0.058 -0.033 -0.038  0.028 -0.551  0.048       
Wealth_levl -0.240  0.049 -0.072 -0.062  0.005  0.011  0.062  0.005 -0.080

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.046
		(0.049)
IG_Identification		0.440***
		(0.047)
OG_Bonds		-0.017
		(0.021)
Empathic_concern		0.140
		(0.074)
Perspective_taking		0.390***
		(0.059)
Age		-0.002
		(0.005)
Female		0.093
		(0.079)
Married		0.030
		(0.099)
Wealth_level		0.002
		(0.002)
Constant	4.700***	0.160
	(0.098)	(0.380)
Observations	1,832	1,166
Log Likelihood	-3,635.000	-1,966.000
Akaike Inf. Crit.	7,276.000	3,957.000
Bayesian Inf. Crit.	7,293.000	4,017.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     3.3     5.0     4.7     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 = 4.70", 
                 xintercept = 4.70, 
                 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 = 2) +
  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        Ghana         4.9
3       Malawi         5.0
4     Pakistan         4.5
5 Sierra Leone         4.7
6     Tanzania         4.6
7       Uganda         5.1
8          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 = 4.70", 
                 xintercept = 4.70, 
                 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")+
# facet_wrap( ~ Country, nrow = 2)+
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  theme_bw()

Section 2. Outcome: Endorsement of BBL

Unconditional means model

Display code

## Varying intercept model with no predictors:

m02<- lmer(formula = 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: 7454

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.302 -0.652  0.244  0.802  1.423 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.115    0.34    
 Residual             3.361    1.83    
Number of obs: 1836, groups:  Country, 8

Fixed effects:
            Estimate Std. Error   df t value     Pr(>|t|)    
(Intercept)     4.89       0.13 7.25    37.5 0.0000000014 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random intercept models

Display code

## Varying intercept models with individual-level predictors:

m03 <- lmer(formula = 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: 3840

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.717 -0.653  0.052  0.639  4.032 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0202   0.142   
 Residual             1.5056   1.227   
Number of obs: 1167, groups:  Country, 7

Fixed effects:
                     Estimate Std. Error         df t value
(Intercept)           0.25685    0.36473  997.00234    0.70
IG_Fusion             0.01866    0.04661 1156.97424    0.40
IG_Identification     0.57283    0.04517 1156.91622   12.68
OG_Bonds             -0.03595    0.01972 1154.35863   -1.82
Empathic_concern      0.15197    0.07095 1156.42051    2.14
Perspective_taking    0.32743    0.05657 1153.76338    5.79
Age                  -0.00821    0.00501 1156.95280   -1.64
Female                0.05017    0.07614 1143.29711    0.66
Married               0.07075    0.09464 1147.90939    0.75
Wealth_level          0.00072    0.00182 1140.69680    0.40
                               Pr(>|t|)    
(Intercept)                       0.481    
IG_Fusion                         0.689    
IG_Identification  < 0.0000000000000002 ***
OG_Bonds                          0.068 .  
Empathic_concern                  0.032 *  
Perspective_taking         0.0000000091 ***
Age                               0.102    
Female                            0.510    
Married                           0.455    
Wealth_level                      0.693    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) IG_Fsn IG_Idn OG_Bnd Empth_ Prspc_ Age    Female Marrid
IG_Fusion   -0.002                                                        
IG_Idntfctn  0.072 -0.876                                                 
OG_Bonds    -0.132 -0.152  0.178                                          
Empthc_cncr -0.662  0.001 -0.085 -0.012                                   
Prspctv_tkn -0.383 -0.095 -0.150 -0.044 -0.151                            
Age         -0.391  0.016 -0.002 -0.025  0.031  0.002                     
Female      -0.113 -0.015  0.002  0.070  0.018  0.006  0.025              
Married      0.152 -0.031  0.051 -0.026 -0.037  0.022 -0.546  0.056       
Wealth_levl -0.236  0.048 -0.074 -0.060  0.003  0.014  0.056  0.000 -0.072

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.019
		(0.047)
IG_Identification		0.570***
		(0.045)
OG_Bonds		-0.036
		(0.020)
Empathic_concern		0.150*
		(0.071)
Perspective_taking		0.330***
		(0.057)
Age		-0.008
		(0.005)
Female		0.050
		(0.076)
Married		0.071
		(0.095)
Wealth_level		0.001
		(0.002)
Constant	4.900***	0.260
	(0.130)	(0.360)
Observations	1,836	1,167
Log Likelihood	-3,727.000	-1,920.000
Akaike Inf. Crit.	7,460.000	3,865.000
Bayesian Inf. Crit.	7,476.000	3,925.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.7     5.3     4.9     6.3     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.90", 
                 xintercept = 4.90, 
                 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 = 2) +
  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        Ghana         5.3
3       Malawi         5.2
4     Pakistan         4.7
5 Sierra Leone         4.6
6     Tanzania         4.6
7       Uganda         5.3
8          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.90", 
                 xintercept = 4.90, 
                 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 = 2) +
  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 = 2) +
  theme_bw()

Faceted plots: Endorse BBL vs Ingroup identification

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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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.046		0.019
		(0.049)		(0.047)
IG_Identification		0.440***		0.570***
		(0.047)		(0.045)
OG_Bonds		-0.017		-0.036
		(0.021)		(0.020)
Empathic_concern		0.140		0.150*
		(0.074)		(0.071)
Perspective_taking		0.390***		0.330***
		(0.059)		(0.057)
Age		-0.002		-0.008
		(0.005)		(0.005)
Female		0.093		0.050
		(0.079)		(0.076)
Married		0.030		0.071
		(0.099)		(0.095)
Wealth_level		0.002		0.001
		(0.002)		(0.002)
Constant	4.700***	0.160	4.900***	0.260
	(0.098)	(0.380)	(0.130)	(0.360)
Observations	1,832	1,166	1,836	1,167
Log Likelihood	-3,635.000	-1,966.000	-3,727.000	-1,920.000
Akaike Inf. Crit.	7,276.000	3,957.000	7,460.000	3,865.000
Bayesian Inf. Crit.	7,293.000	4,017.000	7,476.000	3,925.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(formula = 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: 7673

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.054 -0.707  0.240  0.861  1.277 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0288   0.17    
 Residual             3.6191   1.90    
Number of obs: 1858, groups:  Country, 8

Fixed effects:
            Estimate Std. Error     df t value     Pr(>|t|)    
(Intercept)   4.7589     0.0772 5.5543    61.6 0.0000000043 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random intercept model

Display code

## Varying intercept models with individual-level predictors:

m05 <- lmer(formula = IG_Fusion~event_positive_affect+event_negative_affect+
              event_episodic_recall+event_shared_perception+event_event_reflection+
              event_transformative_indiv+event_transformative_group+Age+Female+
              Married+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_event_reflection + event_transformative_indiv +  
    event_transformative_group + Age + Female + Married + Wealth_level +  
    (1 | Country)
   Data: ds

REML criterion at convergence: 5553

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.497 -0.548  0.016  0.578  4.966 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.00313  0.0559  
 Residual             1.29803  1.1393  
Number of obs: 1770, groups:  Country, 8

Fixed effects:
                              Estimate  Std. Error          df t value
(Intercept)                   0.437327    0.143822  539.363698    3.04
event_positive_affect         0.093461    0.016919 1755.998285    5.52
event_negative_affect         0.003450    0.013394 1733.038854    0.26
event_episodic_recall         0.216597    0.024896 1753.974017    8.70
event_shared_perception       0.191615    0.023995 1757.994771    7.99
event_event_reflection        0.177442    0.026735 1756.554396    6.64
event_transformative_indiv    0.189043    0.024945 1757.919218    7.58
event_transformative_group    0.070452    0.024412 1757.564340    2.89
Age                          -0.001229    0.003015 1020.892466   -0.41
Female                       -0.023615    0.056514 1700.580254   -0.42
Married                       0.014442    0.068679 1364.358175    0.21
Wealth_level                 -0.000826    0.001326 1463.451683   -0.62
                                       Pr(>|t|)    
(Intercept)                              0.0025 ** 
event_positive_affect        0.0000000381181644 ***
event_negative_affect                    0.7968    
event_episodic_recall      < 0.0000000000000002 ***
event_shared_perception      0.0000000000000025 ***
event_event_reflection       0.0000000000425074 ***
event_transformative_indiv   0.0000000000000563 ***
event_transformative_group               0.0040 ** 
Age                                      0.6837    
Female                                   0.6761    
Married                                  0.8335    
Wealth_level                             0.5332    
---
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_v_
evnt_pstv_f      -0.221                                            
evnt_ngtv_f      -0.336  0.419                                     
evnt_psdc_r      -0.015 -0.147     0.020                           
evnt_shrd_p      -0.072  0.002    -0.080  -0.293                   
evnt_vnt_rf      -0.040 -0.089    -0.134  -0.288    -0.139         
evnt_trnsfrmtv_n -0.042 -0.143    -0.044  -0.246    -0.064  -0.309 
evnt_trnsfrmtv_g -0.062 -0.219    -0.033  -0.061    -0.346  -0.143 
Age              -0.578  0.021     0.007  -0.022     0.047  -0.007 
Female           -0.122 -0.009     0.026   0.003    -0.026   0.031 
Married           0.191 -0.033    -0.068   0.012    -0.039   0.021 
Wealth_levl      -0.401  0.026    -0.006  -0.034    -0.025  -0.023 
                 evnt_trnsfrmtv_n evnt_trnsfrmtv_g Age    Female Marrid
evnt_pstv_f                                                            
evnt_ngtv_f                                                            
evnt_psdc_r                                                            
evnt_shrd_p                                                            
evnt_vnt_rf                                                            
evnt_trnsfrmtv_n                                                       
evnt_trnsfrmtv_g -0.201                                                
Age              -0.018            0.011                               
Female           -0.068            0.056           -0.066              
Married           0.042           -0.001           -0.490  0.082       
Wealth_levl       0.009            0.015            0.035 -0.001 -0.080

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.093***
		(0.017)
event_negative_affect		0.003
		(0.013)
event_episodic_recall		0.220***
		(0.025)
event_shared_perception		0.190***
		(0.024)
event_event_reflection		0.180***
		(0.027)
event_transformative_indiv		0.190***
		(0.025)
event_transformative_group		0.070**
		(0.024)
Age		-0.001
		(0.003)
Female		-0.024
		(0.057)
Married		0.014
		(0.069)
Wealth_level		-0.001
		(0.001)
Constant	4.800***	0.440**
	(0.077)	(0.140)
Observations	1,858	1,770
Log Likelihood	-3,836.000	-2,776.000
Akaike Inf. Crit.	7,679.000	5,581.000
Bayesian Inf. Crit.	7,695.000	5,657.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     3.3     5.3     4.8     6.3     7.0 

Display code

ggplot(data = df01, 
       aes(x = IG_Fusion)) +
  geom_histogram(color = "black",
                 bins = 20)+
  xlim(1, 7)+
  geom_textvline(label = "Mean = 4.80", 
                 xintercept = 4.80, 
                 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 = 2) +
  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        Ghana       4.9
3       Malawi       4.9
4     Pakistan       4.7
5 Sierra Leone       4.4
6     Tanzania       4.5
7       Uganda       5.0
8          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 = 4.80", 
                 xintercept = 4.80, 
                 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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  theme_bw()

Faceted plots: Ingroup fusion vs Reflection of event

Display code

ggplot(data = ds, 
       aes(y = IG_Fusion,
           x = event_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 = 2) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_Fusion,
           x = event_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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  theme_bw()

Section 5. Outcome: Ingroup identification

Unconditional means model

Display code

## Varying intercept model with no predictors:

m06 <- lmer(formula = 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: 7862

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.096 -0.666  0.295  0.861  1.194 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0715   0.267   
 Residual             3.9989   2.000   
Number of obs: 1858, groups:  Country, 8

Fixed effects:
            Estimate Std. Error    df t value      Pr(>|t|)    
(Intercept)    4.905      0.109 6.865    45.2 0.00000000094 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random intercept model

Display code

## Varying intercept models with individual-level predictors:

m07 <- lmer(formula = IG_Identification~event_positive_affect+event_negative_affect+
              event_episodic_recall+event_shared_perception+event_event_reflection+
              event_transformative_indiv+event_transformative_group+Age+Female+
              Married+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_event_reflection +  
    event_transformative_indiv + event_transformative_group +  
    Age + Female + Married + Wealth_level + (1 | Country)
   Data: ds

REML criterion at convergence: 5565

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.178 -0.547  0.005  0.582  5.015 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.00     0.00    
 Residual             1.31     1.15    
Number of obs: 1769, groups:  Country, 8

Fixed effects:
                               Estimate   Std. Error           df t value
(Intercept)                   0.5236078    0.1418265 1757.0000000    3.69
event_positive_affect         0.0638384    0.0169954 1757.0000000    3.76
event_negative_affect        -0.0396379    0.0134238 1757.0000000   -2.95
event_episodic_recall         0.3083694    0.0249813 1757.0000000   12.34
event_shared_perception       0.2376408    0.0240968 1757.0000000    9.86
event_event_reflection        0.1380038    0.0268873 1757.0000000    5.13
event_transformative_indiv    0.2090299    0.0250857 1757.0000000    8.33
event_transformative_group    0.0343042    0.0245497 1757.0000000    1.40
Age                          -0.0041910    0.0029920 1757.0000000   -1.40
Female                        0.0442387    0.0565810 1757.0000000    0.78
Married                      -0.0811617    0.0683785 1757.0000000   -1.19
Wealth_level                 -0.0000456    0.0013210 1757.0000000   -0.03
                                       Pr(>|t|)    
(Intercept)                             0.00023 ***
event_positive_affect                   0.00018 ***
event_negative_affect                   0.00319 ** 
event_episodic_recall      < 0.0000000000000002 ***
event_shared_perception    < 0.0000000000000002 ***
event_event_reflection               0.00000032 ***
event_transformative_indiv < 0.0000000000000002 ***
event_transformative_group              0.16249    
Age                                     0.16147    
Female                                  0.43440    
Married                                 0.23541    
Wealth_level                            0.97247    
---
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_v_
evnt_pstv_f      -0.219                                            
evnt_ngtv_f      -0.338  0.418                                     
evnt_psdc_r      -0.018 -0.145     0.025                           
evnt_shrd_p      -0.073  0.002    -0.079  -0.295                   
evnt_vnt_rf      -0.041 -0.089    -0.134  -0.289    -0.139         
evnt_trnsfrmtv_n -0.045 -0.141    -0.043  -0.247    -0.065  -0.310 
evnt_trnsfrmtv_g -0.061 -0.222    -0.035  -0.059    -0.346  -0.143 
Age              -0.579  0.017     0.002  -0.020     0.046  -0.005 
Female           -0.124 -0.009     0.029   0.001    -0.024   0.030 
Married           0.199 -0.037    -0.074   0.018    -0.036   0.015 
Wealth_levl      -0.409  0.023    -0.011  -0.028    -0.024  -0.027 
                 evnt_trnsfrmtv_n evnt_trnsfrmtv_g Age    Female Marrid
evnt_pstv_f                                                            
evnt_ngtv_f                                                            
evnt_psdc_r                                                            
evnt_shrd_p                                                            
evnt_vnt_rf                                                            
evnt_trnsfrmtv_n                                                       
evnt_trnsfrmtv_g -0.201                                                
Age              -0.018            0.011                               
Female           -0.068            0.057           -0.064              
Married           0.047           -0.004           -0.496  0.081       
Wealth_levl       0.012            0.013            0.037  0.001 -0.096
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

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.064***
		(0.017)
event_negative_affect		-0.040**
		(0.013)
event_episodic_recall		0.310***
		(0.025)
event_shared_perception		0.240***
		(0.024)
event_event_reflection		0.140***
		(0.027)
event_transformative_indiv		0.210***
		(0.025)
event_transformative_group		0.034
		(0.025)
Age		-0.004
		(0.003)
Female		0.044
		(0.057)
Married		-0.081
		(0.068)
Wealth_level		-0.00005
		(0.001)
Constant	4.900***	0.520***
	(0.110)	(0.140)
Observations	1,858	1,769
Log Likelihood	-3,931.000	-2,783.000
Akaike Inf. Crit.	7,868.000	5,593.000
Bayesian Inf. Crit.	7,885.000	5,670.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     3.3     5.7     4.9     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 = 4.90", 
                 xintercept = 4.90, 
                 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 = 2) +
  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        Ghana               5.1
3       Malawi               5.1
4     Pakistan               4.7
5 Sierra Leone               4.6
6     Tanzania               4.6
7       Uganda               5.2
8          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 = 4.90", 
                 xintercept = 4.90, 
                 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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  theme_bw()

Faceted plots: Ingroup identification vs Reflection of event

Display code

ggplot(data = ds, 
       aes(y = IG_Identification,
           x = event_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 = 2) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = IG_Identification,
           x = event_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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  theme_bw()

Section 6. Outcome: Outgroup bonds

Unconditional means model

Display code

## Varying intercept model with no predictors:

m08 <- lmer(formula = 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: 5170

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.410 -0.872 -0.124  0.707  2.213 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.143    0.378   
 Residual             3.457    1.859   
Number of obs: 1264, groups:  Country, 7

Fixed effects:
            Estimate Std. Error    df t value   Pr(>|t|)    
(Intercept)    3.206      0.157 5.968    20.5 0.00000094 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random intercept model

Display code

## Varying intercept models with individual-level predictors:

m09 <- lmer(formula = OG_Bonds~event_positive_affect+event_negative_affect+
              event_episodic_recall+event_shared_perception+event_event_reflection+
              event_transformative_indiv+event_transformative_group+Age+Female+
              Married+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_event_reflection + event_transformative_indiv +  
    event_transformative_group + Age + Female + Married + Wealth_level +  
    (1 | Country)
   Data: ds

REML criterion at convergence: 4872

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.017 -0.798 -0.141  0.719  2.584 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.0631   0.251   
 Residual             3.2169   1.794   
Number of obs: 1202, groups:  Country, 7

Fixed effects:
                             Estimate Std. Error         df t value
(Intercept)                   1.96903    0.31396  283.09844    6.27
event_positive_affect         0.12053    0.03207 1187.85090    3.76
event_negative_affect         0.20672    0.02575 1189.90357    8.03
event_episodic_recall        -0.23331    0.04855 1189.76587   -4.81
event_shared_perception      -0.02160    0.04596 1187.60700   -0.47
event_event_reflection        0.05813    0.05220 1187.39683    1.11
event_transformative_indiv    0.03973    0.04876 1188.23402    0.81
event_transformative_group    0.09274    0.04766 1188.67470    1.95
Age                           0.00501    0.00722 1189.44892    0.69
Female                       -0.29844    0.11020 1182.05445   -2.71
Married                       0.08372    0.13798 1187.12864    0.61
Wealth_level                  0.00573    0.00261 1182.63474    2.19
                                     Pr(>|t|)    
(Intercept)                0.0000000013316381 ***
event_positive_affect                 0.00018 ***
event_negative_affect      0.0000000000000024 ***
event_episodic_recall      0.0000017373880534 ***
event_shared_perception               0.63846    
event_event_reflection                0.26569    
event_transformative_indiv            0.41531    
event_transformative_group            0.05190 .  
Age                                   0.48801    
Female                                0.00687 ** 
Married                               0.54415    
Wealth_level                          0.02853 *  
---
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_v_
evnt_pstv_f      -0.216                                            
evnt_ngtv_f      -0.308  0.430                                     
evnt_psdc_r      -0.015 -0.152    -0.008                           
evnt_shrd_p      -0.049 -0.008    -0.095  -0.282                   
evnt_vnt_rf      -0.060 -0.073    -0.113  -0.282    -0.175         
evnt_trnsfrmtv_n -0.026 -0.146    -0.017  -0.243    -0.037  -0.310 
evnt_trnsfrmtv_g -0.039 -0.196    -0.043  -0.101    -0.314  -0.127 
Age              -0.628  0.044     0.021   0.007     0.027   0.012 
Female           -0.156  0.017     0.032   0.029    -0.061   0.041 
Married           0.252 -0.035    -0.076  -0.011    -0.025   0.010 
Wealth_levl      -0.383  0.026     0.009  -0.038    -0.045  -0.005 
                 evnt_trnsfrmtv_n evnt_trnsfrmtv_g Age    Female Marrid
evnt_pstv_f                                                            
evnt_ngtv_f                                                            
evnt_psdc_r                                                            
evnt_shrd_p                                                            
evnt_vnt_rf                                                            
evnt_trnsfrmtv_n                                                       
evnt_trnsfrmtv_g -0.236                                                
Age              -0.040           -0.001                               
Female           -0.076            0.032            0.029              
Married           0.051            0.014           -0.554  0.053       
Wealth_levl       0.015            0.008            0.048  0.018 -0.056

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.120***
		(0.032)
event_negative_affect		0.210***
		(0.026)
event_episodic_recall		-0.230***
		(0.049)
event_shared_perception		-0.022
		(0.046)
event_event_reflection		0.058
		(0.052)
event_transformative_indiv		0.040
		(0.049)
event_transformative_group		0.093
		(0.048)
Age		0.005
		(0.007)
Female		-0.300**
		(0.110)
Married		0.084
		(0.140)
Wealth_level		0.006*
		(0.003)
Constant	3.200***	2.000***
	(0.160)	(0.310)
Observations	1,264	1,202
Log Likelihood	-2,585.000	-2,436.000
Akaike Inf. Crit.	5,176.000	4,900.000
Bayesian Inf. Crit.	5,191.000	4,972.000
Note:	*p<0.05; **p<0.01; ***p<0.001

Histogram: Outgroup bonds

Display code

df01 <- ds %>% drop_na(OG_Bonds)

summary(df01$OG_Bonds)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0     1.5     3.0     3.3     4.8     7.0 

Display code

ggplot(data = df01, 
       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 = 2) +
  theme_bw()

Display code

tbl02 <- aggregate(df01$OG_Bonds, 
                    by=list(df01$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        Ghana      2.9
3       Malawi      2.8
4     Pakistan      3.6
5 Sierra Leone      2.8
6     Tanzania      3.6
7       Uganda      2.9

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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  theme_bw()

Faceted plots: Outgroup bonds vs Reflection of event

Display code

ggplot(data = ds, 
       aes(y = OG_Bonds,
           x = event_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 = 2) +
  theme_bw()

Display code

ggplot(data = ds, 
       aes(y = OG_Bonds,
           x = event_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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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 = 2) +
  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.093***		0.064***		0.120***
		(0.017)		(0.017)		(0.032)
event_negative_affect		0.003		-0.040**		0.210***
		(0.013)		(0.013)		(0.026)
event_episodic_recall		0.220***		0.310***		-0.230***
		(0.025)		(0.025)		(0.049)
event_shared_perception		0.190***		0.240***		-0.022
		(0.024)		(0.024)		(0.046)
event_event_reflection		0.180***		0.140***		0.058
		(0.027)		(0.027)		(0.052)
event_transformative_indiv		0.190***		0.210***		0.040
		(0.025)		(0.025)		(0.049)
event_transformative_group		0.070**		0.034		0.093
		(0.024)		(0.025)		(0.048)
Age		-0.001		-0.004		0.005
		(0.003)		(0.003)		(0.007)
Female		-0.024		0.044		-0.300**
		(0.057)		(0.057)		(0.110)
Married		0.014		-0.081		0.084
		(0.069)		(0.068)		(0.140)
Wealth_level		-0.001		-0.00005		0.006*
		(0.001)		(0.001)		(0.003)
Constant	4.800***	0.440**	4.900***	0.520***	3.200***	2.000***
	(0.077)	(0.140)	(0.110)	(0.140)	(0.160)	(0.310)
Observations	1,858	1,770	1,858	1,769	1,264	1,202
Log Likelihood	-3,836.000	-2,776.000	-3,931.000	-2,783.000	-2,585.000	-2,436.000
Akaike Inf. Crit.	7,679.000	5,581.000	7,868.000	5,593.000	5,176.000	4,900.000
Bayesian Inf. Crit.	7,695.000	5,657.000	7,885.000	5,670.000	5,191.000	4,972.000
Note:	*p<0.05; **p<0.01; ***p<0.001