---
title: Visualizing Uncertain Tropical Cyclone Predictions using Representative Samples
  from Ensembles of Forecast Tracks
author: 'Le Liu, Lace M. K. Padilla, Sarah Creem-Regehr, Donald House '
date: 'Updated: June 28th, 2018'
output:
  html_document:
    code_folding: hide
    fig_height: 6
    fig_width: 8
    highlight: pygments
    theme: yeti
    toc: yes
    toc_float: yes
  pdf_document:
    toc: yes
  word_document:
    toc: yes
always_allow_html: yes
---

<hr/>
<H4>- <a href="http://lacepadilla.com/Downloads/publications/vispathsample_2018/Analysis.html"target="_blank">Click here </a> to view an online version of this analysis.</h4> 

<H4>- Download the data files and code to generate the analysis <a href="http://lacepadilla.com/Downloads/publications/vispathsample_2018/vispathsample_2018.zip"target="_blank">here.</a> </h4>

<h4>- <a href="https://bitbucket.org/leliu712/hurricane-paths-ensemble-sampling.git"target="_blank">Link to the Bitbucket repository</a> that contains the research codes and pseudocodes of our technique.

<hr/>

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, echo=FALSE, warning=FALSE, message=FALSE)


library(lme4)
library(ggplot2)
library(ggthemes)
source("Functions.R")
library(car)
library(plyr)
library(tidyverse)
library(kableExtra)
library(MuMIn)

```


```{r processing, include=FALSE}


#read in data files
total <- read.csv("Data/total.csv", check.names=FALSE)
confidence <- read.csv("Data/confidence.csv", check.names=FALSE)
hurricane <- read.csv("Data/medianAndMean.csv", check.names=FALSE)


#dummy coding time
total$TimeD <- ifelse(total$Time == "24", 0, ifelse(total$Time == "48", 1, 9))
confidence$TimeD <- ifelse(confidence$Time == "24", 0, ifelse(confidence$Time == "48", 1, 9))

#creating codes for distance
#change the location values to actual kilometer values.
total$Distance[total$Location == -178] <- 1.78
total$Distance[total$Location == -11] <- 1.1 
total$Distance[total$Location == -6] <- .6
total$Distance[total$Location == -2] <- .2
total$Distance[total$Location == 2] <- .2
total$Distance[total$Location == 6] <- .6
total$Distance[total$Location == 11] <- 1.1
total$Distance[total$Location == 178] <- 1.78

confidence$Distance[confidence$Location == -178] <- 1.78
confidence$Distance[confidence$Location == -11] <- 1.1 
confidence$Distance[confidence$Location == -6] <- .6
confidence$Distance[confidence$Location == -2] <- .2
confidence$Distance[confidence$Location == 2] <- .2
confidence$Distance[confidence$Location == 6] <- .6
confidence$Distance[confidence$Location == 11] <- 1.1
confidence$Distance[confidence$Location == 178] <- 1.78

#dummy coding Vis Type (Note that one trial in the 15-track condition was coded as 16. Therefore 16 is also included.) 
total$Visualization  <- total$Vis
total$Vis <- ifelse(total$Visualization == "7", "0seven", 
                         ifelse(total$Visualization == "15" & total$Annotation== "0", "15",
                                ifelse(total$Visualization == "16" & total$Annotation == "0", "15",
                                  ifelse(total$Visualization == "63", "63", 
                                       ifelse(total$Annotation == "1", "Full", 999)))))


#dummy coding Vis Type for the confidence ratings (Note that one trial in the 63-track condition was coded as 62.) 

confidence$Visualization  <- confidence$Vis
confidence$Vis <- ifelse(confidence$Visualization == "7", "0seven", 
                         ifelse(confidence$Visualization == "15" & total$Annotation == "0", "15", 
                                ifelse(confidence$Visualization == "63", "63",
                                       ifelse(confidence$Visualization == "62", "63",
                                          ifelse(confidence$Annotation == "1", "Full", 999)))))



#subsetting the data to look at the influence of annotation
totalAn <- subset(total, Vis == "15" | Vis == "Full")
confidenceAn <- subset(confidence, Vis == "15" | Vis == "Full")

#subsetting the data to remove the annotated conditions. 
totalNum <- subset(total, Vis == "0seven" | Vis == "15" | Vis == "63")
confidenceNum <- subset(confidence, Vis == "0seven" | Vis == "15" | Vis == "63")

#merging the data file that includes the metrics on hurricane size and intensity with the main data files
totalNum<- merge(totalNum,hurricane,by=c("Hurricane","Time"))
totalAn<- merge(totalAn,hurricane,by=c("Hurricane","Time"))

#getting male and female count and age
#sex 1 = male, 2 = female 
table<- total %>% 
  group_by(Sex) %>% 
    summarise(length(unique(ResponseID)))

#age
mean(total$Age)
sd(total$Age)
```


<h1>Participants </h1>
7 tracks: 50<br>
15 tracks: 51<br>
63 tracks: 50<br>
Full Vis: 50<br>

Males: 49<br>
Female: 152<br>
Age: Mean = 23.53, SD = 6.6<br>

<hr />

<h1>Results: H1 </h1>
VisNumD0 = 7 tracks (referent)<br>
VisNumD1 = 15 tracks<br>
VisNumD2 = 63 tracks<br>


```{r damage number plot, echo=FALSE, message=FALSE}

#making the plot theme the correct size
theme_set(theme_bw(base_size = 13))


#coding Vis Type for easy interpretation
totalNum$Visualization <- ifelse(totalNum$Vis == "0seven", 7,
                                 ifelse(totalNum$Vis == "15", 15,
                                         ifelse(totalNum$Vis == "63", 63, 999)))

#coding Vis Type for easy interpretation
total$Visualization <- ifelse(total$Vis == "0seven", 7,
                                 ifelse(total$Vis == "15", "15 Unannotated",
                                         ifelse(total$Vis == "63", 63, 
                                                ifelse(total$Vis == "Full", "15 Annotated", 999))))

#changing the order of the factors
total$Visualization <- factor(total$Visualization, levels =c("7","15 Unannotated", "63", "15 Annotated"))

#getting error bars with the Morey 2008 correction
SummaryDamage<- summarySEwithin(total, "value", betweenvars="Visualization", withinvars=c("Distance", "Time"), idvar="ResponseID")

#plot
p<- SummaryDamage %>%
  ggplot(aes(x=Distance, y=value_norm, group=interaction(Visualization, Time)))+
  geom_errorbar(aes(ymin=value_norm-ci, ymax=value_norm+ci), width=0.1)+
  geom_line(aes(color=Time), size = .8)+
  scale_colour_manual(values = c("#d3661c", "#3c8dba"))+
  geom_point(aes(),size = .5) +
  xlab("\nDistance") +
  ylab("Damage\n")+
  scale_y_continuous(breaks=c(1, 2, 3, 4, 5, 6, 7), limits = c(1.1,6.9))

#plotting 7, 15, 63 next to each other   
p <- p + facet_wrap( ~ Visualization, nrow = 1)+
  theme(legend.position="bottom", strip.background = element_rect(colour="#FFFFFF", fill="#FFFFFF"))+
  coord_fixed(ratio = 1/1)
p

ggsave("test.pdf",width = 12, height = 5, p)

```


```{r H1 damage, echo=FALSE, message=FALSE}

#making Vis Type and Time into factors for analysis. 
totalNum$Vis<-as.factor(totalNum$Vis)
totalNum$Time<-as.factor(totalNum$Time)


#multilevel model equation 
m.sem = lmer(value ~ Distance*Time*Vis+ (1|ResponseID), data=totalNum, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
#use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(12, bold = T, color = "white", background = "#d3661c")
```

Confidence intervals using Wald method
```{r damage CI, echo=FALSE, message=FALSE}
#confidence intervals 
fm1W <- confint(m.sem, method="Wald")
fm1W
```

Effect size
```{r damage r2, echo=FALSE, message=FALSE}
#effect size for the model
r.squaredGLMM(m.sem)

```

The results of this analysis revealed a significant three-way interaction between distance, time point, and the 7 vs. 63 track display, b = 0.589, t = 5.814, p = 0.000, 95 % CI [0.39, 0.78] (highlighted in orange). To break down the three-way interaction, the same equation as above was computed but on each of the visualization types separately. 

<hr/>

<h2>Post-Hoc: 63-track</h2> 

```{r 63 model, echo=FALSE, message=FALSE}

#subsetting the data to evaluate each visualization type
seven<- subset(totalNum, Vis=="0seven")
fifteen<- subset(totalNum, Vis=="15")
sixtythree<- subset(totalNum, Vis=="63")


#multilevel model equation for the 48hr timepoint
m.sem = lmer(value ~ Distance*Time + (1|ResponseID), data=sixtythree, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(4, bold = T, color = "white", background = "#d3661c")
```
  
Confidence intervals using Wald method
```{r 63 track CI, echo=FALSE, message=FALSE}
#confidence intervals 
fm1W <- confint(m.sem, method="Wald")
fm1W
```

Effect size
```{r 63 track r2, echo=FALSE, message=FALSE}
#effect size for the model
r.squaredGLMM(m.sem)

```

There is a significant between interaction time point and distance (highlighted in orange). To break down the interaction, the same equation as above was computed but on each time point separately.  

<hr/>
<h3>Post-Hoc: 63 track, 24-hr</h3>  

```{r 63 model 24 hr, echo=FALSE, message=FALSE}

sixtythree24<- subset(sixtythree, Time=="24")

#multilevel model equation for the 48hr timepoint
m.sem = lmer(value ~ Distance + (1|ResponseID), data=sixtythree24, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
  #row_spec(4, bold = T, color = "white", background = "#d3661c")

```
Confidence intervals using Wald method
```{r 63 24hr track CI, echo=FALSE, message=FALSE}
#confidence intervals 
fm1W <- confint(m.sem, method="Wald")
fm1W
```

Effect size
```{r 63 24hr track r2, echo=FALSE, message=FALSE}
#effect size for the model
r.squaredGLMM(m.sem)

```
The 63 track display at 24 hours has an overall slope of -2.45, meaning that for every one unit change in distance, damage ratings decrease by 2.45 on average on the Likert scale from 1-7, conditional R-squared = .55. 

<h3>Post-Hoc: 63 track, 48-hr </h3> 

```{r 63 model 48 hr, echo=FALSE, message=FALSE}

sixtythree48<- subset(sixtythree, Time=="48")

#multilevel model equation for the 48hr timepoint
m.sem = lmer(value ~ Distance + (1|ResponseID), data=sixtythree48, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
  #row_spec(4, bold = T, color = "white", background = "#d3661c")
```

Confidence intervals using Wald method
```{r 63 48hr track CI, echo=FALSE, message=FALSE}
#confidence intervals 
fm1W <- confint(m.sem, method="Wald")
fm1W
```

Effect size
```{r 63 48hr track r2, echo=FALSE, message=FALSE}
#effect size for the model
r.squaredGLMM(m.sem)

```
Whereas, for the 63 track display at 48 hours, the average slope is -1.77. The original time point and distance interaction (highlighted in orange) indicates that there is a significant difference between the 24- and 48-hr slopes of .68  (b = 0.67, t = 9.18, p ¡ 0.000, 95 % CI [0.53, 0.82], conditional R-squared for the model = 0.527).

<hr/>
<h2>Post-Hoc: 15 track </h2>

```{r 15 model, echo=FALSE, message=FALSE}

#multilevel model equation for the 48hr timepoint
  m.sem = lmer(value ~ Distance*Time + (1|ResponseID), data=fifteen, REML=FALSE)
  
  coefs <- data.frame(coef(summary(m.sem)))
  # use normal distribution to approximate p-value
  coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
  knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
    row_spec(4, bold = T, color = "white", background = "#d3661c")
```
  
Confidence intervals using Wald method
```{r 15 track CI, echo=FALSE, message=FALSE}
#confidence intervals 
fm1W <- confint(m.sem, method="Wald")
fm1W
```

Effect size
```{r 15 track r2, echo=FALSE, message=FALSE}
#effect size for the model
r.squaredGLMM(m.sem)

```

  
```{r 15 slopes for 24 and 48, include=FALSE}
####### getting the slopes for the 7-track 24hr and 48hr
fifteen24<- subset(fifteen, Time=="24")
fifteen48<- subset(fifteen, Time=="48")

#24
m.sem = lmer(value ~ Distance+ (1|ResponseID), data=fifteen24, REML=FALSE)
coefs <- data.frame(coef(summary(m.sem)))
coefs

#48
m.sem = lmer(value ~ Distance+ (1|ResponseID), data=fifteen48, REML=FALSE)
coefs <- data.frame(coef(summary(m.sem)))
coefs

```

There is not a significant between interaction time point and distance (highlighted in orange).

<hr/>
<h2>Post-Hoc: 7 track</h2> 
```{r 7 track model, echo=FALSE, message=FALSE}
#multilevel model equation for the 7 track display
m.sem = lmer(value ~ Distance*Time + (1|ResponseID), data=seven, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(4, bold = T, color = "white", background = "#d3661c")

```

Confidence intervals using Wald method
```{r 7 track CI, echo=FALSE, message=FALSE}
#confidence intervals 
fm1W <- confint(m.sem, method="Wald")
fm1W
```

Effect size
```{r 7 track r2, echo=FALSE, message=FALSE}
#effect size for the model
r.squaredGLMM(m.sem)

```


```{r 7 slopes for 24 and 48, include=FALSE}
# getting the slopes for the 7-track 24hr and 48hr
seven24<- subset(seven, Time=="24")
seven48<- subset(seven, Time=="48")

#24
m.sem = lmer(value ~ Distance+ (1|ResponseID), data=seven24, REML=FALSE)
coefs <- data.frame(coef(summary(m.sem)))
coefs

#48
m.sem = lmer(value ~ Distance+ (1|ResponseID), data=seven48, REML=FALSE)
coefs <- data.frame(coef(summary(m.sem)))
coefs
```

There is not a significant between interaction time point and distance (highlighted in orange). <br/>

In sum, the three-way interaction indicates an interaction between time point and distance for the 63 track display but not the 7 track display, which supports H1.   



<hr/>

<h1>Results: H2</h1> 
VisFull: 15 unannotated is the referant
```{r damage annotation, echo=FALSE }

#making Visualization a factor
totalAn$Vis<-as.factor(totalAn$Vis)


#multilevel model equation comparing 15 annotated to 15 unannotated 
m.sem = lmer(value ~ Distance*TimeD*Vis +  (1|ResponseID), data=totalAn, REML=FALSE)
coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(4, bold = T, color = "white", background = "#3c8dba")%>%
  row_spec(6, bold = T, color = "white", background = "#d3661c")

```

Confidence intervals using Wald method
```{r h2 CI, echo=FALSE, message=FALSE}
#confidence intervals 
fm1W <- confint(m.sem, method="Wald")
fm1W
```

Effect size
```{r h2 r2, echo=FALSE, message=FALSE}
#effect size for the model
r.squaredGLMM(m.sem)
```

There is a significant interaction between distance and VisFull (15 annotated vs. 15 unannotated) (highlighted in orange). To break this down, we recenter distance around the farthest time point (1.78).

<hr/>
<h2> Post-Hoc: re-centering</h2> 
```{r damage annotation recenter, echo=FALSE }

#recentering 
totalAn$DistanceCenter1.78 <- totalAn$Distance - 1.78

#multilevel model equation comparing 15 annotated to 15 unannotated with  distance recentered
m.sem = lmer(value ~ DistanceCenter1.78*Vis +  (1|ResponseID), data=totalAn, REML=FALSE)
coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(3, bold = T, color = "white", background = "#3c8dba")
 
#making factors for the plots
totalAn$Time<-as.factor(totalAn$Time)
totalAn$Visualization<-totalAn$Vis
totalAn$Visualization<-as.factor(totalAn$Visualization)


#getting error bars with the Morey 2008 correction
SummaryDamage<- summarySEwithin(totalAn, "value", betweenvars="Vis", withinvars=c("Distance", "Time"), idvar="ResponseID")

#plot
p<- SummaryDamage %>%
  ggplot(aes(x=Distance, y=value_norm, group=interaction(Vis, Time)))+
  geom_errorbar(aes(ymin=value_norm-ci, ymax=value_norm+ci), width=0.1)+
  geom_line(aes(color=Time), size = .8)+
  scale_colour_manual(values = c("#d3661c", "#3c8dba"))+
  geom_point() +
  xlab("\nDistance") +
  ylab("Damage\n")+
  scale_y_continuous(breaks=c(1, 2, 3, 4, 5, 6, 7))

#plotting 7, 15, 63 next to each other   
p <- p + facet_wrap( ~ Vis)+
  theme(legend.position="bottom", strip.background = element_rect(colour="#FFFFFF", fill="#FFFFFF"))
p

```
<br/>
This analysis revealed that at the farthest distance there is not a significant difference between the two Vis types (highlighted in blue). Overall, the 15 track elicited greater damage ratings but this difference diminished at the farthest time point. <br/><br/>

<hr/>
<h3>Results: Size and Intensity for 15 annotated </h3>

```{r 15 track size and intensity breakdown}

#subsetting to look at just 15 annotated 
totalAnFull <- subset(totalAn, Vis == "Full")

#multilevel model equation evaluating the impact of size and intensity on the 15 track annotated display. 
full.model = lmer(value ~ Distance*Time + SizeMean + IntenistyMean +  (1|ResponseID), data=totalAn, REML=FALSE)
coefs <- data.frame(coef(summary(full.model)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(4, bold = T, color = "white", background = "#3c8dba")%>%
  row_spec(5, bold = T, color = "white", background = "#3c8dba")
```

Confidence intervals using Wald method
```{r size and intensity CI, echo=FALSE, message=FALSE}
#confidence intervals 
fm1W <- confint(m.sem, method="Wald")
fm1W
```

Effect size
```{r size and intensity r2, echo=FALSE, message=FALSE}
#effect size for the model
r.squaredGLMM(m.sem)
```


As predicted, both size and intensity account for a significant proportion of variance in damage ratings, such that as both size and intensity increase damage also increases. 

<hr/>
<h3>15 annotated</h3>
```{r size and intensity for 15 annotated, echo=FALSE }

#getting average damage ratings per person at a given size for 15 annotated
totalAnSum<- totalAn%>%
  filter(Annotation=="1")%>%
    dplyr::group_by(ResponseID, SizeMean)%>%
      dplyr::summarise(value = mean(value))



#size plot for 15 annotated 
sizeplot.pub <-  totalAnSum %>%
  ggplot(aes(x=SizeMean, y=value)) + 
  geom_smooth(method="lm", se=TRUE, size=2) +
  geom_jitter(width = .9, alpha=.3)+
  xlab("\nMean Size (nmi)") +
  ylab("Damage\n") + 
  theme(legend.position="bottom") + 
  scale_y_continuous(breaks=c(1, 2, 3, 4, 5, 6, 7)) +
  expand_limits(y=1)+expand_limits(y=7) +
    theme(panel.border = element_blank()) +
    theme(axis.line = element_line(color = 'black'))
sizeplot.pub




#getting average damage ratings per person at a given intenisty  15 annotated
totalAnSumIn<- totalAn%>%
  filter(Annotation=="1")%>%
    dplyr::group_by(ResponseID, IntenistyMean)%>%
     dplyr::summarise(value = mean(value))

#Intenisty plot for 15 annotated 
sizeplot.pub <-  totalAnSumIn %>%
  ggplot(aes(x=IntenistyMean, y=value)) + 
  geom_smooth(method="lm", se=TRUE, size=2) +
  geom_jitter(width = 2.3, alpha=.3)+
  xlab("\nMean Intensity (kt)") +
  ylab("Damage\n") + 
  theme(legend.position="bottom") + 
  scale_x_continuous(breaks=c(50, 60, 70, 80, 90, 100, 110)) +
  scale_y_continuous(breaks=c(1, 2, 3, 4, 5, 6, 7)) +
  expand_limits(y=1)+expand_limits(y=7) +
    theme(panel.border = element_blank()) +
    theme(axis.line = element_line(color = 'black'))
sizeplot.pub
```

<h3>15 unannotated</h3>
```{r size and intensity for 15 unannotated, echo=FALSE }
#getting average damage ratings per person at a given size for 15 unannotated
FifteenUnSize<- totalNum%>%
  filter(Vis=="15")%>%
    dplyr::group_by(ResponseID, SizeMean)%>%
      dplyr::summarise(value = mean(value))

#getting average damage ratings per person at a given size for 15 annotated
totalUnSum<- totalAn%>%
  filter(Annotation=="0")%>%
    dplyr::group_by(ResponseID, SizeMean)%>%
      dplyr::summarise(value = mean(value))


#Size plot for 15 unannotated 
sizeplot.pub <-  FifteenUnSize %>%
  ggplot(aes(x=SizeMean, y=value)) + 
  geom_smooth(method="lm", se=TRUE, size=2) +
  geom_jitter(width = .9, alpha=.3)+
  xlab("\nMean Size (nmi)") +
  ylab("Damage\n") + 
  theme(legend.position="bottom") + 
  scale_y_continuous(breaks=c(1, 2, 3, 4, 5, 6, 7)) +
  expand_limits(y=1)+expand_limits(y=7) +
    theme(panel.border = element_blank()) +
    theme(axis.line = element_line(color = 'black'))
sizeplot.pub


#getting average damage ratings per person at a given intenisty  15 unannotated
FifteenUnIntenisty<- totalNum%>%
  filter(Vis=="15")%>%
  dplyr::group_by(ResponseID, IntenistyMean)%>%
    dplyr::summarise(value = mean(value))


#Intenisty plot for 15 unannotated 
sizeplot.pub <-  FifteenUnIntenisty %>%
  ggplot(aes(x=IntenistyMean, y=value)) + 
  geom_smooth(method="lm", se=TRUE, size=2) +
  geom_jitter(width = .9, alpha=.3)+
  xlab("\nMean Size (nmi)") +
  ylab("Damage\n") + 
  theme(legend.position="bottom") + 
  scale_y_continuous(breaks=c(1, 2, 3, 4, 5, 6, 7)) +
  expand_limits(y=1)+expand_limits(y=7) +
    theme(panel.border = element_blank()) +
    theme(axis.line = element_line(color = 'black'))
sizeplot.pub

```


<hr/>

<h2>Results: 15 annotated vs. 63</h2>  

VisFull: 63 is the referant
```{r damage annotation vs 63, echo=FALSE }

#subsetting to look at 63 vs. 15 annotated 
totalAnvs63 <- subset(total, Vis == "Full" | Vis == "63")

#multilevel model equation 
m.sem = lmer(value ~ Distance*TimeD*Vis + (1|ResponseID), data=totalAnvs63, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(8, bold = T, color = "white", background = "#d3661c")

#confidence intervals 
fm1W <- confint(m.sem, method="Wald")
fm1W

#effect size for the model
r.squaredGLMM(full.model)


```

There is a significant 3-way interaction between distance, timepoint, and vis-type (highlighted orange). Following the same procedures as before, we ran the same model on the 24hr and 48hr time points, separately.

<h3>Post-Hoc: 15 annotated </h3>
```{r 24hr annotation vs 63, echo=FALSE, message=FALSE}

#subsetting the data to look at the 24 and 48 hours timepoints separately
totalAnvs63_24<- subset(totalAnvs63, Time=="24")
totalAnvs63_48<- subset(totalAnvs63, Time=="48")

#multilevel model equation for 24hr timepoint
m.sem = lmer(value ~ Distance*Vis + (1|ResponseID), data=totalAnvs63_24, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(3, bold = T, color = "white", background = "#3c8dba")%>%
  row_spec(4, bold = T, color = "white", background = "#d3661c")



```

<h3>Post-Hoc: 15 annotated </h3>
```{r 15 annotation a, echo=FALSE, message=FALSE}


#multilevel model equation for 24hr timepoint
m.sem = lmer(value ~ Distance*Time + (1|ResponseID), data=totalAnFull, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
  #row_spec(3, bold = T, color = "white", background = "#3c8dba")%>%
  #row_spec(4, bold = T, color = "white", background = "#d3661c")


```

<h3>Post-Hoc: 63 track </h3>
```{r 15 annotation b, echo=FALSE, message=FALSE}


#multilevel model equation for 24hr timepoint
m.sem = lmer(value ~ Distance*Time + (1|ResponseID), data=sixtythree, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
  #row_spec(3, bold = T, color = "white", background = "#3c8dba")%>%
  #row_spec(4, bold = T, color = "white", background = "#d3661c")

```

The analysis of the 24hr time point revealed a significant 2-way interaction between distance and vis-type (highlighted in orange). To break this down, we recentered distance around the farthest time point (1.78). 

<h4>Post-Hoc: 24hr re-centering</h4>

```{r 24hr 63 vs full re-center, echo=FALSE, message=FALSE}

#recentering
totalAnvs63_24$DistanceCenter1.78 <- totalAnvs63_24$Distance - 1.78

#multilevel model equation for 24hr timepoint and recentered distance 
m.sem = lmer(value ~ DistanceCenter1.78*Vis + (1|ResponseID), data=totalAnvs63_24, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(3, bold = T, color = "white", background = "#3c8dba")


```

This analysis revealed that there was still an effect of vis-type (highlighted in blue) but in the opposite direction.  Meaning that, at the center of the storm, participants believed that the 15 track annotated display would receive less damage than the 63 track  display. But at the farthest time point (1.78) participants rated the 15 track annotated display as receiving significantly more damage.  

<h3>Post-Hoc: 48hr </h3>

```{r 48hr 63 vs full, echo=FALSE, message=FALSE}

#multilevel model equation for 24hr timepoint and recentered distance 
m.sem = lmer(value ~ Distance*Vis + (1|ResponseID), data=totalAnvs63_48, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(3, bold = T, color = "white", background = "#3c8dba")%>%
  row_spec(4, bold = T, color = "white", background = "#d3661c")


```

Similar to the 24hr timepoint, there was a significant interaction between distance and vis-type (orange).

<h4>Post-Hoc: 48hr re-centering</h4>

```{r 48hr 63 vs full re-center, echo=FALSE, message=FALSE}

#recentering 
totalAnvs63_48$DistanceCenter1.78 <- totalAnvs63_48$Distance - 1.78

#multilevel model equation for 48hr timepoint and recentered distance 
m.sem = lmer(value ~ DistanceCenter1.78*Vis + (1|ResponseID), data=totalAnvs63_48, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(3, bold = T, color = "white", background = "#3c8dba")


```

When we recentered the 48hr time point around the farthest probe location, there was no longer a significant difference between the 15 track annotated display and the 63 track display (also highlighted in blue). 


<hr/>

<h1>Confidence</h1>
```{r Confidence plot, echo=FALSE, message=FALSE}

#making the plot theme the correct size
theme_set(theme_bw(base_size = 13))

#coding Vis Type for easy interpretation
confidence$Visualization <- ifelse(confidence$Vis == "0seven", 7,
                                 ifelse(confidence$Vis == "15", "15 Unannotated",
                                         ifelse(confidence$Vis == "63", 63, 
                                                ifelse(confidence$Vis == "Full", "15 Annotated", 999))))

#changing the order of the factors
confidence$Visualization <- factor(confidence$Visualization, levels =c("7","15 Unannotated", "63", "15 Annotated"))

#getting error bars with the Morey 2008 correction
SummaryConfidence<- summarySEwithin(confidence, "Confidence", betweenvars="Visualization", withinvars=c("Distance", "Time"), idvar="ResponseID")

#plot
p<- SummaryConfidence %>%
  ggplot(aes(x=Distance, y=Confidence_norm, group=interaction(Visualization, Time)))+
  geom_errorbar(aes(ymin=Confidence_norm-ci, ymax=Confidence_norm+ci), width=0.1)+
  geom_line(aes(color=Time), size = .8)+
  scale_colour_manual(values = c("#d3661c", "#3c8dba"))+
  geom_point(aes(),size = .5) +
  xlab("\nDistance") +
  ylab("Confidence\n")+
  scale_y_continuous(breaks=c(1, 2, 3, 4, 5, 6, 7), limits = c(1.1,6.9))

#plotting 7, 15, 63 next to each other   
p <- p + facet_wrap( ~ Visualization, nrow = 1)+
  theme(legend.position="bottom", strip.background = element_rect(colour="#FFFFFF", fill="#FFFFFF"))+
  coord_fixed(ratio = 1/1)
p


```

<h2>Confidence: 7, 15, 63 </h2>
Vis: 7 is the referant
```{r confidnec number, echo=FALSE}

#making Vis type a factor 
confidenceNum$Vis<-as.factor(confidenceNum$Vis)

#multilevel model equation 
m.sem = lmer(Confidence ~ Distance+TimeD+Vis + (1|ResponseID), data=confidenceNum, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

```

```{r confidence number plot, echo=FALSE}
#making the plot theme the correct size
theme_set(theme_bw(base_size = 24))

#coding Vis Type for easy interpretation 
confidenceNum$Visualization <- ifelse(confidenceNum$Vis == "0seven", 7, ifelse(confidenceNum$Vis == "15", 15, ifelse(confidenceNum$Vis == "63", 63, 9)))

#making Visualization and Time a factor
confidenceNum$Visualization<-as.factor(confidenceNum$Visualization)
confidenceNum$Time<-as.factor(confidenceNum$Time)


```


<hr/>
<h2>Confidence: 15 annotated vs. 15 unannotated </h2>

```{r confience annotation, echo=FALSE}

#multilevel model equation 
m.sem = lmer(Confidence ~ Distance+ Time+ Vis + (1|ResponseID), data=confidenceAn, REML=FALSE)

coefs <- data.frame(coef(summary(m.sem)))
# use normal distribution to approximate p-value
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

```
<hr/>

<h1>Follow-up Questions </h1>
<em>1. We are very interested in how you made your decisions about the level of damage that the oil platform would incur. Please describe in as much detail as possible your decision-making process. This will greatly help us learn about how people make decisions with hurricane visualizations.</em>
<hr/>

<em> 2. Do you have experience with hurricane forecasts? </em><br/>

```{r post-exp quesitons Forecasts,  echo=FALSE, }
#renaming
renameing<- total
renameing$Forecasts <- ifelse(renameing$Forecasts == "1", "Yes", ifelse(renameing$Forecasts == "2", "No", 9))

#1=yes, 2=no 
table<- renameing %>% 
  group_by(Forecasts) %>% 
    summarise(length(unique(ResponseID)))

#print table 
knitr::kable(table, "html")%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")


```
<hr/>
<em> 3. Have you lived or do you live in an area that experiences hurricane threats?</em><br/>

```{r post-exp quesitons Threats,  echo=FALSE, }

#1=yes, 2=no <br/>

#renaming
renameing$Threats <- ifelse(renameing$Threats == "1", "Yes", ifelse(renameing$Threats == "2", "No", 9))

#getting counts
table<- renameing %>% 
  group_by(Threats) %>% 
    summarise(length(unique(ResponseID)))

#print table 
knitr::kable(table, "html")%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")


```
<hr/>
<em> 4. The display shows the hurricane getting large over time.</em> <br/>


```{r post-exp quesitons Larger,  echo=FALSE,fig.width=4.5, fig.height=4 }
#3=False, 1=True<br/>

#renaming
renameing$Vis <- ifelse(renameing$Vis == "0seven", "7", renameing$Vis)
renameing$Larger <- ifelse(renameing$Larger == "3", "False", ifelse(renameing$Larger == "1", "True", 9))

#getting counts 
table<- renameing %>% 
  group_by(Vis, Larger) %>% 
    summarise(length(unique(ResponseID)))

#print table
knitr::kable(table, "html")%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "float_left")

#plotting data
ggplot(data=renameing, aes(Vis))+
  ylab("percent of responses")+
  geom_bar(aes(fill=as.factor(Larger)), position="fill")+
  scale_fill_discrete(name = "")+
  theme(text = element_text(size=12))


              
```

Doing a follow-up regression analysis using the 15 track annotated as the referant, we find that participants viewing the 15 track annotated visualization were more likely to believe that the display showed the hurricane getting larger over time than those viewing the 15 track unannotated display. 


```{r post-exp quesitons Larger reggresion,  echo=FALSE,fig.width=4.5, fig.height=4 }

#getting counts
LargeComp<- total %>%
    group_by(ResponseID, Vis) %>%
    summarize(mean(Larger))

#renaming
LargeComp$Vis <- ifelse(LargeComp$Vis == "0seven", "7",
                                  ifelse(LargeComp$Vis == "Full", "0Full", LargeComp$Vis))

#renaming
LargeComp$Larger <- as.factor(ifelse(LargeComp$`mean(Larger)` == "3", "False", ifelse(LargeComp$`mean(Larger)` == "1", "True", 9)))

#regression analysis
fit<- glm(Larger~Vis, data=LargeComp, family="binomial")

coefs <- data.frame(coef(summary(fit)))
knitr::kable(coefs, "html")%>%   kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

```


<hr/>
<em> 5. The display indicates that the forecasters are less certain about the path of the hurricane as time passes.</em> <br/>


```{r post-exp quesitons Variable,  echo=FALSE, fig.width=4.5, fig.height=4}
#3 = False, 2 = True<br/>

#renaming
renameing$Variable <- ifelse(renameing$Variable == "3", "False", ifelse(renameing$Variable == "2", "True", 9))

#getting counts
table<- renameing %>% 
  group_by(Vis, Variable) %>% 
    summarise(length(unique(ResponseID)))

#print table
knitr::kable(table, "html")%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "float_left")

#make plot
ggplot(data=renameing, aes(Vis))+
  ylab("percent of responses")+
  geom_bar(aes(fill=as.factor(Variable)), position="fill")+
  scale_fill_discrete(name = "")+
  theme(text = element_text(size=12))
  

```
<hr/>
<em> 6. Areas on the map not covered by the visualization will not be hit by the hurricane.</em><br/>


```{r post-exp quesitons Hit,  echo=FALSE, fig.width=4.5, fig.height=4 }
#True=1 , False=2

#renaming
renameing$Hit <- ifelse(renameing$Hit == "1", "True", ifelse(renameing$Hit == "2", "False", 9))

#get counts
table<- renameing %>% 
  group_by(Vis, Hit) %>% 
    summarise(length(unique(ResponseID)))

#print table 
knitr::kable(table, "html")%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "float_left")

#plot data
ggplot(data=renameing, aes(Vis))+
  ylab("percent of responses")+
  geom_bar(aes(fill=as.factor(Hit)), position="fill")+
  scale_fill_discrete(name = "")+
  theme(text = element_text(size=12))

```
<hr/>
<em>7. Did you recognize the areas in the map? </em><br/>

```{r post-exp quesitons Q277,  echo=FALSE, }
#1=yes, 2=no

#renaming
renameing$Recognize<-renameing$Q277
renameing$Recognize <- ifelse(renameing$Recognize == "1", "Yes", ifelse(renameing$Recognize == "2", "No", 9))

#getting counts
table<- renameing %>% 
  group_by(Recognize) %>% 
    summarise(length(unique(ResponseID)))

#print table
knitr::kable(table, "html")%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")


```