rm(list=ls()) ########################################################################################### # Authors: # # Title: Perception of Perception Project # # Journal: # # R Replication Code # ########################################################################################### library(dplyr) library(magrittr) library(lattice) library(sandwich) # for robust standard errors library(stargazer) library(foreign) library(gridExtra) library(psych) library(ggplot2) library(corrplot) library(sjmisc) library(GPArotation) library(nFactors) library(dplyr) library(psy) library(mediation) library(tidyverse) library(haven) library(multiwayvcov) library(lmtest) library(AER) library(ivpack) library(plm) library(lmerTest) ### Load US data #setwd("C:/Users/18735/Documents/Projects/parallel survey/analysis 10182020") dat1 <- read.csv("Parallel+Survey+2020+Wave+1+-+English_April+13,+2021_10.26.csv", encoding="UTF-8", stringsAsFactors=FALSE, na.strings=c(""," ","NA")) dat1_text <- read.csv("Parallel+Survey+2020+Wave+1+-+English_April+13,+2021_10.27_text.csv", encoding="UTF-8", stringsAsFactors=FALSE, na.strings=c(""," ","NA")) ### load chinese data dat2 <- read.csv("CFGU+--+Wave+6+--+Parallel+Chinese_April+13,+2021_10.29.csv", encoding="UTF-8", stringsAsFactors=FALSE, na.strings=c(""," ","NA")) dat2_text <- read.csv("CFGU+--+Wave+6+--+Parallel+Chinese_April+13,+2021_10.30_text.csv", encoding="UTF-8", stringsAsFactors=FALSE, na.strings=c(""," ","NA")) ### get lists of variables q.us<-dat1[c(1,2),] %>% t() q.cn<-dat2_text[c(1,2),] %>% t() ### clean the dataset. Just run all codes in the section dat1<-dat1[c(-1,-2),] dat1_text<-dat1_text[c(-1,-2),] dat2<-dat2[c(-1,-2),] dat2_text<-dat2_text[c(-1,-2),] dat1$date<-as.Date(dat1$StartDate) dat1<-dat1[dat1$date>=as.Date("2020-08-12"),] dat1_text$date<-as.Date(dat1_text$StartDate) dat1_text<-dat1_text[dat1_text$date>=as.Date("2020-08-12"),] dat2$date<-as.Date(dat2$StartDate) dat2<-dat2[dat2$date>=as.Date("2020-08-12"),] dat2_text$date<-as.Date(dat2_text$StartDate) dat2_text<-dat2_text[dat2_text$date>=as.Date("2020-08-12"),] #get those who passed attention checks data.us <- dat1 %>% filter(Q13_135 == 7 & age==age2) #data.us <-data.us[is.na(data.us$immi_1)==FALSE,] data.us.text <- dat1_text %>% filter(Q13_135 == 7 & age==age2) #data.us.text <-data.us.text[is.na(data.us.text$immi_1)==FALSE,] data.ch <- dat2 %>% filter(Q13_13 == 7 & age==age2) #data.ch <-data.ch[is.na(data.ch$media_2_9)==FALSE,] data.ch.text <- dat2_text %>% filter(Q13_13 == 7 & age==age2) #data.ch.text <-data.ch.text[is.na(data.ch.text$media_2_9)==FALSE,] ############################################################################################################## #### US Survey ############################################################################################### #### Perception ############################################################################################## ### Ddependent variable: Favorability (Q3, Q5, Q7, Q9, Q119) How favorable do you feel toward China? # 1= very favorable 6= very unfavorable data.us$favor_ch[data.us$Q3==1 | data.us$Q5==1 | data.us$Q7==1 | data.us$Q9==1 | data.us$Q119==1] <- 6 # very favorable data.us$favor_ch[data.us$Q3==2 | data.us$Q5==5 | data.us$Q7==5 | data.us$Q9==5 | data.us$Q119==4] <- 5 # somewhat favorable data.us$favor_ch[data.us$Q3==3 | data.us$Q5==6 | data.us$Q7==6 | data.us$Q9==6 | data.us$Q119==5] <- 4 # slightly favorable data.us$favor_ch[data.us$Q3==4 | data.us$Q5==7 | data.us$Q7==7 | data.us$Q9==7 | data.us$Q119==6] <- 3 # slightly unfavorable data.us$favor_ch[data.us$Q3==5 | data.us$Q5==8 | data.us$Q7==8 | data.us$Q9==8 | data.us$Q119==7] <- 2 # somewhat unfavorable data.us$favor_ch[data.us$Q3==6 | data.us$Q5==9 | data.us$Q7==9 | data.us$Q9==9 | data.us$Q119==8] <- 1 # very unfavorable hist(data.us$favor_ch) # Treatment groups data.us$treat_sd<-ifelse(is.na(data.us$Q3)==FALSE,0, ifelse(is.na(data.us$Q5)==FALSE,1, ifelse(is.na(data.us$Q7)==FALSE,2, ifelse(is.na(data.us$Q9)==FALSE,3, ifelse(is.na(data.us$Q119)==FALSE,4,NA))))) ## this is the distribution of favorability towards China among US respondents library(ggplot2) library(scales) png(file="p_antiforeign_us_control.png",width=800,height=1000, res=200,pointsize = 6,family = "SimHei", type='cairo') ggplot(data.us[data.us$treat_sd==0,], aes(x = favor_ch)) + geom_bar(aes(y = (..count..)/sum(..count..))) + geom_text(aes(y = ((..count..)/sum(..count..)), label = scales::percent((..count..)/sum(..count..))), stat = "count", vjust = -0.25) + scale_x_continuous(breaks=1:6, labels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", "Slightly favorable","Somewhat favorable", "Very favorable"), guide = guide_axis(angle = 45))+ scale_y_continuous(labels = scales::percent, limits=c(0,0.3)) + theme_bw()+ labs(y = "Percent", x = "") dev.off() ############################################################################################################## #### China Survey ############################################################################################### #### Perception ############################################################################################## ### Favorability (Q3, Q5, Q7, Q9, Q119) How favorable do you feel toward U.S.? #6= very favorable, 1=very unfavorable data.ch$favor_us[data.ch$Q3==1 | data.ch$Q5==1 | data.ch$Q7==1 | data.ch$Q9==1 | data.ch$Q119==1] <- 6 # very favorable data.ch$favor_us[data.ch$Q3==7 | data.ch$Q5==2 | data.ch$Q7==5 | data.ch$Q9==5 | data.ch$Q119==5] <- 5 # somewhat favorable data.ch$favor_us[data.ch$Q3==8 | data.ch$Q5==3 | data.ch$Q7==6 | data.ch$Q9==6 | data.ch$Q119==6] <- 4 # slightly favorable data.ch$favor_us[data.ch$Q3== 9 | data.ch$Q5==4 | data.ch$Q7==7 | data.ch$Q9==7 | data.ch$Q119==7] <- 3 # slightly unfavorable data.ch$favor_us[data.ch$Q3==10 | data.ch$Q5==5 | data.ch$Q7==8 | data.ch$Q9==8 | data.ch$Q119==8] <- 2 # somewhat unfavorable data.ch$favor_us[data.ch$Q3==11 | data.ch$Q5==6 | data.ch$Q7==9 | data.ch$Q9==9 | data.ch$Q119==9] <- 1 # very unfavorable hist(data.ch$favor_us) # Treatment groups data.ch$treat_sd<-ifelse(is.na(data.ch$Q3)==FALSE,0, ifelse(is.na(data.ch$Q5)==FALSE,1, ifelse(is.na(data.ch$Q7)==FALSE,2, ifelse(is.na(data.ch$Q9)==FALSE,3, ifelse(is.na(data.ch$Q119)==FALSE,4,NA))))) ### the plot below shows the distribution of favoritism toward the US among Chinese respondents library(ggplot2) library(scales) png(file="p_antiforeign_cn_control.png",width=800,height=1000, res=200,pointsize = 6,family = "SimHei", type='cairo') ggplot(data.ch[data.ch$treat_sd==0,], aes(x = favor_us)) + geom_bar(aes(y = (..count..)/sum(..count..))) + geom_text(aes(y = ((..count..)/sum(..count..)), label = scales::percent((..count..)/sum(..count..))), stat = "count", vjust = -0.25) + scale_x_continuous(breaks=1:6, labels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", "Slightly favorable","Somewhat favorable", "Very favorable"), guide = guide_axis(angle = 45))+ scale_y_continuous(labels = scales::percent, limits=c(0,0.3)) + theme_bw()+ labs(y = "Percent", x = "") dev.off() # Blog visualization 1 ##### # favor_ch <- data.us[data.us$treat_sd==0,] %>% # group_by(favor=favor_ch) %>% # mutate(favor = case_when( # favor == 1 ~factor("Very unfavorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # favor == 2 ~factor("Somewhat unfavorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # favor == 3 ~factor("Slightly unfavorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # favor == 4 ~factor("Slightly favorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # favor == 5 ~factor("Somewhat favorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # favor == 6 ~factor("Very favorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # )) %>% # summarise(count = n()) %>% # ungroup() %>% # mutate(perc = round(count/sum(count)*100,2), country = 'CH') # favor_us <- data.ch[data.ch$treat_sd==0,] %>% # group_by(favor=favor_us) %>% # mutate(favor = case_when( # favor == 1 ~factor("Very unfavorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # favor == 2 ~factor("Somewhat unfavorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # favor == 3 ~factor("Slightly unfavorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # favor == 4 ~factor("Slightly favorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # favor == 5 ~factor("Somewhat favorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # favor == 6 ~factor("Very favorable",levels=c("Very unfavorable","Somewhat unfavorable","Slightly unfavorable", # "Slightly favorable","Somewhat favorable", # "Very favorable")), # )) %>% # summarise(count = n()) %>% # ungroup() %>% # mutate(perc = round(count/sum(count)*100,2),country = 'US') # bind_rows(favor_us, favor_ch) %>% # write_csv('export.csv') ###################### covariates ######################## ############## Nationalism ############### # 1=strongly agree--> more national pride, 5=strongly disagree-->less national pride data.us$nationlism1<-6-as.numeric(data.us$Q19_1) data.us$nationlism2<-6-as.numeric(data.us$Q19_2) data.us$nationlism3<-6-as.numeric(data.us$Q19_3) ############## Military Assertivness 20_1/2/3 ########### data.us$ma3<-6-as.numeric(data.us$Q20_3) data.us$ma2<-as.numeric(data.us$Q20_2) data.us$ma1<-6-as.numeric(data.us$Q20_1) data.us$ma <- ((data.us$ma1) + (data.us$ma2) + (data.us$ma3))*(1/3) ############## Demographics ############# # gender 1=male, 0=female data.us$gender<-ifelse(data.us$gen==1,1, ifelse(data.us$gen==2,0,NA)) table(data.us$gender) # year of birth data.us$age<-as.character(data.us$age) data.us$age[which(data.us$age=="10/26/1976")]<-"1967" data.us$yrbirth<-as.numeric(data.us$age) data.us$yrbirth<-ifelse(!nchar(data.us$yrbirth)==4,NA,data.us$yrbirth) table(data.us$yrbirth) # Region # 1=Northeast,2=Midwest,3=South,4=West table(data.us$region) table(data.us.text$region) # education #original code: 1= None, 2=Some high school, 3=High school graduate, #4=Some college,5=Four year college degree,6=Graduate or professional degree. data.us$education<-data.us$edu table(data.us$education) # knowledge about politics/ir data.us$knowledge.house.speaker<-ifelse(data.us$know_1 %in% c(2),1,0) data.us$knowledge.senate<-ifelse(data.us$know_2 %in% c(2),1,0) data.us$knowledge.un<-ifelse(data.us$know_3 %in% c(4),1,0) table(data.us$knowledge.un) data.us$know<-data.us$knowledge.un+ data.us$knowledge.senate +data.us$knowledge.house.speaker mean(data.us$knowledge.un) mean(data.us$knowledge.house.speaker) mean(data.us$knowledge.senate) # income # 1= less than 1000 yuan, 9=more than 50000yuan, 10=refuse to answer data.us$income<-ifelse(data.us$Q146==11,NA,data.us$Q146) table(data.us$income) ################# covariates China ################# ############## Nationalism ############### # 1=strongly agree--> more national pride, 5=strongly disagree-->less national pride data.ch$nationlism1<-6-as.numeric(data.ch$Q19_1) data.ch$nationlism2<-6-as.numeric(data.ch$Q19_2) data.ch$nationlism3<-6-as.numeric(data.ch$Q19_3) data.ch$nationalism<-(data.ch$nationlism1+data.ch$nationlism2+data.ch$nationlism3)/3 # Military Assertivness 20_1/2/3 data.ch$ma3<-6-as.numeric(data.ch$Q20_3) data.ch$ma2<-as.numeric(data.ch$Q20_2) data.ch$ma1<-6-as.numeric(data.ch$Q20_1) data.ch$ma <- ((data.ch$ma1) + (data.ch$ma2) + (data.ch$ma3))*(1/3) ############## Demographics ############# # gender 1=male, 0=female data.ch$gender<-ifelse(data.ch$gen==1,1,0) table(data.ch$gender) # year of birth data.ch$yrbirth<-as.numeric(data.ch$age) table(data.ch$yrbirth) # live in city 1=city 2=rural # all live in cities table(data.ch$city) # education #original code: 1= elementary school, 11=senior high school, 2=junior high school #3= vocational high school, 4=Associate degree, 5= master's degree, 6=doctor degree, 8=colledge degree table(data.ch$edu) table(data.ch.text$edu) #recode them as 1= elementary school, 2=junior high school,3= vocational high school, #4=senior high school,5=Associate degree, 6=colledge degree,7= master's degree, 8=doctor degree data.ch$education<-ifelse(data.ch$edu==11,4, ifelse(data.ch$edu==4,5, ifelse(data.ch$edu==8,6, ifelse(data.ch$edu==5,7, ifelse(data.ch$edu==6,8,data.ch$edu))))) table(data.ch$education) # knowledge about politics/ir data.ch$knowledge.un<-ifelse(data.ch$know_2 %in% c(4),1,0) data.ch$knowledge.npc<-ifelse(data.ch$know_3 %in% c(1),1,0) data.ch$knowledge.cpc<-ifelse(data.ch$know_1 %in% c(2),1,0) table(data.ch$knowledge.un) data.ch$know<-data.ch$knowledge.un+ data.ch$knowledge.npc +data.ch$knowledge.cpc mean(data.ch$knowledge.cpc) mean(data.ch$knowledge.npc) mean(data.ch$knowledge.un) # income # 1= less than 1000 yuan, 9=more than 50000yuan, 10=refuse to answer data.ch$income<-ifelse(data.ch$Q213==10,NA,data.ch$Q213) table(data.ch$income) #### main results table ########################### data.ch$income<-as.numeric(as.character(data.ch$income)) data.ch$education<-as.numeric(as.character(data.ch$education)) data.us$income<-as.numeric(as.character(data.us$income)) data.us$education<-as.numeric(as.character(data.us$education)) data.ch$treat_sd_factor <- as.factor(as.numeric(data.ch$treat_sd)) data.us$treat_sd_factor <- as.factor(as.numeric(data.us$treat_sd)) summary(m.ch<-lm(favor_us~treat_sd_factor+ gender+ yrbirth+ education+ income,data=data.ch)) summary(m.us<-lm(favor_ch~treat_sd_factor+ gender+ yrbirth+ education+ income,data=data.us)) library("lmtest") library("sandwich") # m.ch<-coeftest(m.ch, vcov = vcovHC(m.ch, type = "HC0")) # m.ch # m.us<-coeftest(m.us, vcov = vcovHC(m.us, type = "HC0")) # m.us ## plot the results library(sjPlot) library(ggplot2) theme_set(theme_sjplot()) a<-plot_model(m.ch, type = "pred", terms = "treat_sd_factor", axis.title = c("Group","Favorability level"), title="") d_plot<-a[[1]] # d_plot$y<-NA # d_plot$y[1]<-mean(data.ch$favor_us[data.ch$treat_sd==0]) # d_plot$y[2]<-mean(data.ch$favor_us[data.ch$treat_sd==1]) # d_plot$y[3]<-mean(data.ch$favor_us[data.ch$treat_sd==2]) # d_plot$y[4]<-mean(data.ch$favor_us[data.ch$treat_sd==3]) # d_plot$y[5]<-mean(data.ch$favor_us[data.ch$treat_sd==4]) # d_plot$x<-as.factor(d_plot$x) # class(d_plot$x) d_plot$x <- factor(d_plot$x,levels = c(0,3,2,1,4)) png(file="p_marginal effects_cn.png",width=600,height=400, res=100,pointsize = 5,family = "Source Sans Pro", type='cairo') ggplot(d_plot, aes(x,predicted))+ geom_point()+ geom_bar(aes(x,predicted,fill=x),alpha=0.7, stat="identity")+ geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.2)+ labs(x="",y="Favorability Ratings (1-6 scale)")+ coord_flip(ylim = c(1,6))+ scale_x_discrete(labels=c("Control", "Negative Reciprocity", "Positive Reciprocity", "Humanization", "Interaction"))+ scale_fill_manual( values = #c("#264653","#2a9d8f","#e9c46a","#f4a261","#e76f51") #c("#1a535c","#4ecdc4","#f7fff7","#ff6b6b","#ffe66d") #c("#247ba0","#70c1b3","#b2dbbf","#f3ffbd","#ff1654") #c("#360568","#5b2a86","#7785ac","#9ac6c5","#a5e6ba") c("#300049","#4a0072","#00695c","#00838f","#ffd600") )+ geom_hline(yintercept=d_plot$predicted[1], linetype="dashed", color = "red")+ geom_text(aes(label=format(round(predicted,2)),nsmall=2),color='white', hjust=2, size=4)+ theme(legend.position = 'None')+ scale_y_continuous(breaks = seq(1,6)) + ggtitle('Chinese Respondents') dev.off() p_us<-plot_model(m.us, type = "pred", terms = "treat_sd_factor", axis.title = c("Group","Favorability level"), title="") d_plot<-p_us[[1]] # d_plot$y<-NA # d_plot$y[1]<-mean(data.us$favor_ch[data.us$treat_sd==0]) # d_plot$y[2]<-mean(data.us$favor_ch[data.us$treat_sd==1]) # d_plot$y[3]<-mean(data.us$favor_ch[data.us$treat_sd==2]) # d_plot$y[4]<-mean(data.us$favor_ch[data.us$treat_sd==3]) # d_plot$y[5]<-mean(data.us$favor_ch[data.us$treat_sd==4]) d_plot$x <- factor(d_plot$x,levels = c(0,3,2,1,4)) png(file="p_marginal effects_us.png",width=600,height=400, res=100,pointsize = 5,family = "Source Sans Pro", type='cairo') ggplot(d_plot, aes(x,predicted))+ geom_point()+ geom_bar(aes(fill=x),alpha=0.7, stat="identity")+ geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.2)+ labs(x="",y="Favorability Ratings (1-6 scale)")+ coord_flip(ylim = c(1,6))+ scale_x_discrete(labels=c("Control", "Negative Reciprocity", "Positive Reciprocity", "Humanization", "Interaction"))+ scale_fill_manual( values = #c("#264653","#2a9d8f","#e9c46a","#f4a261","#e76f51") #c("#1a535c","#4ecdc4","#f7fff7","#ff6b6b","#ffe66d") #c("#247ba0","#70c1b3","#b2dbbf","#f3ffbd","#ff1654") #c("#360568","#5b2a86","#7785ac","#9ac6c5","#a5e6ba") c("#300049","#4a0072","#00695c","#00838f","#ffd600") )+ geom_hline(yintercept=d_plot$predicted[1], linetype="dashed", color = "red")+ geom_text(aes(label=format(round(predicted,2)),nsmall=2),color='white', hjust=2, size=4)+ theme(legend.position = 'None')+ scale_y_continuous(breaks = seq(1,6)) + ggtitle('U.S. Respondents') dev.off()