################################## 
## Christopher Claassen. Group Entitlement, Anger and Participation in Intergroup Violence
## R code
##################################



library(weights)
library(psych)
library(car)
library(polycor)
library(Amelia)
library(nnet)
library(mediation)
library(survey)
library(mitools)
library(ltm)

# set working directory (change dirname to your directory)
setwd("dirname")

## raw data 
alex <- read.csv("participation_intergroup_violence_data.csv", header=TRUE)



######### Initial recodes



# National ID
alex$Idimp <- Recode(alex$Idimp, "1=4;2=3;3=2;9=2;4=1;8=NA")
alex$Iddiff <- Recode(alex$Iddiff, "9=3;8=NA")
alex$Idinsult <- Recode(alex$Idinsult, "1=5;2=4;3=3;9=3;4=2;5=1;8=NA")

# Support for violence
alex$ViolNorm1 <- Recode(alex$ViolNorm1, "1=5;2=4;3=3;9=3;4=2;5=1;8=NA")
alex$ViolNorm2 <- Recode(alex$ViolNorm2, "9=3;8=NA")

# Environmental blame
alex$PeerBlameJ <- Recode(alex$PeerBlameJ, "1=4;2=3;3=2;4=1;8=NA;9=2")
alex$PeerBlameH <- Recode(alex$PeerBlameH, "1=4;2=3;3=2;4=1;8=NA;9=2")
alex$PeerBlameD <- Recode(alex$PeerBlameD, "1=4;2=3;3=2;4=1;8=NA;9=2")
alex$NoMeeting <- Recode(alex$NoMeeting, "2:3=2;4:5=3;6:100=4;8=NA;9=0")
alex$LeadBlameJ <- Recode(alex$LeadBlameJ, "1=4;2=3;3=2;4=1;8=NA;9=2")
alex$LeadBlameJ <- ifelse(alex$NoMeeting == 0, 1, alex$LeadBlameJ)
alex$LeadBlameH <- Recode(alex$LeadBlameH, "1=4;2=3;3=2;4=1;8=NA;9=2")
alex$LeadBlameH <- ifelse(alex$NoMeeting == 0, 1, alex$LeadBlameH)
alex$LeadBlameD <- Recode(alex$LeadBlameD, "1=4;2=3;3=2;4=1;8=NA;9=2")
alex$LeadBlameD <- ifelse(alex$NoMeeting == 0, 1, alex$LeadBlameD)

# Emotion: jobs
alex$EmotJob1 <- Recode(as.factor(alex$EmotJob1), "1='happy'; 2='proud'; 3='angry'; 
                        4='irritated'; 5='jealous'; 6='worried'; 7='ashamed'; 8='disappointed'; 
                        9='no feeling'; 88=NA; 99='no feeling'")
alex$EmotJob1 <- relevel(alex$EmotJob1, "no feeling")
alex$DegrEmotJ <- Recode(alex$DegrEmotJ, "1=4;2=3;3=2;4=1;8=NA;9=1")
alex$DegrEmotJ <- ifelse(alex$EmotJob1 == 9, 0, alex$DegrEmotJ)
alex$EmotJob2 <- Recode(as.factor(alex$EmotJob2), "1='happy'; 2='proud'; 3='angry'; 
                        4='irritated'; 5='jealous'; 6='worried'; 7='ashamed'; 8='disappointed'; 
                        9='no feeling'; 88=NA; 99='no feeling'")
alex$EmotJob2 <- relevel(alex$EmotJob2, "no feeling")

# Eotion: housing
alex$EmotHous1 <- Recode(as.factor(alex$EmotHous1), "1='happy'; 2='proud'; 3='angry'; 
                         4='irritated'; 5='jealous'; 6='worried'; 7='ashamed'; 8='disappointed'; 
                         9='no feeling'; 88=NA; 99='no feeling'")
alex$EmotHous1 <- relevel(alex$EmotHous1, "no feeling")
alex$DegrEmotH <- Recode(alex$DegrEmotH, "1=4;2=3;3=2;4=1;8=NA;9=1")
alex$DegrEmotH <- ifelse(alex$EmotHous1 == 9, 0, alex$DegrEmotH)
alex$EmotHous2 <- Recode(as.factor(alex$EmotHous2), "1='happy'; 2='proud'; 3='angry'; 
                         4='irritated'; 5='jealous'; 6='worried'; 7='ashamed'; 8='disappointed'; 
                         9='no feeling'; 88=NA; 99='no feeling'")
alex$EmotHous2 <- relevel(alex$EmotHous2, "no feeling")

# Group entitlement and endowment 
alex$SAPosNow <- Recode(alex$SAPosNow, "8=NA;9=NA")
alex$SAPosNorm <- Recode(alex$SAPosNorm, "8=NA;9=NA")
alex$ForPosNow <- Recode(alex$ForPosNow, "8=NA;9=NA")
alex$ForPosNorm <- Recode(alex$ForPosNorm, "8=NA;9=NA")

# Emotion: group entitlement
alex$EmotPos1 <- Recode(as.factor(alex$EmotPos1), "1='happy'; 2='proud'; 3='angry'; 
                        4='irritated'; 5='jealous'; 6='worried'; 7='ashamed'; 8='disappointed'; 
                        9='no feeling'; 88=NA; 99='no feeling'")
alex$EmotPos1 <- relevel(alex$EmotPos1, "no feeling")
alex$DegrEmotPos <- Recode(alex$DegrEmotPos, "1=4;2=3;3=2;4=1;8=NA;9=1")
alex$DegrEmotPos <- ifelse(alex$EmotPos1 == 9, 0, alex$DegrEmotPos)
alex$EmotPos2 <- Recode(as.factor(alex$EmotPos2), "1='happy'; 2='proud'; 3='angry'; 
                        4='irritated'; 5='jealous'; 6='worried'; 7='ashamed'; 8='disappointed'; 
                        9='no feeling'; 88=NA; 99='no feeling'")
alex$EmotPos2 <- relevel(alex$EmotPos2, "no feeling")

# Participation in 08
alex$LiveAlex08 <- Recode(alex$LiveAlex08, "2=0; 8=NA; 9=0")
alex$ParticToyi <- Recode(alex$ParticToyi, "3=0; 2=1; 1=1; NA=0")
alex$ParticIntim <- Recode(alex$ParticIntim, "3=0; 2=1; 1=1; NA=0")
alex$ParticLoot <- Recode(alex$ParticLoot, "3=0; 2=1; 1=1; NA=0")
alex$ParticHarm <- Recode(alex$ParticHarm, "3=0; 2=1; 1=1; NA=0")
alex$ParticDestr <- Recode(alex$ParticDestr, "3=0; 2=1; 1=1; NA=0")
alex$ParticHelp <- Recode(alex$ParticHelp, "3=0; 2=-1; 1=-1; NA=-1")

# Behavioural intentions
alex$IntendToyi <- Recode(alex$IntendToyi, "1=4;2=3;3=2;4=1;8=NA;9=2")
alex$IntendToyi2 <- Recode(alex$IntendToyi2, "1=4;2=3;3=2;4=1;8=NA;9=2")
alex$IntendHit <- Recode(alex$IntendHit, "1=4;2=3;3=2;4=1;8=NA;9=2")
alex$IntendDestr <- Recode(alex$IntendDestr, "1=4;2=3;3=2;4=1;8=NA;9=2")
alex$IntendHelp <- Recode(alex$IntendHelp, "8=NA;9=3")

# Demographics
alex$HousingType <- Recode(as.factor(alex$HousingType), "1='shack'; 2='back shack'; 3='flat'; 
                           4='hostel'; 5='house'; 6='back room'")
alex$AttendCh <- Recode(alex$AttendCh, "1=4;2=3;3=2;4=1;8=NA;9=2")
alex$Married <- Recode(as.factor(alex$Married), "1='married'; 2='relationship'; 3='single'; 
                       4='div/wid'; 8=NA; 9=NA")
alex$YrsAlex <- ifelse(alex$YrsAlex==77,  alex$Age, alex$YrsAlex)
alex$Working <- Recode(as.factor(alex$Working), "1='empl-ft'; 2='empl-pt'; 3='empl-self'; 
                       4='nw-hwf'; 5='nw-stu'; 6='nw-ret'; 7='unemp-lfw'; 8='unemp-nlfw'; 88=NA; 
                       99=NA")
alex$JobR <- Recode(as.factor(alex$JobR), "1='never worked';8=NA")
alex$Trader <- ifelse(alex$Job == "shopkeeper", 1, 0)
alex$Lang <- Recode(as.factor(alex$Lang), "1='zulu'; 2='xhosa'; 3='sotho'; 4='pedi'; 5='tsonga'; 
                    6='venda'; 7='swazi'; 8='tswana'; 9='ndebele'; 10='english'; 11='other'; 88=NA; 
                    99=NA")
alex$Party <- Recode(as.factor(alex$Party), "1='da'; 2='ifp'; 3='anc'; 4='pac'; 5='cope'; 
                     6='udm'; 7='other'; 8='none'; 88=NA; 99='none'")
alex$Education <- Recode(as.factor(alex$Education), "1='ps'; 2='some hs'; 3='matric'; 
                         4='uni deg'; 5='tech dip'; 6='tc dip'; 7='other dip'; 8='artis cert'; 
                         88=NA; 99=NA")
alex$Ownership <- Recode(as.factor(alex$Ownership), "1='own'; 2='rent'; 3='live'; 88=NA; 99=NA")
alex$HaveElect <- Recode(alex$HaveElect, "2=1; 1=0; 8=NA; 9=NA")
alex$HaveTV <- Recode(alex$HaveTV, "2=1; 1=0; 8=NA; 9=0")
alex$HaveCell <- Recode(alex$HaveCell, "2=1; 1=0; 8=NA; 9=0")
alex$HaveFridge <- Recode(alex$HaveFridge, "2=1; 1=0; 8=NA; 9=0")
alex$HaveMicro <- Recode(alex$HaveMicro, "2=1; 1=0; 8=NA; 9=0")
alex$HaveComp <- Recode(alex$HaveComp, "2=1; 1=0; 8=NA; 9=0")
alex$HaveCar <- Recode(alex$HaveCar, "2=1; 1=0; 8=NA; 9=0")



######## Imputation


a.out <- amelia(alex, m = 5,
                idvars = c("ID", "Hhsize", "HousingType", "JobR", "Block"),
                noms = c("EmotJob1", "EmotJob2", "EmotHous1", "EmotHous2", "EmotPos1", 
                         "EmotPos2", "LiveAlex08", "Married", "Working", "Lang", "Party",
                         "Education", "Ownership", "ParticToyi", "ParticIntim", "ParticLoot",
                         "ParticHelp", "ParticHarm", "ParticDestr", "Gender"),
                ords = c("Idimp", "Iddiff", "Idinsult", "ViolNorm1", "ViolNorm2", "PeerBlameJ",
                         "PeerBlameH", "PeerBlameD", "LeadBlameJ", "LeadBlameH", "LeadBlameD",
                         "DegrEmotJ", "DegrEmotH", "SAPosNow", "SAPosNorm", "ForPosNow", 
                         "ForPosNorm", "DegrEmotPos", "IntendToyi", "IntendToyi2",
                         "IntendHelp", "IntendHit", "IntendDestr", "AttendCh", "NoMeeting", 
                         "HaveElect", "HaveTV", "HaveCell", "HaveFridge", "HaveMicro", 
                         "HaveComp", "HaveCar"), empri = 5)

write.amelia(a.out, separate = TRUE, file.stem = "alex.imp.", format = "csv", 
             row.names = FALSE)




########## Load and recode imputed datasets



# imputed datasets
alex1 <- read.csv("alex.imp.1.csv", header=TRUE)
alex2 <- read.csv("alex.imp.2.csv", header=TRUE)
alex3 <- read.csv("alex.imp.3.csv", header=TRUE)
alex4 <- read.csv("alex.imp.4.csv", header=TRUE)
alex5 <- read.csv("alex.imp.5.csv", header=TRUE)

# list of imputed datasets
implist <- list(alex1, alex2, alex3, alex4, alex5)



#### Recode and extract data for Mplus

for (i in 1:5) {
  
  # Intend trichotomous
  
  implist[[i]]$IntToyi2.tri <- Recode(implist[[i]]$IntendToyi2, "1=0; 2=1; 3=1; 4=2")
  implist[[i]]$IntToyi.tri <- Recode(implist[[i]]$IntendToyi, "1=0; 2=1; 3=1; 4=2")
  implist[[i]]$IntHelp.tri <- Recode(implist[[i]]$IntendHelp, "1=0; 2=1; 3=1; 4=2")
  implist[[i]]$IntHit.tri <- Recode(implist[[i]]$IntendHit, "1=0; 2=1; 3=1; 4=2")
  implist[[i]]$IntDestr.tri <- Recode(implist[[i]]$IntendDestr, "1=0; 2=1; 3=1; 4=2")
  
  # recode emotions
  
  # anger
  implist[[i]]$AngH <- ifelse(implist[[i]]$EmotHous1 == "angry" | 
                                implist[[i]]$EmotHous2 == "angry", 1, 0)
  implist[[i]]$AngJ <- ifelse(implist[[i]]$EmotJob1 == "angry" | 
                                implist[[i]]$EmotJob2 == "angry", 1, 0)
  implist[[i]]$AngP <- ifelse(implist[[i]]$EmotPos1 == "angry" | 
                                implist[[i]]$EmotPos2 == "angry", 1, 0)
  implist[[i]]$IrrH <- ifelse(implist[[i]]$EmotHous1 == "irritated" | 
                                implist[[i]]$EmotHous2 == "irritated", 1, 0)
  implist[[i]]$IrrJ <- ifelse(implist[[i]]$EmotJob1 == "irritated" | 
                                implist[[i]]$EmotJob2 == "irritated", 1, 0)
  implist[[i]]$IrrP <- ifelse(implist[[i]]$EmotPos1 == "irritated" | 
                                implist[[i]]$EmotPos2 == "irritated", 1, 0)
  implist[[i]]$temp1 <- ifelse(implist[[i]]$AngH == 1 | implist[[i]]$IrrH == 1, 1, 0)
  implist[[i]]$temp2 <- ifelse(implist[[i]]$AngJ == 1 | implist[[i]]$IrrJ == 1, 1, 0)
  implist[[i]]$temp3 <- ifelse(implist[[i]]$AngP == 1 | implist[[i]]$IrrP == 1, 1, 0)
  implist[[i]]$AngerCnt <- rowSums(implist[[i]][,c("temp1", "temp2", "temp3")])
  
  # anxiety
  implist[[i]]$AnxH <- ifelse(implist[[i]]$EmotHous1 == "worried" | 
                                implist[[i]]$EmotHous2 == "worried", 1, 0)
  implist[[i]]$AnxJ <- ifelse(implist[[i]]$EmotJob1 == "worried" | 
                                implist[[i]]$EmotJob2 == "worried", 1, 0)
  implist[[i]]$AnxP <- ifelse(implist[[i]]$EmotPos1 == "worried" | 
                                implist[[i]]$EmotPos2 == "worried", 1, 0)
  implist[[i]]$AnxCnt <- rowSums(implist[[i]][,c("AnxH", "AnxJ", "AnxP")])
  
  # shame
  implist[[i]]$ShaH <- ifelse(implist[[i]]$EmotHous1 == "ashamed" | 
                                implist[[i]]$EmotHous2 == "ashamed", 1, 0)
  implist[[i]]$ShaJ <- ifelse(implist[[i]]$EmotJob1 == "ashamed" | 
                                implist[[i]]$EmotJob2 == "ashamed", 1, 0)
  implist[[i]]$ShaP <- ifelse(implist[[i]]$EmotPos1 == "ashamed" | 
                                implist[[i]]$EmotPos2 == "ashamed", 1, 0)
  implist[[i]]$ShameCnt <- rowSums(implist[[i]][,c("ShaH", "ShaJ", "ShaP")])
  
  # jealousy
  implist[[i]]$JeaH <- ifelse(implist[[i]]$EmotHous1 == "jealous" | 
                                implist[[i]]$EmotHous2 == "jealous", 1, 0)
  implist[[i]]$JeaJ <- ifelse(implist[[i]]$EmotJob1 == "jealous" | 
                                implist[[i]]$EmotJob2 == "jealous", 1, 0)
  implist[[i]]$JeaP <- ifelse(implist[[i]]$EmotPos1 == "jealous" | 
                                implist[[i]]$EmotPos2 == "jealous", 1, 0)
  implist[[i]]$JealCnt <- rowSums(implist[[i]][,c("JeaH", "JeaJ", "JeaP")])

  # sadness
  implist[[i]]$SadH <- ifelse(implist[[i]]$EmotHous1 == "disappointed" | 
                                implist[[i]]$EmotHous2 == "disappointed", 1, 0)
  implist[[i]]$SadJ <- ifelse(implist[[i]]$EmotJob1 == "disappointed" | 
                                implist[[i]]$EmotJob2 == "disappointed", 1, 0)
  implist[[i]]$SadP <- ifelse(implist[[i]]$EmotPos1 == "disappointed" | 
                                implist[[i]]$EmotPos2 == "disappointed", 1, 0)
  implist[[i]]$SadCnt <- rowSums(implist[[i]][,c("SadH", "SadJ", "SadP")])
  
  # no emotion
  implist[[i]]$NonH <- ifelse(implist[[i]]$EmotHous1 == "no feeling" | 
                                implist[[i]]$EmotHous2 == "no feeling", 1, 0)
  implist[[i]]$NonJ <- ifelse(implist[[i]]$EmotJob1 == "no feeling" | 
                                implist[[i]]$EmotJob2 == "no feeling", 1, 0)
  implist[[i]]$NonP <- ifelse(implist[[i]]$EmotPos1 == "no feeling" | 
                                implist[[i]]$EmotPos2 == "no feeling", 1, 0)
  implist[[i]]$NoEmotCnt <- rowSums(implist[[i]][,c("NonH", "NonJ", "NonP")])
  
  # happiness
  implist[[i]]$HapH <- ifelse(implist[[i]]$EmotHous1 == "happy" | 
                                implist[[i]]$EmotHous2 == "happy", 1, 0)
  implist[[i]]$HapJ <- ifelse(implist[[i]]$EmotJob1 == "happy" | 
                                implist[[i]]$EmotJob2 == "happy", 1, 0)
  implist[[i]]$HapP <- ifelse(implist[[i]]$EmotPos1 == "happy" | 
                                implist[[i]]$EmotPos2 == "happy", 1, 0)
  implist[[i]]$ProH <- ifelse(implist[[i]]$EmotHous1 == "proud" | 
                                implist[[i]]$EmotHous2 == "proud", 1, 0)
  implist[[i]]$ProJ <- ifelse(implist[[i]]$EmotJob1 == "proud" | 
                                implist[[i]]$EmotJob2 == "proud", 1, 0)
  implist[[i]]$ProP <- ifelse(implist[[i]]$EmotPos1 == "proud" | 
                                implist[[i]]$EmotPos2 == "proud", 1, 0)
  implist[[i]]$temp1 <- ifelse(implist[[i]]$HapH == 1 | implist[[i]]$ProH == 1, 1, 0)
  implist[[i]]$temp2 <- ifelse(implist[[i]]$HapJ == 1 | implist[[i]]$ProJ == 1, 1, 0)
  implist[[i]]$temp3 <- ifelse(implist[[i]]$HapP == 1 | implist[[i]]$ProP == 1, 1, 0)
  implist[[i]]$HappyCnt <- rowSums(implist[[i]][,c("temp1", "temp2", "temp3")])
  
  # strata and clusters
  implist[[i]]$Strata.Orig <-
    ifelse(implist[[i]]$Area=="flats" & implist[[i]]$Gender==1, "fl.m",
    ifelse(implist[[i]]$Area=="flats" & implist[[i]]$Gender==2, "fl.f",
    ifelse(implist[[i]]$Area=="hostels" & implist[[i]]$Gender==1, "ho.m",
    ifelse(implist[[i]]$Area=="hostels" & implist[[i]]$Gender==2, "ho.f",
    ifelse(implist[[i]]$Area=="shacks&houses" & implist[[i]]$Gender==1, "sh.m",
    ifelse(implist[[i]]$Area=="shacks&houses" & implist[[i]]$Gender==2, "sh.f", NA))))))
  
  implist[[i]]$Clusters <- as.matrix(implist[[i]]$Block)
  
  implist[[i]]$Sector <- as.factor(
    ifelse(implist[[i]]$Block == "F1", "Sector 2",
    ifelse(implist[[i]]$Block == "H1", "Sector 2",
    ifelse(implist[[i]]$Block == "O1", "Sector 2",
    ifelse(implist[[i]]$Block == "S1", "Sector 5",
    ifelse(implist[[i]]$Block == "S10", "Sector 2",
    ifelse(implist[[i]]$Block == "S13", "Sector 2",
    ifelse(implist[[i]]$Block == "S14", "Sector 6",
    ifelse(implist[[i]]$Block == "S15", "Sector 3",
    ifelse(implist[[i]]$Block == "S16", "Sector 4",
    ifelse(implist[[i]]$Block == "S17", "Sector 3",
    ifelse(implist[[i]]$Block == "S18", "Sector 6",
    ifelse(implist[[i]]$Block == "S19", "Sector 7",
    ifelse(implist[[i]]$Block == "S2", "Sector 5",
    ifelse(implist[[i]]$Block == "S20", "Sector 7",
    ifelse(implist[[i]]$Block == "S21", "Sector 4",
    ifelse(implist[[i]]$Block == "S22", "Sector 4",
    ifelse(implist[[i]]$Block == "S23", "Sector 5",
    ifelse(implist[[i]]$Block == "S4", "Sector 5",
    ifelse(implist[[i]]$Block == "S8", "Sector 2",
    ifelse(implist[[i]]$Block == "S9", "Sector 2", NA)))))))))))))))))))))
  
}


#### Construct sampling weights

adat <- implist[[1]]
adat$w2 <- sqrt(adat$Hhsize) * (dim(adat)[1] / sum(sqrt(adat$Hhsize)))
w2 <- adat$w2

w3 <- ifelse(adat$Strata.Orig=="fl.f", sum(subset(adat, Strata.Orig=="fl.f")$w2), 
      ifelse(adat$Strata.Orig=="fl.m", sum(subset(adat, Strata.Orig=="fl.m")$w2), 
      ifelse(adat$Strata.Orig=="ho.f", sum(subset(adat, Strata.Orig=="ho.f")$w2), 
      ifelse(adat$Strata.Orig=="ho.m", sum(subset(adat, Strata.Orig=="ho.m")$w2), 
      ifelse(adat$Strata.Orig=="sh.f", sum(subset(adat, Strata.Orig=="sh.f")$w2), 
      ifelse(adat$Strata.Orig=="sh.m", sum(subset(adat, Strata.Orig=="sh.m")$w2), NA ))))))

w4 <- ifelse(adat$Strata=="fl.f", .0105, 
      ifelse(adat$Strata=="fl.m", .0141, 
      ifelse(adat$Strata=="ho.f", .0172, 
      ifelse(adat$Strata=="ho.m", .0370, 
      ifelse(adat$Strata=="sh.f", .4438, 
      ifelse(adat$Strata=="sh.m", .4774, NA ))))))

w5 <- w2 / w3 * w4 * dim(adat)[1]

# adding weights to datasets
for(i in 1:5) {
  implist[[i]]$SampWt <- w5
}


#### Extract for Mplus

mplus.data <- list()
for (i in 1:5) {
  mplus.data[[i]] <- implist[[i]][,c("ID","AngH","AngJ","AngP","IrrH","IrrJ","IrrP","JeaH","JeaJ",
                                     "JeaP","ShaH","ShaJ","ShaP","SadH","SadJ","SadP","AnxH","AnxJ",
                                     "AnxP","HapH","HapJ","HapP","ProH","ProJ","ProP","NonH","NonJ",
                                     "NonP","EmotDH","EmotDJ","EmotDP","IntDestr.tri","IntHelp.tri",
                                     "IntHit.tri","IntToyi.tri","IntToyi2.tri","SampWt","Clusters")]
  mplus.data[[i]]$Clusters <- as.numeric(mplus.data[[i]]$Clusters)
}
write.table(mplus.data[[1]], file = "alex.mplus.1.dat", sep = "\t", row.names = FALSE, 
            col.names = FALSE)
write.table(mplus.data[[2]], file = "alex.mplus.2.dat", sep = "\t", row.names = FALSE, 
            col.names = FALSE)
write.table(mplus.data[[3]], file = "alex.mplus.3.dat", sep = "\t", row.names = FALSE, 
            col.names = FALSE)
write.table(mplus.data[[4]], file = "alex.mplus.4.dat", sep = "\t", row.names = FALSE, 
            col.names = FALSE)
write.table(mplus.data[[5]], file = "alex.mplus.5.dat", sep = "\t", row.names = FALSE, 
            col.names = FALSE)



#### IRT model for SES

data.ses <- implist[[1]][,c("HaveTV","HaveCell","HaveFridge","HaveMicro","HaveComp","HaveCar")]
ses.fs <- factor.scores(ltm(data.ses ~ z1))$score.dat[,-c(7,8,10)]
colnames(ses.fs)[7] <- "SESIRT"
implist[[1]] <- merge(implist[[1]], ses.fs)

data.ses <- implist[[2]][,c("HaveTV","HaveCell","HaveFridge","HaveMicro","HaveComp","HaveCar")]
ses.fs <- factor.scores(ltm(data.ses ~ z1))$score.dat[,-c(7,8,10)]
colnames(ses.fs)[7] <- "SESIRT"
implist[[2]] <- merge(implist[[2]], ses.fs)

data.ses <- implist[[3]][,c("HaveTV","HaveCell","HaveFridge","HaveMicro","HaveComp","HaveCar")]
ses.fs <- factor.scores(ltm(data.ses ~ z1))$score.dat[,-c(7,8,10)]
colnames(ses.fs)[7] <- "SESIRT"
implist[[3]] <- merge(implist[[3]], ses.fs)

data.ses <- implist[[4]][,c("HaveTV","HaveCell","HaveFridge","HaveMicro","HaveComp","HaveCar")]
ses.fs <- factor.scores(ltm(data.ses ~ z1))$score.dat[,-c(7,8,10)]
colnames(ses.fs)[7] <- "SESIRT"
implist[[4]] <- merge(implist[[4]], ses.fs)

data.ses <- implist[[5]][,c("HaveTV","HaveCell","HaveFridge","HaveMicro","HaveComp","HaveCar")]
ses.fs <- factor.scores(ltm(data.ses ~ z1))$score.dat[,-c(7,8,10)]
colnames(ses.fs)[7] <- "SESIRT"
implist[[5]] <- merge(implist[[5]], ses.fs)



## Run the IRT models for participation intentions and anger in MPlus 
## Load latent variable estimates from Mplus

anger.data <- list()
anger.data[[1]] <- read.table("anger.mplus.1.dat", header = FALSE)[,8:9]
anger.data[[2]] <- read.table("anger.mplus.2.dat", header = FALSE)[,8:9]
anger.data[[3]] <- read.table("anger.mplus.3.dat", header = FALSE)[,8:9]
anger.data[[4]] <- read.table("anger.mplus.4.dat", header = FALSE)[,8:9]
anger.data[[5]] <- read.table("anger.mplus.5.dat", header = FALSE)[,8:9]

int.data <- list()
int.data[[1]] <- read.table("int.mplus.1.dat", header = FALSE)[,6:7]
int.data[[2]] <- read.table("int.mplus.2.dat", header = FALSE)[,6:7]
int.data[[3]] <- read.table("int.mplus.3.dat", header = FALSE)[,6:7]
int.data[[4]] <- read.table("int.mplus.4.dat", header = FALSE)[,6:7]
int.data[[5]] <- read.table("int.mplus.5.dat", header = FALSE)[,6:7]

for (i in 1:5) {
  implist[[i]]$AngerLT <- anger.data[[i]][,2]
  implist[[i]]$IntendLT <- int.data[[i]][,2]
}


#### More recodes

# create 0-1 function
zero.one <- function (varname){
  (varname - min(varname)) / (max(varname) - min(varname))
}

for (i in 1:5) {
  
  implist[[i]]$AngerLT <- zero.one(implist[[i]]$AngerLT)
  implist[[i]]$AngerStd <- stdz(implist[[i]]$AngerLT, weight = implist[[i]]$SampWt)
  implist[[i]]$JealCnt <- zero.one(implist[[i]]$JealCnt)
  implist[[i]]$JealStd <- stdz(implist[[i]]$JealCnt, weight = implist[[i]]$SampWt)
  implist[[i]]$AnxCnt <- zero.one(implist[[i]]$AnxCnt)
  implist[[i]]$AnxStd <- stdz(implist[[i]]$AnxCnt, weight = implist[[i]]$SampWt)
  implist[[i]]$SadCnt <- zero.one(implist[[i]]$SadCnt)
  implist[[i]]$SadStd <- stdz(implist[[i]]$SadCnt, weight = implist[[i]]$SampWt)
  implist[[i]]$ShameCnt <- zero.one(implist[[i]]$ShameCnt)
  implist[[i]]$ShameStd <- stdz(implist[[i]]$ShameCnt, weight = implist[[i]]$SampWt)
  implist[[i]]$NoEmotCnt <- zero.one(implist[[i]]$NoEmotCnt)
  implist[[i]]$NoEmotStd <- stdz(implist[[i]]$NoEmotCnt, weight = implist[[i]]$SampWt)
  implist[[i]]$HappyCnt <- zero.one(implist[[i]]$HappyCnt)
  implist[[i]]$HappyStd <- stdz(implist[[i]]$HappyCnt, weight = implist[[i]]$SampWt)
  
  implist[[i]]$PercIneq <- zero.one(implist[[i]]$ForPosNow - implist[[i]]$SAPosNow)
  implist[[i]]$EntitViol <- zero.one(implist[[i]]$SAPosNorm - implist[[i]]$SAPosNow) 
  implist[[i]]$NormGrpPos <- zero.one(implist[[i]]$SAPosNorm - implist[[i]]$ForPosNorm) 
  implist[[i]]$GrpPosViol <- zero.one(implist[[i]]$PercIneq + implist[[i]]$NormGrpPos) 
  implist[[i]]$IngViolEnt <- zero.one(implist[[i]]$SAPosNorm - implist[[i]]$SAPosNow)
  implist[[i]]$OutgViolEnt <- zero.one(implist[[i]]$ForPosNow - implist[[i]]$ForPosNorm)
  
  implist[[i]]$SAPosNowStd <- stdz(implist[[i]]$SAPosNow, weight = implist[[i]]$SampWt)
  implist[[i]]$SAPosNormStd <- stdz(implist[[i]]$SAPosNorm, weight = implist[[i]]$SampWt)
  implist[[i]]$ForPosNowStd <- stdz(implist[[i]]$ForPosNow, weight = implist[[i]]$SampWt)
  implist[[i]]$ForPosNormStd <- stdz(implist[[i]]$ForPosNorm, weight = implist[[i]]$SampWt)
  implist[[i]]$PercIneqStd <- stdz(implist[[i]]$PercIneq, weight = implist[[i]]$SampWt)
  implist[[i]]$DeservIneqStd <- stdz(implist[[i]]$NormGrpPos, weight = implist[[i]]$SampWt)
  implist[[i]]$IngViolEntStd <- stdz(implist[[i]]$IngViolEnt, weight = implist[[i]]$SampWt)
  implist[[i]]$OutgViolEntStd <- stdz(implist[[i]]$OutgViolEnt, weight = implist[[i]]$SampWt)
  implist[[i]]$GEVStd <- stdz(implist[[i]]$GrpPosViol, weight = implist[[i]]$SampWt)
  
  implist[[i]]$ViolNorm <- zero.one(implist[[i]]$ViolNorm1 + implist[[i]]$ViolNorm2)
  implist[[i]]$NatID <- rowSums(implist[[i]][,c("Idimp", "Iddiff", "Idinsult")])
  implist[[i]]$NatID <- zero.one(implist[[i]]$NatID)
  implist[[i]]$PeerBlame <- rowSums(implist[[i]][,c("PeerBlameJ", "PeerBlameH", "PeerBlameD")])
  implist[[i]]$PeerBlame <- zero.one(implist[[i]]$PeerBlame)
  implist[[i]]$IntendLT <- zero.one(implist[[i]]$IntendLT)
  implist[[i]]$IntendStd <- stdz(implist[[i]]$IntendLT, weight = implist[[i]]$SampWt)
  implist[[i]]$Educ.Scal <- Recode(implist[[i]]$Education, "'ps'=1; 'some hs'=2; 'matric'=3; 
                                  'artis cert'=3; else=4")
  implist[[i]]$Educ.Scal <- zero.one(implist[[i]]$Educ.Scal)
  implist[[i]]$Age <- round(implist[[i]]$Age, 0)
  implist[[i]]$YrsAlex <- round(implist[[i]]$YrsAlex, 0)
  implist[[i]]$Age <- zero.one(implist[[i]]$Age)
  implist[[i]]$YrsAlex <- zero.one(implist[[i]]$YrsAlex)
  implist[[i]]$NoMeeting.Di <- Recode(implist[[i]]$NoMeeting, "0=0; 1:50=1")
  implist[[i]]$NoMeeting <- zero.one(implist[[i]]$NoMeeting)
  implist[[i]]$SESIRT <- zero.one(implist[[i]]$SESIRT)
  implist[[i]]$Partic08 <- ifelse( 
    (implist[[i]]$ParticIntim == 1 | implist[[i]]$ParticLoot == 1 | 
     implist[[i]]$ParticHarm == 1 | implist[[i]]$ParticDestr == 1), 2, 
    ifelse(implist[[i]]$ParticToyi == 1, 1, 0))
  implist[[i]]$Partic08.any <- ifelse(implist[[i]]$Partic08 > 0, 1, 0)
  implist[[i]]$Male <- Recode(implist[[i]]$Gender, "1=1;2=0;9=NA")
  implist[[i]]$Working.tri <- Recode(implist[[i]]$Working, "'empl-ft'='empl'; 'empl-pt'='empl'; 
                                    'empl-self'='empl'; 'nw-hwf'='nw'; 'nw-stu'='nw'; 
                                    'nw-ret'='nw'; 'unemp-lfw'='unemp'; 'unemp-nlfw'='nw'; 
                                    88=NA; 99=NA")
  implist[[i]]$Unemp <- ifelse(implist[[i]]$Working.tri == "unemp", 1, 0)
  implist[[i]]$NILF <- ifelse(implist[[i]]$Working.tri == "nw", 1, 0)
  implist[[i]]$Party.tri <- Recode(implist[[i]]$Party, "'cope'='other'; 'da'='other'; 'ifp'='other'; 
                                  'pac'='other'")
  
}





################# Data analysis





  
#### Set up survey designs

## 2-level cluster designs

for (i in 1:5){
  implist[[i]]$nclust <- 120
  implist[[i]]$npeop <- 193940 / implist[[i]]$nclust
}
options("survey.lonely.psu" = "adjust")

# imputed datasets
impdata <- imputationList(implist)
designs <- svydesign(ids = ~Clusters + ID, fpc = ~nclust + npeop, data = impdata, weights = ~SampWt)



##### Regression Models

## Step 1: Participation intentions

# Model for paper 
form1 <-  (IntendLT ~ AngerLT + GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab1 <- summary(MIcombine(with(designs, svyglm(form1))))[1:2]
round(cbind(tab1[1:2], (tab1[1]/tab1[2])), 2)

# With partic08
form1.2 <-  (IntendLT ~ AngerLT + GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + Partic08.any + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab1.2 <- summary(MIcombine(with(designs, svyglm(form1.2))))[1:2]
round(cbind(tab1.2[1:2], (tab1.2[1]/tab1.2[2])), 2)

# With anger + other emotions
form1.3 <-  (IntendLT ~ AngerLT + JealCnt + AnxCnt + SadCnt + ShameCnt + GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab1.3 <- summary(MIcombine(with(designs, svyglm(form1.3))))[1:2]
round(cbind(tab1.3[1:2], (tab1.3[1]/tab1.3[2])), 2)

# Model with everything but anger
form1.4 <- (IntendLT ~ GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab1.4 <- summary(MIcombine(with(designs, svyglm(form1.4))))[1:2]
round(cbind(tab1.4[1:2], (tab1.4[1]/tab1.4[2])), 2)

# GEV disagg-4
form1.5 <-  (IntendLT ~ AngerLT + SAPosNow + SAPosNorm + ForPosNow + ForPosNorm + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab1.5 <- summary(MIcombine(with(designs, svyglm(form1.5))))[1:2]
round(cbind(tab1.5[1:2], (tab1.5[1]/tab1.5[2])), 2)

# Without meeting attendance
form1.6 <-  (IntendLT ~ AngerLT + GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab1.6 <- summary(MIcombine(with(designs, svyglm(form1.6))))[1:2]
round(cbind(tab1.6[1:2], (tab1.6[1]/tab1.6[2])), 2)

# Wthout peer blame
form1.7 <-  (IntendLT ~ AngerLT + GrpPosViol + ViolNorm + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab1.7 <- summary(MIcombine(with(designs, svyglm(form1.7))))[1:2]
round(cbind(tab1.7[1:2], (tab1.7[1]/tab1.7[2])), 2)

# Interaction anger and violnorm
form1.9 <-  (IntendLT ~ AngerLT + GrpPosViol + ViolNorm + AngerLT:ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab1.9 <- summary(MIcombine(with(designs, svyglm(form1.9))))[1:2]
round(cbind(tab1.9[1:2], (tab1.9[1]/tab1.9[2])), 2)

# interaction with meetings
form1.10 <-  (IntendLT ~ AngerLT + GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector + AngerLT:NoMeeting)
tab1.10 <- summary(MIcombine(with(designs,  svyglm(form1.10))))
round(cbind(tab1.10[1:2], (tab1.10[1]/tab1.10[2])), 4)

# interaction with meetings - dichot
form1.11 <-  (IntendLT ~ AngerLT + GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting.Di + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector + AngerLT:NoMeeting.Di)
tab1.11 <- summary(MIcombine(with(designs,  svyglm(form1.11))))
round(cbind(tab1.11[1:2], (tab1.11[1]/tab1.11[2])), 4)

# interaction with previous participation
form1.12 <-  (IntendLT ~ AngerLT + GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Partic08.any + Sector + AngerLT:Partic08.any)
tab1.12 <- summary(MIcombine(with(designs,  svyglm(form1.12))))
round(cbind(tab1.12[1:2], (tab1.12[1]/tab1.12[2])), 4)

# with student, unemp, NILF
form1.14 <-  (IntendLT ~ AngerLT + GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + I(recode(Working,"'nw-stu'=1;else=0")) + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab1.14 <- summary(MIcombine(with(designs, svyglm(form1.14))))[1:2]
round(cbind(tab1.14[1:2], (tab1.14[1]/tab1.14[2])), 2)


## Anger ~ 

# Model for paper: without intlang and partic08
form2 <- (AngerLT ~ GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab2 <- summary(MIcombine(with(designs, svyglm(form2))))[1:2]
round(cbind(tab2[1:2], (tab2[1]/tab2[2])), 2)

# Just GEV
form2.5 <-  (AngerLT ~ GrpPosViol)
tab2.5 <- summary(MIcombine(with(designs, svyglm(form2.5))))[1:2]
round(cbind(tab2.5[1:2], (tab2.5[1]/tab2.5[2])), 2)

# GEV disagg-4
form2.6 <-  (AngerLT ~ SAPosNow + SAPosNorm + ForPosNow + ForPosNorm + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Sector + Educ.Scal)
tab2.6 <- summary(MIcombine(with(designs, svyglm(form2.6))))[1:2]
round(cbind(tab2.6[1:2], (tab2.6[1]/tab2.6[2])), 2)

# without peerblame
form2.7 <-  (AngerLT ~ GrpPosViol + ViolNorm + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab2.7 <- summary(MIcombine(with(designs, svyglm(form2.7))))[1:2]
round(cbind(tab2.7[1:2], (tab2.7[1]/tab2.7[2])), 2)

# without meetings
form2.8 <-  (AngerLT ~ GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab2.8 <- summary(MIcombine(with(designs, svyglm(form2.8))))[1:2]
round(cbind(tab2.8[1:2], (tab2.8[1]/tab2.8[2])), 2)

# moderated mediation?
form2.9 <-  (AngerLT ~ GrpPosViol * ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab2.9 <- summary(MIcombine(with(designs, svyglm(form2.9))))[1:2]
round(cbind(tab2.9[1:2], (tab2.9[1]/tab2.9[2])), 2)

# with partic08 
form2.10 <-  (AngerLT ~ GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + Partic08.any + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
tab2.10 <- summary(MIcombine(with(designs, svyglm(form2.10))))[1:2]
round(cbind(tab2.10[1:2], (tab2.10[1]/tab2.10[2])), 2)




#### Figures

## Intention by emotion figure 

# regressions
results <- rbind(
  mod.ang <- summary(MIcombine(with(designs, svyglm(IntendStd ~ AngerStd))))[2,1:4],
  mod.jea <- summary(MIcombine(with(designs, svyglm(IntendStd ~ JealStd))))[2,1:4],
  mod.sha <- summary(MIcombine(with(designs, svyglm(IntendStd ~ ShameStd))))[2,1:4],
  mod.anx <- summary(MIcombine(with(designs, svyglm(IntendStd ~ AnxStd))))[2,1:4],
  mod.sad <- summary(MIcombine(with(designs, svyglm(IntendStd ~ SadStd))))[2,1:4],
  mod.non <- summary(MIcombine(with(designs, svyglm(IntendStd ~ NoEmotStd))))[2,1:4],
  mod.hap <- summary(MIcombine(with(designs, svyglm(IntendStd ~ HappyStd))))[2,1:4]
)

mod.names <- c("Anger", "Jealousy", "Shame", "Anxiety", "Sadness", "No Emotion", "Happiness")

# figure
x.axis <- 1:7
pdf("intxemotreg.pdf", height = 4.5, width = 6.5)
par(mfrow = c(1, 1), mar = c(3.2, 3.8, 0.5, 0.5))
plot(results[,1], ylim = c(-1, 1), pch = 19, axes = FALSE, xlab = "", ylab = "")
segments(0.8, seq(-1, 1, by = 0.25), 7.2, seq(-1, 1, by = 0.25), col = grey(0.4), lty = 3)
segments(0.8, 0, 7.2, 0, lty = 2, col = grey(0.4))
points(results[,1], pch = 19)
segments(x.axis, results[,3], x.axis, results[,4], lwd = 1.5)
axis(2, at = seq(-1, 1, by = 0.25), las = 1, tick = TRUE, cex.axis = 0.75, mgp = c(1, .7, 0))
axis(1, at = x.axis, label = mod.names, las = 0, tick = TRUE, mgp = c(1, .4, 0), cex.axis = 0.7)
mtext("Bivariate correlation with participation intentions", side = 2, line = 2.5, cex = 0.8)
mtext("Emotional reactions to the outgroup", side = 1, line = 1.5, cex = 0.8)
dev.off()



## Robustness interaction models

form.int1 <-  (IntendLT ~ AngerLT * I(ifelse(Partic08 > 0, 1, 0)) + GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
summary(mod.int1 <- MIcombine(with(designs, svyglm(form.int1))))

# marginal effects
(me1.0 <- coef(mod.int1)[2])
(se1.0 <- sqrt(vcov(mod.int1)[2,2]))
up1.0 <- me1.0 + 1.96 * se1.0
lo1.0 <- me1.0 - 1.96 * se1.0
(me1.1 <- coef(mod.int1)[2] + coef(mod.int1)[27] * 1)
(se1.1 <- sqrt(vcov(mod.int1)[2, 2] + (1^2) * vcov(mod.int1)[27, 27] + 2 * 1 * vcov(mod.int1)[2, 27]))
up1.1 <- me1.1 + 1.96 * se1.1
lo1.1 <- me1.1 - 1.96 * se1.1

# NoMeeting
form.int2 <-  (IntendLT ~ AngerLT * I(ifelse(NoMeeting > 0, 1, 0)) + GrpPosViol + ViolNorm + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
summary(mod.int2 <- MIcombine(with(designs, svyglm(form.int2))))

# Marginal effects
(me2.0 <- coef(mod.int2)[2])
(se2.0 <- sqrt(vcov(mod.int2)[2, 2]))
up2.0 <- me2.0 + 1.96 * se2.0
lo2.0 <- me2.0 - 1.96 * se2.0
(me2.1 <- coef(mod.int2)[2] + coef(mod.int2)[26] * 1)
(se2.1 <- sqrt(vcov(mod.int2)[2, 2] + (1^2) * vcov(mod.int2)[26, 26] + 2 * 1 * vcov(mod.int2)[2, 26]))
up2.1 <- me2.1 + 1.96 * se2.1
lo2.1 <- me2.1 - 1.96 * se2.1

# graph
pdf("robust.pdf", height = 3.5, width = 6.5)
xdim <- c(0, 1)
xlabs1 <- c("No", "Yes")
xlabs2 <- c("None", "One or more")
par(mfrow = c(1, 2), mar=c(2.5, 3.25, 0.5, 0.25))
plot(xdim, c(me1.0, me1.1), type = "p", ylim = c(0,0.8), axes = FALSE, mgp = c(1.3, 0, 0), 
     ylab = "", yaxt = "n", xaxt = "n", pch = 16, xlab = "Whether Participated in 2008 Attacks", 
     cex.lab = 0.8, xlim = c(-0.3, 1.3))
axis(2, at = seq(-0.2, 1, by = .2), las = 1, tick = TRUE, cex.axis = 0.75, mgp = c(1.8, 0.7, 0))
axis(1, at = xdim, label = xlabs1, las = 0, tick = FALSE, mgp = c(1.8, 0, 0), cex.axis = 0.75)
mtext("Marginal Effect of Intergroup Anger", side = 2, line = 2, at = 0.4, cex = 0.8)
abline(h = seq(0, 0.8, by = 0.2), lty = 3, lwd = .75)
lines(c(xdim[1], xdim[1]), c(up1.0, lo1.0), lty = 1, lwd = 1.3)
lines(c(xdim[2], xdim[2]), c(up1.1, lo1.1), lty = 1, lwd = 1.3)
plot(xdim, c(me2.0, me2.1), type = "p", ylim = c(0,0.8), mgp = c(1.3, 0, 0), ylab = "", yaxt = "n", 
     xaxt = "n", cex.lab = 0.8, axes = FALSE, pch = 16, xlab = "Number of Community Meetings Attended",
     xlim = c(-0.3, 1.3))
axis(2, at = seq(-0.2, 1, by = .2), las = 1, labels = NULL, tick = TRUE, cex.axis = 0.75)
axis(1, at = xdim, label = xlabs2, las = 0, tick = FALSE, mgp = c(1, 0, 1), cex.axis = 0.75)
abline(h = seq(0, 0.8, by = 0.2), lty = 3, lwd = .75)
lines(c(xdim[1], xdim[1]), c(up2.0, lo2.0), lty = 1, lwd = 1.3)
lines(c(xdim[2], xdim[2]), c(up2.1, lo2.1), lty = 1, lwd = 1.3)
dev.off()



## Interaction between anger and norms of violence

form.int3 <-  (IntendLT ~ AngerLT * ViolNorm + GrpPosViol + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
summary(mod.int3 <- MIcombine(with(designs, svyglm(form.int3))))

form.ang <-  (AngerLT ~ ViolNorm + GrpPosViol + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
summary(mod.ang <- MIcombine(with(designs, svyglm(form.ang))))

# marginal effects of anger
x.val <- seq(0, 1, .01)
me3 <- coef(mod.int3)[2] + coef(mod.int3)[26] * x.val
se3 <- sqrt(vcov(mod.int3)[2,2] + x.val^2 * vcov(mod.int3)[26,26] + 2 * x.val * vcov(mod.int3)[2,26])
ci3 <- matrix(NA, ncol = 2, nrow = 101)
ci3[,1] <- me3 - 1.96 * se3
ci3[,2] <- me3 + 1.96 * se3  
# marginal effects of gev
me4 <- coef(mod.int3)[2] * coef(mod.ang)[3] + coef(mod.int3)[26] * coef(mod.ang)[3] * x.val
se4 <- sqrt(vcov(mod.int3)[2,2] + x.val^2 * vcov(mod.int3)[26,26] + 2 * x.val * vcov(mod.int3)[2,26])
ci4 <- matrix(NA, ncol = 2, nrow = 101)
ci4[,1] <- me4 - 1.96 * se4
ci4[,2] <- me4 + 1.96 * se4  

## graph
pdf("angxnorms.pdf", height = 4, width = 6.5)
mat.layout <- matrix(c(1, 2), nrow = 2, ncol = 1)
layout(mat.layout, heights = c(.70, .28))
# me of anger
par(mar = c(0.2, 3.25, 0, 0.5))
plot(x.val, me3, type = "l", ylim = c(0, 1.2), axes = FALSE, mgp = c(1.3, 0, 0), ylab = "", 
     yaxt = "n", xaxt = "n", xlab = "", cex.lab = 0.75, xaxs = "i", yaxs = "i")
axis(2, at = seq(0, 1, by = .25), las = 1, tick = TRUE, cex.axis = 0.7, mgp = c(2.2, 0.7, 0))
mtext("Marginal Effect of Intergroup Anger", side = 2, line = 2, at = 0.55, cex = 0.75)
abline(h = seq(0, 1, by = .25), lty = 3, lwd = .75)
lines(x.val, me3, lwd = 2)
lines(x.val, ci3[,1], col = grey(0.4), lwd = 1.5)
lines(x.val, ci3[,2], col = grey(0.4), lwd = 1.5)
abline(h = 0)
# dist of violnorm
par(mar = c(3, 3.25, 0.2, 0.5))
plot(density(alex1$ViolNorm, from = 0, to = 1, weights = alex1$SampWt),
     ylab = "", main = "", yaxt = "n", xaxt = "n", cex.lab = 0.8, xlab = "",
     lwd = 1.5, xaxs = "i", yaxs = "i")
polygon(x = c(density(alex1$ViolNorm, from = 0, to = 1, weights = alex1$SampWt)$x, 0), 
        y = c(density(alex1$ViolNorm, from = 0, to = 1, weights = alex1$SampWt)$y, 0), 
        col = grey(0.6))
axis(1, at = c(0, 1), las = 0, tick = TRUE, mgp = c(1.5, .35, 0), cex.axis = 0.7)
mtext("Normative support for violence", side = 1, line = 1, cex = 0.75)
dev.off()



## Anger by group entitlements and group endowments

# regressions
results <- rbind(
  mod1 <- summary(MIcombine(with(designs, svyglm(AngerStd ~ SAPosNowStd))))[2,1:4],
  mod2 <- summary(MIcombine(with(designs, svyglm(AngerStd ~ ForPosNowStd))))[2,1:4],
  mod3 <- summary(MIcombine(with(designs, svyglm(AngerStd ~ SAPosNormStd))))[2,1:4],
  mod4 <- summary(MIcombine(with(designs, svyglm(AngerStd ~ ForPosNormStd))))[2,1:4],
  mod5 <- summary(MIcombine(with(designs, svyglm(AngerStd ~ PercIneqStd))))[2,1:4],
  mod6 <- summary(MIcombine(with(designs, svyglm(AngerStd ~ DeservIneqStd))))[2,1:4],
  mod7 <- summary(MIcombine(with(designs, svyglm(AngerStd ~ IngViolEntStd))))[2,1:4],
  mod8 <- summary(MIcombine(with(designs, svyglm(AngerStd ~ OutgViolEntStd))))[2,1:4],
  mod9 <- summary(MIcombine(with(designs, svyglm(AngerStd ~ GEVStd))))[2,1:4]
)

mod.names <- c("Ingrp\nEnd\n", "Outgrp\nEnd\n", "Ingrp\nEnt\n", "Outgrp\nEnt\n", "Perc\nIneq\n", "Deserv\nIneq\n", "Ingrp\nViol\nEnt", "Outgrp\nViol\nEnt", "Viol\nGrp\nEnt")

# figure
x.axis <- 1:9
pdf("angxgrpend.pdf", height = 4.5, width = 6.5)
par(mfrow = c(1, 1), mar = c(3.8, 3.8, 0.5, 0.5))
plot(results[,1], ylim = c(-1, 1), pch = 19, axes = FALSE, xlab = "", ylab = "")
segments(0.8, seq(-1, 1, by = 0.25), 9.2, seq(-1, 1, by = 0.25), col = grey(0.4), lty = 3)
segments(0.8, 0, 9.2, 0, lty = 2, col = grey(0.4))
points(results[,1], pch = 19)
segments(x.axis, results[,3], x.axis, results[,4], lwd = 1.5)
axis(2, at = seq(-1, 1, by = 0.25), las = 1, tick = TRUE, cex.axis = 0.75, mgp = c(1, .7, 0))
axis(1, at = x.axis, label = mod.names, las = 0, tick = TRUE, mgp = c(1, 1.8, 0), cex.axis = 0.7)
mtext("Bivariate correlation with intergroup anger", side = 2, line = 2.8, cex = 0.8)
mtext("Group evaluations", side = 1, line = 2.5, cex = 0.8)
dev.off()

# significance of differences
cov.wt(implist[[1]][, c("AngerLT", "GrpPosViol", "IngViolEnt", "OutgViolEnt", "PercIneq", "NormGrpPos")], wt = implist[[1]]$SampWt, cor = TRUE)$cor
r.test(n = 497, r12 = 0.3220029, r13 = 0.3086733, r23 = 0.9034460)
r.test(n = 497, r12 = 0.3220029, r13 = 0.2141947, r23 = 0.7502595)
r.test(n = 497, r12 = 0.3220029, r13 = 0.17972951, r23 = 0.70000546)
r.test(n = 497, r12 = 0.3220029, r13 = 0.2895338, r23 = 0.7732075)




#### Mediation analysis and graphs

form1 <- (IntendLT ~ AngerLT + ViolNorm + GrpPosViol + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)
form2 <-  (AngerLT ~ ViolNorm + GrpPosViol + PeerBlame + LeadBlame + NatID + NILF + Unemp + Trader + InformHome + NoMeeting + SESIRT + Party.tri + Age + I(Age^2) + YrsAlex + Male + Educ.Scal + Sector)

slm.y <- with(designs, svyglm(form1))
slm.m <- with(designs, svyglm(form2))

med.out <- matrix(NA, nrow = 200, ncol = 5)
dir.out <- matrix(NA, nrow = 200, ncol = 5)
prop.out <- matrix(NA, nrow = 200, ncol = 5)
med.mods <- list()
for(i in 1:5){
  med.mods[[i]] <- mediate(slm.m[[i]], slm.y[[i]], sims = 200, boot = FALSE, treat = "GrpPosViol",
                           mediator = "AngerLT")
  med.out[,i] <- med.mods[[i]]$d1.sims
  dir.out[,i] <- med.mods[[i]]$z1.sims
  prop.out[,i] <- med.mods[[i]]$n1.sims
}

res1 <- apply(cbind(as.vector(med.out), as.vector(dir.out), as.vector(prop.out)), 2, mean)
res2 <- apply(cbind(as.vector(med.out), as.vector(dir.out), as.vector(prop.out)), 2, sd)
res3 <- apply(cbind(as.vector(med.out), as.vector(dir.out), as.vector(prop.out)), 2, quantile, 
              c(.025, .975))
res <- data.frame(Est = res1, SE = res2, CIlow = res3[1,], CIupp = res3[2,])
row.names(res) <- c("Mediated effect", "Direct effect", "Prop. mediated")
res