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.

Data

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

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 **

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.

ggplot(data=d.q1, aes(y=Response, x=SubQuestion)) + theme_bw() +
  geom_col()
## Warning: Removed 330 rows containing missing values (position_stack).

Question 2

Overview figures and summary table.

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))
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 6 x 5
##   SubQuestion  mean median  mode    sd
##   <chr>       <dbl>  <dbl> <dbl> <dbl>
## 1 A            3.93    4       2  1.87
## 2 B            4.30    5       5  1.57
## 3 C            3.73    4       5  1.62
## 4 D            3.14    3       1  1.95
## 5 E            4.51    5       5  1.47
## 6 H            4.22    4.5     6  1.70

Question 3

Overview figures and summary table.

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))
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 6 x 5
##   names                            mean median  mode    sd
##   <fct>                           <dbl>  <dbl> <dbl> <dbl>
## 1 thermally modified pine          3.99      4     6  1.72
## 2 thermally m. oak + vacuum + wax  3.43      3     4  1.54
## 3 untreated oak                    3.22      3     2  1.45
## 4 acetylated pine + coating        3.88      4     1  2.08
## 5 thermally modified oak + vacuum  3.82      4     3  1.47
## 6 furfurylated pine                3.96      4     5  1.68

Question 4

Overview figures and summary table, basi

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))
## `summarise()` ungrouping output (override with `.groups` argument)
d.q4.sum
## # A tibble: 6 x 5
##   names                            mean median    sd  mode
##   <fct>                           <dbl>  <dbl> <dbl> <dbl>
## 1 thermally modified pine          4.61    5    1.59     5
## 2 thermally m. oak + vacuum + wax  3.03    3    1.39     3
## 3 untreated oak                    3.15    3    1.52     2
## 4 acetylated pine + coating        4.28    4.5  2.04     7
## 5 thermally modified oak + vacuum  3.70    4    1.49     3
## 6 furfurylated pine                3.5     3    1.58     3

Question 5

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))
## `summarise()` ungrouping output (override with `.groups` argument)
d.q5q4 <- d.q4.sum %>% mutate(rankq5 = d.q5.sum$rank)
d.q5q4 
## # A tibble: 6 x 6
##   names                            mean median    sd  mode rankq5
##   <fct>                           <dbl>  <dbl> <dbl> <dbl>  <dbl>
## 1 thermally modified pine          4.61    5    1.59     5      1
## 2 thermally m. oak + vacuum + wax  3.03    3    1.39     3      5
## 3 untreated oak                    3.15    3    1.52     2      6
## 4 acetylated pine + coating        4.28    4.5  2.04     7      3
## 5 thermally modified oak + vacuum  3.70    4    1.49     3      2
## 6 furfurylated pine                3.5     3    1.58     3      4
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

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))
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 6 x 4
##   names                            mean median    sd
##   <fct>                           <dbl>  <dbl> <dbl>
## 1 thermally modified pine          4.12      4  1.74
## 2 thermally m. oak + vacuum + wax  4.05      4  1.51
## 3 untreated oak                    4.09      4  1.72
## 4 acetylated pine + coating        3.41      3  2.05
## 5 thermally modified oak + vacuum  4.47      5  1.42
## 6 furfurylated pine                3.15      3  1.64

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.

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
## # A tibble: 6 x 7
##   names          median.diff median.weath median.fresh mu.diff mu.fresh mu.weath
##   <fct>                <dbl>        <dbl>        <dbl>   <dbl>    <dbl>    <dbl>
## 1 thermally mod…           0            4          4    0.0541     3.93     3.99
## 2 thermally m. …          -1            3          5   -0.865      4.30     3.43
## 3 untreated oak            0            3          4   -0.514      3.73     3.22
## 4 acetylated pi…           0            4          3    0.743      3.14     3.88
## 5 thermally mod…          -1            4          5   -0.689      4.51     3.82
## 6 furfurylated …           0            4          4.5 -0.257      4.22     3.96
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

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))
## # A tibble: 6 x 4
##   names                           median.fresh median.both median.diff
##   <fct>                                  <dbl>       <dbl>       <dbl>
## 1 thermally modified pine                  4           5             1
## 2 thermally m. oak + vacuum + wax          5           3            -1
## 3 untreated oak                            4           3             0
## 4 acetylated pine + coating                3           4.5           1
## 5 thermally modified oak + vacuum          5           4            -1
## 6 furfurylated pine                        4.5         3            -1
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.

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)
## Call:
## polr(formula = Resp.ord ~ names + Gender + age_group + Region, 
##     data = d.q4, Hess = TRUE)
## 
## Coefficients:
##                                         Value Std. Error t value
## namesfurfurylated pine               -0.90166     0.3073 -2.9342
## namesthermally m. oak + vacuum + wax -1.39635     0.3088 -4.5218
## namesthermally modified oak + vacuum -0.65926     0.3037 -2.1711
## namesthermally modified pine          0.28888     0.3057  0.9449
## namesuntreated oak                   -1.26998     0.3109 -4.0850
## GenderM                              -0.26101     0.1749 -1.4920
## GenderNot Specified                  -0.08230     0.7601 -0.1083
## age_group26-37                        0.21893     0.2173  1.0075
## age_group38-50                        0.05590     0.2357  0.2371
## age_group51-64                       -0.17428     0.3190 -0.5463
## Regionsuburban                        0.06916     0.2089  0.3310
## Regionurban                           0.28510     0.2437  1.1698
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -2.9764  0.3433    -8.6704
## 2|3 -1.7960  0.3208    -5.5992
## 3|4 -0.7393  0.3119    -2.3700
## 4|5  0.0823  0.3109     0.2648
## 5|6  0.9932  0.3151     3.1518
## 6|7  2.2277  0.3447     6.4636
## 
## Residual Deviance: 1615.88 
## AIC: 1651.88
q4.or.red <- polr(Resp.ord ~ names, 
              data=d.q4, Hess=TRUE)

summary(q4.or.red)
## Call:
## polr(formula = Resp.ord ~ names, data = d.q4, Hess = TRUE)
## 
## Coefficients:
##                                        Value Std. Error t value
## namesfurfurylated pine               -0.8908     0.3068 -2.9035
## namesthermally m. oak + vacuum + wax -1.3676     0.3075 -4.4470
## namesthermally modified oak + vacuum -0.6474     0.3035 -2.1333
## namesthermally modified pine          0.2997     0.3049  0.9829
## namesuntreated oak                   -1.2642     0.3095 -4.0848
## 
## Intercepts:
##     Value    Std. Error t value 
## 1|2  -2.9528   0.2737   -10.7903
## 2|3  -1.7770   0.2470    -7.1945
## 3|4  -0.7317   0.2367    -3.0909
## 4|5   0.0798   0.2338     0.3412
## 5|6   0.9826   0.2396     4.1015
## 6|7   2.2120   0.2782     7.9498
## 
## Residual Deviance: 1621.659 
## AIC: 1643.659
anova(q4.or.red, q4.or.full)
## Likelihood ratio tests of ordinal regression models
## 
## Response: Resp.ord
##                                 Model Resid. df Resid. Dev   Test    Df
## 1                               names       433   1621.659             
## 2 names + Gender + age_group + Region       426   1615.880 1 vs 2     7
##   LR stat.   Pr(Chi)
## 1                   
## 2 5.778438 0.5658452
brant::brant(q4.or.red) 
## ---------------------------------------------------------------------------- 
## Test for                 X2  df  probability 
## ---------------------------------------------------------------------------- 
## Omnibus                      34.37   25  0.1
## namesfurfurylated pine           12.27   5   0.03
## namesthermally m. oak + vacuum + wax 11.56   5   0.04
## namesthermally modified oak + vacuum 12.5    5   0.03
## namesthermally modified pine         8.18    5   0.15
## namesuntreated oak               8.21    5   0.14
## ---------------------------------------------------------------------------- 
## 
## H0: Parallel Regression Assumption holds
##                                             X2 df probability
## Omnibus                              34.369051 25  0.10025200
## namesfurfurylated pine               12.267245  5  0.03130383
## namesthermally m. oak + vacuum + wax 11.561787  5  0.04131155
## namesthermally modified oak + vacuum 12.502439  5  0.02851547
## namesthermally modified pine          8.179242  5  0.14662999
## namesuntreated oak                    8.210728  5  0.14499801
#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()
Coefficient Value Std..Error t.value p.value
namesfurfurylated pine -0.891 0.307 -2.903 0.004
namesthermally m. oak + vacuum + wax -1.368 0.308 -4.447 0.000
namesthermally modified oak + vacuum -0.647 0.303 -2.133 0.033
namesthermally modified pine 0.300 0.305 0.983 0.326
namesuntreated oak -1.264 0.309 -4.085 0.000
1|2 -2.953 0.274 -10.790 0.000
2|3 -1.777 0.247 -7.195 0.000
3|4 -0.732 0.237 -3.091 0.002
4|5 0.080 0.234 0.341 0.733
5|6 0.983 0.240 4.102 0.000
6|7 2.212 0.278 7.950 0.000
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
##            1          2         3         4         5          6          7
## 1 0.03723490 0.07414753 0.1514172 0.1824280 0.2191408 0.20690657 0.12872503
## 2 0.17005089 0.22900487 0.2547542 0.1557764 0.1033574 0.05992539 0.02713092
## 3 0.15595051 0.21858527 0.2555122 0.1630855 0.1112331 0.06563541 0.02999794
## 4 0.04960282 0.09507528 0.1801356 0.1951154 0.2076855 0.17370825 0.09867718
## 5 0.09067426 0.15356727 0.2346873 0.1952531 0.1619817 0.10963959 0.05419683
## 6 0.11283749 0.17905428 0.2477786 0.1855510 0.1416190 0.09016619 0.04299345
##                          Material
## 1         thermally modified pine
## 2 thermally m. oak + vacuum + wax
## 3                   untreated oak
## 4       acetylated pine + coating
## 5 thermally modified oak + vacuum
## 6               furfurylated pine
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.

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
## Likelihood ratio tests of ordinal regression models
## 
## Response: Resp.ord
##                                                                Model Resid. df
## 1                                         Colour + Gloss + Roughness       435
## 2 Region + Gender + country + age_group + Colour + Gloss + Roughness       427
##   Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1   1647.698                                
## 2   1641.949 1 vs 2     8 5.748359 0.6753956
summary(q4a.or.red)
## Call:
## polr(formula = Resp.ord ~ Colour + Gloss + Roughness, data = d.q4a, 
##     Hess = TRUE)
## 
## Coefficients:
##              Value Std. Error t value
## Colour    -0.03142   0.012854 -2.4445
## Gloss      0.03867   0.041280  0.9367
## Roughness -0.01509   0.009126 -1.6532
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -2.8468  0.2849    -9.9925
## 2|3 -1.6967  0.2600    -6.5246
## 3|4 -0.6906  0.2505    -2.7569
## 4|5  0.0813  0.2476     0.3284
## 5|6  0.9476  0.2533     3.7402
## 6|7  2.1593  0.2922     7.3891
## 
## Residual Deviance: 1647.698 
## AIC: 1665.698
q4a.pr <- profile(q4a.or.red)
plot(q4a.pr)

confint(q4a.pr)
##                 2.5 %      97.5 %
## Colour    -0.05669326 -0.00627657
## Gloss     -0.04214281  0.11976549
## Roughness -0.03302866  0.00277097
pairs(q4a.pr)

summary(colour_change)
##     Sample              Colour           Gloss          Roughness     
##  Length:6           Min.   : 1.937   Min.   :-2.640   Min.   : 0.150  
##  Class :character   1st Qu.: 7.704   1st Qu.:-0.955   1st Qu.: 4.694  
##  Mode  :character   Median :13.439   Median :-0.410   Median :12.537  
##                     Mean   :13.013   Mean   : 0.590   Mean   :14.726  
##                     3rd Qu.:18.255   3rd Qu.: 0.900   3rd Qu.:23.732  
##                     Max.   :23.613   Max.   : 6.800   Max.   :33.587  
##     names          
##  Length:6          
##  Class :character  
##  Mode  :character  
##                    
##                    
## 
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()
Colour, Gloss, & Roughness are changes
Coefficient Value Std..Error t.value p.value
Colour -0.0314 0.0129 -2.4445 0.0145
Gloss 0.0387 0.0413 0.9367 0.3489
Roughness -0.0151 0.0091 -1.6532 0.0983
1|2 -2.8468 0.2849 -9.9925 0.0000
2|3 -1.6967 0.2600 -6.5246 0.0000
3|4 -0.6906 0.2505 -2.7569 0.0058
4|5 0.0813 0.2476 0.3284 0.7426
5|6 0.9476 0.2533 3.7402 0.0002
6|7 2.1593 0.2922 7.3891 0.0000
#write_csv(q4a.or.st, "data/q4_cgr_ord.csv")

Difference between Q4 and Q6

Wilcoxon tests.

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()
estimate statistic p.value conf.low conf.high method alternative Name p.adjusted
-0.9999 446.0 0.0240 -1.0000 0.0000 Wilcoxon signed rank test with continuity correction two.sided thermally modified pine 0.0288
-1.5000 157.0 0.0000 -2.0000 -1.0001 Wilcoxon signed rank test with continuity correction two.sided acetylated pine + coating 0.0001
1.0000 1341.0 0.0001 0.5001 1.5000 Wilcoxon signed rank test with continuity correction two.sided thermally modified oak + vacuum 0.0002
1.5000 1092.5 0.0000 1.0000 2.0000 Wilcoxon signed rank test with continuity correction two.sided untreated oak 0.0001
1.5000 1191.0 0.0000 1.0000 2.0000 Wilcoxon signed rank test with continuity correction two.sided thermally m. oak + vacuum + wax 0.0000
-0.5000 721.0 0.0677 -1.0000 0.0000 Wilcoxon signed rank test with continuity correction two.sided furfurylated pine 0.0677
#write_csv(d.q6q4.out, "data/Q6-Q4_out.csv")

Exploratory q3 to q4

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))
## # A tibble: 6 x 4
##   names                           median.weath median.both median.diff
##   <chr>                                  <dbl>       <dbl>       <dbl>
## 1 acetylated pine + coating                  4         4.5           0
## 2 furfurylated pine                          4         3             0
## 3 thermally m. oak + vacuum + wax            3         3             0
## 4 thermally modified oak + vacuum            4         4             0
## 5 thermally modified pine                    4         5             0
## 6 untreated oak                              3         3             0
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)

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
## 
##  Kendall's rank correlation tau
## 
## data:  d.q4q5$Q4 and d.q4q5$Q5
## z = -14.287, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##        tau 
## -0.5259663

Plots

Plots for presentations and publications

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"))
## `summarise()` ungrouping output (override with `.groups` argument)
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)
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"))
## `summarise()` ungrouping output (override with `.groups` argument)
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)
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"))
## `summarise()` ungrouping output (override with `.groups` argument)
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)
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"))
## `summarise()` ungrouping output (override with `.groups` argument)
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.

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"))
## `summarise()` ungrouping output (override with `.groups` argument)
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)
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

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.1.0 fs_1.4.2         forcats_0.5.0    stringr_1.4.0   
##  [5] dplyr_1.0.0      purrr_0.3.4      readr_1.3.1      tidyr_1.1.0     
##  [9] tibble_3.0.3     ggplot2_3.3.2    tidyverse_1.3.0  MASS_7.3-51.6   
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.0   xfun_0.15          lattice_0.20-41    haven_2.3.1       
##  [5] colorspace_1.4-1   vctrs_0.3.2        generics_0.0.2     viridisLite_0.3.0 
##  [9] htmltools_0.5.0    yaml_2.2.1         utf8_1.1.4         blob_1.2.1        
## [13] rlang_0.4.7        pillar_1.4.6       glue_1.4.1         withr_2.2.0       
## [17] DBI_1.1.0          RColorBrewer_1.1-2 dbplyr_1.4.4       modelr_0.1.8      
## [21] readxl_1.3.1       lifecycle_0.2.0    munsell_0.5.0      gtable_0.3.0      
## [25] cellranger_1.1.0   rvest_0.3.5        evaluate_0.14      labeling_0.3      
## [29] knitr_1.29         brant_0.2-0        fansi_0.4.1        highr_0.8         
## [33] broom_0.7.0        Rcpp_1.0.5         scales_1.1.1       backports_1.1.8   
## [37] webshot_0.5.2      jsonlite_1.7.0     farver_2.0.3       hms_0.5.3         
## [41] digest_0.6.25      stringi_1.4.6      grid_4.0.2         cli_2.0.2         
## [45] tools_4.0.2        magrittr_1.5       crayon_1.3.4       pkgconfig_2.0.3   
## [49] Matrix_1.2-18      ellipsis_0.3.1     xml2_1.3.2         reprex_0.3.0      
## [53] lubridate_1.7.9    assertthat_0.2.1   rmarkdown_2.3      httr_1.4.2        
## [57] rstudioapi_0.11    R6_2.4.1           compiler_4.0.2

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.

References

print(citation(), bibtex=FALSE)
## 
## To cite R in publications use:
## 
##   R Core Team (2020). R: A language and environment for statistical
##   computing. R Foundation for Statistical Computing, Vienna, Austria.
##   URL https://www.R-project.org/.
## 
## We have invested a lot of time and effort in creating R, please cite it
## when using it for data analysis. See also 'citation("pkgname")' for
## citing R packages.
print(citation("MASS"), bibtex=FALSE)
## 
## To cite the MASS package in publications use:
## 
##   Venables, W. N. & Ripley, B. D. (2002) Modern Applied Statistics with
##   S. Fourth Edition. Springer, New York. ISBN 0-387-95457-0
print(citation("tidyverse"), bibtex=FALSE)
## 
##   Wickham et al., (2019). Welcome to the tidyverse. Journal of Open
##   Source Software, 4(43), 1686, https://doi.org/10.21105/joss.01686
print(citation("kableExtra"), bibtex=FALSE)
## 
## To cite package 'kableExtra' in publications use:
## 
##   Hao Zhu (2019). kableExtra: Construct Complex Table with 'kable' and
##   Pipe Syntax. R package version 1.1.0.
##   https://CRAN.R-project.org/package=kableExtra
print(citation("fs"), bibtex=FALSE)
## 
## To cite package 'fs' in publications use:
## 
##   Jim Hester and Hadley Wickham (2020). fs: Cross-Platform File System
##   Operations Based on 'libuv'. R package version 1.4.2.
##   https://CRAN.R-project.org/package=fs
print(citation("stringr"), bibtex=FALSE)
## 
## To cite package 'stringr' in publications use:
## 
##   Hadley Wickham (2019). stringr: Simple, Consistent Wrappers for
##   Common String Operations. R package version 1.4.0.
##   https://CRAN.R-project.org/package=stringr
print(citation("brant"), bibtex=FALSE)
## 
## To cite package 'brant' in publications use:
## 
##   Benjamin Schlegel and Marco Steenbergen (2018). brant: Test for
##   Parallel Regression Assumption. R package version 0.2-0.
##   https://CRAN.R-project.org/package=brant
## 
## ATTENTION: This citation information has been auto-generated from the
## package DESCRIPTION file and may need manual editing, see
## 'help("citation")'.