library(foreign)
library(lme4)
library(margins)
library(dotwhisker)
library(gridExtra)
library(tibble)
library(ggplot2)

data<-read.csv("EJPRReplicationData.csv")

#mepid: MEP ID
#igid: IG ID
#MSid: ID the MS of the MEP
#NPid: National party of the MEP
#HQid: HQ country of the interest group
#meptoig: MEP follows IG on Twitter
#igtomep: IG follows MEP on Twitter
#meetings: Number of EC meetings by the IG
#passes: Number of EP passes by the IG
#type: Type of interest group
#shared_MS: HQ of iG is in the same country as the MS of MEP
#Elsystem: Electoral system in the MS of the MEP
#lrgen: Left-Right ideology of MEP
#prominence: Prominence of IG in relation to the specific MEP
#eu: IG is an EU-level organization
#nofriends_log: Number of Twitter accounts followed by MEP, logged
#followers_ig_log: Number of Twitter accounts following IG, logged
#chair: MEP is a (vice-)chair of an EP committee
#participation.sqr: Share of plenary sessions attended by MEP, squared root
#reportcounts.seq: Number of legislative and non-legislative reports by MEP, squared root
#numberterms: Number of terms MEP served in EP
#cee: MEP from a CEE member state
#SharedLanguage: Estimate of MEP speaking the language of the country of IG's HQ

####################################################
###############Descriptive analyses#################
data$mepid<-as.factor(data$mepid)
data$igid<-as.factor(data$igid)

##Median number of stakeholder friends
median(table(data$mepid[data$meptoig==1]))

#Median number of "friends" by MEP
median(exp(data$nofriends_log)[duplicated(data$mepid)=="FALSE"])

#Number of IGs followed by MEPs
table(table(data$mepid[data$meptoig==1]))

##Tie reciprocation
table(data$meptoig,data$igtomep)

##Median number of IG friends by MEPs
tmp<-table(data$mepid[data$meptoig==1])
tmp<-as.numeric(tmp)
summary(tmp)

##Distribution of MEPs following a group
tmp<-table(data$igid[data$meptoig==1])
tmp1<-as.numeric(tmp)
summary(as.numeric(tmp))
length(tmp1[tmp1==0])
length(tmp1[tmp1>28])
length(tmp1[tmp1>56])
####################################################

####################################################
############Convergent validity (Figure 1)##########
mepfollowers<-as.data.frame(table(data$igid[data$meptoig==1]))
names(mepfollowers)<-c("igid","meptoig")
mepfollowers$sqrmeptoig<-sqrt(mepfollowers$meptoig)

###Info on meetings and passes
igs<-data[duplicated(data$igid)=="FALSE",]
igs<-igs[c("igid","meetings","passes")]
igs$sqrmeetings<-sqrt(igs$meetings)
igs$sqrpasses<-sqrt(igs$passes)
mepfollowers<-merge(mepfollowers,igs,by="igid",all.x="TRUE")

####################################################
####################################################
scat1<-ggplot(mepfollowers, aes(x=sqrmeptoig, y=sqrpasses) ) +
  geom_count() + 
  labs(colour="N", x="Square root (Number of MEPs recognizing the IG)",
       y="Square root (Count of EP passes)")+
  stat_smooth(method = "lm", col = "black")+
  theme_bw()+
  theme(text = element_text(size=15))+
  geom_point(size=3)+
  ylim(0,8)
scat2<-ggplot(mepfollowers, aes(x=sqrmeptoig, y=sqrmeetings) ) +
  geom_count() + 
  labs(colour="N", x="Square root (Number of MEPs recognizing the IG)",
       y="Square root (Count of meetings with EC officials")+
  theme(text = element_text(size=10))+
  stat_smooth(method = "lm", col = "black")+
  theme_bw()+
  theme(text = element_text(size=15))+
  geom_point(size=3)+
  ylim(0,10.5)+scale_y_continuous(breaks=seq(0,10,2))
png("validity.png",res=600,height=5000,width=7500)
grid.arrange(scat1,scat2,ncol=2)
dev.off()
####################################################
####################################################

############Inferential analyses####################
data<-data[complete.cases(data$meptoig,data$shared_MS,data$Elsystem,
                          data$lrgen,data$prominence,
                          data$eu,
                          data$nofriends_log,
                          data$followers_ig_log,
                          data$chair,
                          data$participation.sqr,
                          data$reportscounts.sqr,
                          data$numberterms,
                          data$cee,
                          data$SharedLanguage),]

####################################################
####################################################
###################Models in Table 1################
mbusiness1<-glmer(meptoig~shared_MS*Elsystem+
                    lrgen*prominence+
                    eu+
                    nofriends_log+
                    followers_ig_log+
                    chair+
                    participation.sqr+
                    reportscounts.sqr+
                    numberterms+
                    cee+
                    SharedLanguage+(1|mepid)+(1|igid)+
                    (1|NPid)+(1|MSid)+(1|HQid),
                  data=data,binomial("logit"),
                  nAGQ=0,
                  control=glmerControl(optimizer = "nloptwrap"),
                  subset=type==2)

mngo1<-update(mbusiness1,.~.,subset=type==4)
mtrade1<-update(mbusiness1,.~.,subset=type==8)
mprof1<-update(mbusiness1,.~.,subset=type==5)
mcons1<-update(mbusiness1,.~.,subset=type==3)
mtt1<-update(mbusiness1,.~.,subset=type==7)
mac1<-update(mbusiness1,.~.,subset=type==1)
msub1<-update(mbusiness1,.~.,subset=type==6)
#################################################################################
#################################################################################

################################################################################
#######################Electoral system effects#################################
################################################################################
bus.el<-summary(margins(mbusiness1,at=list(Elsystem=c("Closed","Flex","Open")),type="response",variables="shared_MS"))
ngo.el<-summary(margins(mngo1,at=list(Elsystem=c("Closed","Flex","Open")),type="response",variables="shared_MS"))                
trade.el<-summary(margins(mtrade1,at=list(Elsystem=c("Closed","Flex","Open")),type="response",variables="shared_MS")) 
prof.el<-summary(margins(mprof1,at=list(Elsystem=c("Closed","Flex","Open")),type="response",variables="shared_MS")) 
cons.el<-summary(margins(mcons1,at=list(Elsystem=c("Closed","Flex","Open")),type="response",variables="shared_MS"))
tt.el<-summary(margins(mtt1,at=list(Elsystem=c("Closed","Flex","Open")),type="response",variables="shared_MS"))
ac.el<-summary(margins(mac1,at=list(Elsystem=c("Closed","Flex","Open")),type="response",variables="shared_MS"))
sub.el<-summary(margins(msub1,at=list(Elsystem=c("Closed","Flex","Open")),type="response",variables="shared_MS"))

bus.el$type<-"Business"
ngo.el$type<-"NGO"
trade.el$type<-"TradeU"
prof.el$type<-"Professional"
cons.el$type<-"Consultancy"
tt.el$type<-"ThinkTank"
ac.el$type<-"Academic"
sub.el$type<-"Subnational"
marg<-rbind(bus.el,ngo.el,trade.el,
            prof.el,cons.el,tt.el,
            ac.el,sub.el)

#################################################################################
#################################################################################
############################Preparing estimates for plotting#####################
marg$Elsystem[marg$Elsystem=="Flex"]<-"Flexible"
marg$Elsystem[marg$Elsystem=="Open"]<-"Open/STV"
marg$std.error<-marg$SE
marg$statistic<-marg$z
marg$p.value<-marg$p
marg$term<-marg$type
marg$estimate<-marg$AME
marg$lb<-marg$lower
marg$ub<-marg$upper
marg$model<-marg$Elsystem
marg$model<-as.factor(marg$model)
marg$term<-as.character(marg$term)
marg$term[marg$term=="TradeU"]<-"TradeUnion"
marg$term<-factor(marg$term,
                  levels=c("Business","NGO","TradeUnion",
                           "Professional","Consultancy",
                           "ThinkTank","Academic",
                           "Subnational"))
marg<-marg[order(marg$term,marg$model),]

marg<-marg[c("term","estimate","std.error",
             "statistic","p.value",
             "model")]
marg1<-marg[marg$term%in%c("Business","NGO","TradeUnion","Professional"),]
marg2<-marg[marg$term%in%c("Consultancy","ThinkTank","Academic","Subnational"),]
marg2$term<-as.factor(as.character(marg2$term))
marg2$term<-factor(marg2$term,
                   levels=c("Consultancy","ThinkTank","Academic","Subnational"))
marg2$model<-as.factor(as.character(marg2$model))
marg2$model<-factor(marg2$model,
                    levels=c("Closed","Flexible","Open/STV"))
marg1<-as_tibble(marg1)
marg2<-as_tibble(marg2)
#################################################################################

########################################Plot#####################################
plot1<-small_multiple(marg1) +
  theme_bw() + 
  scale_y_continuous(limits=c(0,0.05),oob = scales::squish)+
  ylab("Marginal Effect") +
  xlab("Electoral system")+
  theme(text = element_text(size=10))+
  theme(legend.position = "none")+
  geom_point(size=4)+
  geom_hline(yintercept = 0, color = "grey60", linetype = 2)
plot2<-small_multiple(marg2) +
  theme_bw() + 
  scale_y_continuous(limits=c(0,0.05),oob = scales::squish)+
  theme(text = element_text(size=10))+
  theme(legend.position = "none")+
  geom_point(size=4)+
  xlab("Electoral system")+
  geom_hline(yintercept = 0, color = "grey60", linetype = 2)
png("SHAREDMS.EFFECT.png",
    height=5000,width=7000,res=600)
grid.arrange(plot1,plot2,ncol=2)
dev.off()
#################################################################################
#################################################################################
#################################################################################


#####################################################################
######################################################################
###################Effect of LR positions (Figure 3)##################
bus.lrgen<-summary(margins(mbusiness1,at=list(qqrmaxcentr=seq(0,2,0.2)),type="response",variables="lrgen"))
ngo.lrgen<-summary(margins(mngo1,at=list(qqrmaxcentr=seq(0,2.2,0.2)),type="response",variables="lrgen"))
trade.lrgen<-summary(margins(mtrade1,at=list(qqrmaxcentr=seq(0,1.8,0.2)),type="response",variables="lrgen"))
prof.lrgen<-summary(margins(mprof1,at=list(qqrmaxcentr=seq(0,1.6,0.2)),type="response",variables="lrgen"))
cons.lrgen<-summary(margins(mcons1,at=list(qqrmaxcentr=seq(0,1.6,0.2)),type="response",variables="lrgen"))
tt.lrgen<-summary(margins(mtt1,at=list(qqrmaxcentr=seq(0,2.0,0.2)),type="response",variables="lrgen"))
ac.lrgen<-summary(margins(mac1,at=list(qqrmaxcentr=seq(0,1.6,0.2)),type="response",variables="lrgen"))
sub.lrgen<-summary(margins(msub1,at=list(qqrmaxcentr=seq(0,1.6,0.2)),type="response",variables="lrgen"))
bus.lrgen$type<-"Business"
ngo.lrgen$type<-"NGO"
trade.lrgen$type<-"TradeU"
prof.lrgen$type<-"Professional"
cons.lrgen$type<-"Consultancy"
tt.lrgen$type<-"ThinkTank"
ac.lrgen$type<-"Academic"
sub.lrgen$type<-"Subnational"
margeffects.lrgen<-rbind(bus.lrgen,ngo.lrgen,trade.lrgen,
                         prof.lrgen,cons.lrgen,tt.lrgen,
                         ac.lrgen,sub.lrgen)

#################################################################################
#################################################################################
############################Preparing estimates for plotting#####################
marg<-margeffects.lrgen
marg$std.error<-marg$SE
marg$statistic<-marg$z
marg$p.value<-marg$p
marg$term<-marg$type
marg$estimate<-marg$AME
marg$lb<-marg$lower
marg$ub<-marg$upper
marg$model<-marg$qqrmaxcentr
marg$model<-as.factor(marg$model)
marg$term<-as.character(marg$term)
marg$term[marg$term=="TradeU"]<-"TradeUnion"
marg$term<-factor(marg$term,
                  levels=c("Business","NGO","TradeUnion",
                           "Professional","Consultancy",
                           "ThinkTank","Academic",
                           "Subnational"))
marg<-marg[order(marg$term,marg$model),]

marg<-marg[c("term","estimate","std.error",
             "statistic","p.value",
             "model")]

marg1<-marg[marg$term%in%c("Business","NGO","TradeUnion","Professional"),]
marg2<-marg[marg$term%in%c("Consultancy","ThinkTank","Academic","Subnational"),]
marg2$term<-as.factor(as.character(marg2$term))
marg2$term<-factor(marg2$term,
                   levels=c("Consultancy","ThinkTank","Academic","Subnational"))
marg2$model<-as.factor(as.character(marg2$model))
marg1<-as_tibble(marg1)
marg2<-as_tibble(marg2)
#################################################################################

########################################Plot#####################################
plot1<-small_multiple(marg1) +
  theme_bw() + 
  scale_y_continuous(limits=c(-0.006,0.002),oob = scales::squish)+
  ylab("Marginal Effect") +
  xlab("Interest group prominence")+
  theme(text = element_text(size=10))+
  theme(legend.position = "none")+
  geom_point(size=2)+
  geom_hline(yintercept = 0, color = "grey60", linetype = 2)
plot2<-small_multiple(marg2) +
  theme_bw() + 
  scale_y_continuous(limits=c(-0.006,0.002),oob = scales::squish)+
  theme(text = element_text(size=10))+
  theme(legend.position = "none")+
  geom_point(size=2)+
  xlab("Interest group prominence")+
  geom_hline(yintercept = 0, color = "grey60", linetype = 2)
png("LR.EFFECT.png",
    height=5000,width=7500,res=600)
grid.arrange(plot1,plot2,ncol=2)
dev.off()
####################################################################
#################################################################################
#################################################################################
#################################################################################



#################################################################################
###################Effect of prominence (Figure 4)###############################
bus.qqrmaxcentr<-summary(margins(mbusiness1,at=list(lrgen=seq(0,10,1)),type="response",variables="qqrmaxcentr"))
ngo.qqrmaxcentr<-summary(margins(mngo1,at=list(lrgen=seq(0,10,1)),type="response",variables="qqrmaxcentr"))
trade.qqrmaxcentr<-summary(margins(mtrade1,at=list(lrgen=seq(0,10,1)),type="response",variables="qqrmaxcentr"))
prof.qqrmaxcentr<-summary(margins(mprof1,at=list(lrgen=seq(0,10,1)),type="response",variables="qqrmaxcentr"))
cons.qqrmaxcentr<-summary(margins(mcons1,at=list(lrgen=seq(0,10,1)),type="response",variables="qqrmaxcentr"))
tt.qqrmaxcentr<-summary(margins(mtt1,at=list(lrgen=seq(0,10,1)),type="response",variables="qqrmaxcentr"))
ac.qqrmaxcentr<-summary(margins(mac1,at=list(lrgen=seq(0,10,1)),type="response",variables="qqrmaxcentr"))
sub.qqrmaxcentr<-summary(margins(msub1,at=list(lrgen=seq(0,10,1)),type="response",variables="qqrmaxcentr"))
bus.qqrmaxcentr$type<-"Business"
ngo.qqrmaxcentr$type<-"NGO"
trade.qqrmaxcentr$type<-"TradeU"
prof.qqrmaxcentr$type<-"Professional"
cons.qqrmaxcentr$type<-"Consultancy"
tt.qqrmaxcentr$type<-"ThinkTank"
ac.qqrmaxcentr$type<-"Academic"
sub.qqrmaxcentr$type<-"Subnational"
margeffects.qqrmaxcentr<-rbind(bus.qqrmaxcentr,ngo.qqrmaxcentr,trade.qqrmaxcentr,
                               prof.qqrmaxcentr,cons.qqrmaxcentr,tt.qqrmaxcentr,
                               ac.qqrmaxcentr,sub.qqrmaxcentr)


#################################################################################
#################################################################################
############################Preparing estimates for plotting#################################
marg<-margeffects.qqrmaxcentr
marg$std.error<-marg$SE
marg$statistic<-marg$z
marg$p.value<-marg$p
marg$term<-marg$type
marg$estimate<-marg$AME
marg$lb<-marg$lower
marg$ub<-marg$upper
marg$model<-marg$lrgen
marg$model<-as.factor(marg$model)
marg$term<-as.character(marg$term)
marg$term[marg$term=="TradeU"]<-"TradeUnion"
marg$term<-factor(marg$term,
                  levels=c("Business","NGO","TradeUnion",
                           "Professional","Consultancy",
                           "ThinkTank","Academic",
                           "Subnational"))
marg<-marg[order(marg$term,marg$model),]

marg<-marg[c("term","estimate","std.error",
             "statistic","p.value",
             "model")]

marg1<-marg[marg$term%in%c("Business","NGO","TradeUnion","Professional"),]
marg2<-marg[marg$term%in%c("Consultancy","ThinkTank","Academic","Subnational"),]
marg2$term<-as.factor(as.character(marg2$term))
marg2$term<-factor(marg2$term,
                   levels=c("Consultancy","ThinkTank","Academic","Subnational"))
marg2$model<-as.factor(as.character(marg2$model))
marg2$model<-factor(marg2$model,
                    levels=c("0","1","2",
                             "3","4",
                             "5","6",
                             "7","8","9","10"))
marg1<-as_tibble(marg1)
marg2<-as_tibble(marg2)

#####################################Plot#################################
plot1<-small_multiple(marg1) +
  theme_bw() + 
  scale_y_continuous(limits=c(0,0.02),oob = scales::squish)+
  ylab("Marginal Effect") +
  xlab("Left-Right Positions of MEP")+
  theme(text = element_text(size=10))+
  theme(legend.position = "none")+
  geom_point(size=2)+
  geom_hline(yintercept = 0, color = "grey60", linetype = 2)
plot2<-small_multiple(marg2) +
  theme_bw() + 
  scale_y_continuous(limits=c(0,0.02),oob = scales::squish)+
  theme(text = element_text(size=10))+
  theme(legend.position = "none")+
  geom_point(size=2)+
  xlab("Left-Right Positions of MEP")+
  geom_hline(yintercept = 0, color = "grey60", linetype = 2)
png("PROMINENCE.EFFECT.png",
    height=5000,width=7500,res=600)
grid.arrange(plot1,plot2,ncol=2)
dev.off()
######################################################################
######################################################################
