---
title: |
  | Data cleaning and analysis for the Master's thesis:
  | DIFFERENCES IN CONSUMER PREFERENCES FOR UNWEATHERED AND WEATHERED WOOD
  | submitted at the University of Primorska
  | Faculty of Mathematics, Natural Sciences, and Information Technologies
author:
  - \*
  - Hana Remešovà, student, UP FAMNIT
  - Michael Burnard, mentor, Assist. prof. University of Primorska
  - \*
date: "11 August 2020"
output:
  html_document: 
    toc: true
    toc_depth: 2
  pdf_document:
    toc: true
    toc_depth: 2
---


# Introduction

This file contains the processing steps associated with a master's student thesis submitted to UP FAMNIT.

The purpose of the thesis was to assess consumer preferences for wood facade materials in their fresh state and after they have weathered for some time.

Interviews were conducted to gather consumer preferences and material properties were measured in the lab.

The necessary data inputs are uploaded along with this file, with the exception of the individual files containing colour measurements that were combined into a single file for further analysis. The combined file is uploaded, instead.

See the masters thesis for more information about the project, the methods, and the data description. Hypotheses are paraphrased in this document.

```{r setup, include=FALSE, message=FALSE, warning=FALSE}
knitr::opts_chunk$set(echo = TRUE, 
                      warning = TRUE, 
                      fig.width=7.5, 
                      fig.height=4)
library(MASS)
library(tidyverse)
library(stringr)
library(fs)
library(kableExtra)

Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

fix.names <- function(colname)
{
  return(str_replace(colname, " ", "."))
}
```

# Data

Import and clean survey data. Add supplemental information about samples (material properties).

```{r data_1, message=FALSE, warning=FALSE}
d <- read_csv("data/remesova_ms_data_v1.csv")
colnames(d)[3:13] <- paste("q1_", colnames(d)[3:13], sep="")
colnames(d)[14:19] <- paste("q2_", colnames(d)[14:19], sep="")
colnames(d)[20:25] <- paste("q3_", colnames(d)[20:25], sep="")
colnames(d)[26:31] <- paste("q4_", colnames(d)[26:31], sep="")
colnames(d)[32:37] <- paste("q5_", str_sub(colnames(d)[32:37], 0,2), sep="")
colnames(d)[38:43] <- paste("q6_", str_sub(colnames(d)[38:43], 0,2), sep="")
colnames(d)[44:48] <- paste("q7_", str_sub(colnames(d)[44:48], 0,2), sep="")
colnames(d)[53:56] <- paste("q10_", str_sub(colnames(d)[53:56], 0,2), sep="")
colnames(d)[49] <- "Gender"

d$Gender <- if_else(d$Gender == "babyYoda", "Not Specified", d$Gender)

d.long <- d %>% 
  pivot_longer(c(50:52), names_to="Region", values_to="Response") %>%
  filter(!is.na(Response)) %>%
  select(-Response)

d.long <- d.long %>% 
  pivot_longer(c(44:48), names_to="Age", values_to="AgeGroup") %>%
  filter(!is.na(AgeGroup)) %>%
  select(-AgeGroup)

d.long <- d.long %>% select(-c(12:13,45:48))



d.long <- d.long %>% pivot_longer(
  c(3:41), names_to="Question", values_to="Response") %>%
  separate(Question, into=c("Question", "SubQuestion"), remove=FALSE)

age <- tibble(col_name = c("q7_a_", "q7_b_", "q7_c_", "q7_d_", "q7_e_"),
              age_group = c("18-25", "26-37", "38-50", "51-64", "65+"))

key <- tibble(code = c("A", "B", "C", "D", "E", "H",
                       "I", "L", "F", "K", "G", "J",
                       "AI", "BL", "CF", "DK", "EG", "HJ"),
              names = rep(c("thermally modified pine",
                            "thermally m. oak + vacuum + wax",
                            "untreated oak",
                            "acetylated pine + coating",
                            "thermally modified oak + vacuum",
                            "furfurylated pine"), times=3))

d.q1 <- d.long %>% filter(Question == "q1")


d.long <- d.long %>% 
  filter(Question != "q1")

d.long <- left_join(d.long, key, by=c("SubQuestion" = "code"))
d.long <- left_join(d.long, age, by=c("Age" = "col_name")) %>%
  mutate(names = fct_relevel(names, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))

colour <- read_csv("data/remesova_colour_summary_v1.csv")
colour <- left_join(colour, key, by=c("Sample" = "code")) %>% 
  mutate(state = rep(c("Fresh", "Weathered"), each=6))
colour_change <- read_csv("data/remesova_colour_change_v1.csv")
colour_change <- left_join(colour_change, key, by = c("Sample" = "code"))
```

## Color data 

Import and merge color data.

** NOT RUN Use combined file: remesova_combined_colour_raw_v1.csv **

```{r data_colour, eval=FALSE}

files <- dir_ls("data/color/", glob="*.txt")

quick_parse <- function(file) {
  tmp <- read_tsv(file, locale=locale(decimal_mark = ","),
                  col_names=c("L", "a", "b")) %>%
    mutate(blrp = str_sub(file, -22,-15)) %>%
    separate(blrp, c("id", "stage"))
}

d.clr <- map_dfr(files, quick_parse)

#write the color data in one file
#write_csv(d.clr, "data/remesova_combined_colour_raw_v1.csv")
```


# Question 1

Overview figure.

```{r question1, fig.width=7.5, fig.height=4}

ggplot(data=d.q1, aes(y=Response, x=SubQuestion)) + theme_bw() +
  geom_col()
```

# Question 2

Overview figures and summary table.

```{r question2, fig.width=7.5, fig.height=4}

d.q2 <- d.long %>% filter(Question == "q2")

ggplot(data=d.q2, aes(y=Response, x=names)) + theme_bw() +
  geom_boxplot() + 
  scale_y_continuous(breaks=seq(1,7,1)) +
  theme(axis.text.x = element_text(angle=70, vjust=0.5))

d.q2 %>% group_by(SubQuestion) %>%
  summarise(mean = mean(Response, na.rm=TRUE),
                   median= median(Response, na.rm=TRUE),
                   mode = Mode(Response),
                   sd = sd(Response, na.rm=TRUE))
```

# Question 3

Overview figures and summary table.

```{r question3, }

d.q3 <- d.long %>% filter(Question == "q3")

ggplot(data=d.q3, aes(y=Response, x=names)) + theme_bw() +
  geom_boxplot() + 
  scale_y_continuous(breaks=seq(1,7,1)) +
  theme(axis.text.x = element_text(angle=70, vjust=0.5))

d.q3 %>% group_by(names) %>%
  summarise(mean = mean(Response, na.rm=TRUE),
                   median= median(Response, na.rm=TRUE),
                   mode = Mode(Response),
                   sd = sd(Response, na.rm=TRUE))
```

# Question 4

Overview figures and summary table, basi

```{r question4}

d.q4 <- d.long %>% filter(Question == "q4") %>%
  mutate(Resp.ord = factor(Response, levels=c(1,2,3,4,5,6,7), ordered=TRUE))

ggplot(data=d.q4, aes(y=Response, x=names)) + theme_bw() +
  geom_boxplot() + 
  scale_y_continuous(breaks=seq(1,7,1)) +
  theme(axis.text.x = element_text(angle=70, vjust=0.5))

d.q4.sum <- d.q4 %>% group_by(names) %>%
  summarise(mean = mean(Response, na.rm=TRUE),
                   median= median(Response, na.rm=TRUE),
                   sd = sd(Response, na.rm=TRUE),
                   mode=Mode(Response))
d.q4.sum

```


# Question 5

```{r question5, fig.width=7.5, height=4}

d.q5 <- d.long %>% filter(Question == "q5")

ggplot(data=d.q5, aes(y=Response, x=names)) + theme_bw() +
  geom_boxplot() + 
  scale_y_continuous(breaks=seq(1,7,1)) +
  theme(axis.text.x = element_text(angle=70, vjust=0.5))

d.q5.sum <- d.q5 %>% group_by(names) %>%
  summarise(mean = mean(Response, na.rm=TRUE),
                   median= median(Response, na.rm=TRUE),
                   sd = sd(Response, na.rm=TRUE)) %>%
  mutate(rank = rank(mean))

d.q5q4 <- d.q4.sum %>% mutate(rankq5 = d.q5.sum$rank)
d.q5q4 
ggplot(data=d.q5q4, aes(y=rankq5, x=mean, xmin=mean-sd, xmax=mean+sd,
                        colour=names, shape=names)) + theme_bw() +
  geom_point(size=4) + 
  geom_errorbarh(height=0.2) +
  scale_y_continuous(limits=c(0.5,6.5), breaks=seq(1:6), expand=c(0,0)) +
  scale_x_continuous(limits=c(0.5,7.5), breaks=seq(1:7), expand=c(0,0)) +
  scale_colour_brewer(palette="Dark2", name="Material") +
  scale_shape_manual(values=c(0:2,4:6), name="Material") +
  labs(y="Rank\n(Lower is better)",
       x="Mean rating\n(Higher is better)")

ggsave("figs/q5q4_rankvrating.png", device=png(), height=4, width=8, units="in", dpi=300)
```

# Question 6



```{r question6}

d.q6 <- d.long %>% filter(Question == "q6")

ggplot(data=d.q6, aes(y=Response, x=names)) + theme_bw() +
  geom_boxplot() + 
  scale_y_continuous(breaks=seq(1,7,1)) +
  theme(axis.text.x = element_text(angle=70, vjust=0.5))

d.q6 %>% group_by(names) %>%
  summarise(mean = mean(Response, na.rm=TRUE),
                   median= median(Response, na.rm=TRUE),
                   sd = sd(Response, na.rm=TRUE))
```

# Hypothesis 1

Ratings of fresh and weathered specimens will differ.

Compared ratings per sample type between unweathered and weathered using paired wilcoxon rank sum test. Ties are problematic.

```{r hyp1, message=FALSE, warning=FALSE}
d.h1 <- full_join(d.q2 %>% select(code, names, Response) %>%
                    rename(Fresh = Response),
                  d.q3 %>% select(code, names, Response) %>%
                    rename(Weathered=Response),
                  by = c("code"= "code", "names"="names")) %>%
  mutate(diff = Weathered - Fresh)

d.h1.sum <- d.h1 %>% group_by(names) %>%
  summarise(median.diff = median(diff),
            median.weath = median(Weathered),
            median.fresh = median(Fresh),
            mu.diff = mean(diff),
            mu.fresh = mean(Fresh),
            mu.weath = mean(Weathered))
d.h1.sum

tmp.o <- d.h1 %>% filter(names=="thermally modified pine")
hyp1.out <- broom::tidy(wilcox.test(tmp.o$Weathered, tmp.o$Fresh, paired=TRUE, conf.int=TRUE)) %>%
  mutate(Name = "thermally modified pine")

acp <- d.h1 %>% filter(names=="acetylated pine + coating")
hyp1.out <- bind_rows(hyp1.out, broom::tidy(wilcox.test(acp$Weathered, acp$Fresh, paired=TRUE, conf.int=TRUE)) %>%
  mutate(Name = "acetylated pine + coating"))

fur.p <- d.h1 %>% filter(names=="furfurylated pine")
hyp1.out <- bind_rows(hyp1.out, broom::tidy(wilcox.test(fur.p$Weathered, fur.p$Fresh, paired=TRUE, conf.int=TRUE)) %>%
  mutate(Name = "furfurylated pine"))

oak <- d.h1 %>% filter(names=="untreated oak")
hyp1.out <- bind_rows(hyp1.out, broom::tidy(wilcox.test(oak$Weathered, oak$Fresh, paired=TRUE, conf.int=TRUE)) %>%
  mutate(Name = "untreated oak"))

tmo.v <- d.h1 %>% filter(names=="thermally modified oak + vacuum")
hyp1.out <- bind_rows(hyp1.out, broom::tidy(wilcox.test(tmo.v$Weathered, tmo.v$Fresh, paired=TRUE, conf.int=TRUE)) %>%
  mutate(Name = "thermally modified oak + vacuum"))

tmo.w <- d.h1 %>% filter(names=="thermally m. oak + vacuum + wax")
hyp1.out <- bind_rows(hyp1.out,broom::tidy(wilcox.test(tmo.w$Weathered, tmo.w$Fresh, paired=TRUE, conf.int=TRUE)) %>%
  mutate(Name = "thermally m. oak + vacuum + wax"))

hyp1.out <- hyp1.out %>% mutate(
  p.adjusted = p.adjust(p.value, method="fdr")
)

#write_csv(hyp1.out, "data/Hyp1_out.csv")
```

# Hypothesis 2

knowing the weathered appearance of a material will change consumers preference for material selection

```{r hyp2, message=FALSE, warning=FALSE}
d.h2 <- full_join(d.q2 %>% select(code, names, Response) %>%
                    rename(Fresh = Response),
                  d.q4 %>% select(code, names, Response) %>%
                    rename(Both=Response),
                  by = c("code"= "code", "names"="names")) %>%
        mutate(diff = Both - Fresh)

d.h2 %>% group_by(names) %>%
  summarise(median.fresh = median(Fresh),
            median.both = median(Both),
            median.diff = median(diff))
tmp.o <- d.h2 %>% filter(names=="thermally modified pine")
hyp2.out <- broom::tidy(wilcox.test(tmp.o$Both, tmp.o$Fresh, paired=TRUE, conf.int=TRUE)) %>%
  mutate(Name = "thermally modified pine")

acp <- d.h2 %>% filter(names=="acetylated pine + coating")
hyp2.out <- bind_rows(hyp2.out, 
                      broom::tidy(wilcox.test(acp$Both, acp$Fresh, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "acetylated pine + coating"))

tmo.v <- d.h2 %>% filter(names=="thermally modified oak + vacuum")
hyp2.out <- bind_rows(hyp2.out, 
                      broom::tidy(wilcox.test(tmo.v$Both, tmo.v$Fresh, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "thermally modified oak + vacuum"))

oak <- d.h2 %>% filter(names=="untreated oak")
hyp2.out <- bind_rows(hyp2.out, 
                      broom::tidy(wilcox.test(oak$Both, oak$Fresh, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "untreated oak"))

tmo.w <- d.h2 %>% filter(names=="thermally m. oak + vacuum + wax")
hyp2.out <- bind_rows(hyp2.out, 
                      broom::tidy(wilcox.test(tmo.w$Both, tmo.w$Fresh, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "thermally m. oak + vacuum + wax"))

fur.p <- d.h2 %>% filter(names=="furfurylated pine")
hyp2.out <- bind_rows(hyp2.out, 
                      broom::tidy(wilcox.test(fur.p$Both, fur.p$Fresh, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "furfurylated pine"))

hyp2.out <- hyp2.out %>% mutate(
  p.adjusted = p.adjust(p.value, method="fdr")
)

#write_csv(hyp2.out, "data/Hyp2_out.csv")

```



# Hypothesis 3

Consumer preference for weathered and unweathered wooden facade materials does not differ based on age, gender, or region.

Use ordered logistic regression. Check assumptions of the model with brant.

```{r hyp3}
d.q4  <- left_join(d.q4, colour_change, by="names")
q4.or.full <- polr(Resp.ord ~ names + Gender + age_group + Region, 
              data=d.q4, Hess=TRUE)

summary(q4.or.full)

q4.or.red <- polr(Resp.ord ~ names, 
              data=d.q4, Hess=TRUE)

summary(q4.or.red)

anova(q4.or.red, q4.or.full)
brant::brant(q4.or.red) 
#omnibuss passes, some factors are problematic. 
#But in the end ok we go with the omnibus.

q4.or.st <- as_tibble(as.matrix(coef(summary(q4.or.red))), rownames=NA,
                      .name_repair = fix.names) %>% 
  mutate(Coefficient = rownames(.),
         p.value = pnorm(abs(t.value),lower.tail = FALSE)*2) %>%
  select(Coefficient, Value, Std..Error, t.value, p.value)


q4.or.st %>% kable(digits=3) %>% kable_styling()


pred.d <- data.frame(names = c("thermally modified pine",
                       "thermally m. oak + vacuum + wax",
                       "untreated oak",
                       "acetylated pine + coating",
                       "thermally modified oak + vacuum",
                       "furfurylated pine"))

hyp3.out <- as.data.frame(predict(q4.or.red, pred.d, type="p")) %>%
  mutate(Material=pred.d$names)

hyp3.out

h3.fd <- hyp3.out %>% mutate(
  Unlikely = `1` + `2` + `3`,
  Neutral = `4`,
  Likely = `5` + `6` + `7`) %>% select(8:11) %>%
  pivot_longer(c(2:4), names_to="Rating", values_to="Probability")
  

ggplot(data=h3.fd, aes(y=Probability, x=Material, fill=Rating)) + theme_bw() +
  geom_col() +
  scale_fill_brewer(palette = "Dark2") +
  coord_flip()

#ggsave("figs/hyp3_probabilitiybars.png", device=png(), height=9, width=16, units="in", dpi=300)
#write_csv(hyp3.out, "data/Hyp3_out.csv")
#write_csv(q4.or.st, "data/Hyp3_ordinalRegression.csv")
```

# Exploratory q6

How does other information (maintenance, price, lifespan) change ratings. 

Doesn't add anything new because of correlation between factors. Raw summary data is still interesting.
```{r q4-exp-ordinal}

d.q4a <- d.long %>% filter(Question == "q4") %>% select(-Age) %>%
  mutate(Resp.ord = factor(Response, levels=c(1,2,3,4,5,6,7), ordered=TRUE))

d.q4a <- full_join(d.q4a, colour_change %>% select(-Sample), by="names")


q4a.or.full  <- polr(Resp.ord ~ Region + Gender + country + age_group + Colour + Gloss + Roughness,
                data=d.q4a, Hess=TRUE)

q4a.or.red <- polr(Resp.ord ~ Colour + Gloss + Roughness,
              data=d.q4a, Hess=TRUE)

anova(q4a.or.full, q4a.or.red) # choose full - significant drop in deviance

summary(q4a.or.red)
q4a.pr <- profile(q4a.or.red)
plot(q4a.pr)
confint(q4a.pr)
pairs(q4a.pr)
summary(colour_change)
hist(colour_change$Colour)
q4a.pred <- tibble(Colour=c(2.5, 10, 20, 7.5,7.5,7.5), 
                   #Gloss=c(0, 0, 0, 0, 0, 0),
                   Gloss=c(-1, -1, -1, -1, -1, -1),
                   Roughness=c(5, 5, 5, 2.5, 10, 20))
q4a.pred <- bind_cols(q4a.pred, as.data.frame(predict(q4a.or.red, q4a.pred, type="p")))
q4a.or.st <- as_tibble(as.matrix(coef(summary(q4a.or.red))), rownames=NA,
                      .name_repair = fix.names) %>% 
  mutate(Coefficient = rownames(.),
         p.value = pnorm(abs(t.value),lower.tail = FALSE)*2) %>%
  select(Coefficient, Value, Std..Error, t.value, p.value)


q4a.or.st %>% 
  kable(digits=4, 
       caption="Colour, Gloss, & Roughness are changes") %>%
  kable_styling()

#write_csv(q4a.or.st, "data/q4_cgr_ord.csv")

```
## Difference between Q4 and Q6

Wilcoxon tests.
```{r q6-q4-exp, message=FALSE, warning=FALSE}

d.q6q4 <- full_join(d.q6 %>% select(code, names, Response), 
                    d.q4 %>% select(code, names, Response),
                    by = c("code"="code", "names"="names"))
names(d.q6q4)[c(3,4)] <- c("Q6", "Q4")

tmp.o <- d.q6q4 %>% filter(names=="thermally modified pine")
d.q6q4.out <- broom::tidy(wilcox.test(tmp.o$Q6, tmp.o$Q4, paired=TRUE, conf.int=TRUE)) %>%
  mutate(Name = "thermally modified pine")

acp <- d.q6q4 %>% filter(names=="acetylated pine + coating")
d.q6q4.out <- bind_rows(d.q6q4.out, 
                      broom::tidy(wilcox.test(acp$Q6, acp$Q4, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "acetylated pine + coating"))

tmo.v <- d.q6q4 %>% filter(names=="thermally modified oak + vacuum")
d.q6q4.out <- bind_rows(d.q6q4.out, 
                      broom::tidy(wilcox.test(tmo.v$Q6, tmo.v$Q4, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "thermally modified oak + vacuum"))

oak <- d.q6q4 %>% filter(names=="untreated oak")
d.q6q4.out <- bind_rows(d.q6q4.out, 
                      broom::tidy(wilcox.test(oak$Q6, oak$Q4, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "untreated oak"))

tmo.w <- d.q6q4 %>% filter(names=="thermally m. oak + vacuum + wax")
d.q6q4.out <- bind_rows(d.q6q4.out, 
                      broom::tidy(wilcox.test(tmo.w$Q6, tmo.w$Q4, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "thermally m. oak + vacuum + wax"))

fur.p <- d.q6q4 %>% filter(names=="furfurylated pine")
d.q6q4.out <- bind_rows(d.q6q4.out, 
                      broom::tidy(wilcox.test(fur.p$Q6, fur.p$Q4, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "furfurylated pine"))

d.q6q4.out <- d.q6q4.out %>% mutate(
  p.adjusted = p.adjust(p.value, method="fdr")
)

d.q6q4.out %>% kable(digits=4) %>%
  kable_styling()

#write_csv(d.q6q4.out, "data/Q6-Q4_out.csv")

```

# Exploratory q3 to q4

```{r expq3toq4, message=FALSE, warning=FALSE}
d.34 <- full_join(d.q3 %>% select(code, names, Response) %>%
                    rename(Weathered = Response),
                  d.q4 %>% select(code, names, Response) %>%
                    rename(Both=Response),
                  by = c("code"= "code", "names"="names")) %>%
        mutate(diff = Both - Weathered)

d.34 %>% group_by(names) %>%
  summarise(median.weath = median(Weathered),
            median.both = median(Both),
            median.diff = median(diff))
tmp.o <- d.34 %>% filter(names=="thermally modified pine")
d34.out <- broom::tidy(wilcox.test(tmp.o$Both, tmp.o$Weathered, paired=TRUE, conf.int=TRUE)) %>%
  mutate(Name = "thermally modified pine")

acp <- d.34 %>% filter(names=="acetylated pine + coating")
d34.out <- bind_rows(d34.out, 
                      broom::tidy(wilcox.test(acp$Both, acp$Weathered, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "acetylated pine + coating"))

tmo.v <- d.34 %>% filter(names=="thermally modified oak + vacuum")
d34.out <- bind_rows(d34.out, 
                      broom::tidy(wilcox.test(tmo.v$Both, tmo.v$Weathered, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "thermally modified oak + vacuum"))

oak <- d.34 %>% filter(names=="untreated oak")
d34.out <- bind_rows(d34.out, 
                      broom::tidy(wilcox.test(oak$Both, oak$Weathered, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "untreated oak"))

tmo.w <- d.34 %>% filter(names=="thermally m. oak + vacuum + wax")
d34.out <- bind_rows(d34.out, 
                      broom::tidy(wilcox.test(tmo.w$Both, tmo.w$Weathered, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "thermally m. oak + vacuum + wax"))

fur.p <- d.34 %>% filter(names=="furfurylated pine")
d34.out <- bind_rows(d34.out, 
                      broom::tidy(wilcox.test(fur.p$Both, fur.p$Weathered, paired=TRUE, conf.int=TRUE)) %>%
                        mutate(Name = "furfurylated pine"))

d34.out <- d34.out %>% mutate(
  p.adjusted = p.adjust(p.value, method="fdr")
)

#write_csv(d34.out, "data/Q3toQ4_out.csv")
```

# Exploratory Q4 to Q5 (rank vs. rating)

```{r}

d.q4q5 <- full_join(d.q4 %>% select(code, names, Response),
                    d.q5 %>% select(code, names, Response),
                    by=c("code"="code", "names"="names"))

names(d.q4q5)[c(3,4)] <- c("Q4", "Q5")

tau <- cor.test(d.q4q5$Q4, d.q4q5$Q5, method="kendall")

tau
```


# Plots

Plots for presentations and publications

```{r dumbel-q6-q4, fig.width=7.5, fig.height=4}

d.q6q4.sum <- d.q6q4 %>% group_by(names) %>%
  summarise(q6.mu = mean(Q6),
            q4.mu = mean(Q4)) %>%
  mutate(names = fct_relevel(names, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))

d.q6q4.long <- d.q6q4.sum %>%
  pivot_longer(cols=c(q6.mu, q4.mu), names_to="Question", values_to="MeanRating") %>%
  mutate(shp = if_else(Question == "q6.mu", "Q6 Rating", "Q4 Rating"),
         names = fct_relevel(names, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))


ggplot(data=d.q6q4.long, aes(x=names, y=MeanRating, colour=names)) + 
  theme_bw() +
  geom_point(size=4, aes(shape=shp)) +
  geom_segment(data=d.q6q4.sum, aes(x=names, xend=names, 
                   y=q4.mu, yend=q6.mu), 
               arrow=arrow(length=unit(0.3, "cm")), 
               size=.5, show.legend=FALSE) +
  scale_colour_brewer(palette="Dark2", name="Material", guide="none") +
  coord_flip() +
  scale_y_continuous(limits=c(1,7), breaks=seq(1:7)) +
  scale_x_discrete(limits=(rev(levels(d.q6q4.long$names)))) +
  scale_shape_manual(values=c(4,1), name="Value") +
  labs(x="Material",
       y="Rating")

#ggsave("figs/q6q4_dumbell.png", device=png(), height=4, width=8, units="in", dpi=300)
```

```{r dumbel-q3-q2, fig.width=7.5, fig.height=4}
d.q3q2 <- full_join(d.q3 %>% select(code, names, Response), 
                    d.q2 %>% select(code, names, Response),
                    by = c("code"="code", "names"="names"))
names(d.q3q2)[c(3,4)] <- c("Q3", "Q2")

d.q3q2.sum <- d.q3q2 %>% group_by(names) %>%
  summarise(q3.mu = mean(Q3),
            q2.mu = mean(Q2)) %>%
  mutate(names = fct_relevel(names, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))

d.q3q2.long <- d.q3q2.sum %>%
  pivot_longer(cols=c(q3.mu, q2.mu), names_to="Question", values_to="MeanRating") %>%
  mutate(shp = if_else(Question == "q3.mu", "Q3 Rating", "Q2 Rating"),
         names = fct_relevel(names, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))


ggplot(data=d.q3q2.long, aes(x=names, y=MeanRating, colour=names)) + 
  theme_bw() +
  geom_point(size=4, aes(shape=shp)) +
  geom_segment(data=d.q3q2.sum, aes(x=names, xend=names, 
                   y=q2.mu, yend=q3.mu), 
               arrow=arrow(length=unit(0.3, "cm")), 
               size=.5, show.legend=FALSE) +
  scale_colour_brewer(palette="Dark2", name="Material", guide="none") +
  coord_flip() +
  scale_y_continuous(limits=c(1,7), breaks=seq(1:7)) +
  scale_x_discrete(limits=(rev(levels(d.q3q2.long$names)))) +
  scale_shape_manual(values=c(4,1), name="Value") +
  labs(x="Material",
       y="Rating")

#ggsave("figs/q3q2_dumbell.png", device=png(), height=4, width=8, units="in", dpi=300)
```

```{r dumbel-q2-q4, fig.width=7.5, fig.height=4}
d.q4q2 <- full_join(d.q4 %>% select(code, names, Response), 
                    d.q2 %>% select(code, names, Response),
                    by = c("code"="code", "names"="names"))
names(d.q4q2)[c(3,4)] <- c("Q4", "Q2")

d.q4q2.sum <- d.q4q2 %>% group_by(names) %>%
  summarise(q4.mu = mean(Q4),
            q2.mu = mean(Q2)) %>%
  mutate(names = fct_relevel(names, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))

d.q4q2.long <- d.q4q2.sum %>%
  pivot_longer(cols=c(q4.mu, q2.mu), names_to="Question", values_to="MeanRating") %>%
  mutate(shp = if_else(Question == "q4.mu", "Q4 Rating", "Q2 Rating"),
         names = fct_relevel(names, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))


ggplot(data=d.q4q2.long, aes(x=names, y=MeanRating, colour=names)) + 
  theme_bw() +
  geom_point(size=4, aes(shape=shp)) +
  geom_segment(data=d.q4q2.sum, aes(x=names, xend=names, 
                   y=q2.mu, yend=q4.mu), 
               arrow=arrow(length=unit(0.3, "cm")), 
               size=.5, show.legend=FALSE) +
  scale_colour_brewer(palette="Dark2", name="Material", guide="none") +
  coord_flip() +
  scale_y_continuous(limits=c(1,7), breaks=seq(1:7)) +
  scale_x_discrete(limits=(rev(levels(d.q4q2.long$names)))) +
  scale_shape_manual(values=c(4,1), name="Value") +
  labs(x="Material",
       y="Rating")

#ggsave("figs/q4q2_dumbell.png", device=png(), height=4, width=8, units="in", dpi=300)
```

```{r dumbel-q2-q3, fig.width=7.5, fig.height=4}
d.q3q4 <- full_join(d.q3 %>% select(code, names, Response), 
                    d.q4 %>% select(code, names, Response),
                    by = c("code"="code", "names"="names"))
names(d.q3q4)[c(3,4)] <- c("Q3", "Q4")

d.q3q4.sum <- d.q3q4 %>% group_by(names) %>%
  summarise(q3.mu = mean(Q3),
            q4.mu = mean(Q4)) %>%
  mutate(names = fct_relevel(names, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))

d.q3q4.long <- d.q3q4.sum %>%
  pivot_longer(cols=c(q3.mu, q4.mu), names_to="Question", values_to="MeanRating") %>%
  mutate(shp = if_else(Question == "q3.mu", "Q3 Rating", "Q4 Rating"),
         names = fct_relevel(names, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))


ggplot(data=d.q3q4.long, aes(x=names, y=MeanRating, colour=names)) + 
  theme_bw() +
  geom_point(size=4, aes(shape=shp)) +
  geom_segment(data=d.q3q4.sum, aes(x=names, xend=names, 
                   y=q3.mu, yend=q4.mu), 
               arrow=arrow(length=unit(0.3, "cm")), 
               size=.5, show.legend=FALSE) +
  scale_colour_brewer(palette="Dark2", name="Material", guide="none") +
  coord_flip() +
  scale_y_continuous(limits=c(1,7), breaks=seq(1:7)) +
  scale_x_discrete(limits=(rev(levels(d.q3q4.long$names)))) +
  scale_shape_manual(values=c(4,1), name="Value") +
  labs(x="Material",
       y="Rating")

#ggsave("figs/q3q4_dumbell.png", device=png(), height=4, width=8, units="in", dpi=300)
```


Be cautious of the plot below, that doesn't show the full x-axis.

```{r dumbel-q234, fig.width=7.5, fig.height=4}
d.234 <- full_join(d.q2 %>% select(code, names, Response), 
                    d.q3 %>% select(code, names, Response),
                    by = c("code"="code", "names"="names"))
names(d.234)[c(3,4)] <- c("Q2", "Q3")

d.234 <- full_join(d.234, 
                    d.q4 %>% select(code, names, Response),
                    by = c("code"="code", "names"="names"))
names(d.234)[c(5)] <- "Q4"

d.234.sum <- d.234 %>% group_by(names) %>%
  summarise(q4.mu = mean(Q4),
            q3.mu = mean(Q3),
            q2.mu = mean(Q2)) %>%
  mutate(names = fct_relevel(names, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))

d.234.long <- d.234.sum %>%
  pivot_longer(cols=c(2:4), names_to="Question", values_to="MeanRating") %>%
  mutate(shp = if_else(Question == "q4.mu", "Q4 Rating", 
               if_else(Question == "q3.mu", "Q3 Rating", "Q2 Rating")),
         names = fct_relevel(names, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))


ggplot(data=d.234.long, aes(x=names, y=MeanRating, colour=names)) + 
  theme_bw() +
  geom_point(size=4, aes(shape=shp)) +
  geom_segment(data=d.234.sum, aes(x=names, xend=names, 
                   y=q2.mu, yend=q3.mu), 
               arrow=arrow(length=unit(0.3, "cm")), 
               size=.5, show.legend=FALSE) +
  geom_segment(data=d.234.sum, aes(x=names, xend=names, 
                   y=q3.mu, yend=q4.mu), 
               arrow=arrow(length=unit(0.3, "cm")), 
               size=.5, show.legend=FALSE) +
  scale_colour_brewer(palette="Dark2", name="Material", guide="none") +
  coord_flip() +
  scale_y_continuous(limits=c(2.5,5.5), breaks=seq(2.5,5.5,.5)) +
  scale_x_discrete(limits=(rev(levels(d.234.long$names)))) +
  scale_shape_manual(values=c(4,1,5), name="Value") +
  labs(x="Material",
       y="Rating")

#ggsave("figs/q234_dumbell.png", device=png(), height=4, width=8, units="in", dpi=300)
```

```{r probabilitybars, fig.width=7.5, fig.height=4}
h3.fd <- hyp3.out %>% mutate(
  Unlikely = `1` + `2` + `3`,
  Neutral = `4`,
  Likely = `5` + `6` + `7`) %>% select(8:11) %>%
  pivot_longer(c(2:4), names_to="Rating", values_to="Probability") %>%
  mutate(Material = fct_relevel(Material, "thermally modified pine",
                          "thermally m. oak + vacuum + wax",
                          "untreated oak",
                          "acetylated pine + coating",
                          "thermally modified oak + vacuum",
                          "furfurylated pine"))
  

ggplot(data=h3.fd, aes(y=Probability, x=Material, fill=Rating)) + theme_bw() +
  geom_col() +
  scale_fill_brewer(palette = "Dark2") +
  scale_x_discrete(limits=(rev(levels(h3.fd$Material)))) +
  coord_flip() +
  theme(text = element_text(size=24))

#ggsave("figs/hyp3_probabilitiybars.png", device=png(), height=9, width=16, units="in", dpi=300)

q4a.fd <- q4a.pred %>% mutate(
  Unlikely = `1` + `2` + `3`,
  Neutral = `4`,
  Likely = `5` + `6` + `7`,
  Properties = str_c("\U00394", "C=", Colour, ", ",
                     "\U00394", "G=", Gloss, ", ",
                     "\U00394", "Ra=", Roughness ),
  Order = 6:1) %>% select(11:15) %>%
  pivot_longer(c(1:3), names_to="Rating", values_to="Probability") %>%
  mutate(Properties = fct_reorder(Properties, Order))
  

ggplot(data=q4a.fd, aes(y=Probability, x=Properties, fill=Rating)) + theme_bw() +
  geom_col() +
  scale_fill_brewer(palette = "Dark2") +
  coord_flip() +
  theme(text = element_text(size=24))

#ggsave("figs/property_probabilitiybars.png", device=png(), height=9, width=16, units="in", dpi=300)
```

# Environment

```{r environment}
sessionInfo()
```

# Acknowledgement

The authors gratefully acknowledge European Commission for funding the InnoRenew project (Grant Agreement N° 739574) and the Slovenian Research Agency (ARRS) ARCHI-BIO project (BI/US-20-054). The authors also acknowledge the CLICKdesign project, which is supported under the umbrella of ERA-NET Cofund ForestValue by the Ministry of Education, Science and Sport (MIZS) - Slovenia; The Ministry of the Environment (YM) - Finland; The Forestry Commissioners (FC) - UK; Research Council of Norway (RCN) - Norway; The French Environment & Energy Management Agency (ADEME) and The French National Research Agency (ANR) - France; The Swedish Research Council for Environment, Agricultural Sciences and Spatial Planning (FORMAS), Swedish Energy Agency (SWEA), Swedish Governmental Agency for Innovation Systems (Vinnova) - Sweden; Federal Ministry of Food and Agriculture (BMEL) and Agency for Renewable Resources (FNR) - Germany. ForestValue has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement N° 773324.

# License

This work is licensed under [CC BY-SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/?ref=chooser-v1).

# References

```{r refs}
print(citation(), bibtex=FALSE)
print(citation("MASS"), bibtex=FALSE)
print(citation("tidyverse"), bibtex=FALSE)
print(citation("kableExtra"), bibtex=FALSE)
print(citation("fs"), bibtex=FALSE)
print(citation("stringr"), bibtex=FALSE)
print(citation("brant"), bibtex=FALSE)
```

