#####################################################################
# Replication file for:                                            #
# Lehrer et al. "Rallying Around the Leader in Times of Crises:    #
# "The Opposing Effects of Perceived Threat and Anxiety"           #
#                                                                  #
# Published in the European Journal of Political Research          #
####################################################################

########################################################################
# Last edited:                                                         #
# 26 July 2024 by Roni Lehrer                                          #
#                                                                      #
# Information on package versions can be found at the end of this file #
########################################################################

## clear working memory
rm(list=ls())

## load packages
library(data.table)
library(stargazer)
library(sandwich)
library(weights)
library(car)
library(lme4)
library(ggplot2)
library(diagis)

#load function weighted.fe.robust
source("functions/weighted_fe_robust.R")

#load data
#to prepare data from original survey files, see "prepare_data.R".
load("prepared_data.Rdata")


############################################
#Summary stats Merkel Approval (Table SI1) #
############################################

#keep weeks in which Merkel approval item was asked
sub <- prepared.data[Woche%in%c(1,2,3,4,5,6,7,8,11,14,16)]

#keep observation only if no missings in key variables
sub <- sub[!is.na(merkel1)&!is.na(threat)&!is.na(fear)&!is.na(wi_rake)]
sub <- sub[!is.na(last7days.cases.bund)&!is.na(hh.incom.prev.month.pc)&!is.na(meets.borders)]

#keep observation only if given respondent provides data at two points in time
sub <- sub[,count:=.N, id_g]
sub <- sub[count>1]

sub <- sub[,list(merkel1, threat, fear,  last7days.cases.bund, hh.incom.prev.month.pc, meets.borders,
                  id_g, Woche, wi_rake)]

summary.data <- sub[,list(merkel1, threat, fear,  last7days.cases.bund, hh.incom.prev.month.pc, meets.borders)]
stargazer(summary.data,
          omit.summary.stat=c("p25", "p75"),
          covariate.labels = c("Merkel Approval",
                               "Perceived Threat",
                               "Anxiety",
                               "COVID-19 Incidence",
                               "HH Income Previous Month",
                               "Policy Congruence: Border Closures"),
          header=FALSE,
          title="Summary Statistics: German Panel Data",
          label="tab:summary_merkel"
)

###########
# Table 1 #
###########
merkel1 <- weighted.fe.robust(merkel1~threat, data=sub, group="id_g", time="Woche", weights="wi_rake")
merkel2 <- weighted.fe.robust(merkel1~
                              threat
                              +last7days.cases.bund
                              +hh.incom.prev.month.pc
                              +meets.borders,
                              data=sub, group="id_g", time="Woche", weights="wi_rake")
merkel3 <- weighted.fe.robust(merkel1~fear, data=sub, group="id_g", time="Woche", weights="wi_rake")
merkel4 <- weighted.fe.robust(merkel1~
                              +fear
                              +last7days.cases.bund
                              +hh.incom.prev.month.pc
                              +meets.borders,
                              data=sub, group="id_g", time="Woche", weights="wi_rake")
merkel5 <- weighted.fe.robust(merkel1~threat+fear, data=sub, group="id_g", time="Woche", weights="wi_rake")
merkel6 <- weighted.fe.robust(merkel1~
                               threat
                              +fear
                              +last7days.cases.bund
                              +hh.incom.prev.month.pc
                              +meets.borders,
                              data=sub, group="id_g", time="Woche", weights="wi_rake")
stargazer(merkel1,
          merkel2,
          merkel3,
          merkel4,
          merkel5,
          merkel6,
          no.space=TRUE,
          omit.stat=c("f", "rsq", "adj.rsq", "ser"),
          header=FALSE,
          title="The Effect of Perceived Threat and Anxiety on Merkel Approval",
          order=c(1, 5, 2, 3, 4),
          covariate.labels = c("Perceived Threat",
                               "Anxiety",
                               "COVID-19 Incidence",
                               "HH Income Previous Month",
                               "Policy Congruence: Border Closures"),
          add.lines = list(c("Individual Fixed effects", "Yes", "Yes", "Yes", "Yes","Yes", "Yes"),
                           c("Number of Respondents",
                             length(unique(sub$id_g)), length(unique(sub$id_g)),
                             length(unique(sub$id_g)), length(unique(sub$id_g)),
                             length(unique(sub$id_g)), length(unique(sub$id_g)))),
          se= list(sqrt(diag(merkel1$clustered.vcov)),
                   sqrt(diag(merkel2$clustered.vcov)),
                   sqrt(diag(merkel3$clustered.vcov)),
                   sqrt(diag(merkel4$clustered.vcov)),
                   sqrt(diag(merkel5$clustered.vcov)),
                   sqrt(diag(merkel6$clustered.vcov))
                   ),
         dep.var.caption  = "",
         dep.var.labels.include = FALSE,
         label="tab:merkel", digits=2
          )


##############
# Table SI 4 #
##############
sub$threat <- scale(sub$threat, center=TRUE, scale=TRUE)
sub$fear <- scale(sub$fear, center=TRUE, scale=TRUE)
sub$last7days.cases.bund <- scale(sub$last7days.cases.bund, center=TRUE, scale=TRUE)
sub$hh.incom.prev.month.pc <- scale(sub$hh.incom.prev.month.pc, center=TRUE, scale=TRUE)
sub$meets.borders <- scale(sub$meets.borders, center=TRUE, scale=TRUE)
sub$merkel1 <- scale(sub$merkel1, center=TRUE, scale=TRUE)
merkel7 <- weighted.fe.robust(merkel1~
                               threat
                              +fear,
                              data=sub, group="id_g", time="Woche", weights="wi_rake")
merkel8 <- weighted.fe.robust(merkel1~
                               threat
                              +fear
                              +last7days.cases.bund
                              +hh.incom.prev.month.pc
                              +meets.borders,
                              data=sub, group="id_g", time="Woche", weights="wi_rake")
stargazer(merkel7,
          merkel8,
          no.space=TRUE,
          omit.stat=c("f", "rsq", "adj.rsq", "ser"),
          header=FALSE,
          title="The Effect of Perceived Threat and Anxiety on Merkel Approval (Standardized Variables)",
          covariate.labels = c("Perceived Threat",
                               "Anxiety",
                               "COVID-19 Incidence",
                               "HH Income Previous Month",
                               "Policy Congruence: Border Closures"),
          add.lines = list(c("Individual Fixed effects", "Yes", "Yes"),
                           c("Number of Respondents",
                             length(unique(sub$id_g)), length(unique(sub$id_g)))),
          se= list(sqrt(diag(merkel7$clustered.vcov)),
                   sqrt(diag(merkel8$clustered.vcov))
                   ),
         dep.var.caption  = "",
         dep.var.labels.include = FALSE,
         label="tab:standardized", digits=2
          )


########################################################################
# Weighted Correltion of Anxiety and Perceived Threat in Demeaned Data #
########################################################################
#Demean the data
demeaned <- sub[,list(id_g, Woche, wi_rake, threat, fear)]
    for (i in names(demeaned)) {
        if (i%in%c("id_g", "Woche", "wi_rake")==FALSE) {
            if (class(demeaned[, get(i)])!="character") {
                demeaned[, (i):=get(i)-mean(get(i)), id_g]
                }
        }
    }
#compute weighted correlation
demeaned <- na.omit(demeaned[,list(threat, fear, wi_rake)])
wtd.cor(na.omit(cbind(demeaned$fear, demeaned$threat)), weight=demeaned$wi_rake)

#############################
# VIF of Model 6 in Table 1 #
#############################
vif(merkel6)

###########
# Table 2 #
###########
dat <- prepared.data
threat1 <- lmer(threat~precondition+age.cat+male+merkel0+(1|id_g)+(1|Woche),
                data=dat, weights=wi_rake)
threat2 <- lmer(threat~precondition+age.cat+male+sat.breg.nov19+(1|id_g)+(1|Woche),
                data=dat, weights=wi_rake)
fear1 <- lmer(fear~precondition+age.cat+male+merkel0+(1|id_g)+(1|Woche),
                data=dat, weights=wi_rake)
fear2 <- lmer(fear~precondition+age.cat+male+sat.breg.nov19+(1|id_g)+(1|Woche),
                data=dat, weights=wi_rake)
stargazer(threat1,
          threat2,
          fear1,
          fear2,
          no.space=TRUE,
          header=FALSE,
          omit.stat=c("aic", "bic"),
          title="Does Merkel cause Threat or Anxiety?",
          covariate.labels = c("Merkel Approval (July 2019)",
                               "Government Approval (November 2020)",
                               "Medical Precondition",
                               "60+ Years",
                               "Male",
                               "Intercept"),
          add.lines = list(c("Respondent Random Effect: Std. Dev.",
                           round(attributes(VarCorr(threat1)$id_g)$stddev,3),
                           round(attributes(VarCorr(threat2)$id_g)$stddev,3),
                           round(attributes(VarCorr(fear1)$id_g)$stddev,3),
                           round(attributes(VarCorr(fear2)$id_g)$stddev,3)),
                           c("MCS Week Random Effect: Std. Dev.",
                           round(attributes(VarCorr(threat1)$Woche)$stddev,3),
                           round(attributes(VarCorr(threat2)$Woche)$stddev,3),
                           round(attributes(VarCorr(fear1)$Woche)$stddev,3),
                           round(attributes(VarCorr(fear2)$Woche)$stddev,3)),
                           c("Number of Respondents",
                             nrow(ranef(threat1)$id_g), nrow(ranef(threat2)$id_g),
                             nrow(ranef(fear1)$id_g), nrow(ranef(fear2)$id_g))),
          column.labels=c("Perceived Threat", "Anxiety"),
          column.separate = c(2, 2),
          dep.var.labels.include = FALSE,
          dep.var.caption  = "",
         order=c(4, 5, 1, 2, 3, 6),
         label="tab:reverse2"
          )



############################
# Figure SI1 and  Figure 2 #
############################
#get demeaned values, weighted and their SEs
dat <- prepared.data
pthreat <- na.omit(dat[,list(Woche, threat, wi_rake, id_g)])
pthreat <- pthreat[,mean:=mean(threat, na.rm=TRUE), id_g]
pthreat <- pthreat[,wmean:=weighted.mean(threat, wi_rake), Woche]
pthreat <- pthreat[,dwmean:=weighted.mean(threat-mean, wi_rake), Woche]
pthreat <- pthreat[,wse:=weighted_se(threat, wi_rake), Woche]
pthreat <- pthreat[,dwse:=weighted_se(threat-mean, wi_rake), Woche]
pthreat <- unique(pthreat[,list(Woche, wmean, wse, dwmean, dwse)])
pthreat$var <- "threat"
pthreat

dat <- prepared.data
pfear <- na.omit(dat[,list(Woche, fear, wi_rake, id_g)])
pfear <- pfear[,mean:=mean(fear, na.rm=TRUE), id_g]
pfear <- pfear[,wmean:=weighted.mean(fear, wi_rake), Woche]
pfear <- pfear[,dwmean:=weighted.mean(fear-mean, wi_rake), Woche]
pfear <- pfear[,wse:=weighted_se(fear, wi_rake), Woche]
pfear <- pfear[,dwse:=weighted_se(fear-mean, wi_rake), Woche]
pfear <- unique(pfear[,list(Woche, wmean, wse, dwmean, dwse)])
pfear$var <- "fear"
pfear

dat <- prepared.data
pmerkel1 <- na.omit(dat[,list(Woche, merkel1, wi_rake, id_g)])
pmerkel1 <- pmerkel1[,mean:=mean(merkel1, na.rm=TRUE), id_g]
pmerkel1 <- pmerkel1[,wmean:=weighted.mean(merkel1, wi_rake)-6, Woche]
pmerkel1 <- pmerkel1[,dwmean:=weighted.mean(merkel1-mean, wi_rake), Woche]
pmerkel1 <- pmerkel1[,wse:=weighted_se(merkel1, wi_rake), Woche]
pmerkel1 <- pmerkel1[,dwse:=weighted_se(merkel1-mean, wi_rake), Woche]
pmerkel1 <- unique(pmerkel1[,list(Woche, wmean, wse, dwmean, dwse)])
pmerkel1$var <- "merkel"
pmerkel1

pdata <- rbindlist(list(pthreat, pfear, pmerkel1))

#Figure SI1
p <- ggplot(data=pdata, aes(x=Woche, group=var))
p <- p + geom_ribbon(aes(ymin=(wmean-(1.96*wse)),
                         ymax=(wmean+(1.96*wse)),
                         xmin=Woche, xmax=Woche), fill="grey75", alpha=.7)
p  <-  p + geom_line(aes(y=wmean, linetype=var))
p <- p + scale_y_continuous("Weighted Mean Anxiety\nWeighted Mean Perceived Threat",
                            sec.axis = sec_axis(~ . +6, name = "Weighted Mean Merkel Approval"))
p <- p + theme_bw()
p <- p + scale_x_continuous("MCS Week", breaks=1:16)
p <- p + theme(legend.position="none",
               text=element_text(size=35))
p <- p + annotate("text", x=16.75, y=pdata[Woche==16&var=="threat", wmean]+.03,
                  label="Perceived\nThreat", size=10)
p <- p + annotate("text", x=16.8, y=pdata[Woche==16&var=="fear", wmean]-.01,
                  label="Anxiety", size=10)
p <- p + annotate("text", x=16.8, y=pdata[Woche==16&var=="merkel", wmean]-.03,
                  label="Merkel\nApproval", size=10)
p

#Figure 2
pdata$dwmean <- ifelse(pdata$var=="merkel", pdata$dwmean/2, pdata$dwmean) #adjust scales for Merkel approval
pdata$dwse <- ifelse(pdata$var=="merkel", pdata$dwse/2, pdata$dwse) #adjust scales for Merkel approval


p <- ggplot(data=pdata, aes(x=Woche, group=var))
p <- p + geom_ribbon(aes(ymin=(dwmean-(1.96*dwse)),
                         ymax=(dwmean+(1.96*dwse)),
                         xmin=Woche, xmax=Woche), fill="grey75", alpha=.7)
p  <-  p + geom_line(aes(y=dwmean, linetype=var))
p <- p + scale_y_continuous("Weighted Mean Deviation from Respondent-Means\nin Perceived Threat and Anxiety",
  sec.axis = sec_axis(~ . * 2, name = "Weighted Mean Deviation from Respondent-Means\nin Merkel Approval"))
p <- p + theme_bw()
p <- p + scale_x_continuous("MCS Week", breaks=1:16)
p <- p + theme(legend.position="none",
               text=element_text(size=34))
p <- p + annotate("text", x=16.5, y=pdata[Woche==16&var=="threat", dwmean]-.05,
                  label="Perceived\nThreat", size=9)
p <- p + geom_curve(aes(x = 15.6,
                          y=-.1,
                          xend = 15.3,
                        yend = -.06),
                    curvature = -0.15,
                      arrow = arrow(length = unit(0.5, "cm")))
p <- p + annotate("text", x=16.7, y=pdata[Woche==16&var=="fear", dwmean]+.05,
                  label="Anxiety", size=9)
p <- p + geom_curve(aes(x = 16,
                          y=.008,
                          xend = 15.3,
                        yend = -.03),
                    curvature = .3,
                      arrow = arrow(length = unit(0.5, "cm")))
p <- p + annotate("text", x=16.7, y=pdata[Woche==16&var=="merkel", dwmean],
                  label="Merkel\nApproval", size=9)
p


###########
# Table 3 #
###########
#restrict dataset to observations as in Table 1 but with Altaier as DV
dat <- prepared.data
sub <- dat[!is.na(altmaier)&!is.na(threat)&!is.na(fear)&!is.na(wi_rake)]
sub <- sub[!is.na(last7days.cases.bund)&!is.na(hh.incom.prev.month.pc)&!is.na(meets.borders)]
sub <- sub[,count:=.N, id_g]
sub <- sub[count>1]
sub <- sub[,list(altmaier, threat, fear,  last7days.cases.bund, hh.incom.prev.month.pc, meets.borders,
                  id_g, Woche, wi_rake)]
n.altmaier <- length(unique(sub$id_g))
altmaier <- weighted.fe.robust(altmaier~
                              threat
                              +fear
                              +last7days.cases.bund
                              +hh.incom.prev.month.pc
                              +meets.borders,
                              data=sub, group="id_g", time="Woche", weights="wi_rake")

#restrict dataset to observations as in Table 1 but with Spahn as DV
dat <- prepared.data
sub <- dat[!is.na(spahn)&!is.na(threat)&!is.na(fear)&!is.na(wi_rake)]
sub <- sub[!is.na(last7days.cases.bund)&!is.na(hh.incom.prev.month.pc)&!is.na(meets.borders)]
sub <- sub[,count:=.N, id_g]
sub <- sub[count>1]
sub <- sub[,list(spahn, threat, fear,  last7days.cases.bund, hh.incom.prev.month.pc, meets.borders,
                  id_g, Woche, wi_rake)]
n.spahn <- length(unique(sub$id_g))
spahn <- weighted.fe.robust(spahn~
                              threat
                              +fear
                              +last7days.cases.bund
                              +hh.incom.prev.month.pc
                              +meets.borders,
                              data=sub, group="id_g", time="Woche", weights="wi_rake")
stargazer(altmaier,
          spahn,
          no.space=TRUE,
          omit.stat=c("f", "rsq", "adj.rsq", "ser"),
          header=FALSE,
          title="The Effect of Perceived Threat and Anxiety on Minister Approval",
          covariate.labels = c("Perceived Threat",
                               "Anxiety",
                               "COVID-19 Incidence",
                               "HH Income Previous Month",
                               "Policy Congruence: Border Closures"),
          add.lines = list(c("Individual Fixed effects", "Yes", "Yes"),
                           c("Number of Respondents",
                             n.altmaier, n.spahn)),
          se= list(sqrt(diag(altmaier$clustered.vcov)),
                   sqrt(diag(spahn$clustered.vcov))
                   ),
         dep.var.caption  = "",
         dep.var.labels = c("Economics Affairs", "Health"),
         label="tab:spahnaltmaier"
          )


###########################
# Table SI2 and Table SI3 #
###########################
#Compute changes since first participation for anxiety,  threat and controls
dat <- prepared.data
setkey(dat, id_g, Woche)
dat <- dat[,min.week:=min(Woche), id_g]
dat <- dat[,min.week.fear:=fear[Woche==min.week], id_g]
dat <- dat[,min.week.threat:=threat[Woche==min.week], id_g]
dat <- dat[,min.week.merkel:=merkel1[Woche==min.week], id_g]
dat <- dat[,min.week.last7days.cases.bund:=last7days.cases.bund[Woche==min.week], id_g]
dat <- dat[,min.week.hh.incom.prev.month.pc:=hh.incom.prev.month.pc[Woche==min.week], id_g]
dat <- dat[,min.week.meets.borders:=meets.borders[Woche==min.week], id_g]

dat <- dat[,lag1.Woche:=shift(Woche, n=1, type="lag"), id_g]
dat <- dat[,lag2.Woche:=shift(Woche, n=2, type="lag"), id_g]
dat <- dat[,c("diff1", "diff2"):=list(Woche-lag1.Woche, Woche-lag2.Woche), id_g]
 dat <- dat[,lag1.threat:=shift(threat, n=1, type="lag"), id_g]
 dat$lag1.threat <- ifelse(dat$diff1>1, NA, dat$lag1.threat)
 dat <- dat[,lag2.threat:=shift(threat, n=2, type="lag"), id_g]
 dat$lag2.threat <- ifelse(dat$diff2>2, NA, dat$lag2.threat)
dat <- dat[,lag1.fear:=shift(fear, n=1, type="lag"), id_g]
dat$lag1.fear <- ifelse(dat$diff1>1, NA, dat$lag1.fear)
dat <- dat[,lag2.fear:=shift(fear, n=2, type="lag"), id_g]
dat$lag2.fear <- ifelse(dat$diff2>2, NA, dat$lag2.fear)
 dat <- dat[,lag1.merkel1:=shift(merkel1, n=1, type="lag"), id_g]
 dat$lag1.merkel1 <- ifelse(dat$diff1>1, NA, dat$lag1.merkel1)
 dat <- dat[,lag2.merkel1:=shift(merkel1, n=2, type="lag"), id_g]
 dat$lag2.merkel1 <- ifelse(dat$diff2>2, NA, dat$lag2.merkel1)

#additional controls
dat <- dat[,lag1.last7days.cases.bund:=shift(last7days.cases.bund, n=1, type="lag"), id_g]
dat$lag1.last7days.cases.bund <- ifelse(dat$diff1>1, NA, dat$lag1.last7days.cases.bund)
dat <- dat[,lag2.last7days.cases.bund:=shift(last7days.cases.bund, n=2, type="lag"), id_g]
dat$lag2.last7days.cases.bund <- ifelse(dat$diff2>2, NA, dat$lag2.last7days.cases.bund)
dat$D2.last7days.cases.bund <- dat$lag1.last7days.cases.bund-dat$lag2.last7days.cases.bund
 dat <- dat[,lag1.hh.incom.prev.month.pc:=shift(hh.incom.prev.month.pc, n=1, type="lag"), id_g]
 dat$lag1.hh.incom.prev.month.pc <- ifelse(dat$diff1>1, NA, dat$lag1.hh.incom.prev.month.pc)
 dat <- dat[,lag2.hh.incom.prev.month.pc:=shift(hh.incom.prev.month.pc, n=2, type="lag"), id_g]
 dat$lag2.hh.incom.prev.month.pc <- ifelse(dat$diff2>2, NA, dat$lag2.hh.incom.prev.month.pc)
 dat$D2.hh.incom.prev.month.pc <- dat$lag1.hh.incom.prev.month.pc-dat$lag2.hh.incom.prev.month.pc
dat <- dat[,lag1.meets.borders:=shift(meets.borders, n=1, type="lag"), id_g]
dat$lag1.meets.borders <- ifelse(dat$diff1>1, NA, dat$lag1.meets.borders)
dat <- dat[,lag2.meets.borders:=shift(meets.borders, n=2, type="lag"), id_g]
dat$lag2.meets.borders <- ifelse(dat$diff2>2, NA, dat$lag2.meets.borders)
dat$D2.meets.borders <- dat$lag1.meets.borders-dat$lag2.meets.borders

dat$D1.merkel1 <- dat$merkel1-dat$min.week.merkel
dat$D2.merkel1 <- dat$lag1.merkel1-dat$min.week.merkel

dat$D1.threat <- dat$threat-dat$min.week.threat
dat$D2.threat <- dat$lag1.threat-dat$min.week.threat

dat$D1.fear <- dat$fear-dat$min.week.fear
dat$D2.fear <- dat$lag1.fear-dat$min.week.fear

dat$D1.last7days.cases.bund <- dat$last7days.cases.bund-dat$min.week.last7days.cases.bund
dat$D2.last7days.cases.bund <- dat$lag1.last7days.cases.bund-dat$min.week.last7days.cases.bund

dat$D1.last7days.cases.bund <- dat$last7days.cases.bund-dat$min.week.last7days.cases.bund
dat$D2.last7days.cases.bund <- dat$lag1.last7days.cases.bund-dat$min.week.last7days.cases.bund

dat$D1.meets.borders <- dat$meets.borders-dat$min.week.meets.borders
dat$D2.meets.borders <- dat$lag1.meets.borders-dat$min.week.meets.borders

#threat as DV
sub <- dat[Woche!=min.week]
sub <- sub[!is.na(D2.last7days.cases.bund)&!is.na(D2.hh.incom.prev.month.pc)&!is.na(D2.meets.borders)]
sub <- sub[!is.na(D2.merkel1)&!is.na(wi_rake)&!is.na(D1.threat)]
FD.threat <- lm(D1.threat~D2.merkel1-1, dat=sub, weights=wi_rake)
FD.threat.vcov <- vcovCL(FD.threat, cluster = sub$id_g)
FD.threat.n <- length(unique(sub$id_g))

FD.threat2 <- lm(D1.threat~D2.merkel1+D2.last7days.cases.bund+D2.hh.incom.prev.month.pc+D2.meets.borders-1,
                 dat=sub, weights=wi_rake)
FD.threat2.vcov <- vcovCL(FD.threat2, cluster = sub$id_g)

#fear as DV
sub <- dat[Woche!=min.week]
sub <- sub[!is.na(D2.last7days.cases.bund)&!is.na(D2.hh.incom.prev.month.pc)&!is.na(D2.meets.borders)]
sub <- sub[!is.na(D2.merkel1)&!is.na(wi_rake)&!is.na(D1.fear)]
FD.fear <- lm(D1.fear~D2.merkel1-1, dat=sub, weights=wi_rake)
FD.fear.vcov <- vcovCL(FD.fear, cluster = sub$id_g)
FD.fear.n <- length(unique(sub$id_g))

FD.fear2 <- lm(D1.fear~D2.merkel1+D2.last7days.cases.bund+D2.hh.incom.prev.month.pc+D2.meets.borders-1,
                 dat=sub, weights=wi_rake)
FD.fear2.vcov <- vcovCL(FD.fear2, cluster = sub$id_g)

#Table SI2
stargazer(FD.threat,
          FD.threat2,
          FD.fear,
          FD.fear2,
          no.space=TRUE,
          omit.stat=c("f", "rsq", "adj.rsq", "ser"),
          header=FALSE,
          title="First Difference Models",
          ## float.env="sidewaystable",
          covariate.labels = c("Change Merkel Support (Lagged)",
                               "Change COVID-19 Incidence (Lagged)",
                               "Change HH Income Previous Month (Lagged)",
                               "Change Policy Congruence: Border Closures (Lagged)"
                               ),
          add.lines = list(c("Number of Respondents",
                             FD.threat.n, FD.threat.n, FD.fear.n, FD.fear.n)),
          se= list(sqrt(diag(FD.threat.vcov, )),
                   sqrt(diag(FD.threat2.vcov,)),
                   sqrt(diag(FD.fear.vcov,   )),
                   sqrt(diag(FD.fear2.vcov,  ))),
          dep.var.labels  = c("Peceived threat",
                              "Anxiety"),
          dep.var.caption  = "",
          label="tab:FD",
          digits=2
          )


#merkel1 as DV
sub <- dat[Woche!=min.week]
sub <- sub[!is.na(D2.last7days.cases.bund)&!is.na(D2.hh.incom.prev.month.pc)&!is.na(D2.meets.borders)]
sub <- sub[!is.na(D1.merkel1)&!is.na(wi_rake)&!is.na(D2.threat)&!is.na(D2.fear)]
FD.merkel <- lm(D1.merkel1~D2.fear+D2.threat-1,
                 dat=sub, weights=wi_rake)
FD.merkel.vcov <- vcovCL(FD.merkel, cluster = sub$id_g)
FD.merkel.n <- length(unique(sub$id_g))

FD.merkel2 <- lm(D1.merkel1~D2.fear+D2.threat+D2.last7days.cases.bund+D2.hh.incom.prev.month.pc+D2.meets.borders-1,
                 dat=sub, weights=wi_rake)
FD.merkel2.vcov <- vcovCL(FD.merkel2, cluster = sub$id_g)

#Table SI3
stargazer(FD.merkel,
          FD.merkel2,
          no.space=TRUE,
          omit.stat=c("f", "rsq", "adj.rsq", "ser"),
          header=FALSE,
          title="First Difference Models",
          order=c(2,1,3,4,5),
          ## float.env="sidewaystable",
          covariate.labels = c("Change Perceived Threat (Lagged)",
                               "Change Anxiety (Lagged)",
                               "Change COVID-19 Incidence (Lagged)",
                               "Change HH Income Previous Month (Lagged)",
                               "Change Policy Congruence: Border Closures (Lagged)"
                               ),
          add.lines = list(c("Number of Respondents",
                             FD.merkel.n, FD.merkel.n)),
          se= list(sqrt(diag(FD.merkel.vcov)),
                   sqrt(diag(FD.merkel2.vcov))),
          dep.var.labels  = c("Merkel approval"),
          dep.var.caption  = "",
          label="tab:FD2",
          digits=2
          )


##############
# Table SI 5 #
##############
dat <- prepared.data
setkey(dat, id_g, Woche)
dat <- dat[,lag.Woche:=c(NA, Woche[-length(Woche)]), id_g]
dat$diff <- dat$Woche-dat$lag.Woche
dat <- dat[,lag.merkel1:=c(NA, merkel1[-length(merkel1)]), id_g]
dat$lag.merkel1 <- ifelse(dat$diff>1, NA, dat$lag.merkel1)
sub <- dat[Woche%in%c(2,3,4,5,6,7,8,9,12,15)] #weeks with lagged Merkel approval
sub <- sub[!is.na(merkel1)&!is.na(lag.merkel1)&!is.na(wi_rake)&!is.na(threat)]
sub <- sub[!is.na(last7days.cases.bund)&!is.na(hh.incom.prev.month.pc)&!is.na(meets.borders)]
sub <- sub[,count:=.N, id_g]
sub <- sub[count>1]
data <- sub
sub <- sub[,list(merkel1, lag.merkel1, threat, fear,
                 last7days.cases.bund, hh.incom.prev.month.pc, meets.borders,
                  id_g, Woche, wi_rake)]
merkel.LDV1 <- weighted.fe.robust(merkel1~lag.merkel1+threat+fear, data=sub, group="id_g", time="Woche", weights="wi_rake")
merkel.LDV2 <- weighted.fe.robust(merkel1~
                               lag.merkel1
                              +threat
                              +fear
                              +last7days.cases.bund
                              +hh.incom.prev.month.pc
                              +meets.borders,
                              data=sub, group="id_g", time="Woche", weights="wi_rake")
stargazer(merkel.LDV1,
          merkel.LDV2,
          no.space=TRUE,
          omit.stat=c("f", "rsq", "adj.rsq", "ser"),
          header=FALSE,
          title="The Effect of Perceived Threat and Anxiety on Merkel Approval (Lagged Dependent Variable Model)",
          covariate.labels = c("Lagged Merkel Approval",
                               "Perceived Threat",
                               "Anxiety",
                               "COVID-19 Incidence",
                               "HH Income Previous Month",
                               "Policy Congruence: Border Closures"),
          add.lines = list(c("Individual Fixed effects", "Yes", "Yes"),
                           c("Number of Respondents",
                             length(unique(sub$id_g)), length(unique(sub$id_g)))),
          se= list(sqrt(diag(merkel.LDV1$clustered.vcov)),
                   sqrt(diag(merkel.LDV2$clustered.vcov))
                   ),
         dep.var.caption  = "",
         dep.var.labels.include = FALSE,
         label="tab:LDV", digits=2
          )



##############
# Table SI 6 #
##############
library(lavaan)

#semTable is currently not available on cran but included in the replication files
#It needs to be installed from soure (see code below).
#You may need to install its dependencies first.

## install.packages("package_not_on_cran/semTable_1.8.tar.gz", repos=NULL, type="source",
##                  dependencies=TRUE)
library(semTable)

dat <- prepared.data
setkey(dat, id_g, Woche)
dat <- dat[,lag.Woche:=c(NA, Woche[-length(Woche)]), id_g]
dat$diff <- dat$Woche-dat$lag.Woche
dat <- dat[,lag.merkel1:=c(NA, merkel1[-length(merkel1)]), id_g]
dat$lag.merkel1 <- ifelse(dat$diff>1, NA, dat$lag.merkel1)
dat <- dat[,lag.fear:=c(NA, fear[-length(fear)]), id_g]
dat$lag.fear <- ifelse(dat$diff>1, NA, dat$lag.fear)
dat <- dat[,lag.threat:=c(NA, threat[-length(threat)]), id_g]
dat$lag.threat <- ifelse(dat$diff>1, NA, dat$lag.threat)
sub <- dat[Woche%in%c(2,3,4,5,6,7,8,9,12,15)] #weeks with lagged Merkel approval
sub <- sub[!is.na(merkel1)&!is.na(lag.merkel1)&!is.na(wi_rake)&!is.na(threat)&!is.na(fear)&!is.na(lag.threat)&!is.na(lag.threat)]
sub <- sub[!is.na(last7days.cases.bund)&!is.na(hh.incom.prev.month.pc)&!is.na(meets.borders)]
sub <- sub[,count:=.N, id_g]
sub <- sub[count>1]
data <- sub
sub <- sub[,list(merkel1, lag.merkel1, threat, fear,
                 lag.threat, lag.fear,
                 last7days.cases.bund, hh.incom.prev.month.pc, meets.borders,
                  id_g, Woche, wi_rake)]
#run model with lavaan
sem.model <-'
    merkel1 ~ AR11*lag.merkel1 + CL21*lag.threat + CL31*lag.fear
    threat ~ AR22*lag.threat + CL12*lag.merkel1 + CL32*lag.fear
    fear ~ AR33*lag.fear + C13*lag.merkel1 + CL23*lag.threat
'
model.fit <- sem(data=sub, model=sem.model)
vlabs <- c ( "merkel1" = "Merkel Approval" , "lag.merkel1" = "Merkel Approval (Lagged)" ,
            "threat" = "Perceived Threat" , "lag.threat" = "Perceived Threat (Lagged)" ,
            "fear" = "Anxiety" , "lag.fear" = "Anxiety (Lagged)")
semtab <- semTable (model.fit, varLabels=vlabs, table.float=TRUE,
                    caption="Cross-lagged panel model", lab="tab_cross", longtable=FALSE)


###########
# Fig SI2 #
###########
rm(list=ls())

#data source:
#https://www.uni-mannheim.de/en/gip/corona-study/methodology/#c192450
dat <- fread("response_rates.csv")
mcs <- dat[Week<20]

#data source:
#https://www.uni-mannheim.de/en/gip/corona-study/methodology/#c373632
aarb <- fread("AARBs.csv")

p2 <- ggplot(data=mcs, aes(y=Total*100, x=Week))
p2 <- p2 + geom_line()
p2 <- p2 + geom_line(linetype="dashed", data=aarb)
p2 <- p2 + theme_bw()
p2 <- p2 + scale_y_continuous("Response Rate (%)\nAARB (%)", limits=c(0, 100), breaks=seq(0,100,10))
p2 <- p2 + scale_x_continuous("MCS Week", limits=c(1, 16), breaks=1:16)
p2


#####################################################################
# Information on R version, packages, and operating system settings #
#####################################################################
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)

## Matrix products: default


## locale:
## [1] LC_COLLATE=German_Germany.utf8  LC_CTYPE=German_Germany.utf8
## [3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
## [5] LC_TIME=German_Germany.utf8

## time zone: Europe/Berlin
## tzcode source: internal

## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base

## other attached packages:
##  [1] semTable_1.8      lavaan_0.6-18     diagis_0.2.3      ggplot2_3.5.1
##  [5] lme4_1.1-35.5     Matrix_1.7-0      car_3.1-2         carData_3.0-5
##  [9] weights_1.0.4     Hmisc_5.1-3       sandwich_3.1-0    stargazer_5.2.3
## [13] data.table_1.15.4

## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1  dplyr_1.1.4       farver_2.1.2      fastmap_1.2.0
##  [5] digest_0.6.36     rpart_4.1.23      lifecycle_1.0.4   cluster_2.1.6
##  [9] survival_3.6-4    gdata_3.0.0       magrittr_2.0.3    compiler_4.4.1
## [13] rlang_1.1.4       tools_4.4.1       utf8_1.2.4        knitr_1.48
## [17] labeling_0.4.3    htmlwidgets_1.6.4 mnormt_2.1.1      plyr_1.8.9
## [21] abind_1.4-5       foreign_0.8-86    withr_3.0.0       purrr_1.0.2
## [25] nnet_7.3-19       grid_4.4.1        stats4_4.4.1      fansi_1.0.6
## [29] jomo_2.7-6        xtable_1.8-4      colorspace_2.1-0  mice_3.16.0
## [33] scales_1.3.0      gtools_3.9.5      iterators_1.0.14  MASS_7.3-60.2
## [37] cli_3.6.3         rmarkdown_2.27    generics_0.1.3    rstudioapi_0.16.0
## [41] minqa_1.2.7       stringr_1.5.1     splines_4.4.1     parallel_4.4.1
## [45] base64enc_0.1-3   vctrs_0.6.5       boot_1.3-30       glmnet_4.1-8
## [49] mitml_0.4-5       Formula_1.2-5     htmlTable_2.4.3   foreach_1.5.2
## [53] tidyr_1.3.1       glue_1.7.0        nloptr_2.1.1      pan_1.9
## [57] codetools_0.2-20  stringi_1.8.4     shape_1.4.6.1     gtable_0.3.5
## [61] quadprog_1.5-8    munsell_0.5.1     stationery_1.1    tibble_3.2.1
## [65] pillar_1.9.0      htmltools_0.5.8.1 R6_2.5.1          evaluate_0.24.0
## [69] pbivnorm_0.6.0    lattice_0.22-6    backports_1.5.0   kutils_1.73
## [73] openxlsx_4.2.5.2  broom_1.0.6       zip_2.3.1         Rcpp_1.0.12
## [77] coda_0.19-4.1     gridExtra_2.3     nlme_3.1-164      checkmate_2.3.1
## [81] xfun_0.45         zoo_1.8-12        pkgconfig_2.0.3
