---
title: "Haptic aesthetics and bodily properties of Ori Gersht's digital art a behavioral and eye tracking study"
author: "N. Langiulli"
output:
  html_notebook:
    df_print: paged
    fig_caption: yes
    fig_height: 8
    fig_width: 8
    number_sections: yes
    toc: yes
    toc_float: yes
---

# Latency of First Fixation -----------------------

Normalized Latency of First Fixtion
Database Import
```{r}
FixLatency <- readr::read_delim("FixLatency.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)
FixLatency <- subset(FixLatency, !(FixLatency$Question == "Baseline"))
FixLatency$AOI <- as.factor(FixLatency$AOI)
FixLatency$Magnitude <- as.factor(FixLatency$Magnitude)
levels(FixLatency$Magnitude) <- c("Detail", "Entire")
print(FixLatency)
```

Latency Definitive Model
```{r}
Lat <- lme4::lmer(LatencyNorm ~  AOI + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixLatency)
```

Latency Overall Statistics
```{r}
AIC(Lat)
```
```{r}
MuMIn::r.squaredGLMM(Lat)
```

Latency Summary & ANOVA
```{r}
summary(Lat)
```
```{r}
aovLat <- aov(LatencyNorm ~ AOI + AOI:Magnitude + Error(RecordingNameII/(Resolution + Magnitude)), data = FixLatency)
summary(aovLat)
```

Latency Robustlmm
```{r}
require(lme4)
require(robustlmm)
robustLat <- robustlmm::rlmer(LatencyNorm ~  AOI + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixLatency, rho.sigma.e = psi2propII(smoothPsi, k = 2.28), rho.sigma.b = psi2propII(smoothPsi,k = 2.28))
```

```{r}
summary(robustLat)
```

Comparing lmer-rlmer models
```{r}
campareLat <- robustlmm::compare(Lat, robustLat, show.rho.functions = FALSE)
```
Get the bootstrapped confidence intervals for the model parameters
```{r}
ciLat <- lme4::confint.merMod(Lat, method = "boot", boot.type = "perc", nsim = 10000, oldNames = F)
```

```{r}
ciLat
```



Latency Post Hoc
```{r}
emmeans::emm_options(pbkrtest.limit = 9203)
emmLat <- emmeans::emmeans(Lat, pairwise ~ AOI:Magnitude, adjust="tukey")
```

```{r}
emmLat
```


## Latency of First Fixation Hierarchical Steps ----

Latency Nested Models
```{r}
null <- nlme::gls(LatencyNorm ~ 1, data = FixLatency)
AOI <- lm(LatencyNorm ~ AOI, data = FixLatency)
AOI.intAOIM <- lm(LatencyNorm ~ AOI + AOI:Magnitude, data = FixLatency)
AOI.intAOIM.rand.intAOIM <- lme4::lmer(LatencyNorm ~ AOI + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixLatency)
```

Latency Nested Models AIC
```{r}
MuMIn::AICc(null)
MuMIn::AICc(AOI)
MuMIn::AICc(AOI.intAOIM)
MuMIn::AICc(AOI.intAOIM.rand.intAOIM)
```

Latency Nested Models BIC
```{r}
BIC(null)
BIC(AOI)
BIC(AOI.intAOIM)
BIC(AOI.intAOIM.rand.intAOIM)
```

Latency Nested Models R squared
```{r}
MuMIn::r.squaredLR(null)
MuMIn::r.squaredGLMM(AOI)
MuMIn::r.squaredGLMM(AOI.intAOIM)
MuMIn::r.squaredGLMM(AOI.intAOIM.rand.intAOIM)
```

Latency Nested Models Likelihood ratio test
```{r}
lmtest::lrtest(null, AOI, AOI.intAOIM, AOI.intAOIM.rand.intAOIM)
```

## Latency of First Fixation Plots -----------------

Interaction AOI * MAgnitude

Create Database with Latency fitted values (representing the effect on the scale of the linear predictor)
```{r}
efLat <- effects::effect("AOI:Magnitude", Lat)
summary(efLat)
efLat <- as.data.frame(efLat)
```
Create Database with Latency means by Subject
```{r}
require(dplyr)
LatmeansSubj <- FixLatency %>% 
  group_by(AOI, Magnitude, RecordingNameII) %>% 
  summarise(RespM = mean(LatencyNorm), Respse =  sjstats::se(LatencyNorm))
print(LatmeansSubj)
```

#### Plot Latency fitted values {.tabset}
```{r}
require(ggplot2)

Latbar <- ggplot(efLat, aes(AOI, fit, fill= Magnitude)) + 
  geom_point() + 
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4) + theme_bw(base_size=12)
Latbar <- ggplot(efLat, aes(AOI, fit, fill=Magnitude)) + geom_bar(stat="identity", position="dodge") +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4, position=position_dodge(width=0.9)) + 
  ylim(0,1.3) +
  labs(x= "AOI", y= "Normalized First Fixation Latency") +
  scale_fill_grey(name = "Magnitude") +
  ggtitle("AOI*Magnitude") +
  theme_classic() +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(Latbar)
```

#### Plot Latency Means by Subject {.tabset}
```{r}
require(ggplot2)
LatsplitSUBJ <- ggplot(LatmeansSubj, aes(AOI, RespM, fill= Magnitude)) + 
    geom_bar(stat="identity", position="dodge") +
    facet_wrap(~ RecordingNameII) + 
    xlab("AOI") + 
    ylab("Normalized First Fixation Latency") +
    scale_fill_grey(name = "Magnitude")
print(LatsplitSUBJ)
```

#### Boxplot Latency {.tabset}
```{r}
Latbox <- ggplot(data = FixLatency, aes(x = Magnitude, y = LatencyNorm)) +
  facet_wrap(. ~ AOI) + 
  geom_jitter(alpha = 0.1, color = "black") +
  geom_boxplot(alpha = 0) +
  theme_light() +
  labs(tag = "A", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized First Fixation Latency") +
  theme(plot.title = element_text(hjust = 0.5))
print(Latbox)

```

#### Violinplot Latency {.tabset}

Function: plot split violins
```{r}
GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, 
                           draw_group = function(self, data, ..., draw_quantiles = NULL) {
                             data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
                             grp <- data[1, "group"]
                             newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
                             newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
                             newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
                             
                             if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
                               stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
                                                                         1))
                               quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
                               aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
                               aesthetics$alpha <- rep(1, nrow(quantiles))
                               both <- cbind(quantiles, aesthetics)
                               quantile_grob <- GeomPath$draw_panel(both, ...)
                               ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
                             }
                             else {
                               ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
                             }
                           })

geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., 
                              draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, 
                              show.legend = NA, inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, 
        position = position, show.legend = show.legend, inherit.aes = inherit.aes, 
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}
```

Violinplot Latency
```{r}

Latsplitviolin <- ggplot(FixLatency, aes(x =Magnitude, y=LatencyNorm, fill=AOI)) + 
  geom_split_violin() +
  scale_fill_grey() +
  theme_light() +
  labs(tag = "A", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized First Fixation Latency") +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(Latsplitviolin)

```


# Total Number of Fixations ----------------------

Normalized Total Number of Fixations
Database Import
```{r}
FixNumber <- readr::read_delim("FixNumber.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)
FixNumber <- subset(FixNumber, !(FixNumber$Question == "Baseline"))
FixNumber$AOI <- as.factor(FixNumber$AOI)
FixNumber$Magnitude <- as.factor(FixNumber$Magnitude)
levels(FixNumber$Magnitude) <- c("Detail", "Entire")
FixNumber$Resolution <- as.factor(FixNumber$Resolution)
print(FixNumber)
```

NFix Definitive Model
```{r}
NFix <- lme4::lmer(NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixNumber)
```

NFix Overall Statistics
```{r}
AIC(NFix)
```
```{r}
MuMIn::r.squaredGLMM(NFix)
```

NFix Summary & ANOVA
```{r}
summary(NFix)
```
```{r}
aovNFix <- aov(NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + Error(RecordingNameII/(Resolution + Magnitude)), data = FixNumber)
summary(aovNFix)
```

NFix Robustlmm
```{r}
require(lme4)
require(robustlmm)
robustNFix <- robustlmm::rlmer(NFixAOIQRMNorm ~  AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixNumber, rho.sigma.e = psi2propII(smoothPsi, k = 2.28), rho.sigma.b = psi2propII(smoothPsi, k = 2.28))
```

```{r}
summary(robustNFix)
```

Comparing lmer-rlmer models
```{r}
campareNFix <- robustlmm::compare(NFix, robustNFix, show.rho.functions = FALSE)
```
Get the bootstrapped confidence intervals for the model parameters
```{r}
ciNFix <- lme4::confint.merMod(NFix, method = "boot", boot.type = "perc", nsim = 10000, oldNames = F)
```

```{r}
ciNFix
```


NFix Post Hoc
```{r}
emmeans::emm_options(pbkrtest.limit = 9278)
emmNFix <- emmeans::emmeans(NFix, pairwise ~ AOI:Magnitude, adjust="tukey") 
```

```{r}
emmNFix
```


## Total Number of Fixations Hierarchical Steps ---

Total Number of Fixations Nested Models
```{r}
null <- nlme::gls(NFixAOIQRMNorm ~ 1, data = FixNumber)
AOI <- lm(NFixAOIQRMNorm ~ AOI, data = FixNumber)
AOIR <- lm(NFixAOIQRMNorm ~ AOI + Resolution, data = FixNumber)
AOIRM <- lm(NFixAOIQRMNorm ~ AOI + Resolution + Magnitude, data = FixNumber)
AOIRM.intAOIM <- lm(NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude, data = FixNumber)
AOIRM.intAOIM.rand.intAOIM <- lme4::lmer(NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixNumber)
```

Total Number of Fixations Nested Models AIC
```{r}
MuMIn::AICc(null)
MuMIn::AICc(AOI)
MuMIn::AICc(AOIR)
MuMIn::AICc(AOIRM)
MuMIn::AICc(AOIRM.intAOIM)
MuMIn::AICc(AOIRM.intAOIM.rand.intAOIM)
```

Total Number of Fixations Nested Models BIC
```{r}
BIC(null)
BIC(AOI)
BIC(AOIR)
BIC(AOIRM)
BIC(AOIRM.intAOIM)
BIC(AOIRM.intAOIM.rand.intAOIM)
```

Total Number of Fixations Nested Models R squared
```{r}
MuMIn::r.squaredLR(null)
MuMIn::r.squaredGLMM(AOI)
MuMIn::r.squaredGLMM(AOIR)
MuMIn::r.squaredGLMM(AOIRM)
MuMIn::r.squaredGLMM(AOIRM.intAOIM)
MuMIn::r.squaredGLMM(AOIRM.intAOIM.rand.intAOIM)
```

Total Number of Fixations Nested Models Likelihood ratio test
```{r}
lmtest::lrtest(null, AOI, AOIR, AOIRM, AOIRM.intAOIM, AOIRM.intAOIM.rand.intAOIM)
```

## Total Number of Fixations Plots ----------------

Interaction AOI * MAgnitude

Create Database with Total Number of Fixations fitted values (representing the effect on the scale of the linear predictor)
```{r}
efNFix <- effects::effect("AOI:Magnitude", NFix)
summary(efNFix)
efNFix <- as.data.frame(efNFix)
```

Create Database with Total Number of Fixations means by Subject
```{r}
require(dplyr)
NFixmeansSubj <- FixNumber %>% 
  group_by(AOI, Magnitude, RecordingNameII) %>% 
  summarise(RespM = mean(NFixAOIQRMNorm), Respse = sjstats::se(NFixAOIQRMNorm))
print(NFixmeansSubj)
```

#### Plot Total Number of Fixations  fitted values {.tabset}
```{r}
require(ggplot2)

NFixbar <- ggplot(efNFix, aes(AOI, fit, fill= Magnitude)) + 
  geom_point() + 
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4) + theme_bw(base_size=12)
NFixbar <- ggplot(efNFix, aes(AOI, fit, fill=Magnitude)) + geom_bar(stat="identity", position="dodge") +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4, position=position_dodge(width=0.9)) + 
  ylim(0,1.3) +
  labs(x= "AOI", y= "Normalized Total Fixations Number") +
  scale_fill_grey(name = "Magnitude") +
  ggtitle("AOI*Magnitude") +
  theme_classic() +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(NFixbar)
```

#### Plot Total Number of Fixations means by Subject {.tabset}
```{r}
require(ggplot2)
NFixsplitSUBJ <- ggplot(NFixmeansSubj, aes(AOI, RespM, fill= Magnitude)) + 
    geom_bar(stat="identity", position="dodge") +
    facet_wrap(~ RecordingNameII) + 
    xlab("AOI") + 
    ylab("Normalized Total Fixations Number") +
    scale_fill_grey(name = "Magnitude", breaks=c("Full", "Detail"), labels=c("Entire", "Detail"))
print(NFixsplitSUBJ)
```

#### Boxplot Total Number of Fixations {.tabset}
```{r}
NFixbox <- ggplot(data = FixNumber, aes(x = Magnitude, y = NFixAOIQRMNorm)) +
  facet_wrap(. ~ AOI) + 
  geom_jitter(alpha = 0.1, color = "black") + 
  geom_boxplot(alpha = 0) + 
  theme_light() + 
  labs(tag = "B", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized Total Fixations Number") +
  theme(plot.title = element_text(hjust = 0.5))
print(NFixbox)

```


#### Violinplot Total Number of Fixations {.tabset}

Function: plot split violins
```{r}
GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, 
                           draw_group = function(self, data, ..., draw_quantiles = NULL) {
                             data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
                             grp <- data[1, "group"]
                             newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
                             newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
                             newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
                             
                             if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
                               stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
                                                                         1))
                               quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
                               aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
                               aesthetics$alpha <- rep(1, nrow(quantiles))
                               both <- cbind(quantiles, aesthetics)
                               quantile_grob <- GeomPath$draw_panel(both, ...)
                               ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
                             }
                             else {
                               ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
                             }
                           })

geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., 
                              draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, 
                              show.legend = NA, inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, 
        position = position, show.legend = show.legend, inherit.aes = inherit.aes, 
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}
```

Violinplot Total Number of Fixations
```{r}
NFixsplitviolin <- ggplot(FixNumber, aes(x=Magnitude, y=NFixAOIQRMNorm, fill=AOI)) + 
  geom_split_violin() +
  scale_fill_grey() +
  theme_light() +
  labs(tag = "B", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized Total Fixations Number") +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(NFixsplitviolin)
```

# Total Fixations Duration -----------------------

Normalized Total Fixations Duration
Database Import
```{r}
FixDuration <- readr::read_delim("FixDuration.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)
FixDuration <- subset(FixDuration, !(FixDuration$Question == "Baseline"))
FixDuration$AOI <- as.factor(FixDuration$AOI)
FixDuration$Magnitude <- as.factor(FixDuration$Magnitude)
levels(FixDuration$Magnitude) <- c("Detail", "Entire")
FixDuration$Resolution <- as.factor(FixDuration$Resolution)
print(FixDuration)
```

Dur Definitive Model
```{r}
Dur <- lme4::lmer(DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixDuration)
```

Dur Overall Statistics
```{r}
AIC(Dur)
```
```{r}
MuMIn::r.squaredGLMM(Dur)
```

Dur Summary & ANOVA
```{r}
summary(Dur)
```
```{r}
aovDur <- aov(DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + Error(RecordingNameII/(Resolution + Magnitude)), data = FixDuration)
summary(aovDur)
```

Dur Robustlmm
```{r}
require(lme4)
require(robustlmm)
robustDur <- robustlmm::rlmer(DurAOIQRMNorm ~  AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixDuration, rho.sigma.e = psi2propII(smoothPsi, k = 2.28), rho.sigma.b = psi2propII(smoothPsi, k = 2.28))
```

```{r}
summary(robustDur)
```

Comparing lmer-rlmer models
```{r}
campareDur <- robustlmm::compare(Dur, robustDur, show.rho.functions = FALSE)
```
Get the bootstrapped confidence intervals for the model parameters
```{r}
ciDur <- lme4::confint.merMod(Dur, method = "boot", boot.type = "perc", nsim = 10000, oldNames = F)
```

```{r}
ciDur
```


Dur Post Hoc
```{r}
emmeans::emm_options(pbkrtest.limit = 9278)
emmDur <- emmeans::emmeans(Dur, pairwise ~ AOI:Magnitude, adjust="tukey")
```

```{r}
emmDur
```


## Total Fixations Duration Hierarchical Steps ----

Total Fixations Duration Nested Models
```{r}
null <- nlme::gls(DurAOIQRMNorm ~ 1, data = FixDuration)
AOI <- lm(DurAOIQRMNorm ~ AOI, data = FixDuration)
AOIR <- lm(DurAOIQRMNorm ~ AOI + Resolution, data = FixDuration)
AOIRM <- lm(DurAOIQRMNorm ~ AOI + Resolution + Magnitude, data = FixDuration)
AOIRM.intAOIM <- lm(DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude, data = FixDuration)
AOIRM.intAOIM.rand.intAOIM <- lme4::lmer(DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixDuration)
```

Total Fixations Duration Nested Models AIC
```{r}
MuMIn::AICc(null)
MuMIn::AICc(AOI)
MuMIn::AICc(AOIR)
MuMIn::AICc(AOIRM)
MuMIn::AICc(AOIRM.intAOIM)
MuMIn::AICc(AOIRM.intAOIM.rand.intAOIM)
```

Total Fixations Duration Nested Models BIC
```{r}
BIC(null)
BIC(AOI)
BIC(AOIR)
BIC(AOIRM)
BIC(AOIRM.intAOIM)
BIC(AOIRM.intAOIM.rand.intAOIM)
```

Total Fixations Duration Nested Models R squared
```{r}
MuMIn::r.squaredLR(null)
MuMIn::r.squaredGLMM(AOI)
MuMIn::r.squaredGLMM(AOIR)
MuMIn::r.squaredGLMM(AOIRM)
MuMIn::r.squaredGLMM(AOIRM.intAOIM)
MuMIn::r.squaredGLMM(AOIRM.intAOIM.rand.intAOIM)
```

Total Fixations Duration Nested Models Likelihood ratio test
```{r}
lmtest::lrtest(null, AOI, AOIR, AOIRM, AOIRM.intAOIM, AOIRM.intAOIM.rand.intAOIM)
```

## Total Fixations Duration Plots -----------------

Create Database with Total Fixations Duration fitted values (representing the effect on the scale of the linear predictor)
Interaction AOI * MAgnitude
```{r}
efDur <- effects::effect("AOI:Magnitude", Dur)
summary(efDur)
efDur <- as.data.frame(efDur)
```

Plot Total Fixations Duration fitted values
```{r}
require(ggplot2)

Durbar <- ggplot(efDur, aes(AOI, fit, fill= Magnitude)) + 
  geom_point() + 
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4) + theme_bw(base_size=12)
Durbar <- ggplot(efDur, aes(AOI, fit, fill=Magnitude)) + geom_bar(stat="identity", position="dodge") +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4, position=position_dodge(width=0.9)) + 
  ylim(0,1.3) +
  labs(x= "AOI", y= "Normalized Total Fixations Duration") +
  scale_fill_grey(name = "Magnitude") +
  ggtitle("AOI*Magnitude") +
  theme_classic() +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(Durbar)
```

Create Database with Total Fixations Duration means by Subject
```{r}
require(dplyr)
DurmeansSubj <- FixDuration %>% 
  group_by(AOI, Magnitude, RecordingNameII) %>% 
  summarise(RespM = mean(DurAOIQRMNorm), Respse = sjstats::se(DurAOIQRMNorm))
print(DurmeansSubj)
```

Plot Total Fixations Duration means by Subject
```{r}
require(ggplot2)
DursplitSUBJ <- ggplot(DurmeansSubj, aes(AOI, RespM, fill= Magnitude)) + 
    geom_bar(stat="identity", position="dodge") +
    facet_wrap(~ RecordingNameII) + 
    xlab("AOI") + 
    ylab("Normalized Total Fixations Duration") +
    scale_fill_grey(name = "Magnitude", breaks=c("Full", "Detail"), labels=c("Entire", "Detail"))
print(DursplitSUBJ)
```

Boxplot Total Fixations Duration
```{r}
Durbox <- ggplot(data = FixDuration, aes(x = Magnitude, y = DurAOIQRMNorm)) +
  facet_wrap(. ~ AOI) + 
  geom_jitter(alpha = 0.1, color = "black") + 
  geom_boxplot(alpha = 0) +
  theme_light() +
  labs(tag = "C", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized Total Fixations Duration") +
  theme(plot.title = element_text(hjust = 0.5))
print(Durbox)
```

Function: plot split violins
```{r}
GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, 
                           draw_group = function(self, data, ..., draw_quantiles = NULL) {
                             data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
                             grp <- data[1, "group"]
                             newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
                             newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
                             newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
                             
                             if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
                               stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
                                                                         1))
                               quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
                               aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
                               aesthetics$alpha <- rep(1, nrow(quantiles))
                               both <- cbind(quantiles, aesthetics)
                               quantile_grob <- GeomPath$draw_panel(both, ...)
                               ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
                             }
                             else {
                               ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
                             }
                           })

geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., 
                              draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, 
                              show.legend = NA, inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, 
        position = position, show.legend = show.legend, inherit.aes = inherit.aes, 
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}
```

Violinplot Total Fixations Duration
```{r}
Dursplitviolin <- ggplot(FixDuration, aes(x=Magnitude, y=DurAOIQRMNorm, fill=AOI)) + 
  geom_split_violin() +
  scale_fill_grey() +
  theme_light() +
  labs(tag = "C", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized Total Fixations Duration") +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(Dursplitviolin)
```

# Behavioral -------------------------------------

Normalized Behavioral Scores
Database Import
```{r}

BehavioralD <- readr::read_delim("BehavioralD.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)
BehavioralD <- subset(BehavioralD, !(BehavioralD$Question == "Baseline"))
BehavioralD <- subset(BehavioralD, !(is.na(BehavioralD$Resp)))
BehavioralD$Magnitude <- as.factor(BehavioralD$Magnitude)
levels(BehavioralD$Magnitude) <- c("Detail", "Entire")
print(BehavioralD)
```

Behavioral Cumulative Link Model for Ordinal Data (Likert)
```{r}
Behavioral <- ordinal::clm(as.factor(Resp) ~ Question + Resolution + Magnitude + Question:Resolution + Question:Magnitude, data = BehavioralD, method = c("Newton"),  threshold="equidistant")

```

Behavioral Overall Statistics
```{r}
AIC(Behavioral)
```

Behavioral Summary & ANOVA
```{r}
summary(Behavioral)
```

```{r}
stats::anova(Behavioral, type="II")
```



Behavioral Post Hoc

Interaction Question * Resolution
```{r}
intQR <- emmeans::emmeans(Behavioral, pairwise ~  Question:Resolution, adjust="tukey")
# contrast                            estimate    SE  df z.ratio p.value
# Beautiful,Detail - Like,Detail      -0.05378 0.116 Inf -0.465  1.0000
# Beautiful,Detail - Movement,Detail  -0.11577 0.117 Inf -0.992  0.9928
# Beautiful,Detail - Proximity,Detail -0.24804 0.114 Inf -2.171  0.4770
# Beautiful,Detail - Touch,Detail      0.46965 0.117 Inf  4.009  0.0025 *
# Like,Detail - Movement,Detail       -0.06198 0.117 Inf -0.531  1.0000
# Like,Detail - Proximity,Detail      -0.19425 0.114 Inf -1.698  0.7971
# Like,Detail - Touch,Detail           0.52344 0.117 Inf  4.461  0.0003 *
# Movement,Detail - Proximity,Detail  -0.13227 0.115 Inf -1.146  0.9799
# Movement,Detail - Touch,Detail       0.58542 0.118 Inf  4.949  <.0001 *
# Proximity,Detail - Touch,Detail      0.71769 0.116 Inf  6.185  <.0001 *
# 
# Beautiful,Detail - Beautiful,Full   -0.05532 0.115 Inf -0.482  1.0000
# Like,Detail - Like,Full              0.05639 0.115 Inf  0.489  1.0000
# Movement,Detail - Movement,Full     -0.05216 0.119 Inf -0.439  1.0000
# Proximity,Detail - Proximity,Full    0.45094 0.112 Inf  4.027  0.0023 *
# Touch,Detail - Touch,Full            0.25190 0.117 Inf  2.150  0.4920
# 
# Beautiful,Full - Like,Full           0.05793 0.115 Inf  0.506  1.0000
# Beautiful,Full - Movement,Full      -0.11261 0.117 Inf -0.960  0.9943
# Beautiful,Full - Proximity,Full      0.25822 0.113 Inf  2.295  0.3923
# Beautiful,Full - Touch,Full          0.77687 0.115 Inf  6.730  <.0001 *
# Like,Full - Movement,Full           -0.17054 0.118 Inf -1.450  0.9113
# Like,Full - Proximity,Full           0.20029 0.113 Inf  1.775  0.7515
# Like,Full - Touch,Full               0.71894 0.116 Inf  6.212  <.0001 *
# Movement,Full - Proximity,Full       0.37082 0.116 Inf  3.207  0.0438
# Movement,Full - Touch,Full           0.88948 0.118 Inf  7.508  <.0001 *
# Proximity,Full - Touch,Full          0.51866 0.113 Inf  4.573  0.0002 *
```




Interaction Question * Magnitude
```{r}
intQM <- emmeans::emmeans(Behavioral, pairwise ~  Question:Magnitude, adjust="tukey")
# contrast                        estimate    SE  df z.ratio p.value
# Beautiful,High - Like,High        0.0413 0.117 Inf  0.353  1.0000
# Beautiful,High - Movement,High    0.1791 0.120 Inf  1.495  0.8946
# Beautiful,High - Proximity,High   0.2671 0.115 Inf  2.318  0.3770
# Beautiful,High - Touch,High       0.7910 0.118 Inf  6.682  <.0001 *
# Like,High - Movement,High         0.1378 0.120 Inf  1.152  0.9791
# Like,High - Proximity,High        0.2258 0.115 Inf  1.963  0.6252
# Like,High - Touch,High            0.7497 0.118 Inf  6.347  <.0001 *
# Movement,High - Proximity,High    0.0880 0.118 Inf  0.747  0.9992
# Movement,High - Touch,High        0.6118 0.121 Inf  5.073  <.0001 *
# Proximity,High - Touch,High       0.5239 0.116 Inf  4.517  0.0003 *
# 
# Beautiful,High - Beautiful,Low    0.9442 0.116 Inf  8.171  <.0001 *
# Like,High - Like,Low              0.8657 0.116 Inf  7.470  <.0001 *
# Movement,High - Movement,Low      0.3576 0.119 Inf  3.001  0.0802
# Proximity,High - Proximity,Low    0.4202 0.112 Inf  3.754  0.0067 *
# Touch,High - Touch,Low            0.6088 0.117 Inf  5.188  <.0001 *
# 
# Beautiful,Low - Like,Low         -0.0372 0.113 Inf -0.328  1.0000
# Beautiful,Low - Movement,Low     -0.4075 0.114 Inf -3.566  0.0134 *
# Beautiful,Low - Proximity,Low    -0.2569 0.112 Inf -2.302  0.3872
# Beautiful,Low - Touch,Low         0.4555 0.114 Inf  3.986  0.0027 *
# Like,Low - Movement,Low          -0.3703 0.115 Inf -3.222  0.0419
# Like,Low - Proximity,Low         -0.2198 0.112 Inf -1.957  0.6292
# Like,Low - Touch,Low              0.4927 0.115 Inf  4.285  0.0008 *
# Movement,Low - Proximity,Low      0.1506 0.113 Inf  1.330  0.9470
# Movement,Low - Touch,Low          0.8631 0.116 Inf  7.429  <.0001 *
# Proximity,Low - Touch,Low         0.7125 0.113 Inf  6.280  <.0001 *
```
## Behavioral Hierarchical Steps ------------------

BehavioralNested Models
```{r}
null <- ordinal::clm(as.factor(Resp) ~ 1, data = BehavioralD, method = c("Newton"),  threshold="equidistant")
Q <- ordinal::clm(as.factor(Resp)  ~ Question, data = BehavioralD, method = c("Newton"),  threshold="equidistant")
QR <- ordinal::clm(as.factor(Resp)  ~ Question + Resolution, data = BehavioralD, method = c("Newton"),  threshold="equidistant")
QRM <- ordinal::clm(as.factor(Resp)  ~ Question + Resolution + Magnitude, data = BehavioralD, method = c("Newton"),  threshold="equidistant")
QRM.intQR <- ordinal::clm(as.factor(Resp)  ~ Question + Resolution + Magnitude + Question:Resolution, data = BehavioralD, method = c("Newton"),  threshold="equidistant")
QRM.intQR.intQM <- ordinal::clm(as.factor(Resp)  ~ Question + Resolution + Magnitude + Question:Resolution +  Question:Magnitude, data = BehavioralD, method = c("Newton"),  threshold="equidistant")

```

Behavioral Nested Models AIC
```{r}
MuMIn::AICc(null)
MuMIn::AICc(Q)
MuMIn::AICc(QR)
MuMIn::AICc(QRM) 
MuMIn::AICc(QRM.intQR)
MuMIn::AICc(QRM.intQR.intQM)
```


Behavioral Nested Models Likelihood ratio test
```{r}
lmtest::lrtest(null, Q, QR, QRM, QRM.intQR, QRM.intQR.intQM)
```

## Behavioral Plots -------------------------------

#### Diverging stacked bar chart Question * Resolution {.tabset}
```{r}
QRLikert <- readr::read_delim("QRLikert.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)

require(sjPlot)
require(sjmisc)

levels_5 <- c("1", "2", "3", "4", "5")

itemsQR <- c("Beautiful-High","Like-High","Movement-High","Proximity-High","Touch-High",
              "Beautiful-Low","Like-Low","Movement-Low","Proximity-Low","Touch-Low")


plot_likert(QRLikert,
            expand.grid = FALSE,
            axis.labels = itemsQR,
            legend.labels = levels_5,
            cat.neutral = 3,
            title = "Question * Resolution",
            legend.title = "Likert Scores",
            show.n = FALSE,
            coord.flip = T) + 
  font_size(labels.y = 13, title = 15, base.theme = legend_style(pos = "bottom", justify = "center", base.theme = theme_sjplot(base_size = 10)))

```

#### Diverging stacked bar chart Question * Magnitude {.tabset}
```{r}
QMLikert <- readr::read_delim("QMLikert.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)

require(sjPlot)
require(sjmisc)

levels_5 <- c("1", "2", "3", "4", "5")

itemsQM <- c("Beautiful-Detail","Like-Detail","Movement-Detail","Proximity-Detail","Touch-Detail",
             "Beautiful-Entire","Like-Entire","Movement-Entire","Proximity-Entire","Touch-Entire")


plot_likert(QMLikert,
            expand.grid = FALSE,
            axis.labels = itemsQM,
            legend.labels = levels_5,
            cat.neutral = 3,
            title = "Question * Magnitude",
            legend.title = "Likert Scores",
            show.n = FALSE,
            coord.flip = T) + 
  font_size(labels.y = 13, title = 15, base.theme = legend_style(pos = "bottom", justify = "center", base.theme = theme_sjplot(base_size = 10)))

```



## Behavioral Corr --------------------------------

Normalized Behavioral Scores means by Subject
Dataset Import
```{r}
Correlations <- readr::read_delim("Correlations.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)
print(Correlations)
```

Correlation Among Questions (Bonferroni's Correction -> 0.05/10 = 0.005)
```{r}
BLc <- cor.test(Correlations$BeautifulMeansNorm, Correlations$LikeMeansNorm, method = c("spearman"))
ifelse(BLc$p.value < 0.005, print(BLc), NA)

BMc <- cor.test(Correlations$BeautifulMeansNorm, Correlations$MovementMeansNorm, method = c("spearman"))
ifelse(BMc$p.value < 0.005, print(BMc), NA)

BPc <- cor.test(Correlations$BeautifulMeansNorm, Correlations$ProximityMeansNorm, method = c("spearman"))
ifelse(BPc$p.value < 0.005, print(BPc), NA)

BTc <- cor.test(Correlations$BeautifulMeansNorm, Correlations$TouchMeansNorm, method = c("spearman"))
ifelse(BTc$p.value < 0.005, print(BTc), NA)

LMc <- cor.test(Correlations$LikeMeansNorm, Correlations$MovementMeansNorm, method = c("spearman"))
ifelse(LMc$p.value < 0.005, print(LMc), NA)

LPc <- cor.test(Correlations$LikeMeansNorm, Correlations$ProximityMeansNorm, method = c("spearman"))
ifelse(LPc$p.value < 0.005, print(LPc), NA)

LTc <- cor.test(Correlations$LikeMeansNorm, Correlations$TouchMeansNorm, method = c("spearman"))
ifelse(LTc$p.value < 0.005, print(LTc), NA)

MPc <- cor.test(Correlations$MovementMeansNorm, Correlations$ProximityMeansNorm, method = c("spearman"))
ifelse(MPc$p.value < 0.005, print(MPc), NA)

MTc <- cor.test(Correlations$MovementMeansNorm, Correlations$TouchMeansNorm, method = c("spearman"))
ifelse(MTc$p.value < 0.005, print(MTc), NA)

PTc <- cor.test(Correlations$ProximityMeansNorm, Correlations$TouchMeansNorm, method = c("spearman"))
ifelse(PTc$p.value < 0.005, print(PTc), NA)
```

Correlationtion Questions * IRI Fantasy Scale (Bonferroni's Correction -> 0.05/5 = 0.01)
```{r}
BIc <- cor.test(Correlations$BeautifulMeansNorm, Correlations$IRIF.z, method = c("spearman"))
ifelse(BIc$p.value < 0.01, print(BIc), NA)

MIc <- cor.test(Correlations$MovementMeansNorm, Correlations$IRIF.z, method = c("spearman"))
ifelse(MIc$p.value < 0.01, print(MIc), NA)

PIc <- cor.test(Correlations$ProximityMeansNorm, Correlations$IRIF.z, method = c("spearman"))
ifelse(PIc$p.value < 0.01, print(PIc), NA)

TIc <- cor.test(Correlations$TouchMeansNorm, Correlations$IRIF.z, method = c("spearman"))
ifelse(TIc$p.value < 0.01, print(TIc), NA)

LIc <- cor.test(Correlations$LikeMeansNorm, Correlations$IRIF.z, method = c("spearman"))
ifelse(LIc$p.value < 0.01, print(LIc), NA)
```


### Plot Stat Significative Correlations

#### Like ~ Beautiful {.tabset}
```{r}
require(ggplot2)
LB <- ggplot(Correlations, aes(x=LikeMeansNorm, y=BeautifulMeansNorm)) +
  geom_point(size=2) +
  geom_smooth(method=lm, se=T, color="black") +
  ylim(0.3, 1.5) +
  xlim(0.5, 1.4) +
  labs(x= "Like (mean score)", y= "Beautiful (mean score)") +
  ggtitle("") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))
print(LB)
```

#### Like ~ Touch {.tabset}
```{r}
require(ggplot2)
LT <- ggplot(Correlations, aes(x=LikeMeansNorm, y=TouchMeansNorm)) +
  geom_point(size=2) +
  geom_smooth(method=lm, se=T, color="black") +
  ylim(0.3, 1.5) +
  xlim(0.5, 1.4) +
  labs(x= "Like (mean score)", y= "Touch (mean score)") +
  ggtitle("") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))
print(LT)
```

#### Beautiful ~ Touch {.tabset}
```{r}
require(ggplot2)
BT <- ggplot(Correlations, aes(x=BeautifulMeansNorm, y=TouchMeansNorm)) +
  geom_point(size=2) +
  geom_smooth(method=lm, se=T, color="black") +
  ylim(0.3, 1.5) +
  xlim(0.5, 1.4) +
  labs(x= "Beautiful (mean score)", y= "Touch (mean score)") +
  ggtitle("") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))
print(BT)
```

# Citations --------------------------------------

Print References

General Packages
```{r}
citation(package ='readr')
citation(package ='sjstats')
citation(package ='MuMIn')
citation(package ='apaTables')
citation(package ='dplyr')
```

Packages to cite for Statistics
```{r}
citation(package ='lme4')
citation(package ='ez')
citation(package ='emmeans')
citation(package = 'robustlmm')
citation(package = 'ordinal')
```

Packages to cite for Plotting
```{r}
citation(package = 'effects')
citation(package ='ggplot2')
```

