#split it into two files
#check which packages are needed in each file
#see if functions can go somewhere else
#add version info
#write read me



####################################################################
# 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:                 #
# 25 July 2024 by Roni Lehrer  #
################################

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

## load packages
## Information on package versions can be found at the end of this file
library(haven)
library(plyr)
library(dplyr)
library(tidyverse)
library(data.table)
library(stargazer)
library(plm)
library(sandwich)
library(readstata13)
library(lme4)
library(lmtest)
library(diagis)
library(weights)
library(wCorr)
library(car)
library(lavaan)
library(tidySEM)

#The semTable is currently not available on Cran. But it is included in the replication files. To install run the following code:
## install.packages("additions/semTable_1.8.tar.gz", repos=NULL, type="source",
##                  dependencies=TRUE)
library(semTable)

weighted.fe.robust <- function(formula, data, group, time, weights) {
    demeaned <- copy(data)
    for (i in names(data)) {
        if (i%in%c(group, time, weights)==FALSE) {
            if (class(demeaned[, get(i)])!="character") {
                demeaned[, (i):=get(i)-mean(get(i)), get(group)]
                }
        }
    }
    if (any(grepl("-1", as.character(formula[[3]]), fixed = TRUE))==FALSE) {
     formula <- as.formula(paste0(formula[[2]], formula[[1]], formula[3], "-1"))
    }
    model <- lm(formula, data=demeaned, weights= get(weights))
    model$clustered.vcov <- vcovCL(model, cluster = as.formula(paste0("~", group)))
    return(model)
}

## Read MCS data
dat <- as.data.table(read_dta("ZA7745_v1-0-0.dta"))
dat <- dat[Teilnahme==1] #keep only row of respondents that provide data in a given MCS week

## Adjust date format
dates <- strsplit(dat$dDatum, ".", TRUE)
dat$date <- as.Date(sapply(dates, function(x) paste0(x[3], "-", x[2], "-", x[1])))
rm(dates)


#####################
# VARIABLE RECODING #
#####################

#Altmaier skalometer
dat$merkel1 <- ifelse(dat$SCPX001<0, NA, dat$SCPX001-1)

#Altmaier skalometer
dat$altmaier <- ifelse(dat$SCPX003<0, NA, dat$SCPX003-1)

#Spahn skalometer
dat$spahn <- ifelse(dat$SCPX002<0, NA, dat$SCPX002-1)

#perceived threat
dat$threat <- ifelse(dat$SCBX003<0, NA, dat$SCBX003)/10

#anxiety index
dat$worried <- ifelse(dat$SCBX009<0, NA, dat$SCBX009)
dat$nervous <- ifelse(dat$SCBX011<0, NA, dat$SCBX011)
dat$fear <- (dat$worried+dat$nervous -2)/6

#gender
dat$male <- ifelse(dat$gender_19<0, NA,
              ifelse(dat$gender_19==1, 1, 0))

#age
dat$age.cat <- ifelse(dat$year_of_birth_cat_19<=0, NA,
               ifelse(dat$year_of_birth_cat_19<=6, 1,0))
dat$age.cat <- factor(dat$age.cat, levels=0:1, labels=c("-60", "60+"))

##income per capita in household in previous month
dat$hh.incom.prev.month <- ifelse(dat$Woche%in%c(1,2), dat$SCDX001, #March
                           ifelse(dat$Woche%in%c(3,4,5,6), dat$SCDX005, #April
                           ifelse(dat$Woche%in%c(7,8,9,10,11), dat$SCDX007, #May
                           ifelse(dat$Woche%in%c(12,13,14,15), dat$SCDX008, dat$SCDX009)))) #June, July
dat$hh.incom.prev.month <-  ifelse(dat$hh.incom.prev.month<0, NA,
                            ifelse(dat$hh.incom.prev.month==1, 150,
                            ifelse(dat$hh.incom.prev.month==2, 400,
                            ifelse(dat$hh.incom.prev.month>2&dat$hh.incom.prev.month<15,
                                (dat$hh.incom.prev.month*500)-500, 7500))))
dat$number_hh_members_19 <- ifelse(dat$number_hh_members_19<0, NA, dat$number_hh_members_19)
dat$hh.incom.prev.month.pc <-  dat$hh.incom.prev.month/dat$number_hh_members_19
dat$hh.incom.prev.month.pc <- dat$hh.incom.prev.month.pc/1000

#medical precondition (was asked only once and needs to be assigned to all observation for a respondent)
dat$precondition <- ifelse(dat$SCTX001<0, NA,
                    ifelse(dat$SCTX001==1, 1, 0))


panel.expand <- function(x) {
  un <- unique(x)
  if (all(is.na(un))) {
    return(rep(NA, length(x)))
  } else {
    return(rep(sort(un), length(x)))
  }
}
dat <- dat[,precondition:=panel.expand(precondition), id_g]

#Are government policy and respondent preference for open/closed borders identical?
dat$demand.borders <- ifelse(dat$SCPX006_b<0, NA, dat$SCPX006_b)
dat$meets.borders <- ifelse(is.na(dat$demand.borders), NA,
                            ifelse((dat$demand.borders==1&dat$date<"2020-06-15")|(dat$demand.borders==0&dat$date>="2020-06-15"),
                                   1, 0))

###############################################
# Obtain information from pre-COVID GIP waves #
###############################################
#Merkel satisfaction from GIP Wave 36 (Juli 2018) (ZA6956)
gip <- as.data.table(read.dta13("ZA6956_v1-0-0.dta"))
gip <- gip[,list(id_g, CE36233)]
dat <- merge(dat, gip, by="id_g", all.x=TRUE)
dat$merkel0 <- ifelse(dat$CE36233<0, NA, dat$CE36233-1) #juli 2018
dat$merkel0 <- dat$merkel0/10
rm(gip)

#Government satisfaction from GIP Wave 44 (November 2019) (ZA7593)
gip <- as.data.table(read.dta13("ZA7593_v1-1-0.dta"))
gip <- gip[,list(id_g, CE44152)]
dat <- merge(dat, gip, by="id_g", all.x=TRUE)
dat$sat.breg.nov19 <- ifelse(dat$CE44152<0, NA, dat$CE44152)
dat$sat.breg.nov19 <- (dat$sat.breg.nov19-1)/10
rm(gip)

#################################################
## Merge with data of COVID-19 cases in Germany #
#################################################
load("rki.daily.pop.Rdata")
rki.daily <- unique(rki.daily[!is.na(last7days.cases.bund),list(date, last7days.cases.bund)])
dat <- merge(dat, rki.daily, by="date")
rm(rki.daily)

manipulated <- dat[,list(date, id_g, Woche, wi_rake, merkel0, merkel1,
                         altmaier, spahn, threat, fear, age.cat, male,
                         hh.incom.prev.month.pc, precondition, meets.borders, sat.breg.nov19,
                         last7days.cases.bund)]

#keep weeks in which Merkel approval item was asked
sub <- manipulated[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 stats Merkel Approval (Table SI1) #
############################################
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 <- manipulated
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 2 and  Figure SI1 #
############################
#get demeaned values, weighted and their SEs
dat <- manipulated
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 <- manipulated
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 <- manipulated
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 SI2
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 <- manipulated
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 <- manipulated
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 <- manipulated
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)
coeftest(FD.threat, vcov=vcovCL(FD.threat, cluster = sub$id_g))
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)
coeftest(FD.threat2, vcov=vcovCL(FD.threat2, cluster = sub$id_g))
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)
coeftest(FD.fear, vcov=vcovCL(FD.fear, cluster = sub$id_g))
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)
coeftest(FD.fear2, vcov=vcovCL(FD.fear2, cluster = sub$id_g))
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)
coeftest(FD.merkel, vcov=vcovCL(FD.merkel, cluster = sub$id_g))
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)
coeftest(FD.merkel2, vcov=vcovCL(FD.merkel2, cluster = sub$id_g))
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
          )



##############
# Figure SI2 #
##############


##############
# Table SI 5 #
##############
dat <- manipulated
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")
coeftest(merkel.LDV1, vcov=merkel.LDV1$clustered.vcov)

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")
coeftest(merkel.LDV2, vcov=merkel.LDV2$clustered.vcov)


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 #
##############
dat <- manipulated
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 sources:
#https://www.uni-mannheim.de/media/Einrichtungen/gip/Response_Rates/German_Internet_Panel_Response_Rates.pdf
#https://www.uni-mannheim.de/en/gip/corona-study/methodology/#c192450

dat <- fread("response_rates.csv")
mcs <- dat[Week<20]

#https://www.uni-mannheim.de/gip/corona-studie/methodik/#c192445
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=c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100))
p2 <- p2 + scale_x_continuous("MCS Week", breaks=1:16)
p2
