library(broom)
library(car)
library(data.table)
library(dplyr)
library(effectsize)
library(effects)
library(emmeans)
library(ggplot2)
library(ggpubr)
library(gghalves)
library(glmmTMB)
library(Hmisc)
library(lmerTest)
library(lsr)
library(MLMusingR)
library(patchwork)
library(plotly)
library(psych)
library(pwr)
library(pwrss)
library(readxl)
library(sjstats)
library(tidyr)
library(tidyverse)
#uploading data and changing type 
X1_8_Pilot_Analysis <- read_excel("C:/Users/lupac/Desktop/1-8 Pilot Analysis.xlsx")
X9_16_Pilot_Analysis <- read_excel("C:/Users/lupac/Desktop/9-16 Pilot Analysis.xlsx")
X17_24_Pilot_Analysis <- read_excel("C:/Users/lupac/Desktop/17-24 Pilot Analysis.xlsx")
X25_32_Pilot_Analysis <- read_excel("C:/Users/lupac/Desktop/25-32 Pilot Analysis.xlsx")
X33_40_Pilot_Analysis <- read_excel("C:/Users/lupac/Desktop/33-40 Pilot Analysis.xlsx")
newdata <- rbind.data.frame(X1_8_Pilot_Analysis,X9_16_Pilot_Analysis, X17_24_Pilot_Analysis, X25_32_Pilot_Analysis, X33_40_Pilot_Analysis)
newdata$ID <- as.factor(newdata$ID)
#transform in factors and change names
newdata$Type <- as.factor(newdata$Type)
newdata$ID <- as.factor(newdata$ID)
newdata$Odour <- as.factor(newdata$Odour)
newdata$Image <- as.factor(newdata$Image)
levels(newdata$Odour)[2] <- "Lavender"
levels(newdata$Odour)[2]
newdata$Type[newdata$Type == "LL"] <- "L"
#LEOS
AOL <- read_excel("C:/Users/lupac/Desktop/Analysis of Leos.xlsx")
TOG <- merge(AOL, newdata, by = c("ID", "Odour"))
AOL$Pleasantness <- as.numeric(AOL$Pleasantness)

#Updates on mixed-effect model
#model baseline
M0 <- lmer(Rating ~ Odour * Valence * Type + (1|ID), data = newdata, REML = FALSE)
#other models
M1 <- lmer(Rating ~ Odour * Valence * Type + (1|ID) + (1|Image), data = newdata, REML = FALSE)
M2 <- lmer(Rating ~ Odour * Valence * Type + (1 + Odour|ID) + (1|Image), data = newdata, REML = FALSE)
M3 <- lmer(Rating ~ Odour * Valence * Type + (1 + Odour|ID) + (1 + Valence|Image), data = newdata, REML = FALSE)

#check singularity
isSingular(M1)
isSingular(M2) #failed
isSingular(M3) #failed

#check best model 
AIC(M0, M1)

#final model check
M1 <- lmer(Rating ~ Odour * Valence * Type + (1|ID) + (1|Image), data = newdata, REML = TRUE)

#check assumption 
#residual normality
qqnorm(residuals(M1))
qqline(residuals(M1))
ks.test(residuals(M1), "pnorm")
hist(residuals(M1))
#homoscendeicity
plot(fitted(M1), residuals(M1))
plot(M1)

#run test
anova(M1)
summary(M1)
confint(M1, method = "profile")
#follow up test
emmeans(M1, pairwise ~ Odour, adjust = "Tukey")
emm_options(lmer.df = "kenward-roger")
emmeans(M1, pairwise ~ Odour | Type, adjust = "None", emm_options, lmer.df = "satterthwaite")
# Effect sizes
es <- eff_size(
  emmeans(M1, ~ Odour | Type),
  sigma = sigma(M1),
  edf = df.residual(M1)
)

es
emmeans(M1, pairwise ~ Odour | Type, adjust = "None")

#check other details of the model
varcomp <- as.data.frame(VarCorr(M1))
var_fixed <- var(as.vector(model.matrix(M1) %*% fixef(M1)))
var_random <- sum(varcomp$vcov)
var_resid <- attr(VarCorr(M1), "sc")^2

R2_marginal <- var_fixed / (var_fixed + var_random + var_resid)
R2_conditional <- (var_fixed + var_random) / (var_fixed + var_random + var_resid)

R2_marginal
R2_conditional

emm_options(lmerTest.limit = 8399)

pc <- pairs(
  emmeans(M1, ~ Valence | Type),
  infer = TRUE, adjust = "none"
)

pc
eta_squared(M1, partial = TRUE)

#updated ANCOVA model
#first check collinearity between emotions
leos <- AOL %>%
  select(Happy, Nostalgic, Disgusted, Sensual, Thirsty, Relaxed, Energised)
cor(leos, method = "pearson")
rcorr(as.matrix(leos), type = "pearson")

#run ancova since no multicollinarity issue
Mcov <- lmer(Rating ~ Odour * Valence * Type +
               Happy + Nostalgic + Disgusted + Relaxed + Sensual + Thirsty + Energised +
               (1|ID) + (1|Image),
             data = TOG,
             REML = FALSE)
isSingular(Mcov)
anova(Mcov)
summary(Mcov)
#effectsize
confint(Mcov, method = "profile")
#R2
varcomp <- as.data.frame(VarCorr(Mcov))
var_fixed <- var(as.vector(model.matrix(Mcov) %*% fixef(Mcov)))
var_random <- sum(varcomp$vcov)
var_resid <- attr(VarCorr(Mcov), "sc")^2

R2_marginal <- var_fixed / (var_fixed + var_random + var_resid)
R2_conditional <- (var_fixed + var_random) / (var_fixed + var_random + var_resid)

R2_marginal
R2_conditional
#eta
eta_squared(Mcov, partial = TRUE)
#post-hoc tests
emmeans(Mcov, pairwise ~ Valence | Type, adjust = "None")
# Effect sizes
es <- eff_size(
  emmeans(Mcov, ~ Valence | Type),
  sigma = sigma(Mcov),
  edf = df.residual(Mcov)
)

es

cor.test(AOL$Happy, AOL$Nostalgic, method = "spearman")

#Plots
#Happiness
meanT <- TOG |>
  group_by(ID, Odour, Happy)|>
  dplyr::summarise(mean = mean(Rating, na.rm=TRUE),
                   .groups = "drop") 

meanT$Odour[meanT$Odour == "Air"] <- "Clean Air"
levels(meanT$Odour)[1]

meanT$Odour <- as.factor(meanT$Odour)
ggplot(
  data = meanT, 
  aes(x = Happy, y = mean, color = Odour)) +
  geom_point(size = 3, alpha = 0.7, shape = 19) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Scatterplot with mean pictures' ratings and happiness by odours",
       x = "Happiness",
       y = "Mean of image pleasantness ratings") +
  scale_color_manual(values=c("Clean Air"="lightblue","Lavender"="purple","Lemon"="#E3B505", "Rose"="pink","Vanilla"="brown"))+ theme(
    plot.title = element_text(size = 13, face = "bold", hjust = 0.5))+
  theme_classic()

ggsave("SH", device = "pdf", width = 3000, height = 2300, units = c("px"))
##
#Nostalgia
meanT <- TOG |>
  group_by(ID, Odour, Nostalgic)|>
  dplyr::summarise(mean = mean(Rating, na.rm=TRUE),
                   .groups = "drop") 

meanT$Odour[meanT$Odour == "Air"] <- "Clean Air"
levels(meanT$Odour)[1]

meanT$Odour <- as.factor(meanT$Odour)
ggplot(
  data = meanT, 
  aes(x = Happy, y = mean, color = Odour)) +
  geom_point(size = 3, alpha = 0.7, shape = 19) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Scatterplot with mean pictures' ratings and nostalgia by odours",
       x = "Nostalgia",
       y = "Mean of image pleasantness ratings") +
  scale_color_manual(values=c("Clean Air"="lightblue","Lavender"="purple","Lemon"="#E3B505", "Rose"="pink","Vanilla"="brown"))+ theme(
    plot.title = element_text(size = 13, face = "bold", hjust = 0.5))+
  theme_classic()

ggsave("SN", device = "pdf", width = 3000, height = 2300, units = c("px"))
##
#Disgust
meanT <- TOG |>
  group_by(ID, Odour, Disgusted)|>
  dplyr::summarise(mean = mean(Rating, na.rm=TRUE),
                   .groups = "drop") 

meanT$Odour <- as.character(meanT$Odour)
levels(meanT$Odour)[1] <- "Clean Air"
levels(meanT$Odour)[1]

meanT$Odour <- as.factor(meanT$Odour)
ggplot(
  data = meanT, 
  aes(x = Disgusted, y = mean, color = Odour)) +
  geom_point(size = 3, alpha = 0.7, shape = 19) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Scatterplot with mean pictures' ratings and disgust by odours",
       x = "Disgust",
       y = "Mean of image pleasantness ratings") +
  scale_color_manual(values=c("Clean Air"="lightblue","Lavender"="purple","Lemon"="#E3B505", "Rose"="pink","Vanilla"="brown"))+ theme(
    plot.title = element_text(size = 13, face = "bold", hjust = 0.5))+
  theme_classic()
ggsave("SD", device = "pdf", width = 3000, height = 2300, units = c("px"))
##
#Thirstiness
meanT <- TOG |>
  group_by(ID, Odour, Thirsty)|>
  dplyr::summarise(mean = mean(Rating, na.rm=TRUE),
                   .groups = "drop") 

meanT$Odour[meanT$Odour == "Air"] <- "Clean Air"
levels(meanT$Odour)[1]

meanT$Odour <- as.factor(meanT$Odour)
ggplot(
  data = meanT, 
  aes(x = Thirsty, y = mean, color = Odour)) +
  geom_point(size = 3, alpha = 0.7, shape = 19) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Scatterplot with mean pictures' ratings and thirstiness by odours",
       x = "Thirstiness",
       y = "Mean of image pleasantness ratings") +
  scale_color_manual(values=c("Clean Air"="lightblue","Lavender"="purple","Lemon"="#E3B505", "Rose"="pink","Vanilla"="brown"))+ theme(
    plot.title = element_text(size = 13, face = "bold", hjust = 0.5))+
  theme_classic()

ggsave("ST", device = "pdf", width = 3000, height = 2300, units = c("px"))
##

#Sensitivity analysis
library(pwr)

pwr.anova.test(
  k         = 5,      # 5 odour conditions
  n         = 40,     # 40 participants (within-subjects, approximated as groups of 45)
  sig.level = 0.05,
  power     = 0.80
)
f <- 0.2472734
eta_p2 <- f^2 / (1 + f^2)
eta_p2

#Complementary analysis
# subsetting to take out lemon
pleasant_data <- subset(newdata, Odour %in% c("Lavender", "Rose", "Vanilla", "Air"))
pleasant_TOG <- subset(TOG, Odour %in% c("Lavender", "Rose", "Vanilla", "Air"))
#now run M1
M1 <- lmer(Rating ~ Odour * Valence * Type + (1|ID) + (1|Image), data = pleasant_data, REML = TRUE)
anova(M1)

emmeans(M1, pairwise ~ Odour)
emmeans(M1, pairwise ~ Type)
emmeans(M1, pairwise ~ Odour | Type, adjust = "Tukey")
emmeans(M1, pairwise ~ Valence | Type, adjust = "Tukey")

#now ancova
Mcov <- lmer(Rating ~ Odour * Valence * Type +
               Happy + Nostalgic + Disgusted + Relaxed + Sensual + Thirsty + Energised +
               (1|ID) + (1|Image),
             data = pleasant_TOG,
             REML = FALSE)
anova(Mcov)
summary(Mcov)
eta_squared(Mcov, partial = TRUE)