#-----------------------Load in libraries and set wd-----------------------------------------------------------------------------------------------------------------------

require(readr)
require(formattable)
require(ggplot2)
require(ggpubr)
require(markovchain)
require(plotly)
require(reshape2)
require(dunn.test)

#set working directory
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#load data for contact frequency analysis
data.freq<-read.csv('data_face_contact.csv')

#for transitional probabilities
data.seq<-read.csv('data_sequence.csv')

#----------------------Figure 1, Males vs. Females--------------------------------------------------

#all face contacts
A<-ggplot(data=data.freq[data.freq$contactarea=="all",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  ggtitle("All")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+theme(legend.position="none")

#mouth contacts
B<-ggplot(data=data.freq[data.freq$contactarea=="mouth",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position="none")+
  ggtitle("Mouth")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)

#eye contacts
C<-ggplot(data=data.freq[data.freq$contactarea=="eyes",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="Sex",labels=c("Female","Male"))+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  ggtitle("Eyes")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+
  theme_pubr()+theme(legend.position="none")

#nose contacts
D<-ggplot(data=data.freq[data.freq$contactarea=="nose",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position="none")+
  ggtitle("Nose")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)

#other contacts
E<-ggplot(data=data.freq[data.freq$contactarea=="other",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position="none")+
  ggtitle("Other")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)

windows()
ggarrange(A,B,C,D,E,nrow=1)

summary(data.freq$frequency[data.freq$gender=="Male" & data.freq$contactarea=="nose"])
sd(data.freq$frequency[data.freq$gender=="Male" & data.freq$contactarea=="nose"])

summary(data.freq$frequency[data.freq$gender=="Female" & data.freq$contactarea=="nose"])
sd(data.freq$frequency[data.freq$gender=="Female" & data.freq$contactarea=="nose"])

summary(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea=="nose"])

summary(data.freq$frequency[data.freq$gender=="Male" & data.freq$contactarea=="nose" & data.freq$age=="adult"])
sd(data.freq$frequency[data.freq$gender=="Male" & data.freq$contactarea=="nose" & data.freq$age=="adult"])

summary(data.freq$frequency[data.freq$gender=="Female" & data.freq$contactarea=="nose" & data.freq$age=="adult"])
sd(data.freq$frequency[data.freq$gender=="Female" & data.freq$contactarea=="nose" & data.freq$age=="adult"])


#-----------------------------Figure S1, interaction between age and gender?----------------------------------------------
plot1<-ggplot(data=data.freq[data.freq$contactarea=="all" & data.freq$age=="adult",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  ggtitle("All: Adults")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+
  theme(legend.position="none")

plot2<-ggplot(data=data.freq[data.freq$contactarea=="all" & data.freq$age=="child",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  ggtitle("All: Children")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+
  theme(legend.position="none")
  
plot3<-ggplot(data=data.freq[data.freq$contactarea=="eyes" & data.freq$age=="adult",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  ggtitle("Eyes: Adults")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",textsize = 6,vjust=0.1,
              map_signif_level=FALSE)+
  theme(legend.position="none")
  
plot4<-ggplot(data=data.freq[data.freq$contactarea=="eyes" & data.freq$age=="child",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  ggtitle("Eyes: Children")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+
  theme(legend.position="none")
  
plot5<-ggplot(data=data.freq[data.freq$contactarea=="mouth" & data.freq$age=="adult",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  ggtitle("Mouth: Adults")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+theme(legend.position="none")
  
plot6<-ggplot(data=data.freq[data.freq$contactarea=="mouth" & data.freq$age=="child",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  ggtitle("Mouth: Children")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+theme(legend.position="none")
  
plot7<-ggplot(data=data.freq[data.freq$contactarea=="nose" & data.freq$age=="adult",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  ggtitle("Nose: Adult")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+theme(legend.position="none")
  
plot8<-ggplot(data=data.freq[data.freq$contactarea=="nose" & data.freq$age=="child",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  ggtitle("Nose: Children")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+theme(legend.position="none")

plot9<-ggplot(data=data.freq[data.freq$contactarea=="other" & data.freq$age=="adult",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  ggtitle("Other: Adults")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+theme(legend.position="none")

plot10<-ggplot(data=data.freq[data.freq$contactarea=="other" & data.freq$age=="child",],aes(x=gender,y=frequency))+
  geom_boxplot(aes(x=gender,y=frequency,fill=gender))+
  geom_jitter(aes(x=gender,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  ggtitle("Other: Children")+
  geom_signif(comparisons=list(c("Female","Male")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+theme(legend.position="none")

windows()
ggarrange(plot1,plot2,plot3,plot4,plot5,ncol=5)
windows()
ggarrange(plot6,plot7,plot8,plot9,plot10,ncol=5)

#---------Figure 2, adult vs. children ---------------------------------------------------------------------------------

#all face contacts
A<-ggplot(data=data.freq[data.freq$contactarea=="all",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
               map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("All")

#mouth contacts
B<-ggplot(data=data.freq[data.freq$contactarea=="mouth",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Mouth")

#eye contacts
C<-ggplot(data=data.freq[data.freq$contactarea=="eyes",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="Age",labels=c("Adult","Child"))+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+theme_pubr()+theme(legend.position = "none")+ggtitle("Eyes")

#nose contacts
D<-ggplot(data=data.freq[data.freq$contactarea=="nose",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Nose")

#other contacts
E<-ggplot(data=data.freq[data.freq$contactarea=="other",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Other")

windows()
ggarrange(A,B,C,D,E,nrow=1)


#---------Figure S2, adult vs. children stratified by eating vs. non-eating ---------------------------------------------------------------------------------

A<-ggplot(data=data.freq[data.freq$contactarea=="all" & data.freq$eatonly=="Eating",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("All: Eating")

#mouth contacts: highly significantly different
B<-ggplot(data=data.freq[data.freq$contactarea=="mouth" & data.freq$eatonly=="Eating",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Mouth: Eating")

#eye contacts: highly significantly different
C<-ggplot(data=data.freq[data.freq$contactarea=="eyes" & data.freq$eatonly=="Eating",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="Age",labels=c("Adult","Child"))+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+theme_pubr()+
  theme(legend.position = "none")+ggtitle("Eyes: Eating")

#nose contacts: highly significantly different
D<-ggplot(data=data.freq[data.freq$contactarea=="nose" & data.freq$eatonly=="Eating",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Nose: Eating")

#other contacts: highly significantly different
E<-ggplot(data=data.freq[data.freq$contactarea=="other" & data.freq$eatonly=="Eating",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Other: Eating")

#all face contacts: very sigificant (both t test and wilcoxon rank sum)
G<-ggplot(data=data.freq[data.freq$contactarea=="all" & data.freq$eatonly=="Non-eating",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("All: Non-eating")

#mouth contacts: highly significantly different
H<-ggplot(data=data.freq[data.freq$contactarea=="mouth" & data.freq$eatonly=="Non-eating",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Mouth: Non-eating")

#eye contacts: highly significantly different
I<-ggplot(data=data.freq[data.freq$contactarea=="eyes" & data.freq$eatonly=="Non-eating",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="Age",labels=c("Adult","Child"))+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+theme_pubr()+
  theme(legend.position = "none")+ggtitle("Eyes: Non-eating")

#nose contacts: highly significantly different
J<-ggplot(data=data.freq[data.freq$contactarea=="nose" & data.freq$eatonly=="Non-eating",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Nose: Non-eating")

#other contacts: highly significantly different
K<-ggplot(data=data.freq[data.freq$contactarea=="other" & data.freq$eatonly=="Non-eating",],aes(x=age,y=frequency))+
  geom_boxplot(aes(x=age,y=frequency,fill=age))+
  geom_jitter(aes(x=age,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="")+
  scale_x_discrete(name="",labels=c("Adult","Child"))+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none")+
  geom_signif(comparisons=list(c("adult","child")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Other: Non-eating")


windows()
ggarrange(A,B,C,D,E,nrow=1,ncol=5)
windows()
ggarrange(G,H,I,J,K,nrow=1,ncol=5)

#----------locations------------------------------------------------------------------------------------------

#all contacts - significant difference
kruskal.test(frequency ~ location, data = data.freq[data.freq$contactarea=="all" & data.freq$age=="adult" & data.freq$eatonly=="Eating",])
kruskal.test(frequency ~ location, data = data.freq[data.freq$contactarea=="all" & data.freq$age=="adult" & data.freq$eatonly=="Non-eating",])

dunn.all<-dunn.test(data.freq$frequency[data.freq$contactarea=="all" & data.freq$eatonly=="Non-eating" & data.freq$age=="adult"],data.freq$location[data.freq$contactarea=="all" & data.freq$eatonly=="Non-eating" & data.freq$age=="adult"],
          method="hochberg")

A<-ggplot(data=data.freq[data.freq$contactarea=="all" & data.freq$age=="adult",])+
  geom_violin(aes(x=location,y=frequency,fill=location))+
  geom_jitter(aes(x=location,y=frequency),position=position_jitter(width = .2,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="Location")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+theme_pubr()+
  theme(axis.text.x=element_text(angle=45),legend.position = "none")

#mouth contacts - significant difference
kruskal.test(frequency ~ location, data = data.freq[data.freq$contactarea=="mouth" & data.freq$age=="adult" & data.freq$eatonly=="Eating",])
kruskal.test(frequency ~ location, data = data.freq[data.freq$contactarea=="mouth" & data.freq$age=="adult" & data.freq$eatonly=="Non-eating",])
dunn.mouth<-dunn.test(data.freq$frequency[data.freq$contactarea=="mouth" & data.freq$age=="adult" & data.freq$eatonly=="Non-eating"],data.freq$location[data.freq$contactarea=="mouth" & data.freq$age=="adult" & data.freq$eatonly=="Non-eating"],
          method="Hochberg")

B<-ggplot(data=data.freq[data.freq$contactarea=="mouth" & data.freq$age=="adult",])+
  geom_violin(aes(x=location,y=frequency,fill=location))+
  geom_jitter(aes(x=location,y=frequency),position=position_jitter(width = .3,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="Location")+
  scale_x_discrete(name="")+theme_pubr()+
  scale_y_continuous(name="Frequency")+
  theme(axis.text.x=element_text(angle=45),legend.position = "none")

#eyes contacts
kruskal.test(frequency ~ location, data = data.freq[data.freq$contactarea=="eyes" & data.freq$age=="adult" & data.freq$eatonly=="Eating",])
kruskal.test(frequency ~ location, data = data.freq[data.freq$contactarea=="eyes" & data.freq$age=="adult" & data.freq$eatonly=="Non-eating",])
dunn.eye<-dunn.test(data.freq$frequency[data.freq$contactarea=="eyes" & data.freq$age=="adult" & data.freq$eatonly=="Non-eating"],data.freq$location[data.freq$contactarea=="eyes" & data.freq$age=="adult" & data.freq$eatonly=="Non-eating"],
                    method="Hochberg")

C<-ggplot(data=data.freq[data.freq$contactarea=="eyes" & data.freq$age=="adult",])+
  geom_violin(aes(x=location,y=frequency,fill=location))+
  geom_jitter(aes(x=location,y=frequency),position=position_jitter(width = .3,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="Location")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+theme_pubr()+
  theme(axis.text.x=element_text(angle=45),legend.position = "none")

#nose contacts - no significant difference
kruskal.test(frequency ~ location, data = data.freq[data.freq$contactarea=="nose" & data.freq$age=="adult" & data.freq$eatonly=="Eating",])
kruskal.test(frequency ~ location, data = data.freq[data.freq$contactarea=="nose" & data.freq$age=="adult" & data.freq$eatonly=="Non-eating",])

D<-ggplot(data=data.freq[data.freq$contactarea=="nose" & data.freq$age=="adult",])+
  geom_violin(aes(x=location,y=frequency,fill=location))+
  geom_jitter(aes(x=location,y=frequency),position=position_jitter(width = .3,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="Location")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+theme_pubr()+
  theme(axis.text.x=element_text(angle=45),legend.position = "none")


#other contacts - significant difference
kruskal.test(frequency ~ location, data = data.freq[data.freq$contactarea=="other" & data.freq$age=="adult" & data.freq$eatonly=="Eating",])
kruskal.test(frequency ~ location, data = data.freq[data.freq$contactarea=="other" & data.freq$age=="adult" & data.freq$eatonly=="Non-eating",])
dunn.other<-dunn.test(data.freq$frequency[data.freq$contactarea=="other" & data.freq$age=="adult" & data.freq$eatonly=="Non-eating"],data.freq$location[data.freq$contactarea=="other" & data.freq$age=="adult" & data.freq$eatonly=="Non-eating"],
          method="Hochberg")

E<-ggplot(data=data.freq[data.freq$contactarea=="other" & data.freq$age=="adult",])+
  geom_violin(aes(x=location,y=frequency,fill=location))+
  geom_jitter(aes(x=location,y=frequency),position=position_jitter(width = .3,height=0.0),alpha=0.5)+
  scale_fill_discrete(name="Location")+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+theme_pubr()+
  theme(axis.text.x=element_text(angle=45),legend.position="none")

A<-annotate_figure(A,top=text_grob("All",face="bold",size=14))
B<-annotate_figure(B,top=text_grob("Mouth",face="bold",size=14))
C<-annotate_figure(C,top=text_grob("Eyes",face="bold",size=14))
D<-annotate_figure(D,top=text_grob("Nose",face="bold",size=14))
E<-annotate_figure(E,top=text_grob("Other",face="bold",size=14))

#--------------------Figure S3, locations---------------------
windows()
ggarrange(A,B,C,D,E,common.legend = TRUE)


#--------------------Figure 4, heat maps for dunns test results--------------------------

#hand-to-head
library(plotly)
handhead.dunn<-matrix(nrow=9,ncol=9)
mouth.dunn<-matrix(nrow=9,ncol=9)
eye.dunn<-matrix(nrow=9,ncol=9)
other.dunn<-matrix(nrow=9,ncol=9)

rownames(handhead.dunn)<-c("airport","bar","church","classroom","food court","museum",
                           "public library","sporting event","university library")
colnames(handhead.dunn)<-c("bar","church","classroom","foood court",
                           "museum","public library","sporting event","university library","zoo")
rownames(mouth.dunn)<-c("airport","bar","church","classroom","food court","museum",
                           "public library","sporting event","university library")
colnames(mouth.dunn)<-c("bar","church","classroom","foood court",
                           "museum","public library","sporting event","university library","zoo")

rownames(eye.dunn)<-c("airport","bar","church","classroom","food court","museum",
                        "public library","sporting event","university library")
colnames(eye.dunn)<-c("bar","church","classroom","foood court",
                        "museum","public library","sporting event","university library","zoo")

rownames(other.dunn)<-c("airport","bar","church","classroom","food court","museum",
                           "public library","sporting event","university library")
colnames(other.dunn)<-c("bar","church","classroom","foood court",
                           "museum","public library","sporting event","university library","zoo")

handhead.dunn[,1]<-c(0.53,0.38,0.0004,1,1,0.99,0.0545,1,0.0198)
handhead.dunn[,2]<-c(NA,0.45,0.48,0.72,0.03,1,1,1,1)
handhead.dunn[,3]<-c(NA,NA,0.54,0.63,0.016,1,1,1,1)
handhead.dunn[,4]<-c(NA,NA,NA,0.049,0,0.16,1,0.0062,1)
handhead.dunn[,5]<-c(NA,NA,NA,NA,0.81,0.96,0.30,1,0.19)
handhead.dunn[,6]<-c(NA,NA,NA,NA,NA,0.08,0.0009,0.83,0.0002)
handhead.dunn[,7]<-c(NA,NA,NA,NA,NA,NA,1,1,1)
handhead.dunn[,8]<-c(NA,NA,NA,NA,NA,NA,NA,0.32,1)
handhead.dunn[,9]<-c(NA,NA,NA,NA,NA,NA,NA,NA,0.14)

mouth.dunn[,1]<-c(0.85,1,0.47,0.88,1,1,1,1,0.41)
mouth.dunn[,2]<-c(NA,1,1,0.08,0.08,1,1,1,1)
mouth.dunn[,3]<-c(NA,NA,1,0.23,0.47,1,1,1,1)
mouth.dunn[,4]<-c(NA,NA,NA,0.06,0.03,1,1,1,0.91)
mouth.dunn[,5]<-c(NA,NA,NA,NA,1,0.12,0.11,0.36,0.05)
mouth.dunn[,6]<-c(NA,NA,NA,NA,NA,0.14,0.12,0.97,0.03)
mouth.dunn[,7]<-c(NA,NA,NA,NA,NA,NA,0.47,1,1)
mouth.dunn[,8]<-c(NA,NA,NA,NA,NA,NA,NA,1,1)
mouth.dunn[,9]<-c(NA,NA,NA,NA,NA,NA,NA,NA,1)

eye.dunn[,1]<-c(1,1,0.48,1,1,1,1,1,1)
eye.dunn[,2]<-c(NA,1,1,1,1,1,1,0.61,1)
eye.dunn[,3]<-c(NA,NA,1,1,1,1,0.5,1,1)
eye.dunn[,4]<-c(NA,NA,NA,1,0.09,0.4,1,0.0055,1)
eye.dunn[,5]<-c(NA,NA,NA,NA,0.98,1,1,1,1)
eye.dunn[,6]<-c(NA,NA,NA,NA,NA,1,1,1,1)
eye.dunn[,7]<-c(NA,NA,NA,NA,NA,NA,1,1,1)
eye.dunn[,8]<-c(NA,NA,NA,NA,NA,NA,NA,1,1)
eye.dunn[,9]<-c(NA,NA,NA,NA,NA,NA,NA,NA,0.47)

other.dunn[,1]<-c(1,0.027,0.0018,1,1,1,0.15,1,1)
other.dunn[,2]<-c(NA,1,0.33,0.94,0.24,1,1,1,1)
other.dunn[,3]<-c(NA,NA,1,1,0.0006,0.16,1,0.52,1)
other.dunn[,4]<-c(NA,NA,NA,1,0,0.016,1,0.08,0.5)
other.dunn[,5]<-c(NA,NA,NA,NA,1,1,1,1,0.49)
other.dunn[,6]<-c(NA,NA,NA,NA,NA,1,0.0058,0.60,0.13)
other.dunn[,7]<-c(NA,NA,NA,NA,NA,NA,0.59,1,1)
other.dunn[,8]<-c(NA,NA,NA,NA,NA,NA,NA,1,1)
other.dunn[,9]<-c(NA,NA,NA,NA,NA,NA,NA,NA,1)

windows()
library(plotly)

f <- list(
  family = "Arial, monospace",
  size = 18,
  color = "black")

A<-plot_ly(z=handhead.dunn,type="heatmap",x=c("airport","bar","church","classroom","food court","museum",
                                           "public library","sporting event","university library"),
                                       y=c("bar","church","classroom","foood court",
                                          "museum","public library","sporting event","university library","zoo"),
        zmin=0,zmax=0.03,reversescale=TRUE,
        colorbar=list(title="Dunn's Test p-values"),showscale=FALSE) %>%
        layout(annotations=list(text="A",y=0.1,x=3.5,font=f,showarrow=FALSE),showlegend=FALSE,xaxis=list(side="top",showgrid=FALSE),yaxis=list(showgrid=FALSE))

B<-plot_ly(z=eye.dunn,type="heatmap",x=c("airport","bar","church","classroom","food court","museum",
                                           "public library","sporting event","university library"),
           y=c("bar","church","classroom","foood court",
               "museum","public library","sporting event","university library","zoo"),
           zmin=0,zmax=0.03,reversescale=TRUE,
           colorbar=list(title="Dunn's Test p-values"),showscale=TRUE) %>%
  layout(annotations=list(text="B",y=0.1,x=3.5,font=f,showarrow=FALSE),xaxis=list(side="top",showgrid=FALSE),yaxis=list(showgrid=FALSE))

C<-plot_ly(z=other.dunn,type="heatmap",x=c("airport","bar","church","classroom","food court","museum",
                                         "public library","sporting event","university library"),
           y=c("bar","church","classroom","foood court",
               "museum","public library","sporting event","university library","zoo"),
           zmin=0,zmax=0.03,reversescale=TRUE,
           colorbar=list(title="Dunn's Test p-values"),showscale=FALSE) %>%
  layout(annotations=list(text="C",y=0.1,x=3.5,font=f,showarrow=FALSE),xaxis=list(side="top",showgrid=FALSE),yaxis=list(showgrid=FALSE))


subplot(A,B,C,nrows=3,margin=0.05)


#------------------------------Figure 3, Eating only vs. no Eating comparison---------------------------------------------------------

#all face contacts:significantly different
A.3<-ggplot(data=data.freq[data.freq$contactarea=="all",],aes(x=eatonly,y=frequency))+
  geom_boxplot(aes(x=eatonly,y=frequency,fill=eatonly))+
  geom_jitter(aes(x=eatonly,y=frequency),position=position_jitter(width = .1,height=0.0),alpha=0.5)+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none",plot.title = element_text(hjust = 0.5))+
  geom_signif(comparisons=list(c("Non-eating","Eating")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("All")
  
  
#mouth contacts: statistically significantly different
B.3<-ggplot(data=data.freq[data.freq$contactarea=="mouth",],aes(x=eatonly,y=frequency))+
  geom_boxplot(aes(x=eatonly,y=frequency,fill=eatonly))+
  geom_jitter(aes(x=eatonly,y=frequency),position=position_jitter(width = .1,height=0.0),alpha=0.5)+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none",plot.title = element_text(hjust = 0.5))+
  geom_signif(comparisons=list(c("Non-eating","Eating")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Mouth")

#nose contacts: statistically significantly different
C.3<-ggplot(data=data.freq[data.freq$contactarea=="nose",],aes(x=eatonly,y=frequency))+
  geom_boxplot(aes(x=eatonly,y=frequency,fill=eatonly))+
  geom_jitter(aes(x=eatonly,y=frequency),position=position_jitter(width = .1,height=0.0),alpha=0.5)+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  scale_fill_discrete(name="Eating Location",labels=c("No","Yes"))+
  geom_signif(comparisons=list(c("Non-eating","Eating")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+theme_pubr()+theme(legend.position="none",plot.title = element_text(hjust = 0.5))+ggtitle("Nose")

#eye contacts: statistically significantly different
D.3<-ggplot(data=data.freq[data.freq$contactarea=="eyes",],aes(x=eatonly,y=frequency))+
  geom_boxplot(aes(x=eatonly,y=frequency,fill=eatonly))+
  geom_jitter(aes(x=eatonly,y=frequency),position=position_jitter(width = .1,height=0.0),alpha=0.5)+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none",plot.title = element_text(hjust = 0.5))+
  geom_signif(comparisons=list(c("Non-eating","Eating")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Eyes")

#other: statistically significantly different
E.3<-ggplot(data=data.freq[data.freq$contactarea=="other",],aes(x=eatonly,y=frequency))+
  geom_boxplot(aes(x=eatonly,y=frequency,fill=eatonly))+
  geom_jitter(aes(x=eatonly,y=frequency),position=position_jitter(width = .1,height=0.0),alpha=0.5)+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Frequency")+
  theme_pubr()+theme(legend.position = "none",plot.title = element_text(hjust = 0.5))+
  geom_signif(comparisons=list(c("Non-eating","Eating")),
              test="wilcox.test",
              map_signif_level=FALSE,textsize = 6,vjust=0.1)+ggtitle("Other")

library(ggpubr)

windows()
ggarrange(A.3,B.3,C.3,D.3,E.3,nrow=1)

median(data.freq$frequency[data.freq$eatonly=="Eating" & data.freq$contactarea=="mouth"])
sd(data.freq$frequency[data.freq$eatonly=="Eating" & data.freq$contactarea=="mouth"])

median(data.freq$frequency[data.freq$eatonly=="Non-eating" & data.freq$contactarea=="mouth"])
sd(data.freq$frequency[data.freq$eatonly=="Non-eating" & data.freq$contactarea=="mouth"])

#------------------------------------------Table 1, summary statistics ----------------------------------------------------------

min<-rep(0,10)
max<-rep(0,10)
mean.save<-rep(0,10)
sd.save<-rep(0,10)
percentile.5<-rep(0,10)
percentile.25<-rep(0,10)
percentile.50<-rep(0,10)
percentile.75<-rep(0,10)
percentile.95<-rep(0,10)
percentile.99<-rep(0,10)
type<-c("mouth","eyes","nose","other","all")

#Children vs. Adults (Eating) - 1st 5 values are for children, next 5 values are for adults in each vector
for(i in 1:5){
  min[i]<-min(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"])
  max[i]<-max(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"])
  mean.save[i]<-mean(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"])
  sd.save[i]<-sd(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"])
  percentile.5[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.05)
  percentile.25[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.25)
  percentile.50[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.50)
  percentile.75[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.75)
  percentile.95[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.95)
  percentile.99[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.99)
  
  min[i+5]<-min(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"])
  max[i+5]<-max(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"])
  mean.save[i+5]<-mean(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"])
  sd.save[i+5]<-sd(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"])
  percentile.5[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.05)
  percentile.25[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.25)
  percentile.50[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.50)
  percentile.75[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.75)
  percentile.95[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.95)
  percentile.99[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Eating"],probs=0.99)
}


#Children vs. Adults (Non-Eating) - 1st 5 values are for children, next 5 values are for adults in each vector
min<-rep(0,10)
max<-rep(0,10)
mean.save<-rep(0,10)
sd.save<-rep(0,10)
percentile.5<-rep(0,10)
percentile.25<-rep(0,10)
percentile.50<-rep(0,10)
percentile.75<-rep(0,10)
percentile.95<-rep(0,10)
percentile.99<-rep(0,10)
type<-c("mouth","eyes","nose","other","all")

for(i in 1:5){
  min[i]<-min(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"])
  max[i]<-max(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"])
  mean.save[i]<-mean(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"])
  sd.save[i]<-sd(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"])
  percentile.5[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.05)
  percentile.25[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.25)
  percentile.50[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.50)
  percentile.75[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.75)
  percentile.95[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.95)
  percentile.99[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.99)
  
  min[i+5]<-min(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"])
  max[i+5]<-max(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"])
  mean.save[i+5]<-mean(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"])
  sd.save[i+5]<-sd(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"])
  percentile.5[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.05)
  percentile.25[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.25)
  percentile.50[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.50)
  percentile.75[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.75)
  percentile.95[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.95)
  percentile.99[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i] & data.freq$eatonly=="Non-eating"],probs=0.99)
}


#Children vs. Adults (All macro-activities) - 1st 5 values are for children, next 5 values are for adults in each vector
min<-rep(0,10)
max<-rep(0,10)
mean.save<-rep(0,10)
sd.save<-rep(0,10)
percentile.5<-rep(0,10)
percentile.25<-rep(0,10)
percentile.50<-rep(0,10)
percentile.75<-rep(0,10)
percentile.95<-rep(0,10)
percentile.99<-rep(0,10)

type<-c("mouth","eyes","nose","other","all")

for(i in 1:5){
  min[i]<-min(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i]])
  max[i]<-max(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i]])
  mean.save[i]<-mean(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i]])
  sd.save[i]<-sd(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i]])
  percentile.5[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i]],probs=0.05)
  percentile.25[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i]],probs=0.25)
  percentile.50[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i]],probs=0.50)
  percentile.75[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i]],probs=0.75)
  percentile.95[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i]],probs=0.95)
  percentile.99[i]<-quantile(data.freq$frequency[data.freq$age=="child" & data.freq$contactarea==type[i]],probs=0.99)
  
  min[i+5]<-min(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i]])
  max[i+5]<-max(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i]])
  mean.save[i+5]<-mean(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i]])
  sd.save[i+5]<-sd(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i]])
  percentile.5[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i]],probs=0.05)
  percentile.25[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i]],probs=0.25)
  percentile.50[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i]],probs=0.50)
  percentile.75[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i]],probs=0.75)
  percentile.95[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i]],probs=0.95)
  percentile.99[i+5]<-quantile(data.freq$frequency[data.freq$age=="adult" & data.freq$contactarea==type[i]],probs=0.99)
}


#---------------------------------Table S2: Transitional Probability Matrices -------------------------------------------------------------
library(markovchain)

list.adult<-list()
list.child<-list()
list.adult.eating<-list()
list.adult.noneating<-list()
list.child.eating<-list()
list.child.noneating<-list()

for(i in 1:max(data.seq$observation)){
  #putting each observation into a list so that transitions between
  #the last observation of someone and the first observation of the next person
  #are not mistakenly counted
  
  list.adult[[i]]<-data.seq$Area.contacted[data.seq$observation==i & data.seq$Age=="adult"]
  list.child[[i]]<-data.seq$Area.contacted[data.seq$observation==i & data.seq$Age=="child"]
  
  list.adult.eating[[i]]<-data.seq$Area.contacted[data.seq$observation==i & data.seq$Age=="adult" & data.seq$eatonly=="yes"]
  list.adult.noneating[[i]]<-data.seq$Area.contacted[data.seq$observation==i & data.seq$Age=="adult" & data.seq$eatonly=="no"]
  
  list.child.eating[[i]]<-data.seq$Area.contacted[data.seq$observation==i & data.seq$Age=="child" & data.seq$eatonly=="yes"]
  list.child.noneating[[i]]<-data.seq$Area.contacted[data.seq$observation==i & data.seq$Age=="child" & data.seq$eatonly=="no"]
}

adults<-markovchainFit(list.adult)
children<-markovchainFit(list.child)

#adults eating vs. noneating
eating.adult<-markovchainFit(list.adult.eating)
noneating.adult<-markovchainFit(list.adult.noneating)

#children eating vs. noneating
eating.child<-markovchainFit(list.child.eating)
noneating.child<-markovchainFit(list.child.noneating)


f <- list(
  family = "Arial, monospace",
  size = 20,
  color = "black")


adults<-melt(adults$estimate@transitionMatrix)
a<-plot_ly(adults, y = ~Var1, x = ~Var2, z = ~value, colors="Blues",colorbar=list(
  title='Probability',showlegend=FALSE
)) %>%
  add_heatmap()%>%
  colorbar(limits = c(0, 1))%>%
  layout(annotations=list(text="A",x=1.5,y=3.7,showarrow=FALSE,font=f),xaxis = list(title="To",showgrid=FALSE,size=15), yaxis = list(title="From",showgrid=FALSE,size=15)) %>%
  hide_colorbar()


children<-melt(children$estimate@transitionMatrix)
b<-plot_ly(children, y = ~Var1, x = ~Var2,colors="Blues", z = ~value, colorbar=list(
  title='Probability',showlegend=FALSE
)) %>%
  add_heatmap()%>%
  colorbar(limits = c(0, 1),x=4)%>%
  layout(annotations=list(text="B",x=1.5,y=3.7,showarrow=FALSE,font=f),xaxis = list(title="To",showgrid=FALSE,size=15), yaxis = list(title="From",showgrid=FALSE,size=15)) %>%
  hide_colorbar()

#-----------------------------

eating.adult<-melt(eating.adult$estimate@transitionMatrix)
noneating.adult<-melt(noneating.adult$estimate@transitionMatrix)


c<-plot_ly(eating.adult, y = ~Var1, x = ~Var2,colors="Blues", z = ~value, colorbar=list(
  title='Probability',showlegend=FALSE
)) %>%
  add_heatmap()%>%
  colorbar(limits = c(0, 1))%>%
  layout(annotations=list(text="C",x=1.5,y=3.7,showarrow=FALSE,font=f),xaxis = list(title="To",showgrid=FALSE), yaxis = list(title="From",showgrid=FALSE)) %>%
  hide_colorbar()


d<-plot_ly(noneating.adult, y = ~Var1, x = ~Var2,colors="Blues", z = ~value, colorbar=list(
  title='Probability',showlegend=TRUE
)) %>%
  add_heatmap()%>%
  colorbar(limits = c(0, 1),x=4)%>%
  layout(annotations=list(text="D",x=1.5,y=3.7,showarrow=FALSE,font=f),xaxis = list(title="To",showgrid=FALSE), yaxis = list(title="From",showgrid=FALSE)) %>%
  hide_colorbar()

#------------------------------------------------------------

eating.child<-melt(eating.child$estimate@transitionMatrix)
noneating.child<-melt(noneating.child$estimate@transitionMatrix)

e<-plot_ly(eating.child, y = ~Var1, x = ~Var2, colors="Blues",z = ~value, colorbar=list(
  title='Probability',showlegend=FALSE
)) %>%
  add_heatmap()%>%
  colorbar(limits = c(0, 1))%>%
  layout(annotations=list(text="E",x=1.5,y=3.7,showarrow=FALSE,font=f),xaxis = list(title="To",showgrid=FALSE), yaxis = list(title="From",showgrid=FALSE)) %>%
  hide_colorbar()

f<-plot_ly(noneating.child, y = ~Var1, x = ~Var2, colors="Blues",z = ~value, colorbar=list(
  title='Probability'
)) %>%
  add_heatmap()%>%
  colorbar(limits = c(0, 1),x=4)%>%
  layout(annotations=list(text="F",x=1.5,y=3.7,showarrow=FALSE,font=f),xaxis = list(title="To",showgrid=FALSE), yaxis = list(title="From",showgrid=FALSE))

#all plots together
subplot(a,b,c,d,e,f,shareX=TRUE,shareY=TRUE,nrows=3)