Data Analysis and Visualization
Tanzania Community Data

Section 1. Demographics

Display code

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

## Install "pacman" package if not installed
# (remove the # symbol from the line below):
# install.packages("pacman")

## Load R packages:
pacman::p_load(data.table, tidyverse, haven, labelled, vtable, 
               psych, scales, weights, clipr, forcats,
               stargazer, ggthemes, ggcharts, geomtextpath,
               corrplot, tm, readxl, arsenal)

## Import dataset:

#ds <- read_excel("~/Desktop/oxford/data/tanzania/Tanzania Data Entry.xlsx")

ds <- fread("~/Desktop/oxford/data/tanzania/Tanzania Data Entry - Sheet1.csv")

Variable: Age

Display code

ds$Age1 <- as.numeric(gsub("\\D", "", ds$Age))
ds$age <- ifelse(ds$Age1 %in% 18:90, ds$Age1, NA)

summary(ds$age)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
     18      27      34      36      43      84      19 

Display code

ds %>% drop_na(age)%>%
ggplot(aes(x = age))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 36.00", 
                 xintercept = 36.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Age", 
       y = "Frequency", 
       title = "Age distribution: Tanzania")+
  theme_bw()

Variable: Gender

Display code

ds$gender <- ifelse(ds$Gender == "Male", "Male",
             ifelse(ds$Gender == "Female", "Female", NA))

lp02 <- ds %>% drop_na(gender) %>%
lollipop_chart(x = gender,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Gender distribution: Tanzania")+
  theme_bw()

lp02

Display code

bp01 <- ds %>% drop_na(gender, age) %>% 
  ggplot(aes(y = age, 
             x = gender))+
geom_boxplot(fill = "grey")+
  labs(y = "Age",
       x = "",
       title = "Age distribution by gender: Tanzania")+
  coord_flip()+
  theme_bw()

bp01

Variable: Province

Display code

ds$Province <- ifelse(ds$Province == "DSM", "Dar Es Salaam",
               ifelse(ds$Province == "2422 DSM", "Dar Es Salaam",
               ifelse(ds$Province == "Morogoro Mjini", "Morogoro", ds$Province)))

prov <- as.data.frame(table(ds$Province))
prov$Var1 <- as.character(prov$Var1)
prov$province <- ifelse(prov$Freq < 3, "Other", prov$Var1)
l1 <- as.list(prov$province)

ds2 <- ds[, c("Province")]

ds2$province <- ifelse(ds2$Province %in% l1, ds2$Province, "Other")
ds2$province <- ifelse(ds2$province == "", NA, ds2$province)
ds2$province <- as.character(ds2$province)

lp03 <- ds2 %>% drop_na(province) %>%
lollipop_chart(x = province,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Province distribution: Tanzania")+
  theme_bw()

lp03

Variable: Household wealth

Display code

ds$hh_wealth <- haven::as_factor(ds$Household_level_of_wealth)

## Convert to sentence case:
ds$hh_wealth <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$hh_wealth, perl = TRUE)

ds$hh_wealth <- factor(ds$hh_wealth, 
                       levels=c("Much lower", "Slightly lower", "Average",
                                "Slightly higher", "Much higher"))

table(ds$hh_wealth)

     Much lower  Slightly lower         Average Slightly higher     Much higher 
             33             119             154              21               4 

Display code

lp04 <- ds %>% drop_na(hh_wealth) %>%
lollipop_chart(x = hh_wealth,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Household wealth distribution: Tanzania")+
  theme_bw()

lp04

Variable: Socioeconomic status

Display code

ds$ses1 <- haven::as_factor(ds$Socioeconomic_status)

## Convert to sentence case:
ds$ses1 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ses1, perl = TRUE)

ds$ses <- factor(ds$ses1, 
                 levels=c("Lower middle", "Middle", "Upper middle", "Upper"))

table(ds$ses)

Lower middle       Middle Upper middle        Upper 
          96          201           28            4 

Display code

lp05 <- ds %>% drop_na(ses) %>%
lollipop_chart(x = ses,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Socioeconomic status: Tanzania")+
  theme_bw()

lp05

Variable: Annual income

Display code

ds$income <- as.numeric(gsub("\\D", "", ds$Annual_income))

## Remove outliers (greater than 5 million):
# ds <- ds[ds$income < 5000000, ]

summary(ds$income)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
       0  1177500  2510000  4110781  5000000 50000000       60 

Display code

ds %>% drop_na(income)%>%
ggplot(aes(x = income))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 80)+
  geom_textvline(label = "Mean = 4.11M", 
                 xintercept = 4110781, 
                 vjust = 1.3, 
                 lwd = 1.05, 
                 linetype = 2)+
  scale_x_continuous(labels = scales::unit_format(unit = "M", 
                                                  scale = 1e-6))+
  labs(x = "Annual income", 
       y = "Frequency", 
       title = "Annual income distribution: Tanzania")+
  theme_bw()

Correlation between income, wealth, and socioeconomic status

Display code

## HH Wealth:
ds$hh_wealth <- haven::as_factor(ds$Household_level_of_wealth)

## Convert to sentence case:
ds$hh_wealth <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$hh_wealth, perl = TRUE)

ds$hh_wealth <- factor(ds$hh_wealth, 
                       levels=c("Much lower", "Slightly lower", "Average",
                                "Slightly higher", "Much higher"))

table(ds$hh_wealth)

     Much lower  Slightly lower         Average Slightly higher     Much higher 
             33             119             154              21               4 

Display code

## SES:
ds$ses1 <- haven::as_factor(ds$Socioeconomic_status)

## Convert to sentence case:
ds$ses1 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ses1, perl = TRUE)

ds$ses <- factor(ds$ses1, 
                 levels=c("Lower middle", "Middle", "Upper middle", "Upper"))

ds$fses <- as.numeric(ds$ses)
ds$fhhw <- as.numeric(ds$hh_wealth)

ds3 <- cbind.data.frame(ds$fses, ds$fhhw, ds$income)
ds3 <- na.omit(ds3)

ds3$Annual_income <- ds3$`ds$income`
ds3$Household_wealth <- ds3$`ds$fhhw`
ds3$SES <- ds3$`ds$fses`

mtx <- cor(ds3[, c(4:6)])

corrplot(mtx, method = "number", number.cex = 1.0,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Variable: Nature of employment

Display code

ds$jobnature <- haven::as_factor(ds$Job_Nature)
ds$jobnature <- as.character(ds$jobnature)

ds$jobnature <- ifelse(ds$jobnature == "", NA, ds$jobnature)

lp06 <- ds %>% drop_na(jobnature) %>%
lollipop_chart(x = jobnature,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Nature of employment: Tanzania")+
  theme_bw()

lp06

Variable: Religious affiliation

Display code

ds$r1 <- ds$Religious_Affiliation

ds$religion <- ifelse(ds$r1 == "Christian – Catholic", "Christian: Catholic",
               ifelse(ds$r1 == "Christian – Other", "Christian: Other",
               ifelse(ds$r1 == "Christian- Pentecost", "Christian: Pentecostal",
               ifelse(ds$r1 == "Christian – Protestant", "Christian: Protestant",
               ifelse(ds$r1 == "Muslim – Sunni", "Muslim: Sunni",
               ifelse(ds$r1 == "Muslim – Shia", "Muslim: Shia",
               ifelse(ds$r1 == "Spiritual not Religious", "Spiritual not religious",
               ifelse(ds$r1 == "Other", "Other", NA))))))))

table(ds$religion)

    Christian: Catholic        Christian: Other  Christian: Pentecostal 
                    145                      37                      33 
  Christian: Protestant            Muslim: Shia           Muslim: Sunni 
                     28                      42                      30 
                  Other Spiritual not religious 
                      5                      14 

Display code

lp07 <- ds %>% drop_na(religion) %>%
lollipop_chart(x = religion,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Religious affiliation: Tanzania")+
  theme_bw()

lp07

Variable: Marital status

Display code

ds$married <- ifelse(ds$Marital_status == "Married", "Married",
              ifelse(ds$Marital_status == "Single", "Unmarried",
              ifelse(ds$Marital_status == "In a relationship (not married)", "Unmarried",
                     "Other")))

table(ds$married)

  Married     Other Unmarried 
      163        58       131 

Display code

bp01 <- ds %>% drop_na(married) %>%
  ggplot(aes(x = fct_rev(fct_infreq(married))))+
  geom_bar(fill = "black")+
  labs(x = "",
       y = "Frequency",
       title = "Marital status")+
  coord_flip()+
  theme_bw()

# bp01

lp08 <- ds %>% drop_na(married) %>%
lollipop_chart(x = married,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Marital status: Tanzania")+
  theme_bw()

lp08

Variable: Number of children

Display code

ds$children <- as.numeric(gsub("\\D", "", ds$Number_of_Children))

summary(ds$children)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    0.0     0.0     2.0     2.4     4.0    17.0      20 

Display code

ds %>% drop_na(children)%>%
ggplot(aes(x = children))+
  geom_bar(color = "black",
                 fill = "gray",
                 width = 0.75)+
  geom_textvline(label = "Mean = 2.40", 
                 xintercept = 2.40, 
                 vjust = 1.3, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Number of children", 
       y = "Frequency", 
       title = "Number of children: Tanzania")+
  theme_bw()

Variable: Ethnicity

Note: I did my best to clean up the variable, but the categories need to be checked by an expert

Display code

eth <- as.data.frame(table(ds$Ethnicity))
eth$Var1 <- as.character(eth$Var1)
eth$ethnicity <- ifelse(eth$Freq < 2, "Other", eth$Var1)
l1 <- as.list(eth$ethnicity)

ds2 <- ds[, c("Ethnicity")]

ds2$ethnicity <- ifelse(ds2$Ethnicity %in% l1, ds2$Ethnicity, "Other")
ds2$ethnicity <- ifelse(ds2$ethnicity == "Masai", "Maasai",
                 ifelse(ds2$ethnicity == "", NA, ds2$ethnicity))

table(ds2$ethnicity)

    Chaga      Haya    Maasai     Mbena    Mchaga     Mdigo     Mfipa     Mgogo 
        2         2       125         3         9         4         3         8 
    Mhaya     Mhehe     Mjita   Mkaguru   Mlugulu   Mluguru  Mmakonde    Mngoni 
        2         5         2         2        11        15         3         8 
 Mnyakywa Mnyamwezi  Mnyaturu Mnyiramba     Mpare   Mpogolo   Mpogoro    Mrangi 
        2         4         6         2         9         3         6         2 
  Msambaa   Msukuma      Muha Muhangaza    Muhaya      Muwa   Mzanaki   Mzaramo 
       11        10         7         2         2         2         2         3 
   Mzigua  Nyakyusa     Other 
        7         2        51 

Display code

lp08 <- ds2 %>% drop_na(ethnicity) %>%
lollipop_chart(x = ethnicity,
               line_color = "black",
               point_color = "black")+
  labs(y = "Frequency",
       x = "",
       title = "Ethnic distribution: Tanzania")+
  theme_bw()

lp08

Variable: Education

Note: This needs a lot of work from someone familiar with Tanzanian education. The spreadsheet to edit can be found here:

Spreadsheet for cleaning: Tanzania

Display code

# ds$id <- ds$`180`
# ds5 <- ds[, c("id", "education")]
# fwrite(ds5, file = "~/Desktop/Tanzania_for_cleaning.csv")

ds$education <- as.character(haven::as_factor(ds$Level_of_education))
table(ds$education)

                              4           4m  4               7   Advance level 
             23               8               1               7               1 
Bachelor Degree     Certificate    Certtificate           Cheti            Chuo 
              1               2               1               6               6 
    Chuo Ckikuu      Chuo kikuu      Chuo Kikuu     Claiasa nne         College 
              1              11               3               1               1 
       Darasa 7        Darasa 8     Darasa la 7  Darasa la Saba      Darsa la 7 
              8               1               1               1               2 
     Darsa La 7   Darsa la saba          Degree  Degree-Shahada   Didato channe 
              3               1              15               1               1 
  Digri/Shahada         Diploma    Elimu Msingi          Form 4          Form 5 
              1              15               1              14               1 
         Form 6       Form Five     Form Five 5       Form Four     Form Four 4 
              2               2               2              73               3 
    Form Four 5         Form IV      Form One 1        Form Six      Form Three 
              1               1               1              14               1 
       Form Two      Form Two 2              IV        Kidato 4    Kidato cha 1 
              6               1               1               3               1 
   Kidato cha 4    Kidato cha 6  Kidato cha nne Kidato Cha seta   Kidato Channe 
              1               2               2               1               1 
  Kidato chaune            La 4            La 7            LA 7          La nne 
              1               3               7               8               4 
     La nne (4)  La nne mkoloni         La saba       La saba 7   Lanne Mkoloni 
              1               1              17               8               2 
         Lasaba          Msingi         O Level   Post Graduate         Primary 
              1               3               1               1               1 
      Secondary        Sekonorw         Shadada         Shahada     Shahada yai 
              6               1               3               5               1 
    Shule Msiry Shule ya msingi        Sokondal      StaShahada         Std VII 
              1               1               1               2               3 
        STD VII         Uzamili              VI         VII Std 
              8               1               1               1 

Section 2. Factor Analysis: Group fusion/identification (Version 1)

Display code

## Four ingroup fusion items:

# I have a deep emotional bond with the [ingroup].
ds$IGF01 <- as.numeric(ds$`Q11.1 Group Bonds`)

# I am strong because of the [ingroup].
ds$IGF02 <- as.numeric(ds$`Q11.2 Group Bonds`)

# I make the [ingroup] strong.  
ds$IGF03 <- as.numeric(ds$`Q11.3 Group Bonds`)

# I am one with the [ingroup]
ds$IGF04 <- as.numeric(ds$`Q11.4 Group Bonds`)


## Four outgroup fusion items:

# I have a deep emotional bond with the [outgroup].
ds$OGF01 <- as.numeric(ds$`Q11.5 Group Bonds`)

# I am strong because of the [outgroup].
ds$OGF02 <- as.numeric(ds$`Q11.6 Group Bonds`)

# I make the [outgroup] strong. 
ds$OGF03 <- as.numeric(ds$`Q11.7 Group Bonds`)

# I am one with the [outgroup].
ds$OGF04 <- as.numeric(ds$`Q11.8 Group Bonds`)


## Four ingroup identification items:

# I identify with the [ingroup].
ds$IGI01 <- as.numeric(ds$`Q11.9 Group Bonds`)

# I have a lot in common with the [ingroup].
ds$IGI02 <- as.numeric(ds$`Q11.10 Group Bonds`)

# I connect with the values of the [ingroup].
ds$IGI03 <- as.numeric(ds$`Q11.11 Group Bonds`)

# I feel a sense of belonging with the [ingroup].
ds$IGI04 <- as.numeric(ds$`Q11.12 Group Bonds`)


## Four outgroup identification items:

# I identify with the [outgroup].   
ds$OGI01 <- as.numeric(ds$`Q11.13 Group Bonds`)

# I have a lot in common with the [outgroup].   
ds$OGI02 <- as.numeric(ds$`Q11.14 Group Bonds`)

# I connect with the values of the [outgroup].  
ds$OGI03 <- as.numeric(ds$`Q11.15 Group Bonds`)

# I feel a sense of belonging with the [outgroup].  
ds$OGI04 <- as.numeric(ds$`Q11.16 Group Bonds`)

Display code

## Bonds dataframe:
bonds <- cbind.data.frame(ds$IGF01, ds$IGF02, ds$IGF03, ds$IGF04,
                          ds$IGI01, ds$IGI02, ds$IGI03, ds$IGI04,
                          ds$OGF01, ds$OGF02, ds$OGF03, ds$OGF04,
                          ds$OGI01, ds$OGI02, ds$OGI03, ds$OGI04)

names(bonds) <- sub('ds\\$', '', names(bonds))

bonds <- na.omit(bonds)
mtx1 <- cor(bonds[, c(1:16)])

Shorthand for item names:
IG = Ingroup, OG = Outgroup
F = Fusion, I = Identification
So for example:
IGF01 = “ingroup fusion: item 1”
OGI04 = “outgroup identification: item 4”, and so on

Display code

corrplot(mtx1, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

The correlation plot above shows the correlation coefficients between the sixteen different items (eight for ingroup fusion/identification and eight for outgroup fusion/identification).

KMO and Bartlett’s test

Display code

KMO(r=cor(bonds))

Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = cor(bonds))
Overall MSA =  0.87
MSA for each item = 
IGF01 IGF02 IGF03 IGF04 IGI01 IGI02 IGI03 IGI04 OGF01 OGF02 OGF03 OGF04 OGI01 
 0.79  0.77  0.82  0.83  0.89  0.78  0.87  0.76  0.87  0.91  0.88  0.91  0.91 
OGI02 OGI03 OGI04 
 0.88  0.88  0.88 

Display code

cortest.bartlett(bonds)

$chisq
[1] 1987

$p.value
[1] 0

$df
[1] 120

Parallel test

Display code

parallel <- fa.parallel(bonds)

Parallel analysis suggests that the number of factors =  2  and the number of components =  2 

Scree plot suggests a two factor structure. We will proceed with promax rotation, which assumes that the items are inter-correlated (that is, not independent from each other).

Two factor model

Display code

fit02 <- factanal(bonds, 2, rotation="promax")
fit02

Call:
factanal(x = bonds, factors = 2, rotation = "promax")

Uniquenesses:
IGF01 IGF02 IGF03 IGF04 IGI01 IGI02 IGI03 IGI04 OGF01 OGF02 OGF03 OGF04 OGI01 
 0.56  0.56  0.70  0.65  0.68  0.85  0.77  0.81  0.39  0.52  0.45  0.58  0.32 
OGI02 OGI03 OGI04 
 0.32  0.35  0.31 

Loadings:
      Factor1 Factor2
IGF01          0.694 
IGF02  0.134   0.705 
IGF03          0.537 
IGF04          0.591 
IGI01          0.534 
IGI02          0.362 
IGI03 -0.149   0.404 
IGI04 -0.107   0.383 
OGF01  0.776         
OGF02  0.721         
OGF03  0.758         
OGF04  0.660         
OGI01  0.841         
OGI02  0.801         
OGI03  0.803         
OGI04  0.835         

               Factor1 Factor2
SS loadings       4.89    2.36
Proportion Var    0.31    0.15
Cumulative Var    0.31    0.45

Factor Correlations:
        Factor1 Factor2
Factor1   1.000   0.376
Factor2   0.376   1.000

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 305 on 89 degrees of freedom.
The p-value is 0.00000000000000000000000016 

We want to split the ingroup/outgroup fusion/identification items into distinct factors. We can examine what the three and four factor solutions looks like.

Three factor model

Display code

fit03 <- factanal(bonds, 3, rotation="promax")
fit03

Call:
factanal(x = bonds, factors = 3, rotation = "promax")

Uniquenesses:
IGF01 IGF02 IGF03 IGF04 IGI01 IGI02 IGI03 IGI04 OGF01 OGF02 OGF03 OGF04 OGI01 
 0.54  0.53  0.70  0.65  0.68  0.85  0.77  0.81  0.30  0.47  0.29  0.54  0.32 
OGI02 OGI03 OGI04 
 0.22  0.28  0.34 

Loadings:
      Factor1 Factor2 Factor3
IGF01  0.165   0.684  -0.110 
IGF02          0.726   0.164 
IGF03          0.532         
IGF04          0.594         
IGI01          0.531         
IGI02          0.362         
IGI03 -0.144   0.402         
IGI04          0.379         
OGF01  0.819                 
OGF02  0.743                 
OGF03  0.892          -0.102 
OGF04  0.679                 
OGI01  0.537           0.446 
OGI02  0.412           0.604 
OGI03  0.447           0.553 
OGI04  0.593           0.352 

               Factor1 Factor2 Factor3
SS loadings       3.55    2.36   1.052
Proportion Var    0.22    0.15   0.066
Cumulative Var    0.22    0.37   0.435

Factor Correlations:
        Factor1 Factor2 Factor3
Factor1   1.000  -0.338   0.476
Factor2  -0.338   1.000  -0.318
Factor3   0.476  -0.318   1.000

Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 206 on 75 degrees of freedom.
The p-value is 0.000000000000034 

Four factor model

Display code

fit04 <- factanal(bonds, 4, rotation="promax")
fit04

Call:
factanal(x = bonds, factors = 4, rotation = "promax")

Uniquenesses:
IGF01 IGF02 IGF03 IGF04 IGI01 IGI02 IGI03 IGI04 OGF01 OGF02 OGF03 OGF04 OGI01 
 0.54  0.55  0.69  0.65  0.67  0.83  0.77  0.79  0.19  0.49  0.32  0.50  0.33 
OGI02 OGI03 OGI04 
 0.26  0.18  0.03 

Loadings:
      Factor1 Factor2 Factor3 Factor4
IGF01  0.164   0.692  -0.132         
IGF02          0.709   0.119         
IGF03          0.537          -0.107 
IGF04          0.595                 
IGI01          0.535                 
IGI02          0.357  -0.103   0.132 
IGI03 -0.124   0.402                 
IGI04          0.401   0.106  -0.125 
OGF01  0.930                  -0.212 
OGF02  0.678                         
OGF03  0.867                         
OGF04  0.620          -0.154   0.288 
OGI01  0.437           0.343   0.232 
OGI02  0.330           0.563         
OGI03  0.297           0.722         
OGI04  0.391                   0.700 

               Factor1 Factor2 Factor3 Factor4
SS loadings       3.07    2.37   1.057   0.737
Proportion Var    0.19    0.15   0.066   0.046
Cumulative Var    0.19    0.34   0.406   0.452

Factor Correlations:
        Factor1 Factor2 Factor3 Factor4
Factor1   1.000  -0.458   0.350   0.559
Factor2  -0.458   1.000  -0.232  -0.467
Factor3   0.350  -0.232   1.000   0.302
Factor4   0.559  -0.467   0.302   1.000

Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 127 on 62 degrees of freedom.
The p-value is 0.0000024 

Section 3. Factor Analysis: Group fusion/identification (Version 2)

The fourth item has been removed from each sub-scale:

Display code

## Remove 4th item from all subscales and retry:

## New dataframe:
bonds <- cbind.data.frame(ds$IGF01, ds$IGF02, ds$IGF03,
                          ds$IGI01, ds$IGI02, ds$IGI03,
                          ds$OGF01, ds$OGF02, ds$OGF03,
                          ds$OGI01, ds$OGI02, ds$OGI03)

names(bonds) <- sub('ds\\$', '', names(bonds))

bonds <- na.omit(bonds)
mtx1 <- cor(bonds[, c(1:12)])

corrplot(mtx1, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

parallel <- fa.parallel(bonds)

Parallel analysis suggests that the number of factors =  2  and the number of components =  2 

Two factor model:

Display code

# Two factor model
fit04 <- factanal(bonds, 2, rotation="promax")
fit04

Call:
factanal(x = bonds, factors = 2, rotation = "promax")

Uniquenesses:
IGF01 IGF02 IGF03 IGI01 IGI02 IGI03 OGF01 OGF02 OGF03 OGI01 OGI02 OGI03 
 0.59  0.61  0.65  0.69  0.81  0.78  0.35  0.53  0.44  0.35  0.32  0.33 

Loadings:
      Factor1 Factor2
IGF01          0.664 
IGF02  0.102   0.654 
IGF03          0.587 
IGI01 -0.108   0.510 
IGI02          0.400 
IGI03 -0.155   0.395 
OGF01  0.810         
OGF02  0.705         
OGF03  0.764         
OGI01  0.819         
OGI02  0.808         
OGI03  0.809         

               Factor1 Factor2
SS loadings       3.77    1.80
Proportion Var    0.31    0.15
Cumulative Var    0.31    0.46

Factor Correlations:
        Factor1 Factor2
Factor1   1.000   0.369
Factor2   0.369   1.000

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 142 on 43 degrees of freedom.
The p-value is 0.0000000000018 

Three factor model

Display code

# Three factor model:
fit05 <- factanal(bonds, 3, rotation="promax")
fit05

Call:
factanal(x = bonds, factors = 3, rotation = "promax")

Uniquenesses:
IGF01 IGF02 IGF03 IGI01 IGI02 IGI03 OGF01 OGF02 OGF03 OGI01 OGI02 OGI03 
 0.58  0.59  0.66  0.69  0.81  0.78  0.25  0.51  0.28  0.35  0.22  0.26 

Loadings:
      Factor1 Factor2 Factor3
IGF01  0.189  -0.104   0.656 
IGF02          0.184   0.671 
IGF03                  0.582 
IGI01                  0.508 
IGI02                  0.404 
IGI03 -0.126           0.396 
OGF01  0.818                 
OGF02  0.632   0.110         
OGF03  0.916                 
OGI01  0.248   0.618         
OGI02          0.851         
OGI03          0.792         

               Factor1 Factor2 Factor3
SS loadings       2.04    1.80    1.80
Proportion Var    0.17    0.15    0.15
Cumulative Var    0.17    0.32    0.47

Factor Correlations:
        Factor1 Factor2 Factor3
Factor1   1.000  -0.329  -0.753
Factor2  -0.329   1.000   0.360
Factor3  -0.753   0.360   1.000

Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 49 on 33 degrees of freedom.
The p-value is 0.036 

Four factor model

Display code

# Four factor model:
fit06 <- factanal(bonds, 4, rotation="promax")
fit06

Call:
factanal(x = bonds, factors = 4, rotation = "promax")

Uniquenesses:
IGF01 IGF02 IGF03 IGI01 IGI02 IGI03 OGF01 OGF02 OGF03 OGI01 OGI02 OGI03 
0.508 0.607 0.650 0.707 0.005 0.783 0.261 0.501 0.254 0.352 0.218 0.262 

Loadings:
      Factor1 Factor2 Factor3 Factor4
IGF01          0.160   0.754  -0.115 
IGF02  0.186  -0.106   0.603         
IGF03                  0.587         
IGI01                  0.453         
IGI02                          1.024 
IGI03         -0.130   0.334         
OGF01  0.118   0.771                 
OGF02  0.153   0.587   0.101         
OGF03          0.939                 
OGI01  0.664   0.190                 
OGI02  0.892                         
OGI03  0.834                         

               Factor1 Factor2 Factor3 Factor4
SS loadings       2.02    1.92    1.61   1.094
Proportion Var    0.17    0.16    0.13   0.091
Cumulative Var    0.17    0.33    0.46   0.554

Factor Correlations:
        Factor1 Factor2 Factor3 Factor4
Factor1   1.000  -0.243  -0.443  -0.295
Factor2  -0.243   1.000   0.309   0.767
Factor3  -0.443   0.309   1.000   0.311
Factor4  -0.295   0.767   0.311   1.000

Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 31 on 24 degrees of freedom.
The p-value is 0.14 

FA basically suggests that ingroup / outgroup fusion/identification are not distinct factors. Only real split is between ingroup (identification + fusion) as one factor and outgroup (identification + fusion) as another factor. However, I’ll replicate the analysis as before - two constructs for ingroup fusion / identification and a single “bonds” construct for outgroup fusion+identification. Here is the reliability and inter-correlations of the three sub-scales:

Reliability
Scale 1: Ingroup fusion

Display code

igfusion <- cbind.data.frame(ds$IGF01, ds$IGF02, ds$IGF03)

names(igfusion) <- sub('ds\\$', '', names(igfusion))

igfusion <- na.omit(igfusion)
mtx2 <- cor(igfusion[, c(1:3)])

corrplot(mtx2, method = "number", number.cex = 0.7,
         col=c("darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

## Reliability:
alpha02 <- psych::alpha(igfusion)
alpha02

Reliability analysis   
Call: psych::alpha(x = igfusion)

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
       0.6       0.6     0.5      0.34 1.5 0.038  6.2 0.88     0.35

    95% confidence boundaries 
         lower alpha upper
Feldt     0.51   0.6  0.67
Duhachek  0.52   0.6  0.67

 Reliability if an item is dropped:
      raw_alpha std.alpha G6(smc) average_r  S/N alpha se var.r med.r
IGF01      0.44      0.46    0.30      0.30 0.84    0.059    NA  0.30
IGF02      0.52      0.52    0.35      0.35 1.08    0.053    NA  0.35
IGF03      0.52      0.53    0.36      0.36 1.13    0.052    NA  0.36

 Item statistics 
        n raw.r std.r r.cor r.drop mean  sd
IGF01 326  0.76  0.76  0.57   0.44  6.2 1.2
IGF02 326  0.69  0.74  0.52   0.40  6.3 1.0
IGF03 326  0.78  0.74  0.51   0.39  6.0 1.3

Non missing response frequency for each item
         1    2    3    4    5    6    7 miss
IGF01 0.01 0.02 0.02 0.01 0.09 0.34 0.51    0
IGF02 0.01 0.01 0.02 0.02 0.10 0.31 0.55    0
IGF03 0.01 0.03 0.06 0.01 0.13 0.29 0.48    0

Reliability
Scale 2: Ingroup identification

Display code

igidentification <- cbind.data.frame(ds$IGI01, ds$IGI02, ds$IGI03)

names(igidentification) <- sub('ds\\$', '', names(igidentification))

igidentification <- na.omit(igidentification)
mtx3 <- cor(igidentification[, c(1:3)])

corrplot(mtx3, method = "number", number.cex = 0.7,
         col=c("darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

## Reliability:
alpha03 <- psych::alpha(igidentification)
alpha03

Reliability analysis   
Call: psych::alpha(x = igidentification)

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
      0.54      0.54    0.44      0.28 1.2 0.044  6.2 0.86     0.28

    95% confidence boundaries 
         lower alpha upper
Feldt     0.44  0.54  0.62
Duhachek  0.45  0.54  0.62

 Reliability if an item is dropped:
      raw_alpha std.alpha G6(smc) average_r  S/N alpha se var.r med.r
IGI01      0.43      0.43    0.28      0.28 0.76    0.063    NA  0.28
IGI02      0.44      0.45    0.29      0.29 0.82    0.061    NA  0.29
IGI03      0.43      0.44    0.28      0.28 0.77    0.063    NA  0.28

 Item statistics 
        n raw.r std.r r.cor r.drop mean  sd
IGI01 325  0.75  0.72  0.48   0.35  6.1 1.3
IGI02 325  0.72  0.72  0.47   0.34  6.1 1.2
IGI03 325  0.69  0.72  0.48   0.35  6.3 1.1

Non missing response frequency for each item
         1    2    3    4    5    6    7 miss
IGI01 0.01 0.02 0.04 0.02 0.09 0.30 0.52    0
IGI02 0.01 0.02 0.04 0.03 0.07 0.36 0.48    0
IGI03 0.01 0.01 0.03 0.03 0.06 0.33 0.54    0

Reliability
Scale 3: Outgroup fusion/identification

Display code

ogbonds <- cbind.data.frame(ds$OGF01, ds$OGF02, ds$OGF03,
                            ds$OGI01, ds$OGI02, ds$OGI03)

names(ogbonds) <- sub('ds\\$', '', names(ogbonds))

ogbonds <- na.omit(ogbonds)
mtx3 <- cor(ogbonds[, c(1:6)])

corrplot(mtx3, method = "number", number.cex = 0.7,
         col=c("darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

## Reliability:
alpha04 <- psych::alpha(ogbonds)
alpha04

Reliability analysis   
Call: psych::alpha(x = ogbonds)

  raw_alpha std.alpha G6(smc) average_r S/N    ase mean  sd median_r
       0.9       0.9     0.9      0.61 9.3 0.0089  3.3 1.7     0.59

    95% confidence boundaries 
         lower alpha upper
Feldt     0.88   0.9  0.92
Duhachek  0.88   0.9  0.92

 Reliability if an item is dropped:
      raw_alpha std.alpha G6(smc) average_r S/N alpha se  var.r med.r
OGF01      0.88      0.88    0.87      0.60 7.4   0.0110 0.0077  0.56
OGF02      0.90      0.90    0.89      0.64 8.7   0.0096 0.0056  0.61
OGF03      0.89      0.89    0.87      0.61 7.9   0.0105 0.0071  0.60
OGI01      0.88      0.88    0.87      0.60 7.5   0.0110 0.0068  0.59
OGI02      0.88      0.88    0.87      0.60 7.5   0.0109 0.0036  0.59
OGI03      0.88      0.88    0.87      0.60 7.4   0.0110 0.0050  0.59

 Item statistics 
        n raw.r std.r r.cor r.drop mean  sd
OGF01 295  0.84  0.84  0.80   0.76  3.5 2.1
OGF02 295  0.76  0.76  0.69   0.65  3.4 2.1
OGF03 295  0.81  0.81  0.76   0.72  3.5 2.1
OGI01 295  0.83  0.84  0.79   0.75  3.1 2.0
OGI02 295  0.83  0.83  0.81   0.75  3.2 2.1
OGI03 295  0.84  0.84  0.81   0.76  3.3 2.1

Non missing response frequency for each item
         1    2    3    4    5    6    7 miss
OGF01 0.19 0.28 0.13 0.04 0.11 0.13 0.13    0
OGF02 0.22 0.24 0.16 0.05 0.13 0.08 0.13    0
OGF03 0.21 0.24 0.15 0.04 0.12 0.14 0.11    0
OGI01 0.26 0.27 0.14 0.03 0.11 0.11 0.08    0
OGI02 0.27 0.24 0.15 0.02 0.10 0.12 0.10    0
OGI03 0.27 0.23 0.13 0.03 0.12 0.12 0.11    0

Section 4. Factor Analysis: BCL/BBL endorsement/experience

Display code

## BCL and BBL items:

# BCL_01:
# Seek out opportunities to bridge social divisions with their opponents, enemies, opposition groups, or other outgroups.   
# Variables: Q6.1, Q7.1

# BCL_02:
# Demonstrate willingness to compromise with their opponents, enemies, opposition groups, or other outgroups. 
# Variables: Q6.2, Q7.2

# BCL_03:
# Try to understand and empathize with their opponents, enemies, opposition groups, or other outgroups.  
# Variables: Q6.3, Q7.3

# BBL_01:
# Represent the interests of the communities and groups that they belong to even at the cost of other groups.
# Variables: Q6.4, Q7.4

# BBL_02:
# Focus on building stronger connections within the communities and groups they belong to rather than building stronger relationships with other groups across boundaries.
# Variables: Q6.5, Q7.5

# BBL_03:
# Try to gain benefits for the communities and groups they belong to even at the expense of other groups.
# Variables: Q6.6, Q7.6

ds$ENDBCL01 <- as.numeric(ds$`Q6.1 Leadership Quality`)
ds$ENDBCL02 <- as.numeric(ds$`Q6.2 Leadership Quality`)
ds$ENDBCL03 <- as.numeric(ds$`Q6.3 Leadership Quality`)
ds$ENDBBL01 <- as.numeric(ds$`Q6.4 Leadership Quality`)
ds$ENDBBL02 <- as.numeric(ds$`Q6.5 Leadership Quality`)
ds$ENDBBL03 <- as.numeric(ds$`Q6.6 Leadership Quality`)

ds$EXPBCL01 <- as.numeric(ds$`Q7.1 Leadership Experience`)
ds$EXPBCL02 <- as.numeric(ds$`Q7.2 Leadership Experience`)
ds$EXPBCL03 <- as.numeric(ds$`Q7.3 Leadership Experience`)
ds$EXPBBL01 <- as.numeric(ds$`Q7.4 Leadership Experience`)
ds$EXPBBL02 <- as.numeric(ds$`Q7.5 Leadership Experience`)
ds$EXPBBL03 <- as.numeric(ds$`Q7.6 Leadership Experience`)

leadership <- cbind.data.frame(ds$ENDBCL01, ds$ENDBCL02, ds$ENDBCL03, 
                               ds$ENDBBL01, ds$ENDBBL02, ds$ENDBBL03,
                               ds$EXPBCL01, ds$EXPBCL02, ds$EXPBCL03,
                               ds$EXPBBL01, ds$EXPBBL02, ds$EXPBBL03)

names(leadership) <- sub('ds\\$', '', names(leadership))

leadership <- na.omit(leadership)
mtx1 <- cor(leadership[, c(1:12)])

Shorthand for item names:
END = Endorse, EXP = Experience
BCL = BCL, BBL = BBL
So for example:
ENDBCL01 = “Endorsement of BCL, item 1”
EXPBBL02 = “Experience of BBL, item 2”, and so on

Display code

corrplot(mtx1, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

The correlation plot above shows the correlation coefficients between the twelve different items (six for endorsement of BCL/BBL, and six for experience of BCL/BBL).

KMO and Bartlett’s test:

Display code

## Kaiser-Meyer-Olkin (KMO) test of factorability
KMO(r=cor(leadership))

Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = cor(leadership))
Overall MSA =  0.77
MSA for each item = 
ENDBCL01 ENDBCL02 ENDBCL03 ENDBBL01 ENDBBL02 ENDBBL03 EXPBCL01 EXPBCL02 
    0.90     0.78     0.79     0.76     0.64     0.67     0.85     0.76 
EXPBCL03 EXPBBL01 EXPBBL02 EXPBBL03 
    0.82     0.72     0.74     0.69 

Display code

## Bartlett's test of sphericity
cortest.bartlett(leadership)

$chisq
[1] 853

$p.value
[1] 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000033

$df
[1] 66

Parallel test

Display code

parallel <- fa.parallel(leadership)

Parallel analysis suggests that the number of factors =  3  and the number of components =  3 

Three factor model

Display code

fit03 <- factanal(leadership, 3, rotation="promax")
fit03

Call:
factanal(x = leadership, factors = 3, rotation = "promax")

Uniquenesses:
ENDBCL01 ENDBCL02 ENDBCL03 ENDBBL01 ENDBBL02 ENDBBL03 EXPBCL01 EXPBCL02 
    0.65     0.41     0.56     0.63     0.40     0.47     0.47     0.22 
EXPBCL03 EXPBBL01 EXPBBL02 EXPBBL03 
    0.44     0.40     0.43     0.42 

Loadings:
         Factor1 Factor2 Factor3
ENDBCL01  0.576  -0.110         
ENDBCL02  0.731  -0.173         
ENDBCL03  0.557           0.284 
ENDBBL01          0.241   0.492 
ENDBBL02         -0.138   0.805 
ENDBBL03 -0.199           0.696 
EXPBCL01  0.733   0.113         
EXPBCL02  0.899          -0.170 
EXPBCL03  0.746   0.189         
EXPBBL01          0.788         
EXPBBL02          0.743         
EXPBBL03          0.761         

               Factor1 Factor2 Factor3
SS loadings       3.13    1.94    1.51
Proportion Var    0.26    0.16    0.13
Cumulative Var    0.26    0.42    0.55

Factor Correlations:
        Factor1 Factor2 Factor3
Factor1  1.0000 -0.0207   0.161
Factor2 -0.0207  1.0000   0.311
Factor3  0.1613  0.3105   1.000

Test of the hypothesis that 3 factors are sufficient.
The chi square statistic is 72 on 33 degrees of freedom.
The p-value is 0.00011 

Four factor model

Display code

fit04 <- factanal(leadership, 4, rotation="promax")
fit04

Call:
factanal(x = leadership, factors = 4, rotation = "promax")

Uniquenesses:
ENDBCL01 ENDBCL02 ENDBCL03 ENDBBL01 ENDBBL02 ENDBBL03 EXPBCL01 EXPBCL02 
   0.642    0.005    0.541    0.628    0.416    0.451    0.417    0.212 
EXPBCL03 EXPBBL01 EXPBBL02 EXPBBL03 
   0.433    0.424    0.431    0.356 

Loadings:
         Factor1 Factor2 Factor3 Factor4
ENDBCL01  0.544  -0.150                 
ENDBCL02  0.261                   0.855 
ENDBCL03  0.334           0.226   0.352 
ENDBBL01          0.198   0.513         
ENDBBL02         -0.157   0.804         
ENDBBL03 -0.151           0.726         
EXPBCL01  0.812                  -0.113 
EXPBCL02  0.870          -0.133         
EXPBCL03  0.754   0.116                 
EXPBBL01  0.112   0.737                 
EXPBBL02          0.737                 
EXPBBL03 -0.131   0.848           0.124 

               Factor1 Factor2 Factor3 Factor4
SS loadings       2.52    1.91    1.52   0.908
Proportion Var    0.21    0.16    0.13   0.076
Cumulative Var    0.21    0.37    0.50   0.572

Factor Correlations:
        Factor1 Factor2 Factor3 Factor4
Factor1   1.000  0.2103  0.5197   0.175
Factor2   0.210  1.0000 -0.0591  -0.354
Factor3   0.520 -0.0591  1.0000   0.147
Factor4   0.175 -0.3541  0.1470   1.000

Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 41 on 24 degrees of freedom.
The p-value is 0.016 

Again, FA structure does not really support a four factor Endorse/Experience BCL/BBL structure. I will still split the sub-scales into those four constructs.

Display code

end_bcl <- cbind.data.frame(ds$ENDBCL01, ds$ENDBCL02, ds$ENDBCL03)
end_bbl <- cbind.data.frame(ds$ENDBBL01, ds$ENDBBL02, ds$ENDBBL03)
exp_bcl <- cbind.data.frame(ds$EXPBCL01, ds$EXPBCL02, ds$EXPBCL03)
exp_bbl <- cbind.data.frame(ds$EXPBBL01, ds$EXPBBL02, ds$EXPBBL03)

names(end_bcl) <- sub('ds\\$', '', names(end_bcl))
names(end_bbl) <- sub('ds\\$', '', names(end_bbl))
names(exp_bcl) <- sub('ds\\$', '', names(exp_bcl))
names(exp_bbl) <- sub('ds\\$', '', names(exp_bbl))

end_bcl <- na.omit(end_bcl)
end_bbl <- na.omit(end_bbl)
exp_bcl <- na.omit(exp_bcl)
exp_bbl <- na.omit(exp_bbl)

Reliability
Scale 1: Endorsement of BCL

Display code

mtx1 <- cor(end_bcl[, c(1:3)])

corrplot(mtx1, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

alph01 <- psych::alpha(end_bcl)
alph01

Reliability analysis   
Call: psych::alpha(x = end_bcl)

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
      0.65      0.65    0.58      0.39 1.9 0.033  5.3 1.3     0.45

    95% confidence boundaries 
         lower alpha upper
Feldt     0.58  0.65  0.71
Duhachek  0.59  0.65  0.72

 Reliability if an item is dropped:
         raw_alpha std.alpha G6(smc) average_r  S/N alpha se var.r med.r
ENDBCL01      0.63      0.63    0.46      0.46 1.69    0.041    NA  0.46
ENDBCL02      0.40      0.41    0.25      0.25 0.68    0.065    NA  0.25
ENDBCL03      0.61      0.62    0.45      0.45 1.63    0.042    NA  0.45

 Item statistics 
           n raw.r std.r r.cor r.drop mean  sd
ENDBCL01 322  0.70  0.74  0.52   0.41  5.8 1.5
ENDBCL02 322  0.83  0.83  0.71   0.57  5.3 1.8
ENDBCL03 322  0.77  0.74  0.53   0.43  4.9 1.9

Non missing response frequency for each item
            1    2    3    4    5    6    7 miss
ENDBCL01 0.04 0.03 0.03 0.01 0.17 0.32 0.40    0
ENDBCL02 0.04 0.11 0.07 0.01 0.19 0.30 0.29    0
ENDBCL03 0.05 0.13 0.08 0.06 0.18 0.27 0.22    0

Reliability
Scale 2: Endorsement of BBL

Display code

mtx2 <- cor(end_bbl[, c(1:3)])

corrplot(mtx2, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

alph02 <- psych::alpha(end_bbl)
alph02

Reliability analysis   
Call: psych::alpha(x = end_bbl)

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
      0.67      0.67    0.58       0.4   2 0.032  5.2 1.4     0.37

    95% confidence boundaries 
         lower alpha upper
Feldt      0.6  0.67  0.72
Duhachek   0.6  0.67  0.73

 Reliability if an item is dropped:
         raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
ENDBBL01      0.64      0.64    0.47      0.47 1.8    0.040    NA  0.47
ENDBBL02      0.53      0.53    0.36      0.36 1.1    0.052    NA  0.36
ENDBBL03      0.54      0.54    0.37      0.37 1.2    0.051    NA  0.37

 Item statistics 
           n raw.r std.r r.cor r.drop mean  sd
ENDBBL01 324  0.75  0.74  0.52   0.42  5.0 1.8
ENDBBL02 324  0.78  0.79  0.63   0.51  5.4 1.7
ENDBBL03 324  0.80  0.79  0.62   0.50  5.2 1.9

Non missing response frequency for each item
            1    2    3    4    5    6    7 miss
ENDBBL01 0.06 0.10 0.09 0.05 0.19 0.31 0.21    0
ENDBBL02 0.04 0.06 0.05 0.03 0.18 0.35 0.29    0
ENDBBL03 0.09 0.07 0.05 0.03 0.16 0.31 0.29    0

Reliability
Scale 3: Experience of BCL

Display code

mtx3 <- cor(exp_bcl[, c(1:3)])

corrplot(mtx3, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

alph03 <- psych::alpha(exp_bcl)
alph03

Reliability analysis   
Call: psych::alpha(x = exp_bcl)

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
      0.78      0.78    0.71      0.54 3.5 0.021  4.7 1.6     0.58

    95% confidence boundaries 
         lower alpha upper
Feldt     0.74  0.78  0.82
Duhachek  0.74  0.78  0.82

 Reliability if an item is dropped:
         raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
EXPBCL01      0.74      0.74    0.59      0.59 2.9    0.028    NA  0.59
EXPBCL02      0.62      0.62    0.45      0.45 1.7    0.042    NA  0.45
EXPBCL03      0.73      0.74    0.58      0.58 2.8    0.029    NA  0.58

 Item statistics 
           n raw.r std.r r.cor r.drop mean  sd
EXPBCL01 330  0.81  0.81  0.66   0.58  4.8 1.9
EXPBCL02 330  0.87  0.87  0.78   0.69  4.8 2.0
EXPBCL03 330  0.82  0.82  0.67   0.59  4.5 2.0

Non missing response frequency for each item
            1    2    3    4    5    6    7 miss
EXPBCL01 0.09 0.08 0.05 0.14 0.19 0.26 0.18    0
EXPBCL02 0.06 0.13 0.09 0.10 0.12 0.27 0.23    0
EXPBCL03 0.08 0.11 0.16 0.09 0.13 0.22 0.20    0

Reliability
Scale 2: Experience of BBL

Display code

mtx4 <- cor(exp_bbl[, c(1:3)])

corrplot(mtx4, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

alph04 <- psych::alpha(exp_bbl)
alph04

Reliability analysis   
Call: psych::alpha(x = exp_bbl)

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
       0.8       0.8    0.73      0.58 4.1 0.025  3.9 1.6     0.58

    95% confidence boundaries 
         lower alpha upper
Feldt     0.75   0.8  0.85
Duhachek  0.75   0.8  0.85

 Reliability if an item is dropped:
         raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
EXPBBL01      0.73      0.73    0.58      0.58 2.8    0.039    NA  0.58
EXPBBL02      0.72      0.72    0.56      0.56 2.6    0.040    NA  0.56
EXPBBL03      0.74      0.74    0.59      0.59 2.8    0.038    NA  0.59

 Item statistics 
           n raw.r std.r r.cor r.drop mean  sd
EXPBBL01 190  0.84  0.85  0.72   0.65  3.9 1.8
EXPBBL02 190  0.86  0.85  0.74   0.66  4.0 2.0
EXPBBL03 190  0.85  0.84  0.72   0.64  3.7 1.9

Non missing response frequency for each item
            1    2    3    4    5    6    7 miss
EXPBBL01 0.08 0.22 0.14 0.19 0.13 0.14 0.09    0
EXPBBL02 0.09 0.22 0.12 0.16 0.08 0.19 0.13    0
EXPBBL03 0.13 0.26 0.08 0.23 0.07 0.12 0.11    0

The four subscales have varying reliability scores.

Section 5. Preliminary regression analysis

Research question:
Do group fusion / identification predict BCL vs BBL?

Regression models

Display code

## Proper (usable) variables in the model:

ds$Age1 <- as.numeric(gsub("\\D", "", ds$Age))
ds$Age <- as.numeric(ifelse(ds$Age1 == "6465", 64.5,
                     ifelse(ds$Age1 == "806", 80.5, 
                     ifelse(ds$Age1 < 18, NA, ds$Age1))))

ds$Female <- ifelse(ds$gender == "Female", 1, 0)
ds$Married <- ifelse(ds$married == "Married", 1, 0)
ds$`SES-` <- ds$ses

ds$Endorse_BCL <- as.numeric((ds$ENDBCL01+ds$ENDBCL02+ds$ENDBCL03)/3)
ds$Endorse_BBL <- as.numeric((ds$ENDBBL01+ds$ENDBBL02+ds$ENDBBL03)/3)
ds$Experience_BCL <- as.numeric((ds$EXPBCL01+ds$EXPBCL02+ds$EXPBCL03)/3)
ds$Experience_BBL <- as.numeric((ds$EXPBBL01+ds$EXPBBL02+ds$EXPBBL03)/3)

ds$IG_Fusion <- as.numeric((ds$IGF01+ds$IGF02+ds$IGF03)/3)
ds$IG_Identification <- as.numeric((ds$IGI01+ds$IGI02+ds$IGI03)/3)
ds$OG_Bonds <- as.numeric((ds$OGF01+ds$OGF02+ds$OGF03+
                             ds$OGI01+ds$OGI02+ds$OGI03)/6)

Display code

## Four regression models predicting endorsement and experience of BCL / BBL:

lm01 <- lm(Endorse_BCL~IG_Fusion+IG_Identification+OG_Bonds+Age+Female+Married+`SES-`,
           data = ds)

lm02 <- lm(Experience_BCL~IG_Fusion+IG_Identification+OG_Bonds+Age+Female+Married+`SES-`, 
           data = ds)

lm03 <- lm(Endorse_BBL~IG_Fusion+IG_Identification+OG_Bonds+Age+Female+Married+`SES-`, 
           data = ds)

lm04 <- lm(Experience_BBL~IG_Fusion+IG_Identification+OG_Bonds+Age+Female+Married+`SES-`, 
           data = ds)

Display code

## Tabulated results:

stargazer(lm01, lm02,
          lm03, lm04, 
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	Endorse_BCL	Experience_BCL	Endorse_BBL	Experience_BBL
	(1)	(2)	(3)	(4)
IG_Fusion	0.290**	0.099	0.060	-0.100
	(0.100)	(0.120)	(0.100)	(0.150)
IG_Identification	-0.008	-0.005	-0.034	-0.340*
	(0.110)	(0.130)	(0.120)	(0.170)
OG_Bonds	0.160**	0.200**	0.068	-0.001
	(0.051)	(0.062)	(0.054)	(0.081)
Age	0.00001	0.00003	0.00001	-0.009
	(0.00003)	(0.00004)	(0.00003)	(0.010)
Female	-0.180	-0.150	-0.160	0.042
	(0.160)	(0.190)	(0.170)	(0.260)
Married	-0.250	-0.520**	0.070	-0.450
	(0.160)	(0.190)	(0.170)	(0.280)
`SES-`Middle	0.300	0.270	-0.200	0.260
	(0.180)	(0.220)	(0.190)	(0.280)
`SES-`Upper middle	0.063	0.280	0.300	1.200
	(0.300)	(0.380)	(0.330)	(0.710)
`SES-`Upper	1.200	1.100	-0.380	1.200
	(0.660)	(0.810)	(0.690)	(0.960)
Constant	3.100***	3.600***	5.000***	7.000***
	(0.760)	(0.950)	(0.800)	(1.200)
Observations	268	274	268	156
R2	0.090	0.081	0.029	0.093
Adjusted R2	0.058	0.050	-0.005	0.037
Residual Std. Error	1.300 (df = 258)	1.600 (df = 264)	1.300 (df = 258)	1.600 (df = 146)
F Statistic	2.800** (df = 9; 258)	2.600** (df = 9; 264)	0.870 (df = 9; 258)	1.700 (df = 9; 146)
Note:	*p<0.05; **p<0.01; ***p<0.001

Section 6. Ingroup/outgroup empathy/hostility

Empathic concern

Display code

## Empathic concern

ds$empathic_concern_01 <- as.numeric(ds$`Q18.1 Empathy`)

ds$empathic_concern_02a <- as.numeric(ds$`Q18.2 Empathy`)
ds$empathic_concern_02 <- as.numeric(8 - ds$empathic_concern_02a)

ds$empathic_concern_03 <- as.numeric(ds$`Q18.3 Empathy`)

ds$empathic_concern <- (ds$empathic_concern_01+ds$empathic_concern_02+
                        ds$empathic_concern_03)/3

summary(ds$empathic_concern)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    2.7     4.7     5.3     5.4     6.3     7.0      27 

Display code

ds %>% drop_na(empathic_concern)%>%
ggplot(aes(x = empathic_concern))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 30)+
  geom_textvline(label = "Mean = 5.40", 
                 xintercept = 5.40, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Empathic concern score", 
       y = "Frequency", 
       title = "Empathic concern: Tanzania")+
  theme_bw()

Perspective taking

Display code

ds$perspective_taking_01 <- as.numeric(ds$`Q18.4 Empathy`)
ds$perspective_taking_02 <- as.numeric(ds$`Q18.5 Empathy`)
ds$perspective_taking_03 <- as.numeric(ds$`Q18.6 Empathy`)
ds$perspective_taking_04 <- as.numeric(ds$`Q18.7 Empathy`)

ds$perspective_taking <- (ds$perspective_taking_01+ds$perspective_taking_02+
                          ds$perspective_taking_03+ds$perspective_taking_04)/4

summary(ds$perspective_taking)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      2       5       6       6       6       7      36 

Display code

ds %>% drop_na(perspective_taking)%>%
ggplot(aes(x = perspective_taking))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 30)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Perspective taking score", 
       y = "Frequency", 
       title = "Perspective taking: Tanzania")+
  theme_bw()

Outgroup cooperation

Display code

ds$og_coop_01 <- as.numeric(ds$`Q12.1 Outgroup`)
ds$og_coop_02 <- as.numeric(ds$`Q12.2 Outgroup`)

ds$og_cooperation <- (ds$og_coop_01+ds$og_coop_02)/2

summary(ds$og_cooperation)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     5.0     6.0     5.4     6.5     7.0      18 

Display code

ds %>% drop_na(og_cooperation)%>%
ggplot(aes(x = og_cooperation))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 30)+
  geom_textvline(label = "Mean = 5.40", 
                 xintercept = 5.40, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Outgroup cooperation score", 
       y = "Frequency", 
       title = "Outgroup cooperation: Tanzania")+
  theme_bw()

Perceived history of discrimination

Display code

ds$history_discrimination <- as.numeric(ds$`Q12.3 Outgroup`)

summary(ds$history_discrimination)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     2.0     3.0     3.6     6.0     7.0      18 

Display code

ds %>% drop_na(history_discrimination)%>%
ggplot(aes(x = history_discrimination))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 20)+
  geom_textvline(label = "Mean = 3.60", 
                 xintercept = 3.60, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Perceived history of discrimination score", 
       y = "Frequency", 
       title = "Perceived history of discrimination: Tanzania")+
  theme_bw()

Outgroup hostility

Display code

ds$og_host_01 <- as.numeric(ds$`Q12.4 Outgroup`)
ds$og_host_02 <- as.numeric(ds$`Q12.5 Outgroup`)

ds$og_hostility <- (ds$og_host_01+ds$og_host_02)/2

summary(ds$og_hostility)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     1.5     3.5     3.6     6.0     7.0      25 

Display code

ds %>% drop_na(og_hostility)%>%
ggplot(aes(x = og_hostility))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 30)+
  geom_textvline(label = "Mean = 3.60", 
                 xintercept = 3.60, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Outgroup hostility score", 
       y = "Frequency", 
       title = "Outgroup hostility: Tanzania")+
  theme_bw()

Willingess to fight outgroup

Display code

ds$fight_outgroup <- as.numeric(ds$`Q12.6 Outgroup`)

summary(ds$fight_outgroup)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     1.0     2.0     3.5     6.0     7.0      16 

Display code

ds %>% drop_na(fight_outgroup)%>%
ggplot(aes(x = fight_outgroup))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 30)+
  geom_textvline(label = "Mean = 3.50", 
                 xintercept = 3.50, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Willingness to fight outgroup score", 
       y = "Frequency", 
       title = "Willingness to fight outgroup: Tanzania")+
  theme_bw()

Section 7. Correlation between ingroup/outgroup empathy/hostility items

Display code

ds3 <- cbind.data.frame(ds$empathic_concern, ds$perspective_taking,
                        ds$og_cooperation, ds$history_discrimination,
                        ds$og_hostility, ds$fight_outgroup)
ds3 <- na.omit(ds3)

names(ds3) <- sub('ds\\$', '', names(ds3))

mtx <- cor(ds3[, c(1:6)])

corrplot(mtx, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Section 8. Regression with new variables

This set of regression models have additional predictor variables compared to the last set. It adds empathic concern, perspective taking, and perceived history of discrimination as additional predictors.

Display code

## Four regression models predicting endorsement and experience of BCL / BBL:

lm01 <- lm(Endorse_BCL~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`, 
           data = ds)

lm02 <- lm(Experience_BCL~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`, 
           data = ds)

lm03 <- lm(Endorse_BBL~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`, 
           data = ds)

lm04 <- lm(Experience_BBL~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`, 
           data = ds)

Display code

## Tabulated results:

stargazer(lm01, lm02,
          lm03, lm04, 
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	Endorse_BCL	Experience_BCL	Endorse_BBL	Experience_BBL
	(1)	(2)	(3)	(4)
IG_Fusion	0.250*	0.110	-0.030	-0.086
	(0.100)	(0.130)	(0.120)	(0.150)
IG_Identification	0.034	0.042	-0.017	-0.330*
	(0.100)	(0.130)	(0.120)	(0.160)
OG_Bonds	0.160***	0.180**	0.044	-0.012
	(0.049)	(0.063)	(0.056)	(0.083)
empathic_concern	-0.270***	-0.340***	-0.040	-0.088
	(0.074)	(0.094)	(0.084)	(0.140)
perspective_taking	0.370***	0.240*	0.250*	-0.099
	(0.095)	(0.120)	(0.110)	(0.160)
history_discrimination	-0.110**	-0.045	0.042	0.210***
	(0.037)	(0.047)	(0.042)	(0.063)
Age	-0.00000	0.00002	0.00001	-0.006
	(0.00003)	(0.00004)	(0.00003)	(0.010)
Female	-0.077	-0.078	-0.180	-0.043
	(0.150)	(0.200)	(0.170)	(0.270)
Married	-0.200	-0.510**	0.046	-0.340
	(0.150)	(0.200)	(0.170)	(0.280)
`SES-`Middle	0.160	0.110	-0.280	0.140
	(0.180)	(0.230)	(0.200)	(0.280)
`SES-`Upper middle	0.076	0.120	0.170	0.990
	(0.300)	(0.390)	(0.350)	(0.740)
`SES-`Upper	0.850	0.810	-0.550	0.940
	(0.620)	(0.800)	(0.700)	(0.920)
Constant	2.900***	4.000***	4.300***	7.200***
	(0.860)	(1.100)	(0.980)	(1.500)
Observations	255	260	255	146
R2	0.200	0.130	0.049	0.200
Adjusted R2	0.160	0.087	0.002	0.130
Residual Std. Error	1.200 (df = 242)	1.500 (df = 247)	1.300 (df = 242)	1.500 (df = 133)
F Statistic	5.000*** (df = 12; 242)	3.100*** (df = 12; 247)	1.000 (df = 12; 242)	2.900** (df = 12; 133)
Note:	*p<0.05; **p<0.01; ***p<0.001

Section 9. Alternative regression models with different outcomes

This set of regression models have different outcome variables which are: outgroup cooperation, outgroup hostility, and willingness to fight outgroup.

Display code

lm01 <- lm(og_cooperation~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`, 
           data = ds)

lm02 <- lm(og_hostility~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`, 
           data = ds)

lm03 <- lm(fight_outgroup~IG_Fusion+IG_Identification+OG_Bonds+empathic_concern+perspective_taking+history_discrimination+Age+Female+Married+`SES-`, 
           data = ds)

Display code

## Tabulated results:

stargazer(lm01, lm02,
          lm03, type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	og_cooperation	og_hostility	fight_outgroup
	(1)	(2)	(3)
IG_Fusion	0.290*	-0.120	0.017
	(0.120)	(0.120)	(0.140)
IG_Identification	0.070	-0.120	-0.070
	(0.120)	(0.120)	(0.150)
OG_Bonds	0.430***	0.027	0.140*
	(0.056)	(0.056)	(0.071)
empathic_concern	0.026	-0.330***	-0.420***
	(0.084)	(0.085)	(0.110)
perspective_taking	0.086	-0.066	-0.320*
	(0.110)	(0.110)	(0.130)
history_discrimination	-0.110*	0.570***	0.550***
	(0.042)	(0.042)	(0.053)
Age	0.00001	0.00003	-0.00002
	(0.00003)	(0.00004)	(0.00004)
Female	0.087	0.032	-0.100
	(0.170)	(0.170)	(0.220)
Married	0.120	0.340	0.270
	(0.170)	(0.170)	(0.220)
`SES-`Middle	-0.390*	0.190	-0.001
	(0.200)	(0.200)	(0.250)
`SES-`Upper middle	-0.580	0.220	0.510
	(0.340)	(0.350)	(0.430)
`SES-`Upper	-1.900**	0.400	0.540
	(0.710)	(0.720)	(0.900)
Constant	1.800	4.700***	5.300***
	(0.970)	(0.980)	(1.200)
Observations	262	260	263
R2	0.260	0.560	0.470
Adjusted R2	0.220	0.530	0.440
Residual Std. Error	1.400 (df = 249)	1.400 (df = 247)	1.700 (df = 250)
F Statistic	7.200*** (df = 12; 249)	26.000*** (df = 12; 247)	18.000*** (df = 12; 250)
Note:	*p<0.05; **p<0.01; ***p<0.001

Section 10. Religious freedom measures

Social perception of religious freedom

Display code

ds$religious_freedom_per_01 <- as.numeric(ds$`Q16.1 Religious Freedom`)
ds$religious_freedom_per_02a <- as.numeric(ds$`Q16.2 Religious Freedom`)
ds$religious_freedom_per_02 <- (8 - ds$religious_freedom_per_02a)
ds$religious_freedom_per_03 <- as.numeric(ds$`Q16.3 Religious Freedom`)
ds$religious_freedom_per_04 <- as.numeric(ds$`Q16.4 Religious Freedom`)
ds$religious_freedom_per_05 <- as.numeric(ds$`Q16.5 Religious Freedom`)
ds$religious_freedom_per_06 <- as.numeric(ds$`Q16.6 Religious Freedom`)
ds$religious_freedom_per_07 <- as.numeric(ds$`Q16.7 Religious Freedom`)
ds$religious_freedom_per_08 <- as.numeric(ds$`Q16.8 Religious Freedom`)

ds$sprf <- (ds$religious_freedom_per_01+ds$religious_freedom_per_02+
            ds$religious_freedom_per_03+ds$religious_freedom_per_04+
            ds$religious_freedom_per_05+ds$religious_freedom_per_06+
            ds$religious_freedom_per_07+ds$religious_freedom_per_08)/8

summary(ds$sprf)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    2.5     5.0     5.5     5.4     5.8     7.0      29 

Display code

ds %>% drop_na(sprf)%>%
ggplot(aes(x = sprf))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 30)+
  geom_textvline(label = "Mean = 5.40", 
                 xintercept = 5.40, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "SPRF score", 
       y = "Frequency", 
       title = "Social perception of religious freedom: Tanzania")+
  theme_bw()

Experience of religious freedom

Display code

ds$religious_freedom_exp_01a <- as.numeric(ds$`Q17.1 Life experience`)
ds$religious_freedom_exp_01 <- (8 - ds$religious_freedom_exp_01a)

ds$religious_freedom_exp_02a <- as.numeric(ds$`Q17.2 Life experience`)
ds$religious_freedom_exp_02 <- (8 - ds$religious_freedom_exp_02a)

ds$religious_freedom_exp_03a <- as.numeric(ds$`Q17.3 Life experience`)
ds$religious_freedom_exp_03 <- (8 - ds$religious_freedom_exp_03a)

ds$religious_freedom_exp_04a <- as.numeric(ds$`Q17.4 Life experience`)
ds$religious_freedom_exp_04 <- (8 - ds$religious_freedom_exp_04a)

ds$religious_freedom_exp_05a <- as.numeric(ds$`Q17.1 Life experience`)
ds$religious_freedom_exp_05 <- (8 - ds$religious_freedom_exp_05a)

ds$religious_freedom_exp_06a <- as.numeric(ds$`Q17.1 Life experience`)
ds$religious_freedom_exp_06 <- (8 - ds$religious_freedom_exp_06a)

ds$exp_religious_freedom <- (ds$religious_freedom_exp_01+ds$religious_freedom_exp_02+
                             ds$religious_freedom_exp_03+ds$religious_freedom_exp_04+
                              ds$religious_freedom_exp_05+ds$religious_freedom_exp_06)/6

summary(ds$exp_religious_freedom)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       5       6       6       7       7      35 

Display code

ds %>% drop_na(exp_religious_freedom)%>%
ggplot(aes(x = exp_religious_freedom))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 30)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Experience of religious freedom score", 
       y = "Frequency", 
       title = "Experience of religious freedom: Tanzania")+
  theme_bw()

Experience of religious freedom based on religious affiliation

Display code

ds$Religion <- ds$religion

ds2 <- ds %>% 
drop_na(Religion, exp_religious_freedom)%>%
 group_by(Religion) %>% 
 summarise(Exp_religious_freedom=mean(exp_religious_freedom),.groups = 'drop') %>%
  as.data.frame()

ds2

                 Religion Exp_religious_freedom
1     Christian: Catholic                   4.9
2        Christian: Other                   6.5
3  Christian: Pentecostal                   6.4
4   Christian: Protestant                   5.7
5            Muslim: Shia                   6.6
6           Muslim: Sunni                   6.5
7                   Other                   6.2
8 Spiritual not religious                   5.8

Display code

ds %>% drop_na(religion, exp_religious_freedom)%>%
ggplot(aes(y = exp_religious_freedom, 
           x = religion))+
  geom_boxplot()+
  labs(x = "", 
       y = "Experience of religious freedom score", 
       title = "Experience of religious freedom: Tanzania")+
  coord_flip()+
  theme_bw()

Section 11. Positive/Negative contact with outgroup

Positive contact with outgroup

Display code

ds$pc01 <- ds$`Never`

## convert to sentence case:
ds$pc01 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$pc01, perl = TRUE)

ds$pc01 <- factor(ds$pc01,
                  levels = c("Never", "Very rarely", "Rarely",
                             "Sometimes", "Often", "Very often", "Always"))

table(ds$pc01)

      Never Very rarely      Rarely   Sometimes       Often  Very often 
         19          23          50          51          45          56 
     Always 
         66 

Display code

# never, very rarely, rarely
# sometimes, often, very often, always

ds$pc02 <- as.numeric(ds$pc01)

summary(ds$pc02)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       3       5       5       6       7      42 

Display code

ds %>% drop_na(pc02)%>%
ggplot(aes(x = pc02))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 30)+
  geom_textvline(label = "Mean = 5.00", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Positive contact score", 
       y = "Frequency", 
       title = "Positive contact with outgroup: Tanzania")+
  theme_bw()

Negative contact with outgroup

Display code

ds$nc01 <- ds$`Q14 Negativ Outgroup`

## convert to sentence case:
ds$nc01 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$nc01, perl = TRUE)

ds$nc01 <- ifelse(ds$nc01 == "Ver rarely", "Very rarely", ds$nc01)

ds$nc01 <- factor(ds$nc01,
                  levels = c("Never", "Very rarely", "Rarely",
                             "Sometimes", "Often", "Very often", "Always"))

table(ds$nc01)

      Never Very rarely      Rarely   Sometimes       Often  Very often 
        112          86          59          38           5           5 
     Always 
          5 

Display code

# never, very rarely, rarely
# sometimes, often, very often, always

ds$nc02 <- as.numeric(ds$nc01)

summary(ds$nc02)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       1       2       2       3       7      42 

Display code

ds %>% drop_na(pc02)%>%
ggplot(aes(x = pc02))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 30)+
  geom_textvline(label = "Mean = 2.00", 
                 xintercept = 2.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Negative contact score", 
       y = "Frequency", 
       title = "Negative contact with outgroup: Tanzania")+
  theme_bw()

Section 12. Inter-group marriage

No data for Tanzania

Section 13. Outgroup affect

Outgroup affect: Scale construction

Display code

## Feel outgroup: negative to positive:

ds$ogaf1 <- as.character(ds$`Q13.1 Outgroup`)

## convert to sentence case:
ds$ogaf1 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf1, perl = TRUE)

ds$og_aff_01 <- factor(ds$ogaf1, 
                       levels=c("Very negative", "Moderately negative",
                                 "A little negative", "Neutral", 
                                 "A little positive", "Moderately positive", 
                                 "Very positive"))


table(ds$og_aff_01)

      Very negative Moderately negative   A little negative             Neutral 
                 25                  27                  41                  36 
  A little positive Moderately positive       Very positive 
                 74                  82                  50 

Display code

## Feel Outgroup: Hostile to friendly:

ds$ogaf2 <- as.character(ds$`Q13.2 Outgroup`)

## convert to sentence case:
ds$ogaf2 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf2, perl = TRUE)

ds$og_aff_02 <- factor(ds$ogaf2, 
                       levels = c("Very hostile", "Moderately hostile", 
                                  "A little hostile", "Neutral", 
                                  "A little friendly", "Moderately friendly",
                                  "Very friendly"))

table(ds$og_aff_02)

       Very hostile  Moderately hostile    A little hostile             Neutral 
                 25                  19                  42                  30 
  A little friendly Moderately friendly       Very friendly 
                 88                  85                  46 

Display code

## Feel Outgroup: Suspicious to trusting:

ds$ogaf3 <- as.character(ds$`Q13.3 Outgroup`)

## convert to sentence case:
ds$ogaf3 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf3, perl = TRUE)

# remove extra spaces:
ds$ogaf3 <- gsub("\\s+"," ",ds$ogaf3, perl = TRUE)

ds$og_aff_03 <- factor(ds$ogaf3, 
                       levels = c("Very suspicious", "Moderately suspicious", 
                                  "A little suspicious", "Neutral",
                                  "A little trusting", "Moderately trusting",
                                  "Very trusting"))
table(ds$og_aff_03)

      Very suspicious Moderately suspicious   A little suspicious 
                   26                    35                    47 
              Neutral     A little trusting   Moderately trusting 
                   28                    77                    89 
        Very trusting 
                   33 

Display code

## Feel Outgroup: Contempt to respect:

ds$ogaf4 <- as.character(ds$`Q13.4 Outgroup`)

## convert to sentence case:
ds$ogaf4 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf4, perl = TRUE)

ds$og_aff_04 <- factor(ds$ogaf4, 
                       levels = c("A lot of contempt", "Moderate contempt",
                                  "A little contempt", "Neutral", 
                                  "A little respect", "Moderate respect",
                                  "A lot of respect"))

table(ds$og_aff_04)

A lot of contempt Moderate contempt A little contempt           Neutral 
               24                21                30                31 
 A little respect  Moderate respect  A lot of respect 
               79                90                59 

Display code

## Feel Outgroup: Concerned to Unconcerned:

ds$ogaf5 <- as.character(ds$`Q13.5 Outgroup`)

## convert to sentence case:
ds$ogaf5 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf5, perl = TRUE)

# remove extra spaces:
ds$ogaf5 <- gsub("\\s+"," ",ds$ogaf5, perl = TRUE)

ds$og_aff_05 <- factor(ds$ogaf5, 
                       levels = c("Very concerned", "Moderately concerned", 
                                  "A little concerned", "Neutral",
                                  "A little unconcerned", "Moderately unconcerned",
                                  "Very unconcerned"))
table(ds$og_aff_05)

        Very concerned   Moderately concerned     A little concerned 
                    37                     67                     62 
               Neutral   A little unconcerned Moderately unconcerned 
                    29                     63                     53 
      Very unconcerned 
                    18 

Display code

## Feel Outgroup: Threatened to Relaxed:

ds$ogaf6 <- as.character(ds$`Q13.6 Outgroup`)

## convert to sentence case:
ds$ogaf6 <- gsub("(\\D)(\\D+)", "\\U\\1\\L\\2", ds$ogaf6, perl = TRUE)

ds$og_aff_06 <- factor(ds$ogaf6, 
                       levels = c("Very threatened", "Moderately threatened", 
                                  "A little threatened", "Neutral", 
                                  "A little relaxed", "Moderately relaxed", 
                                  "Very relaxed"))
table(ds$og_aff_06)

      Very threatened Moderately threatened   A little threatened 
                   25                    33                    35 
              Neutral      A little relaxed    Moderately relaxed 
                   53                    53                    89 
         Very relaxed 
                   47 

Outgroup affect: Factor analysis

Display code

ds$NegativeToPositive <- as.numeric(ds$og_aff_01) 
ds$HostileToFriendly <- as.numeric(ds$og_aff_02) 
ds$SuspiciousToTrusting <- as.numeric(ds$og_aff_03) 
ds$ContemptToRespect <- as.numeric(ds$og_aff_04) 
ds$ConcernedToUnconcerned <- as.numeric(ds$og_aff_05) 
ds$ThreatenedToRelaxed <- as.numeric(ds$og_aff_06) 

og_affect <- cbind.data.frame(ds$NegativeToPositive, ds$HostileToFriendly,
                              ds$SuspiciousToTrusting, ds$ContemptToRespect,
                              ds$ConcernedToUnconcerned, ds$ThreatenedToRelaxed)

names(og_affect) <- sub('ds\\$', '', names(og_affect))

og_affect <- na.omit(og_affect)

mtx1 <- cor(og_affect[, c(1:6)])

corrplot(mtx1, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

## Kaiser-Meyer-Olkin (KMO) test of factorability
KMO(r=cor(og_affect))

Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = cor(og_affect))
Overall MSA =  0.9
MSA for each item = 
    NegativeToPositive      HostileToFriendly   SuspiciousToTrusting 
                  0.91                   0.89                   0.91 
     ContemptToRespect ConcernedToUnconcerned    ThreatenedToRelaxed 
                  0.89                   0.48                   0.93 

Display code

## Bartlett's test of sphericity
cortest.bartlett(og_affect)

$chisq
[1] 1442

$p.value
[1] 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000013

$df
[1] 15

Display code

## Parallel test
parallel <- fa.parallel(og_affect)

Parallel analysis suggests that the number of factors =  1  and the number of components =  1 

Display code

## Reliability test
psych::alpha(og_affect)

Reliability analysis   
Call: psych::alpha(x = og_affect)

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
      0.87      0.87    0.89      0.54 6.9 0.011  4.5 1.4     0.72

    95% confidence boundaries 
         lower alpha upper
Feldt     0.85  0.87  0.89
Duhachek  0.85  0.87  0.89

 Reliability if an item is dropped:
                       raw_alpha std.alpha G6(smc) average_r  S/N alpha se
NegativeToPositive          0.83      0.83    0.85      0.50  4.9   0.0151
HostileToFriendly           0.82      0.82    0.85      0.48  4.6   0.0158
SuspiciousToTrusting        0.82      0.83    0.85      0.49  4.8   0.0155
ContemptToRespect           0.83      0.83    0.85      0.49  4.8   0.0154
ConcernedToUnconcerned      0.94      0.94    0.93      0.76 15.6   0.0053
ThreatenedToRelaxed         0.84      0.84    0.87      0.51  5.2   0.0146
                        var.r med.r
NegativeToPositive     0.1134  0.70
HostileToFriendly      0.1156  0.69
SuspiciousToTrusting   0.1238  0.70
ContemptToRespect      0.1100  0.69
ConcernedToUnconcerned 0.0023  0.77
ThreatenedToRelaxed    0.1346  0.77

 Item statistics 
                         n raw.r std.r r.cor r.drop mean  sd
NegativeToPositive     326  0.87  0.87  0.86   0.80  4.7 1.8
HostileToFriendly      326  0.90  0.90  0.90   0.85  4.7 1.7
SuspiciousToTrusting   326  0.88  0.89  0.87   0.82  4.5 1.8
ContemptToRespect      326  0.88  0.88  0.88   0.82  4.9 1.8
ConcernedToUnconcerned 326  0.32  0.31  0.12   0.11  3.7 1.8
ThreatenedToRelaxed    326  0.84  0.84  0.80   0.76  4.6 1.8

Non missing response frequency for each item
                          1    2    3    4    5    6    7 miss
NegativeToPositive     0.07 0.08 0.12 0.10 0.23 0.25 0.15    0
HostileToFriendly      0.07 0.06 0.13 0.08 0.27 0.25 0.14    0
SuspiciousToTrusting   0.07 0.10 0.14 0.08 0.23 0.27 0.10    0
ContemptToRespect      0.07 0.06 0.09 0.09 0.24 0.27 0.17    0
ConcernedToUnconcerned 0.11 0.21 0.19 0.09 0.19 0.16 0.06    0
ThreatenedToRelaxed    0.07 0.10 0.10 0.16 0.16 0.27 0.14    0

Based on the output of factor analysis as well as reliability, it seems like one item (“concerned to unconcerned”) is problematic. Here is the scale without that item:

Outgroup affect (minus problematic item)

Display code

og_affect <- cbind.data.frame(ds$NegativeToPositive, ds$HostileToFriendly,
                              ds$SuspiciousToTrusting, ds$ContemptToRespect,
                              ds$ThreatenedToRelaxed)

names(og_affect) <- sub('ds\\$', '', names(og_affect))

og_affect <- na.omit(og_affect)
mtx1 <- cor(og_affect[, c(1:5)])

corrplot(mtx1, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

## Reliability analysis:
psych::alpha(og_affect)

Reliability analysis   
Call: psych::alpha(x = og_affect)

  raw_alpha std.alpha G6(smc) average_r S/N    ase mean  sd median_r
      0.94      0.94    0.92      0.75  15 0.0055  4.7 1.6     0.75

    95% confidence boundaries 
         lower alpha upper
Feldt     0.93  0.94  0.95
Duhachek  0.93  0.94  0.95

 Reliability if an item is dropped:
                     raw_alpha std.alpha G6(smc) average_r S/N alpha se   var.r
NegativeToPositive        0.92      0.92    0.90      0.74  12   0.0072 0.00183
HostileToFriendly         0.91      0.92    0.89      0.73  11   0.0077 0.00148
SuspiciousToTrusting      0.92      0.92    0.90      0.75  12   0.0071 0.00198
ContemptToRespect         0.92      0.92    0.90      0.75  12   0.0070 0.00281
ThreatenedToRelaxed       0.93      0.93    0.91      0.78  14   0.0060 0.00054
                     med.r
NegativeToPositive    0.73
HostileToFriendly     0.73
SuspiciousToTrusting  0.75
ContemptToRespect     0.75
ThreatenedToRelaxed   0.78

 Item statistics 
                       n raw.r std.r r.cor r.drop mean  sd
NegativeToPositive   332  0.90  0.90  0.87   0.84  4.7 1.8
HostileToFriendly    332  0.92  0.92  0.90   0.87  4.7 1.7
SuspiciousToTrusting 332  0.90  0.90  0.87   0.84  4.5 1.8
ContemptToRespect    332  0.89  0.89  0.86   0.83  4.9 1.8
ThreatenedToRelaxed  332  0.86  0.85  0.79   0.77  4.6 1.8

Non missing response frequency for each item
                        1    2    3    4    5    6    7 miss
NegativeToPositive   0.08 0.08 0.12 0.11 0.22 0.24 0.15    0
HostileToFriendly    0.08 0.06 0.13 0.08 0.27 0.25 0.14    0
SuspiciousToTrusting 0.08 0.10 0.14 0.08 0.23 0.27 0.10    0
ContemptToRespect    0.07 0.06 0.09 0.09 0.24 0.27 0.18    0
ThreatenedToRelaxed  0.08 0.10 0.10 0.16 0.16 0.27 0.14    0

Based on the above output for factor analysis and reliability analysis, the “outgroup affect” scale is highly reliable if we drop the problematic item (“concerned to unconcerned”). Here is the visualization for the scale after dropping the problematic item:

Display code

ds$og_affect <- (ds$NegativeToPositive+ds$HostileToFriendly+
                 ds$SuspiciousToTrusting+ds$ContemptToRespect+
                 ds$ThreatenedToRelaxed)/5

summary(ds$og_affect)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     3.6     5.0     4.7     6.0     7.0      20 

Display code

ds %>% drop_na(og_affect)%>%
ggplot(aes(x = og_affect))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 4.70", 
                 xintercept = 4.70, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Outgroup affect score", 
       y = "Frequency", 
       title = "Outgroup affect: Tanzania")+
  theme_bw()

Section 14. Updated regression models

This set of regression models have additional predictors: frequency of positive and negative contact, and interaction between perspective taking and perceived history of discrimination

Display code

ds$freq_positive_contact <- ds$pc02
ds$freq_negative_contact <- ds$nc02
ds$perspectiveXdiscrimination <- ds$perspective_taking*ds$history_discrimination

## Four regression models predicting endorsement and experience of BCL / BBL:

lm01 <- lm(Endorse_BCL~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`, 
           data = ds)

lm02 <- lm(Experience_BCL~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`, 
           data = ds)

lm03 <- lm(Endorse_BBL~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`, 
           data = ds)

lm04 <- lm(Experience_BBL~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`, 
           data = ds)

Display code

## Tabulated results:

stargazer(lm01, lm02,
          lm03, lm04, 
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	Endorse_BCL	Experience_BCL	Endorse_BBL	Experience_BBL
	(1)	(2)	(3)	(4)
IG_Fusion	0.200	0.029	-0.024	-0.035
	(0.110)	(0.130)	(0.120)	(0.150)
IG_Identification	0.010	-0.004	-0.060	-0.380*
	(0.110)	(0.130)	(0.120)	(0.160)
OG_Bonds	0.160**	0.120	0.018	-0.019
	(0.053)	(0.066)	(0.061)	(0.085)
freq_positive_contact	0.094	0.240***	0.110*	-0.036
	(0.048)	(0.059)	(0.054)	(0.089)
freq_negative_contact	0.016	0.015	0.140	0.270*
	(0.065)	(0.081)	(0.074)	(0.100)
empathic_concern	-0.320***	-0.340***	-0.042	-0.089
	(0.081)	(0.100)	(0.093)	(0.140)
perspective_taking	0.450*	0.540*	0.051	-0.024
	(0.180)	(0.230)	(0.210)	(0.290)
history_discrimination	0.092	0.580	-0.220	0.190
	(0.250)	(0.310)	(0.280)	(0.370)
perspectiveXdiscrimination	-0.033	-0.100	0.047	0.001
	(0.043)	(0.053)	(0.049)	(0.063)
Age	-0.00000	0.00003	0.00002	-0.007
	(0.00003)	(0.00004)	(0.00004)	(0.010)
Female	-0.067	-0.044	-0.200	-0.096
	(0.160)	(0.200)	(0.180)	(0.260)
Married	-0.180	-0.440*	0.017	-0.320
	(0.160)	(0.200)	(0.180)	(0.270)
`SES-`Middle	0.150	0.120	-0.270	0.140
	(0.180)	(0.230)	(0.210)	(0.280)
`SES-`Upper middle	0.140	0.240	0.160	0.890
	(0.320)	(0.400)	(0.370)	(0.730)
`SES-`Upper	0.870	0.940	-0.550	0.730
	(0.630)	(0.790)	(0.720)	(0.910)
Constant	2.600*	1.800	4.900***	6.500***
	(1.200)	(1.500)	(1.400)	(1.800)
Observations	237	243	237	144
R2	0.220	0.190	0.073	0.260
Adjusted R2	0.170	0.130	0.010	0.170
Residual Std. Error	1.200 (df = 221)	1.500 (df = 227)	1.400 (df = 221)	1.500 (df = 128)
F Statistic	4.100*** (df = 15; 221)	3.500*** (df = 15; 227)	1.200 (df = 15; 221)	2.900*** (df = 15; 128)
Note:	*p<0.05; **p<0.01; ***p<0.001

Section 15. Updated regression models (alternate DVs)

Same predictors as section 14 above, but with different outcomes. This set of regression models has an additional outcome variable: affect towards outgroup.

Display code

## Four regression models:

lm01 <- lm(og_cooperation~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`, 
           data = ds)

lm02 <- lm(og_hostility~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`, 
           data = ds)

lm03 <- lm(fight_outgroup~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`, 
           data = ds)

lm04 <- lm(og_affect~IG_Fusion+IG_Identification+OG_Bonds+freq_positive_contact+freq_negative_contact+empathic_concern+perspective_taking+history_discrimination+perspectiveXdiscrimination+Age+Female+Married+`SES-`, 
           data = ds)

Display code

## Tabulated results:

stargazer(lm01, lm02,
          lm03, lm04, 
          type = "html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "table1.html")

Display code

htmltools::includeHTML("table1.html")

	Dependent variable:
	og_cooperation	og_hostility	fight_outgroup	og_affect
	(1)	(2)	(3)	(4)
IG_Fusion	0.230	0.003	0.140	-0.140
	(0.120)	(0.110)	(0.150)	(0.100)
IG_Identification	-0.014	-0.079	-0.057	0.150
	(0.120)	(0.110)	(0.150)	(0.100)
OG_Bonds	0.380***	-0.019	0.140	0.120*
	(0.057)	(0.056)	(0.074)	(0.052)
freq_positive_contact	0.310***	-0.140**	-0.200**	0.550***
	(0.051)	(0.050)	(0.066)	(0.046)
freq_negative_contact	0.160*	-0.034	-0.065	0.100
	(0.070)	(0.068)	(0.091)	(0.063)
empathic_concern	-0.011	-0.200*	-0.350**	-0.140
	(0.087)	(0.086)	(0.110)	(0.079)
perspective_taking	0.062	-0.240	-0.460	0.150
	(0.200)	(0.200)	(0.260)	(0.180)
history_discrimination	0.047	0.240	0.190	-0.130
	(0.270)	(0.260)	(0.340)	(0.240)
perspectiveXdiscrimination	-0.019	0.063	0.063	-0.009
	(0.046)	(0.045)	(0.059)	(0.041)
Age	0.00003	0.00003	-0.00002	-0.00001
	(0.00003)	(0.00003)	(0.00004)	(0.00003)
Female	0.150	0.041	-0.100	0.015
	(0.170)	(0.170)	(0.220)	(0.160)
Married	0.170	0.089	0.150	-0.064
	(0.170)	(0.170)	(0.230)	(0.160)
`SES-`Middle	-0.390*	0.066	-0.082	0.038
	(0.190)	(0.190)	(0.250)	(0.180)
`SES-`Upper middle	-0.560	0.230	0.380	-0.078
	(0.340)	(0.340)	(0.450)	(0.310)
`SES-`Upper	-1.800*	0.310	0.470	-0.770
	(0.690)	(0.670)	(0.890)	(0.620)
Constant	1.200	4.900***	6.000***	2.100
	(1.300)	(1.300)	(1.700)	(1.200)
Observations	244	242	245	243
R2	0.350	0.600	0.480	0.510
Adjusted R2	0.310	0.570	0.450	0.480
Residual Std. Error	1.300 (df = 228)	1.300 (df = 226)	1.700 (df = 229)	1.200 (df = 227)
F Statistic	8.200*** (df = 15; 228)	23.000*** (df = 15; 226)	14.000*** (df = 15; 229)	16.000*** (df = 15; 227)
Note:	*p<0.05; **p<0.01; ***p<0.001

Section 16. Imagistic Experience: Factor analysis

Display code

ds$event_negative_affect <- as.numeric(ds$`Q9.1 Feelings about Event`)
table(ds$event_negative_affect)

  1   2   3   4   5   6   7 
 60  28   8   4  11  63 158 

Display code

ds$event_positive_affect_01 <- (8 - ds$event_negative_affect)
table(ds$event_positive_affect_01)

  1   2   3   4   5   6   7 
158  63  11   4   8  28  60 

Display code

ds$event_positive_affect_02 <- as.numeric(ds$`Q9.2 Feelings about Event`)
table(ds$event_positive_affect_02)

  1   2   3   4   5   6   7 
125  38  10   5  30  31  91 

Display code

ds$event_positive_affect <- ds$event_positive_affect_02

ds$event_episodic_recall_01 <- as.numeric(ds$`Q9.3 Feelings about Event`)
table(ds$event_episodic_recall_01)

  1   2   3   4   5   6   7 
  6   9  15  11  27 121 136 

Display code

ds$event_episodic_recall_02 <- as.numeric(ds$`Q9.4 Feelings about Event`)
table(ds$event_episodic_recall_02)

  1   2   3   4   5   6   7 
  3   9  15  10  20 123 152 

Display code

ds$event_shared_perception_01 <- as.numeric(ds$`Q9.5 Feelings about Event`)
table(ds$event_shared_perception_01)

 1  2  3  4  5  6  7 
17 28 33  7 66 81 99 

Display code

ds$event_shared_perception_02 <- as.numeric(ds$`Q9.6 Feelings about Event`)
table(ds$event_shared_perception_02)

  1   2   3   4   5   6   7 
  5  23  40   9  57  94 104 

Display code

ds$event_event_reflection_01 <- as.numeric(ds$`Q9.9 Feelings about Event`)
table(ds$event_event_reflection_01)

  1   2   3   4   5   6   7 
  7  26  30   8  55  95 103 

Display code

ds$event_event_reflection_02 <- as.numeric(ds$`Q9.10 Feelings about Event`)
table(ds$event_event_reflection_02)

 1  2  3  4  5  6  7 
18 39 37  9 66 88 75 

Display code

ds$event_transformative_indiv_01 <- as.numeric(ds$`Q9.7 Feelings about Event`)

table(ds$event_transformative_indiv_01)

  1   2   3   4   5   6   7 
 28  36  21   2  34 103 106 

Display code

ds$event_transformative_indiv_02 <- as.numeric(ds$`Q9.8 Feelings about Event`)
table(ds$event_transformative_indiv_02)

  1   2   3   4   5   6   7 
 23  43  39   8  53  56 104 

Display code

ds$event_transformative_group_01 <- as.numeric(ds$`Q9.11 Feelings about Event`)
table(ds$event_transformative_group_01)

 1  2  3  4  5  6  7 
21 41 36  7 62 87 78 

Display code

ds$event_transformative_group_02 <- as.numeric(ds$`Q9.12 Feelings about Event`)
table(ds$event_transformative_group_02)

 1  2  3  4  5  6  7 
55 55 50 10 44 61 57 

Display code

imagistic <- cbind.data.frame(ds$event_negative_affect, ds$event_positive_affect, 
                              ds$event_episodic_recall_01, ds$event_episodic_recall_02,
                              ds$event_shared_perception_01, ds$event_shared_perception_02,
                              ds$event_event_reflection_01, ds$event_event_reflection_02,
                              ds$event_transformative_indiv_01, ds$event_transformative_indiv_02,
                              ds$event_transformative_group_01, ds$event_transformative_group_02)

imagistic <- na.omit(imagistic)

names(imagistic) <- sub('ds\\$event\\_', '', names(imagistic))

mtx <- cor(imagistic[, c(1:12)])

corrplot(mtx, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Imagistic experience: Factor analysis

Display code

## Kaiser-Meyer-Olkin (KMO) test of factorability
KMO(r=cor(imagistic))

Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = cor(imagistic))
Overall MSA =  0.58
MSA for each item = 
        negative_affect         positive_affect      episodic_recall_01 
                   0.40                    0.49                    0.62 
     episodic_recall_02    shared_perception_01    shared_perception_02 
                   0.61                    0.59                    0.53 
    event_reflection_01     event_reflection_02 transformative_indiv_01 
                   0.69                    0.65                    0.40 
transformative_indiv_02 transformative_group_01 transformative_group_02 
                   0.70                    0.67                    0.59 

Display code

## Bartlett's test of sphericity
cortest.bartlett(imagistic)

$chisq
[1] 707

$p.value
[1] 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000049

$df
[1] 66

Display code

## Parallel test
parallel <- fa.parallel(imagistic)

Parallel analysis suggests that the number of factors =  5  and the number of components =  5 

Scree plot suggests there needs to be five factors, but theoretically we want six factors. Here is how the five factors break down:

Display code

fit05 <- factanal(imagistic, 5, rotation="promax")
fit05

Call:
factanal(x = imagistic, factors = 5, rotation = "promax")

Uniquenesses:
        negative_affect         positive_affect      episodic_recall_01 
                  0.005                   0.535                   0.518 
     episodic_recall_02    shared_perception_01    shared_perception_02 
                  0.439                   0.683                   0.350 
    event_reflection_01     event_reflection_02 transformative_indiv_01 
                  0.612                   0.562                   0.725 
transformative_indiv_02 transformative_group_01 transformative_group_02 
                  0.619                   0.502                   0.344 

Loadings:
                        Factor1 Factor2 Factor3 Factor4 Factor5
negative_affect                  1.008                         
positive_affect          0.529  -0.377                         
episodic_recall_01      -0.159                   0.655   0.181 
episodic_recall_02                               0.735   0.119 
shared_perception_01     0.218   0.186                   0.408 
shared_perception_02     0.128          -0.170   0.223   0.838 
event_reflection_01     -0.115  -0.109   0.589   0.171         
event_reflection_02      0.135           0.660                 
transformative_indiv_01  0.265   0.233  -0.163   0.259  -0.259 
transformative_indiv_02  0.158           0.546   0.140  -0.242 
transformative_group_01  0.693                           0.102 
transformative_group_02  0.717           0.196  -0.235   0.132 

               Factor1 Factor2 Factor3 Factor4 Factor5
SS loadings       1.49    1.27    1.20    1.20    1.08
Proportion Var    0.12    0.11    0.10    0.10    0.09
Cumulative Var    0.12    0.23    0.33    0.43    0.52

Factor Correlations:
        Factor1  Factor2 Factor3 Factor4 Factor5
Factor1  1.0000  0.02162  0.0333 0.12311   0.274
Factor2  0.0216  1.00000  0.0833 0.00517  -0.111
Factor3  0.0333  0.08325  1.0000 0.10781  -0.154
Factor4  0.1231  0.00517  0.1078 1.00000   0.326
Factor5  0.2742 -0.11136 -0.1535 0.32626   1.000

Test of the hypothesis that 5 factors are sufficient.
The chi square statistic is 31 on 16 degrees of freedom.
The p-value is 0.013 

The factor structure above makes sense but is a little incoherent. Here is what happens if we “force” six factors:

Display code

fit06 <- factanal(imagistic, 6, rotation="promax")
fit06

Call:
factanal(x = imagistic, factors = 6, rotation = "promax")

Uniquenesses:
        negative_affect         positive_affect      episodic_recall_01 
                  0.005                   0.519                   0.434 
     episodic_recall_02    shared_perception_01    shared_perception_02 
                  0.494                   0.456                   0.534 
    event_reflection_01     event_reflection_02 transformative_indiv_01 
                  0.623                   0.512                   0.672 
transformative_indiv_02 transformative_group_01 transformative_group_02 
                  0.632                   0.561                   0.005 

Loadings:
                        Factor1 Factor2 Factor3 Factor4 Factor5 Factor6
negative_affect          1.071                           0.141   0.273 
positive_affect         -0.441                                   0.337 
episodic_recall_01                       0.776          -0.116         
episodic_recall_02              -0.166   0.649                   0.234 
shared_perception_01     0.158  -0.182  -0.148   0.157   0.845         
shared_perception_02             0.169   0.318  -0.195   0.443  -0.294 
event_reflection_01                      0.210   0.510          -0.222 
event_reflection_02                              0.691   0.164         
transformative_indiv_01  0.216           0.132          -0.122   0.644 
transformative_indiv_02          0.103   0.108   0.511           0.134 
transformative_group_01          0.290                   0.224   0.298 
transformative_group_02          1.106                  -0.134         

               Factor1 Factor2 Factor3 Factor4 Factor5 Factor6
SS loadings       1.44    1.43    1.23    1.08   1.068   0.916
Proportion Var    0.12    0.12    0.10    0.09   0.089   0.076
Cumulative Var    0.12    0.24    0.34    0.43   0.520   0.597

Factor Correlations:
        Factor1   Factor2   Factor3 Factor4 Factor5 Factor6
Factor1  1.0000 -0.242759  0.038077  0.6215  0.2100  0.4698
Factor2 -0.2428  1.000000  0.000127 -0.1633  0.1441 -0.3482
Factor3  0.0381  0.000127  1.000000 -0.2292 -0.3078  0.0408
Factor4  0.6215 -0.163294 -0.229237  1.0000  0.0998  0.1936
Factor5  0.2100  0.144112 -0.307798  0.0998  1.0000  0.1300
Factor6  0.4698 -0.348185  0.040808  0.1936  0.1300  1.0000

Test of the hypothesis that 6 factors are sufficient.
The chi square statistic is 17 on 9 degrees of freedom.
The p-value is 0.051 

FA solution still seems incoherent.

Imagistic experience: Reliability

Display code

imagistic <- cbind.data.frame(ds$event_negative_affect, ds$event_positive_affect, 
                              ds$event_episodic_recall_01, ds$event_episodic_recall_02,
                              ds$event_shared_perception_01, ds$event_shared_perception_02,
                              ds$event_event_reflection_01, ds$event_event_reflection_02,
                              ds$event_transformative_indiv_01, ds$event_transformative_indiv_02,
                              ds$event_transformative_group_01, ds$event_transformative_group_02)


## Reliability:
alph01 <- psych::alpha(imagistic)
alph01

Reliability analysis   
Call: psych::alpha(x = imagistic)

  raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
       0.6      0.62     0.7      0.12 1.6 0.032    5 0.86     0.11

    95% confidence boundaries 
         lower alpha upper
Feldt     0.53   0.6  0.66
Duhachek  0.54   0.6  0.66

 Reliability if an item is dropped:
                                 raw_alpha std.alpha G6(smc) average_r S/N
ds$event_negative_affect              0.61      0.62    0.68      0.13 1.7
ds$event_positive_affect              0.63      0.63    0.69      0.14 1.7
ds$event_episodic_recall_01           0.58      0.60    0.67      0.12 1.5
ds$event_episodic_recall_02           0.58      0.59    0.67      0.12 1.5
ds$event_shared_perception_01         0.56      0.59    0.67      0.11 1.4
ds$event_shared_perception_02         0.59      0.61    0.68      0.13 1.6
ds$event_event_reflection_01          0.58      0.60    0.67      0.12 1.5
ds$event_event_reflection_02          0.55      0.58    0.66      0.11 1.4
ds$event_transformative_indiv_01      0.60      0.62    0.69      0.13 1.7
ds$event_transformative_indiv_02      0.55      0.58    0.66      0.11 1.4
ds$event_transformative_group_01      0.54      0.57    0.65      0.11 1.3
ds$event_transformative_group_02      0.53      0.58    0.64      0.11 1.4
                                 alpha se var.r med.r
ds$event_negative_affect            0.030 0.021 0.119
ds$event_positive_affect            0.029 0.019 0.144
ds$event_episodic_recall_01         0.033 0.022 0.136
ds$event_episodic_recall_02         0.033 0.024 0.104
ds$event_shared_perception_01       0.035 0.024 0.100
ds$event_shared_perception_02       0.032 0.023 0.119
ds$event_event_reflection_01        0.033 0.024 0.136
ds$event_event_reflection_02        0.036 0.024 0.100
ds$event_transformative_indiv_01    0.032 0.024 0.144
ds$event_transformative_indiv_02    0.036 0.024 0.085
ds$event_transformative_group_01    0.037 0.022 0.100
ds$event_transformative_group_02    0.037 0.021 0.100

 Item statistics 
                                   n raw.r std.r r.cor r.drop mean  sd
ds$event_negative_affect         332  0.35  0.31  0.22  0.106  5.1 2.4
ds$event_positive_affect         330  0.33  0.27  0.17  0.077  3.7 2.6
ds$event_episodic_recall_01      325  0.35  0.43  0.36  0.209  5.9 1.4
ds$event_episodic_recall_02      332  0.39  0.46  0.39  0.258  6.0 1.3
ds$event_shared_perception_01    331  0.49  0.49  0.42  0.328  5.2 1.9
ds$event_shared_perception_02    332  0.33  0.37  0.28  0.171  5.4 1.7
ds$event_event_reflection_01     324  0.42  0.46  0.37  0.254  5.4 1.7
ds$event_event_reflection_02     332  0.52  0.54  0.48  0.383  4.9 1.9
ds$event_transformative_indiv_01 330  0.35  0.31  0.19  0.149  5.2 2.0
ds$event_transformative_indiv_02 326  0.54  0.53  0.47  0.371  4.9 2.1
ds$event_transformative_group_01 332  0.58  0.55  0.51  0.429  4.9 1.9
ds$event_transformative_group_02 332  0.60  0.54  0.52  0.436  4.0 2.2

Non missing response frequency for each item
                                    1    2    3    4    5    6    7 miss
ds$event_negative_affect         0.18 0.08 0.02 0.01 0.03 0.19 0.48 0.06
ds$event_positive_affect         0.38 0.12 0.03 0.02 0.09 0.09 0.28 0.06
ds$event_episodic_recall_01      0.02 0.03 0.05 0.03 0.08 0.37 0.42 0.08
ds$event_episodic_recall_02      0.01 0.03 0.05 0.03 0.06 0.37 0.46 0.06
ds$event_shared_perception_01    0.05 0.08 0.10 0.02 0.20 0.24 0.30 0.06
ds$event_shared_perception_02    0.02 0.07 0.12 0.03 0.17 0.28 0.31 0.06
ds$event_event_reflection_01     0.02 0.08 0.09 0.02 0.17 0.29 0.32 0.08
ds$event_event_reflection_02     0.05 0.12 0.11 0.03 0.20 0.27 0.23 0.06
ds$event_transformative_indiv_01 0.08 0.11 0.06 0.01 0.10 0.31 0.32 0.06
ds$event_transformative_indiv_02 0.07 0.13 0.12 0.02 0.16 0.17 0.32 0.07
ds$event_transformative_group_01 0.06 0.12 0.11 0.02 0.19 0.26 0.23 0.06
ds$event_transformative_group_02 0.17 0.17 0.15 0.03 0.13 0.18 0.17 0.06

Imagistic experience: Visualization

If we conceptualize “imagistic experience” as a single measure that is the sum of the six scales, here is the visualization of the measure:

Display code

ds$imagistic <- (((ds$event_positive_affect)+
                 (ds$event_negative_affect)+ 
                ((ds$event_episodic_recall_01+ds$event_episodic_recall_02)/2)+
                ((ds$event_shared_perception_01+ds$event_shared_perception_02)/2)+
                ((ds$event_event_reflection_01+ds$event_event_reflection_02)/2)+
                ((ds$event_transformative_indiv_01+ds$event_transformative_indiv_02)/2)+
                ((ds$event_transformative_group_01+ds$event_transformative_group_02)/2)))

summary(ds$imagistic)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
     12      32      36      35      38      49      53 

Display code

ds %>% drop_na(imagistic)%>%
ggplot(aes(x = imagistic))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 35.00", 
                 xintercept = 35.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Imagistic experience", 
       y = "Frequency", 
       title = "Imagistic experience: Tanzania")+
  theme_bw()

Section 17. Imagistic experience: sub-scales

Below are the individual Cronbach’s alpha for each of the seven subscales, followed by their individual visualizations, followed by the correlation between the seven sub-scales

Negative affect about event

Display code

## Positive affect:
negative_affect <- cbind.data.frame(ds$event_negative_affect)
negative_affect <- na.omit(negative_affect)

## Visualization:
ds$event_negative_affect <- (ds$event_negative_affect)
summary(ds$event_negative_affect)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     2.0     6.0     5.1     7.0     7.0      20 

Display code

ds %>% drop_na(event_negative_affect)%>%
ggplot(aes(x = event_negative_affect))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 5.10", 
                 xintercept = 5.10, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Negative affect about event", 
       y = "Frequency", 
       title = "Negative affect about event: Tanzania")+
  theme_bw()

Positive affect about event

Display code

## Positive affect:
positive_affect <- cbind.data.frame(ds$event_positive_affect)
positive_affect <- na.omit(positive_affect)

## Visualization:
ds$event_positive_affect <- (ds$event_positive_affect)
summary(ds$event_positive_affect)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     1.0     3.0     3.7     7.0     7.0      22 

Display code

ds %>% drop_na(event_positive_affect)%>%
ggplot(aes(x = event_positive_affect))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 3.70", 
                 xintercept = 3.70, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Positive affect about event", 
       y = "Frequency", 
       title = "Positive affect about event: Tanzania")+
  theme_bw()

Episodic recall of event

Display code

## Episodic recall:
episodic_recall <- cbind.data.frame(ds$event_episodic_recall_01, ds$event_episodic_recall_02)
episodic_recall <- na.omit(episodic_recall)

names(episodic_recall) <- sub('ds\\$event\\_', '', names(episodic_recall))

mtx <- cor(episodic_recall[, c(1:2)])

## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

## Reliability:
alph01 <- psych::alpha(episodic_recall)
summary(alph01)

Reliability analysis   
 raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
      0.65      0.65    0.49      0.49 1.9 0.038    6 1.1     0.49

Display code

## Visualization:
ds$event_episodic_recall <- ((ds$event_episodic_recall_01+ds$event_episodic_recall_02)/2)
summary(ds$event_episodic_recall)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     5.5     6.5     6.0     7.0     7.0      29 

Display code

ds %>% drop_na(event_episodic_recall)%>%
ggplot(aes(x = event_episodic_recall))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 6.00", 
                 xintercept = 6.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Episodic recall of event", 
       y = "Frequency", 
       title = "Episodic recall of event: Tanzania")+
  theme_bw()

Shared perception of event

Display code

## Shared perception:
shared_perception <- cbind.data.frame(ds$event_shared_perception_01, ds$event_shared_perception_02)
shared_perception <- na.omit(shared_perception)

names(shared_perception) <- sub('ds\\$event\\_', '', names(shared_perception))

mtx <- cor(shared_perception[, c(1:2)])

## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

## Reliability:
alph01 <- psych::alpha(shared_perception)
summary(alph01)

Reliability analysis   
 raw_alpha std.alpha G6(smc) average_r  S/N   ase mean  sd median_r
      0.49      0.49    0.33      0.33 0.98 0.056  5.3 1.4     0.33

Display code

## Visualization:
ds$event_shared_perception <- ((ds$event_shared_perception_01+ds$event_shared_perception_02)/2)
summary(ds$event_shared_perception)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     4.5     5.5     5.3     6.5     7.0      23 

Display code

ds %>% drop_na(event_shared_perception)%>%
ggplot(aes(x = event_shared_perception))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 5.30", 
                 xintercept = 5.30, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Shared perception of event", 
       y = "Frequency", 
       title = "Shared perception of event: Tanzania")+
  theme_bw()

Reflection of event

Display code

## Reflection:
event_reflection <- cbind.data.frame(ds$event_event_reflection_01, ds$event_event_reflection_02)
event_reflection <- na.omit(event_reflection)

names(event_reflection) <- sub('ds\\$event\\_', '', names(event_reflection))

mtx <- cor(event_reflection[, c(1:2)])

## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

## Reliability:
alph01 <- psych::alpha(event_reflection)
summary(alph01)

Reliability analysis   
 raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
      0.53      0.53    0.36      0.36 1.1 0.053  5.2 1.5     0.36

Display code

## Visualization:
ds$event_event_reflection <- ((ds$event_event_reflection_01+ds$event_event_reflection_02)/2)
summary(ds$event_event_reflection)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     4.0     5.5     5.2     6.5     7.0      30 

Display code

ds %>% drop_na(event_event_reflection)%>%
ggplot(aes(x = event_event_reflection))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 5.20", 
                 xintercept = 5.20, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Reflection of event", 
       y = "Frequency", 
       title = "Reflection of event: Tanzania")+
  theme_bw()

Transformative event for individual

Display code

## Reflection:
event_transformative_indiv <- cbind.data.frame(ds$event_transformative_indiv_01, ds$event_transformative_indiv_02)

event_transformative_indiv <- na.omit(event_transformative_indiv)

names(event_transformative_indiv) <- sub('ds\\$event\\_', '', names(event_transformative_indiv))

mtx <- cor(event_transformative_indiv[, c(1:2)])

## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

## Reliability:
psych::alpha(event_transformative_indiv)

Reliability analysis   
Call: psych::alpha(x = event_transformative_indiv)

  raw_alpha std.alpha G6(smc) average_r  S/N   ase mean  sd median_r
      0.18      0.18     0.1       0.1 0.22 0.091    5 1.5      0.1

    95% confidence boundaries 
         lower alpha upper
Feldt    -0.02  0.18  0.34
Duhachek  0.00  0.18  0.36

 Reliability if an item is dropped:
                        raw_alpha std.alpha G6(smc) average_r  S/N alpha se
transformative_indiv_01     0.099       0.1    0.01       0.1 0.11       NA
transformative_indiv_02     0.101       0.1    0.01       0.1 0.11       NA
                        var.r med.r
transformative_indiv_01     0   0.1
transformative_indiv_02     0   0.1

 Item statistics 
                          n raw.r std.r r.cor r.drop mean  sd
transformative_indiv_01 323  0.74  0.74  0.23    0.1  5.1 2.0
transformative_indiv_02 323  0.74  0.74  0.23    0.1  4.9 2.1

Non missing response frequency for each item
                           1    2    3    4    5    6    7 miss
transformative_indiv_01 0.08 0.11 0.07 0.01 0.10 0.31 0.32    0
transformative_indiv_02 0.07 0.13 0.12 0.02 0.16 0.17 0.32    0

Display code

## Visualization:
ds$event_transformative_indiv <- ((ds$event_transformative_indiv_01+ds$event_transformative_indiv_02)/2)
summary(ds$event_transformative_indiv)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      1       4       5       5       6       7      29 

Display code

ds %>% drop_na(event_transformative_indiv)%>%
ggplot(aes(x = event_transformative_indiv))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 5.00", 
                 xintercept = 5.00, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Transformative event score", 
       y = "Frequency", 
       title = "Transformative event for individual: Tanzania")+
  theme_bw()

Transformative event for group

Display code

event_transformative_group <- cbind.data.frame(ds$event_transformative_group_01, ds$event_transformative_group_02)
event_transformative_group <- na.omit(event_transformative_group)

names(event_transformative_group) <- sub('ds\\$event\\_', '', names(event_transformative_group))

mtx <- cor(event_transformative_group[, c(1:2)])

## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Display code

## Reliability:
alph01 <- psych::alpha(event_transformative_group)
summary(alph01)

Reliability analysis   
 raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
       0.7       0.7    0.54      0.54 2.3 0.033  4.5 1.8     0.54

Display code

## Visualization:
ds$event_transformative_group <- ((ds$event_transformative_group_01+ds$event_transformative_group_02)/2)
summary(ds$event_transformative_group)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    1.0     3.0     4.5     4.5     6.0     7.0      21 

Display code

ds %>% drop_na(event_transformative_group)%>%
ggplot(aes(x = event_transformative_group))+
  geom_histogram(color = "black",
                 fill = "gray",
                 bins = 50)+
  geom_textvline(label = "Mean = 4.50", 
                 xintercept = 4.50, 
                 vjust = 1.1, 
                 lwd = 1.05, 
                 linetype = 2)+
  labs(x = "Transformative event score", 
       y = "Frequency", 
       title = "Transformative event for group: Tanzania")+
  theme_bw()

Correlation between seven subscales

Display code

imagistic_subscales <- cbind.data.frame(ds$event_positive_affect, ds$event_negative_affect,
                                        ds$event_episodic_recall, 
                                        ds$event_shared_perception, ds$event_event_reflection, 
                                        ds$event_transformative_indiv, ds$event_transformative_group)


imagistic_subscales <- na.omit(imagistic_subscales)

names(imagistic_subscales) <- sub('ds\\$event\\_', '', names(imagistic_subscales))

mtx <- cor(imagistic_subscales[, c(1:7)])

## Correlation plot:
corrplot(mtx, method = "number", number.cex = 0.7,
         col=c("white", "darkred", "red",
               "darkgrey", "blue", "darkblue"))

Section 18. Regression models: Relationship between imagistic experience and IG/OG fusion/identification

Display code

## Four regression models predicting IG/OG Fusion/Identification:

lm01 <- lm(IG_Fusion~event_positive_affect+event_negative_affect+event_episodic_recall+event_shared_perception+event_event_reflection+
event_transformative_indiv+event_transformative_group+Age+Female+Married+`SES-`,
           data = ds)

lm02 <- lm(IG_Identification~event_positive_affect+event_negative_affect+event_episodic_recall+event_shared_perception+event_event_reflection+
event_transformative_indiv+event_transformative_group+Age+Female+Married+`SES-`,
           data = ds)

lm03 <- lm(OG_Bonds~event_positive_affect+event_negative_affect+event_episodic_recall+event_shared_perception+event_event_reflection+
event_transformative_indiv+event_transformative_group+Age+Female+Married+`SES-`,
           data = ds)

Display code

## Tabulated results:

stargazer(lm01, lm02,
          lm03, 
          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)
event_positive_affect	-0.043	-0.092***	0.210***
	(0.025)	(0.024)	(0.048)
event_negative_affect	-0.001	-0.057*	0.270***
	(0.026)	(0.025)	(0.050)
event_episodic_recall	0.140**	0.120*	-0.180*
	(0.048)	(0.046)	(0.089)
event_shared_perception	-0.012	-0.001	-0.033
	(0.039)	(0.038)	(0.074)
event_event_reflection	0.012	-0.027	0.100
	(0.039)	(0.037)	(0.076)
event_transformative_indiv	0.017	0.039	0.026
	(0.040)	(0.039)	(0.077)
event_transformative_group	-0.025	-0.063	0.071
	(0.035)	(0.034)	(0.067)
Age	0.00000	0.00001	0.00003
	(0.00002)	(0.00002)	(0.00004)
Female	0.055	0.060	0.190
	(0.100)	(0.098)	(0.200)
Married	-0.004	-0.025	0.320
	(0.100)	(0.100)	(0.200)
`SES-`Middle	-0.120	0.012	0.420
	(0.120)	(0.120)	(0.240)
`SES-`Upper middle	-0.064	-0.025	0.280
	(0.210)	(0.200)	(0.400)
`SES-`Upper	-0.370	0.640	-0.650
	(0.510)	(0.490)	(0.930)
Constant	5.600***	6.300***	0.910
	(0.370)	(0.350)	(0.690)
Observations	281	281	256
R2	0.077	0.150	0.230
Adjusted R2	0.032	0.110	0.190
Residual Std. Error	0.850 (df = 267)	0.810 (df = 267)	1.500 (df = 242)
F Statistic	1.700 (df = 13; 267)	3.600*** (df = 13; 267)	5.700*** (df = 13; 242)
Note:	*p<0.05; **p<0.01; ***p<0.001