###This script is for replicating analyses
###Run in R statistical software, version 4.0.3
library(dplyr)
library(tidyr)
library(irr)
library(tidyverse)
library(gplots)
library(RColorBrewer)
library(psych)
library(gridExtra)

all_data<-read.csv("alldata.csv")%>%mutate(X=NULL)
dem<-read.csv("dem.csv")%>%mutate(X=NULL)
origmatrix<-read.csv("origmatrix.csv")%>%mutate(X=NULL)
flwpmatrix<-read.csv("flwpmatrix.csv")%>%mutate(X=NULL)
w1.labels<-read.csv("w1labels.csv")%>%mutate(X=NULL)

###Individual questions table
att_diff_table <- all_data %>%
  gather(key, val, -pernum) %>%
  extract(key, c('question', 'wave'), '(.*)_(orig|flwp)') %>%
  spread(wave, val) %>%
  group_by(question) %>%
  summarize(
    'Mean, Wave 1'=mean(orig, na.rm=T),
    'Mean, Wave 2'=mean(flwp, na.rm=T),
    'Mean absolute difference'=mean(abs(orig-flwp), na.rm=T),
    'Percent exactly matching'=mean(orig == flwp, na.rm=T) * 100,
    'Percent within 1'=mean(abs(orig - flwp) <= 1, na.rm=T) * 100,
    'ICC'=icc(na.omit(cbind(orig, flwp)), model="twoway", type="agreement", unit="single")$value
  ) %>%
  merge(w1.labels, by.x='question', by.y='variable', all.x=T, all.y=F) %>%
  mutate(across(ends_with('difference') | starts_with('Mean') | starts_with('ICC'), ~round(.x, 2))) %>%
  mutate(across(starts_with('Percent'), ~paste0(round(.x), '%'))) %>%
  select(-question) %>%
  relocate(label)
View(att_diff_table)

###Heatmaps of individual questions
##Society is overreacting to the coronavirus
socmapscores<-as.data.frame(cbind(flwpmatrix[,"att_covid_4"], origmatrix[,"att_covid_overreact"]))
colnames(socmapscores)<-c("Wave_2", "Wave_1")
socmapscores<-mutate(socmapscores, Wave_1=(Wave_1+3), Wave_2=(Wave_2+3))
socmapscores<-socmapscores%>%group_by(Wave_2, Wave_1)%>%tally()
socmaptransition<-matrix(data=socmapscores$n,nrow=5,ncol=5,byrow=TRUE)

par(oma=c(2,2,2,2))
heatmap.2(socmaptransition, labRow=c("5","4","3","2","1"), labCol=c("1","2","3","4","5"),
          Rowv=NA, Colv=NA, dendrogram="none", symm = FALSE, revC=TRUE, xlab="Wave 1 Response", ylab="Wave 2 Response",
          main="''Society is overreacting to\nthe coronavirus''", col=colorRampPalette(brewer.pal(4, "Blues"))(25),
          cexRow=2.5, cexCol=2.5, cex.lab=2.5, srtRow=0, srtCol=0, key=FALSE, trace="none", hline="none", vline="none",
          lhei=c(1,5), lwid=c(0.2,5), margins=c(6,6.25))

##Videocalling is a good alternative to in-person business meetings
vidmapscores<-as.data.frame(cbind(flwpmatrix[,"att_tech_1"], origmatrix[,"att_tech_videobusiness"]))
colnames(vidmapscores)<-c("Wave_2", "Wave_1")
vidmapscores<-mutate(vidmapscores, Wave_1=(Wave_1+3), Wave_2=(Wave_2+3))
vidmapscores<-vidmapscores%>%group_by(Wave_2, Wave_1)%>%tally()
vidmaptransition<-matrix(data=vidmapscores$n,nrow=5,ncol=5,byrow=TRUE)

par(oma=c(2,2,2,2))
heatmap.2(vidmaptransition, labRow=c("5","4","3","2","1"), labCol=c("1","2","3","4","5"),
          Rowv=NA, Colv=NA, dendrogram="none", symm = FALSE, revC=TRUE, xlab="Wave 1 Response", ylab="Wave 2 Response",
          main="''Videocalling is a good alternative\nto in-person business meetings''", col=colorRampPalette(brewer.pal(4, "Blues"))(25),
          cexRow=2.5, cexCol=2.5, cex.lab=2.5, srtRow=0, srtCol=0, key=FALSE, trace="none", hline="none", vline="none",
          lhei=c(1,5), lwid=c(0.2,5), margins=c(6,6.25))

###Wave 1 factor analysis
##Factor solution
faorig<-fa(origmatrix,nfactors=6,rotate="varimax", SMC=TRUE)
##Factor loadings
loadingsorig<-print(faorig$loadings, cutoff=0.3, sort=TRUE)
##Factor scores
#weights.thurstone() calculates weights by multiplying the matrix of correlations between items and factors by the inverse correlaiton matrix of items
#weights.thurstone() uses the method of calculating weights that is used in the psych package
weights.thurstone <- function (efa, mtx) {
  mtx.cor <- cor(mtx)
  w <- solve(mtx.cor, efa$loadings)
  return(w)
}
#scores.thurstone() calculates scores by multiplying the matrix of weights by the matrix of observed scores
#scores.thurstone() standardizes both waves of attitudinal data using the first wave of scores to allow for direct comparability
scores.thurstone <- function (mtx, w) {
  return(scale(mtx, center=apply(origmatrix, 2, mean), scale=apply(origmatrix, 2, sd)) %*% w)
}
#Calculating scores using the above functions
weights <- weights.thurstone(faorig, origmatrix)
origscores.origfa <- scores.thurstone(origmatrix, weights)
flwpscores.origfa <- scores.thurstone(flwpmatrix, weights)
#Standardizing scores so that a 1-point change is equal to a 1-standard deviation change for Wave 1; 
#Wave 2 was standardized using Wave 1, therefore, each Wave 2 factor has a mean that is slightly off from 0 and a standard deviation that is slightly off from 1
center<-apply(origscores.origfa, 2, mean)
scale<-apply(origscores.origfa, 2, sd)
origscores.origfa<-scale(origscores.origfa, center=center, scale=scale)
flwpscores.origfa<-scale(flwpscores.origfa, center=center, scale=scale)
#Means and standard deviations for factors in Wave 1 and Wave 2
apply(origscores.origfa, 2, mean)
apply(flwpscores.origfa, 2, mean)
apply(origscores.origfa, 2, sd)
apply(flwpscores.origfa, 2, sd)
#Difference in scores between waves
change.origfa<-as.data.frame(flwpscores.origfa-origscores.origfa)

###Results table for Wave 1 factor analysis
##Mean change in score
meansorig<-c(mean(change.origfa$MR1), mean(change.origfa$MR2), mean(change.origfa$MR3), mean(change.origfa$MR4), mean(change.origfa$MR5), mean(change.origfa$MR6))
##Mean absolute change in score
meansabsorig<-c(mean(abs(change.origfa$MR1)), mean(abs(change.origfa$MR2)), mean(abs(change.origfa$MR3)), mean(abs(change.origfa$MR4)), mean(abs(change.origfa$MR5)), mean(abs(change.origfa$MR6)))
##Percent within one SD of score
smallchangeorig<-c(sum(abs(change.origfa$MR1)<1), sum(abs(change.origfa$MR2)<1), sum(abs(change.origfa$MR3)<1), sum(abs(change.origfa$MR4)<1), sum(abs(change.origfa$MR5)<1), sum(abs(change.origfa$MR6)<1))/nrow(change.origfa)
##Intraclass correlation coefficient
ICC.orig<-function(col){
  values<-cbind(as.data.frame(origscores.origfa)[,col], as.data.frame(flwpscores.origfa)[,col])
  icc<-icc(values, model="twoway", type="agreement", unit="single")
  print(icc$value)
}
ICCorig<-c(ICC.orig(1), ICC.orig(3), ICC.orig(2), ICC.orig(4), ICC.orig(5), ICC.orig(6))
##Table of results
resultsorig<-data.frame(meansorig, meansabsorig, smallchangeorig, ICCorig, row.names=c("Covid-19 concerned", "Pro-videoconferencing", "Environmentalist city lover", "Anti-in-person-shopping", "Anti-working from home", "Home-oriented"))%>%rename("Mean Change in Score"=meansorig, "Mean Absolute Change in Score"=meansabsorig, "Percent within one SD"=smallchangeorig, "ICC"=ICCorig)

###Wave 2 factor analysis
##Factor solution
faflwp<-fa(flwpmatrix,nfactors=6,rotate="varimax", SMC=TRUE)
##Factor loadings
loadingsflwp<-print(faflwp$loadings, cutoff=0.3, sort=TRUE)
##Factor scores
#weights.thurstone() calculates weights by multiplying the matrix of correlations between items and factors by the inverse correlaiton matrix of items
#weights.thurstone() uses the method of calculating weights that is used in the psych package
weights.thurstone <- function (efa, mtx) {
  mtx.cor <- cor(mtx)
  w <- solve(mtx.cor, efa$loadings)
  return(w)
}
#scores.thurstone.flwp() calculates scores by multiplying the matrix of weights by the matrix of observed scores
#scores.thurstone.flwp() standardizes both waves of attitudinal data using the second wave of scores to allow for direct comparability
scores.thurstone.flwp <- function (mtx, w) {
  return(scale(mtx, center=apply(flwpmatrix, 2, mean), scale=apply(flwpmatrix, 2, sd)) %*% w)
}
#Calculating scores using the above functions
weights.flwp <- weights.thurstone(faflwp, flwpmatrix)
origscores.flwpfa <- scores.thurstone.flwp(origmatrix, weights.flwp)
flwpscores.flwpfa <- scores.thurstone.flwp(flwpmatrix, weights.flwp)
#Standardizing scores so that a 1-point change is equal to a 1-standard deviation change for Wave 2; 
#Wave 1 was standardized using Wave 2, therefore, each Wave 1 factor has a mean that is slightly off from 0 and a standard deviation that is slightly off from 1
center.flwp<-apply(flwpscores.flwpfa, 2, mean)
scale.flwp<-apply(flwpscores.flwpfa, 2, sd)
origscores.flwpfa<-scale(origscores.flwpfa, center=center.flwp, scale=scale.flwp)
flwpscores.flwpfa<-scale(flwpscores.flwpfa, center=center.flwp, scale=scale.flwp)
#Means and standard deviations for factors in Wave 1 and Wave 2
apply(origscores.flwpfa, 2, sd)
apply(flwpscores.flwpfa, 2, sd)
apply(origscores.flwpfa, 2, mean)
apply(flwpscores.flwpfa, 2, mean)
#Difference in scores between waves
change.flwpfa<-as.data.frame(flwpscores.flwpfa-origscores.flwpfa)

###Results table for Wave 2 factor analysis
##Mean change in score
meansflwp<-c(mean(change.flwpfa$MR1), mean(change.flwpfa$MR5), mean(change.flwpfa$MR4), mean(change.flwpfa$MR2), mean(change.flwpfa$MR3), mean(change.flwpfa$MR6))
##Mean absolute change in score
meansabsflwp<-c(mean(abs(change.flwpfa$MR1)), mean(abs(change.flwpfa$MR5)), mean(abs(change.flwpfa$MR4)), mean(abs(change.flwpfa$MR2)), mean(abs(change.flwpfa$MR3)), mean(abs(change.flwpfa$MR6)))
##Percent within one SD of score
smallchangeflwp<-c(sum(abs(change.flwpfa$MR1)<1), sum(abs(change.flwpfa$MR5)<1), sum(abs(change.flwpfa$MR4)<1), sum(abs(change.flwpfa$MR2)<1), sum(abs(change.flwpfa$MR3)<1), sum(abs(change.flwpfa$MR6)<1))/nrow(change.flwpfa)
##Intraclass correlation coefficient
ICC.flwp<-function(col){
  values<-cbind(as.data.frame(origscores.flwpfa)[,col], as.data.frame(flwpscores.flwpfa)[,col])
  icc<-icc(values, model="twoway", type="agreement", unit="single")
  print(icc$value)
}
ICCflwp<-c(ICC.flwp(1), ICC.flwp(3), ICC.flwp(2), ICC.flwp(4), ICC.flwp(5), ICC.flwp(6))
##Table of results
resultsflwp<-data.frame(meansflwp, meansabsflwp, smallchangeflwp, ICCflwp, row.names=c("Covid-19 concerned", "Pro-videoconferencing", "Environmentalist city lover", "Anti-in-person-shopping", "Anti-working from home", "Home-oriented"))%>%rename("Mean Change in Score"=meansflwp, "Mean Absolute Change in Score"=meansabsflwp, "Percent within one SD"=smallchangeflwp, "ICC"=ICCflwp)

###Histograms for Wave 1 factor analysis
##historig() creates proprely scaled and labeled histograms for displaying data on the change in attitudinal scores between waves 
historig<-function(df, col, title){
  hist<-ggplot(data=df, aes(x=col))+
    geom_histogram(breaks=seq(-4.25,4.25,by=0.25), color="black", fill="gray")+ggtitle(title)+
    xlab("Change in Score")+ylab("Number of Respondents")+theme_classic()+scale_x_continuous(n.breaks=9)+scale_y_continuous(n.breaks=8)+
    theme(plot.title=element_text(lineheight=1, hjust=0.5, vjust=5, face="bold", size=14), 
          axis.title.x =element_text(size=12, lineheight=1, vjust=-2),
          axis.title.y =element_text(size=12, lineheight=1, vjust=4),
          axis.text = element_text(size=12, color="black"),
          plot.margin=unit(c(1,1,1,1), "cm"))
  print(hist)
}
#Create histograms for each of the 6 factors
histMR1orig<-historig(change.origfa, change.origfa$MR1, "Change in 'Covid-19 Concerned' Factor \nScore Between Waves 1 and 2")
histMR2orig<-historig(change.origfa, change.origfa$MR2, "Change in 'Pro-Videoconferencing' Factor \nScore Between Waves 1 and 2")
histMR3orig<-historig(change.origfa, change.origfa$MR3, "Change in 'Environmentalist City Lover' Factor \nScore Between Waves 1 and 2")
histMR4orig<-historig(change.origfa, change.origfa$MR4, "Change in 'Anti-In Person Shopping' Factor \nScore Between Waves 1 and 2")
histMR5orig<-historig(change.origfa, change.origfa$MR5, "Change in 'Anti-Working from Home' Factor \nScore Between Waves 1 and 2")
histMR6orig<-historig(change.origfa, change.origfa$MR6, "Change in 'Home-Oriented' Factor \nScore Between Waves 1 and 2")
#Combine into a single figure and save
historiggrid<-grid.arrange(histMR1orig, histMR2orig, histMR3orig, histMR4orig, histMR5orig, histMR6orig, ncol=2)
historiggrid
ggsave("Historiggrid.png", plot=historiggrid, device="png", width=12, height=18, units="in")

###Scatterplots for Wave 1 factor analyis
##scatterorig() creates proprely scaled and labeled scatterplots for displaying data on the correlation between attitudinal scores between waves 
scatterorig<-function(col, ICC, title){
  data<-as.data.frame(cbind(origscores.origfa[,col], flwpscores.origfa[,col]))
  scatter<-ggplot(data=data, aes(x=origscores.origfa[,col],
                                  y=flwpscores.origfa[,col]))+
    geom_point(size=1, pch=19, alpha=0.2)+
    xlim(-3.5,3)+ylim(-4,3.25)+geom_abline(slope=1, intercept=0, color="red")+
    geom_text(x=-3, y=2.75, label=ICC)+xlab("Factor Score (Wave 1 Data)")+ylab("Factor Score (Wave 2 Data)")+
    ggtitle(title)+theme_bw()+theme(plot.title = element_text(hjust=0.5,size = 16,face="bold"),
          axis.title = element_text(size=16),
          axis.text = element_text(size=10, color="black"))
  print(scatter)
}
#Create scatterplots for each of the 6 factors
scatterMR1<-scatterorig(1, "ICC=0.83", "Correlation Between 'Covid-19 Concerned' \nScores from Waves 1 and 2 (Wave 1 factor analysis)")
scatterMR2<-scatterorig(3, "ICC=0.62", "Correlation Between 'Pro-Videoconferencing' \nScores from Waves 1 and 2 (Wave 1 factor analysis)")
scatterMR3<-scatterorig(2, "ICC=0.74", "Correlation Between 'Environmentalist City Lover' \nScores from Waves 1 and 2 (Wave 1 factor analysis)")
scatterMR4<-scatterorig(4, "ICC=0.76", "Correlation Between 'Anti-In-Person Shopping' \nScores from Waves 1 and 2 (Wave 1 factor analysis)")
scatterMR5<-scatterorig(5, "ICC=0.67", "Correlation Between 'Anti-Working From Home' \nScores from Waves 1 and 2 (Wave 1 factor analysis)")
scatterMR6<-scatterorig(6, "ICC=0.74", "Correlation Between 'Home-Oriented' \nScores from Waves 1 and 2 (Wave 1 factor analysis)")
#Combine into a single figure and save
scattergrid<-grid.arrange(scatterMR1, scatterMR2, scatterMR3, scatterMR4, scatterMR5, scatterMR6, ncol=2)
scattergrid
ggsave("Scattergrid.png", plot=scattergrid, device="png", width=12, height=18, units="in")

###Results tables by subgroups of respondents
##Create subgroups of demographics
young<-filter(dem, age<35)
midage<-filter(dem, age>34 & age<65)
old<-filter(dem, age>64)
female<-filter(dem, gender=="Female")
male<-filter(dem, gender=="Male")
highered<-filter(dem, education=="Completed graduate degree(s)"|education=="Bachelor's degree(s) or some graduate school")
nohighered<-filter(dem, education=="Some grade/high school"|education=="Completed high school or GED"| education=="Some college or technical school")
longgap<-filter(dem, time>135)
shortgap<-filter(dem, time<135)
urban<-filter(dem, dens>854.4)
suburban<-filter(dem, dens<=854.4 & dens>=39.4)
rural<-filter(dem, dens<39.4)
early<-filter(dem, start<as.POSIXct("2020-07-01 00:00:00"))
late<-filter(dem, start>=as.POSIXct("2020-07-01 00:00:00"))
qualtrics<-filter(dem, source == "qualtrics_online_panel")
cs<-filter(dem, source == "asu_convenience_sample")
asuemail<-filter(dem, source == "asu_email_deployment")
uicemail<-filter(dem, source == "uic_email_deployment")
##Create key for turning Likert-scaled responses into numeric values
att_to_int <- c(
  "Strongly disagree"=-2,
  "Somewhat disagree"=-1,
  "Neutral"=0,
  "Somewhat agree"=1,
  "Strongly agree"=2,
  "Variable not available in datasource"=NA,
  "Seen but unanswered"=NA,
  "Question not displayed to respondent"=NA
)
##Re-create analysis for demographic groups
resultstable<-function(df){
  #Create separate matrices for each wave of data
  orig<-as.matrix(df %>% select(matches('_orig$')) %>% mutate(across(everything(), ~recode(.x, !!!att_to_int))) %>% rename_with(~str_replace(.x, '_orig$', '')))
  flwp<-as.matrix(df %>% select(matches('_flwp$')) %>% mutate(across(everything(), ~recode(.x, !!!att_to_int))) %>% rename_with(~str_replace(.x, '_flwp$', '')) %>% mutate(att_wfh_3=NULL))
  #Score subset using Wave 1 factor analysis
  weights<- weights.thurstone(faorig, orig)
  origscore<-(scale(orig, center=apply(orig, 2, mean), scale=apply(orig, 2, sd))%*%weights)
  flwpscore<-(scale(flwp, center=apply(orig, 2, mean), scale=apply(orig, 2, sd))%*%weights)
  change<-as.data.frame(flwpscore-origscore)
  #Build results table
  means<-c(mean(change$MR1), mean(change$MR2), mean(change$MR3), mean(change$MR4), mean(change$MR5), mean(change$MR6))
  meansabs<-c(mean(abs(change$MR1)), mean(abs(change$MR2)), mean(abs(change$MR3)), mean(abs(change$MR4)), mean(abs(change$MR5)), mean(abs(change$MR6)))
  smallchange<-c(sum(abs(change$MR1)<1), sum(abs(change$MR2)<1), sum(abs(change$MR3)<1), sum(abs(change$MR4)<1), sum(abs(change$MR5)<1), sum(abs(change$MR6)<1))/nrow(change)
  ICC.dem<-function(col){
    values<-cbind(as.data.frame(origscore)[,col], as.data.frame(flwpscore)[,col])
    icc<-icc(values, model="twoway", type="agreement", unit="single")
    print(icc$value)
  }
  ICC<-c(ICC.dem(1), ICC.dem(3), ICC.dem(2), ICC.dem(4), ICC.dem(5), ICC.dem(6))
  results<-data.frame(means, meansabs, smallchange, ICC, row.names=c("Covid-19 concerned", "Pro-videoconferencing", "Environmentalist city lover", "Anti-in-person-shopping", "Anti-working from home", "Home-oriented"))%>%rename("Mean Change in Score"=means, "Mean Absolute Change in Score"=meansabs, "Percent within one SD"=smallchange, "ICC"=ICC)
  print(results)
}
youngtable<-resultstable(young)
midagetable<-resultstable(midage)
oldtable<-resultstable(old)
femaletable<-resultstable(female)
maletable<-resultstable(male)
higheredtable<-resultstable(highered)
nohigheredtable<-resultstable(nohighered)
longgaptable<-resultstable(longgap)
shortgaptable<-resultstable(shortgap)
urbantable<-resultstable(urban)
suburbantable<-resultstable(suburban)
ruraltable<-resultstable(rural)
earlytable<-resultstable(early)
latetable<-resultstable(late)
qualtricstable<-resultstable(qualtrics)
cstable<-resultstable(cs)
asuemailtable<-resultstable(asuemail)
uicemailtable<-resultstable(uicemail)