#Siin koodis on kahe erineva autori koos kokku pandud. Esimesest osast leiad koodi, millega on tehtud
#eelkatse analüüs ning mõisteliste tunnuste hindamise katse mudeldamise osa. Teisest faili osast leiad
#koodi, mille järgi on loodud mõisteliste tunnuste hindamise katsel põhinev korrelatsiooniplot, heatmap 
#ning klasteranalüüs.


library(ggplot2)
library(standardize)
library(tidyr)
library(reshape2)
library(dplyr)
library(ggpubr)
library(wesanderson)
library(qqplotr)
library(stringr)
library(lme4)
library(sjmisc) 
library(effects)
library(sjstats)
library(sjPlot)
library(redres)
library(olsrr)
library(performance)
library(MuMIn)
library(colorBlindness)
library(Hmisc)
library(corrplot)
library(ggcorrplot)
library(psych)
library(MuMIn)

#---EELKATSE---

hinnangud = read.csv2(choose.files()) #hinnangukatse.csv

#pikka formaati:
hinnangud_pikk=melt(hinnangud, id.vars=c("HINDAJA"))

#hindaja järgi standardiseerimine
hinnangud_pikk$skaleeritud = scale_by(value~HINDAJA, hinnangud_pikk)
standardiseeritud_keskmised=aggregate(hinnangud_pikk$skaleeritud, list(hinnangud_pikk$variable), FUN=mean)

#keskmine hinnang igale dimensioonile
standardiseeritud_keskmised=aggregate(hinnangud_pikk$skaleeritud, list(hinnangud_pikk$variable), FUN=mean)
write.csv2(standardiseeritud_keskmised, "standardiseeritud_keskmised.csv")

#hinnangu järgi järjekora
standardiseeritud_keskmised_jrk=standardiseeritud_keskmised[order(standardiseeritud_keskmised$x),]
write.csv2(standardiseeritud_keskmised_jrk, "standardiseeritud_järjestus.csv")

#---LAUSETE HINDAMISE KATSE---

#---data sisse + puhastus + pikka formaati---

seisma<-read.csv(choose.files(), header = T, encoding = "UTF-8") #hindamiskatse_toorandmed.csv
seisma_puhas<-seisma[,c(1,14:222)]
seisma_demograafia<-seisma[,c(1,8:11,223)]

seisma_puhas[seisma_puhas == "Ei ole seotud1"] <- 1
seisma_puhas[seisma_puhas == "On tugevalt seotud7"] <- 7

seisma_puhas<-seisma_puhas %>%
  select(matches(c("id","SQ"))) #võta ainult õiged tulbad
write.csv(seisma_puhas, "seisma_puhas.csv")

seisma_pikk<-melt(seisma_puhas, id.vars=c("id"))
library(data.table)
seisma_pikk$variable <- str_replace_all(seisma_pikk$variable, '.SQ001.', '') #puhasta data
detach("package:data.table", unload=TRUE)

seisma_pikk_eraldi <-
  tidyr::separate(
    data = seisma_pikk,
    col = variable,
    sep = "v",
    into = c("lause","tunnus"),
    remove = FALSE)   #olulised tunnused eraldi tulpadesse

seisma_pikk_eraldi$tunnus[seisma_pikk_eraldi$tunnus == "1"] <- "vert"   #nimeta mugavalt ümber
seisma_pikk_eraldi$tunnus[seisma_pikk_eraldi$tunnus == "2tasa"] <- "tasa"
seisma_pikk_eraldi$tunnus[seisma_pikk_eraldi$tunnus == "3kont"] <- "kont"
seisma_pikk_eraldi$tunnus[seisma_pikk_eraldi$tunnus == "4"] <- "voim"
seisma_pikk_eraldi$tunnus[seisma_pikk_eraldi$tunnus == "5liik"] <- "liik"
  
seisma <- subset(seisma_pikk_eraldi, select = -variable )
write.csv(seisma, "seisma_pikk_data.csv")

---#lausenumbrite parandus---
  
seisma$lause[seisma$lause == "L01"] <- "S01"
seisma$lause[seisma$lause == "L04"] <- "S11"
seisma$lause[seisma$lause == "L11"] <- "S15"
seisma$lause[seisma$lause == "L02"] <- "S04"
seisma$lause[seisma$lause == "L05"] <- "S21"
seisma$lause[seisma$lause == "L09"] <- "S20"
seisma$lause[seisma$lause == "L03"] <- "S03"
seisma$lause[seisma$lause == "L08"] <- "S07"
seisma$lause[seisma$lause == "L10"] <- "S17"
seisma$lause[seisma$lause == "L07"] <- "S02"
seisma$lause[seisma$lause == "L12"] <- "S16"
seisma$lause[seisma$lause == "L17"] <- "S08"
seisma$lause[seisma$lause == "L15"] <- "S12"
seisma$lause[seisma$lause == "L13"] <- "S13"
seisma$lause[seisma$lause == "L14"] <- "S09"
seisma$lause[seisma$lause == "L16"] <- "S14"
seisma$lause[seisma$lause == "L19"] <- "S06"
seisma$lause[seisma$lause == "L06"] <- "S18"
seisma$lause[seisma$lause == "L18"] <- "S05"
seisma$lause[seisma$lause == "L21"] <- "S19"
seisma$lause[seisma$lause == "L20"] <- "S10"


#lisa rühmade numbrid (meie klasterdus)
seisma<-seisma[, 1:4]
seisma <- seisma %>%
  mutate(rühm = case_when(lause == "S08" | lause == "S19" ~ "R1",
                          lause == "S17" | lause == "S03" | lause == "S12" | lause == "S06" ~ "R3",
                          lause == "S15" | lause == "S21" | lause == "S01" | lause == "S11" ~ "R4",
                          lause == "S05"| lause == "S18"|lause == "S09"| lause == "S04" | lause == "S07" ~ "R5", 
                          TRUE ~ "R2"))
write.csv(seisma, "seisma_meie_rühmadega.csv")

seisma$value=as.numeric(seisma$value)

#----DEMOGRAAFIA----

seisma_demograafia=rename(seisma_demograafia, keel=G03Q08)
seisma_demograafia=rename(seisma_demograafia, vanus=G03Q09)
seisma_demograafia=rename(seisma_demograafia, sugu=G03Q10)
seisma_demograafia=rename(seisma_demograafia, haridus=G03Q11)
seisma_demograafia=rename(seisma_demograafia, seade=G04Q14)

keskmine_vanus=mean(seisma_demograafia$vanus)
keskmine_vanus_med=median(seisma_demograafia$vanus)


count(seisma_demograafia, haridus)
count(seisma_demograafia, sugu)
count(seisma_demograafia, keel)
count(seisma_demograafia, seade)


#tegemise aeg
seisma_toor<-read.csv(choose.files(), h=T, encoding="UTF-8") #hinnangukatse_toorandmed.csv
ajatempel<-seisma_toor[,c(1,6:7)]

ajatempel <- ajatempel %>%
  mutate(kestus = difftime(startdate, datestamp), units = "mins")

ajatempel[c('kestus', 'min')] <- str_split_fixed(ajatempel$kestus, ' ', 2)
str(ajatempel)
ajatempel$kestus=as.numeric(ajatempel$kestus)
mean(ajatempel$kestus)

#ajatempel$startdate=as.POSIXct(ajatempel$startdate)
#ajatempel$datestamp=as.POSIXct(ajatempel$datestamp)
keskmine_kestus<-median(ajatempel$kestus)


#-----MUDELDAMINE------

#subsetin analüüsiks
vertikaalsus <- seisma[seisma$tunnus == "vert", ]
tasakaal <- seisma[seisma$tunnus == "tasa", ]
võimelisus <- seisma[seisma$tunnus == "voim", ]
liikumisvõime <- seisma[seisma$tunnus == "liik", ]
kontakt <- seisma[seisma$tunnus == "kont", ]

#boxplotid ja mudelid
#vertikaalsus
vertikaalsus.means<-aggregate(value ~ id + rühm, data = vertikaalsus, FUN= "mean" )
vertikaalsus.means$rühm<-as.character(vertikaalsus.means$rühm)

vertikaalsus_plot <- ggplot(vertikaalsus.means, aes(x = rühm, y=value,ylab="Hinnang", xlab="Rühm")) + geom_boxplot(outlier.shape = NA, aes(fill = rühm)) + theme_minimal() + theme(legend.position ="top") + 
  scale_fill_brewer(palette = "RdBu") + ggtitle("Vertikaalsus") + labs(x="Rühm", y="Hinnang")

vertikaalsus_plot
cvdPlot(vertikaalsus_plot)
#kui tahad individuaalseid observationeid juurde siis lisa  + geom_jitter(width = 0.2)

vertikaalsus.model<-lmer(value ~ rühm + (1|id) + (1|lause), data=vertikaalsus,REML=FALSE)
vertikaalsus.null<-lmer(value ~ (1|id) + (1|lause), data=vertikaalsus,REML=FALSE)#taking out the verb factor
anova(vertikaalsus.model,vertikaalsus.null) #Likelihood ratio test, comparing the full model to the model without the fixed effect. If the resulting p-value is significant, than chances the model would be the same with and without the fixed effect are very low. I.e. the effect is significant if the p-value is significant.

check_convergence(vertikaalsus.model,tolerance = 0.001)  #kui TRUE siis konvergeerub
check_convergence(vertikaalsus.null, tolerance = 0.001)

r.squaredGLMM(vertikaalsus.model) #R-ruudus mõõdikute saamiseks


#tasakaal
tasakaal.means=aggregate(value ~ id + rühm, data = tasakaal, FUN= "mean" )

tasakaal_plot <- ggplot(tasakaal.means, aes(x = rühm, y=value)) + geom_boxplot(outlier.shape = NA, aes(fill = rühm)) + theme_minimal() + theme(legend.position ="top") + 
  scale_fill_brewer(palette = "RdBu") + ggtitle("Tasakaal") + labs(x="Rühm", y="Hinnang")
tasakaal_plot

tasakaal.model=lmer(value ~ rühm + (1|id) + (1|lause), data=tasakaal,REML=FALSE)
tasakaal.null=lmer(value ~ (1|id) + (1|lause), data=tasakaal,REML=FALSE)#taking out the verb factor
anova(tasakaal.model,tasakaal.null) #Likelihood ratio test, comparing the full model to the model without the fixed effect. If the resulting p-value is significant, than chances the model would be the same with and without the fixed effect are very low. I.e. the effect is significant if the p-value is significant.

r.squaredGLMM(tasakaal.model)

check_convergence(tasakaal.model,tolerance = 0.001)
check_convergence(tasakaal.null, tolerance = 0.001)


#võimelisus
võimelisus.means=aggregate(value ~ id + rühm, data = võimelisus, FUN= "mean" )

võimelisus_plot <- ggplot(võimelisus.means, aes(x = rühm, y=value)) + geom_boxplot(outlier.shape = NA, aes(fill = rühm)) + theme_minimal() + theme(legend.position ="top") + 
  scale_fill_brewer(palette = "RdBu") + ggtitle("Võimelisus") + labs(x="Rühm", y="Hinnang")
võimelisus_plot

võimelisus.model=lmer(value ~ rühm + (1|id) + (1|lause), data=võimelisus,REML=FALSE)
võimelisus.null=lmer(value ~ (1|id) + (1|lause), data=võimelisus,REML=FALSE)#taking out the verb factor
anova(võimelisus.model,võimelisus.null) #Likelihood ratio test, comparing the full model to the model without the fixed effect. If the resulting p-value is significant, than chances the model would be the same with and without the fixed effect are very low. I.e. the effect is significant if the p-value is significant.

r.squaredGLMM(võimelisus.model)

check_convergence(võimelisus.model,tolerance = 0.001)
check_convergence(võimelisus.null, tolerance = 0.001)


#liikumisvõime
liikumisvõime.means=aggregate(value ~ id + rühm, data = liikumisvõime, FUN= "mean" )

liikumisvõime_plot <- ggplot(liikumisvõime.means, aes(x = rühm, y=value)) + geom_boxplot(outlier.shape = NA, aes(fill = rühm)) + theme_minimal() + theme(legend.position ="top") + 
  scale_fill_brewer(palette = "RdBu") + ggtitle("Liikumisvõime") + labs(x="Rühm", y="Hinnang")
liikumisvõime_plot

liikumisvõime.model=lmer(value ~ rühm + (1|id) + (1|lause), data=liikumisvõime,REML=FALSE)
liikumisvõime.null=lmer(value ~ (1|id) + (1|lause), data=liikumisvõime,REML=FALSE)#taking out the verb factor
anova(liikumisvõime.model,liikumisvõime.null) #Likelihood ratio test, comparing the full model to the model without the fixed effect. If the resulting p-value is significant, than chances the model would be the same with and without the fixed effect are very low. I.e. the effect is significant if the p-value is significant.

r.squaredGLMM(liikumisvõime.model)

check_convergence(liikumisvõime.model,tolerance = 0.001)
check_convergence(liikumisvõime.null, tolerance = 0.001)


#kontakt
kontakt.means=aggregate(value ~ id + rühm, data = kontakt, FUN= "mean" )

kontakt_plot <- ggplot(kontakt.means, aes(x = rühm, y=value)) + geom_boxplot(outlier.shape = NA, aes(fill = rühm)) + theme_minimal() + theme(legend.position ="top") + 
  scale_fill_brewer(palette = "RdBu") + ggtitle("Kontakt") + labs(x="Rühm", y="Hinnang")
kontakt_plot

kontakt.model=lmer(value ~ rühm + (1|id) + (1|lause), data=kontakt,REML=FALSE)
kontakt.null=lmer(value ~ (1|id) + (1|lause), data=kontakt,REML=FALSE)#taking out the verb factor
anova(kontakt.model,kontakt.null) #Likelihood ratio test, comparing the full model to the model without the fixed effect. If the resulting p-value is significant, than chances the model would be the same with and without the fixed effect are very low. I.e. the effect is significant if the p-value is significant.

r.squaredGLMM(kontakt.model)

check_convergence(kontakt.model,tolerance = 0.001)
check_convergence(kontakt.null, tolerance = 0.001)

#multiple boxplots on one

ggarrange(vertikaalsus_plot, tasakaal_plot, võimelisus_plot, liikumisvõime_plot, kontakt_plot, 
          ncol = 3, nrow = 2)


###---HEATMAP, KLASTERANALÜÜS, KORRELATSIOONIPLOT---

#########################################
### Moos ja buss: seisma (ESUKA 2023) ###
#########################################

# Artikkel "Kas moos ja buss seisavad endiselt? seisma-verbi pol?seemia ja seismise kehaline kogemus"
# Autorid Ann Veismann, Mariann Proos ja Piia Taremaa
# Ajakiri: ...

library(dplyr)
library(reshape2)
library(Hmisc)
library(standardize)

citation("pvclust")

seisma = read.delim("seisma_pikk300922.txt", header = T, encoding = "UTF-8")
seisma <- seisma %>% mutate_if(is.character,as.factor)
nrow(seisma)

seisma_pikk = seisma
seisma_pikk$skaleeritud = scale_by(value~ID, seisma_pikk)
#standardiseeritud_keskmised=aggregate(seisma_pikk$skaleeritud, list(seisma_pikk$lause, seisma_pikk$tunnus, seisma_pikk$r?hm), FUN=mean)

seisma = seisma_pikk

head(seisma)
str(seisma)
summary(seisma)
colnames(seisma)
nrow(seisma)
seisma$tunnus <- recode(seisma$tunnus, 'kont'='kontakt')
seisma$tunnus <- recode(seisma$tunnus, 'liik'='liikumisv?ime')
seisma$tunnus <- recode(seisma$tunnus, 'tasa'='tasakaal')
seisma$tunnus <- recode(seisma$tunnus, 'voim'='v?imelisus')
seisma$tunnus <- recode(seisma$tunnus, 'vert'='vertikaalsus')
#write.table(data_wide, file = "seisma_standardiseeritud_pikk.txt", row.names = FALSE, sep = "\t")

part <- group_by(seisma, lause, tunnus)
head(part)

(part_mean3 <- summarise(part, m_value = mean(skaleeritud)))
v = melt(part_mean3, id.vars=c("lause", "tunnus"))
head(v)
v2 = v
v2 = select(v, c(lause, tunnus, value)) # pikas formaadis andmed nii, nagu vaja
head(v2)

data_wide <- dcast(v2, lause ~ tunnus, value.var="value") # laias formaadis andmed nii, nagu vaja
head(data_wide)
nrow(data_wide)

write.table(data_wide, file = "seisma_vahefail.txt", row.names = FALSE, sep = "\t")
seisma_tabel = read.delim("seisma_vahefail.txt", header = TRUE, sep = "\t") 
seisma_tabel <- seisma_tabel %>% mutate_if(is.character,as.factor)
str(seisma_tabel)
nrow(seisma_tabel)


### HEATMAP: joonis 1 ###

# ?petus siit: https://rpubs.com/melike/heatmapTable

co=melt(seisma_tabel)
head(co)

library(ggplot2)
library(scales) # for muted function
ggplot(co, aes(x=factor(variable, level=c('vertikaalsus', 'tasakaal', 'kontakt', 'v?imelisus', 'liikumisv?ime')), y=lause)) + # x and y axes => Var1 and Var2
  geom_tile(aes(fill = value)) + # background colours are mapped according to the value column
  geom_text(aes(fill = value, label = round(value, 2))) + # write the values
  scale_fill_gradient(low = "white", high = "gray50") + # determine the colour
  theme(panel.grid.major.x=element_blank(), #no gridlines
        panel.grid.minor.x=element_blank(), 
        panel.grid.major.y=element_blank(), 
        panel.grid.minor.y=element_blank(),
        panel.background=element_rect(fill="white"), # background=white
        axis.text.x = element_text(angle=40, hjust = 1,vjust=1,size = 14,face = "bold"),
        plot.title = element_text(size=20,face="bold"),
        axis.text.y = element_text(size = 12,face = "bold")) + 
  ggtitle("") +
  theme(legend.title=element_text(face="bold", size=14)) + 
  scale_x_discrete(name="") +
  scale_y_discrete(name="", limits=rev) +
  labs(fill="Keskmised \nhinnangud")


### KLASTERDUS: joonis 2 ###

head(seisma_tabel)
colnames(seisma_tabel)
row.names(seisma_tabel) = seisma_tabel[, 1] # esimene, j?rjekorranumbritega veerg maha, et joonisele tuleks verbid, mitte numbrid
seisma_tabel = seisma_tabel[, -1]
head(seisma_tabel)
write.table(seisma_tabel, file = "vahefail.txt", row.names = T, sep = "\t")
seisma_tabel = read.delim("vahefail.txt", sep = "\t") 

# HCA selle j?rgi: https://corpling.hypotheses.org/2675

library(pvclust)
mat = as.data.frame.matrix(seisma_tabel)
head(mat)
mat = t(mat)
fit = pvclust(mat, method.hclust="ward.D", method.dist="euclidean") # oli "canberra"
fit
plot(fit, hang = -1, cex = 1.2, main="", lwd=1,lty=1, asp = 2:1, reverse = T)


### KORRELATSIOONID: joonis 3 ###

# Poollai andmestik (?hel real unikaalne kombinatsioon katseisikust ja lausest), kood Mariannilt:

xseisma = read.delim("seisma_pikk300922.txt", header = T, encoding = "UTF-8")
xseisma <- xseisma %>% mutate_if(is.character,as.factor)
xseisma$skaleeritud = scale_by(value~ID, xseisma)
xxseisma<-xseisma[,c(1,3:4,7)]
nrow(xxseisma)
head(xxseisma)

wide_xseisma <- xxseisma %>%
  pivot_wider(names_from = tunnus, values_from = skaleeritud)
head(wide_xseisma)
nrow(wide_xseisma)
xtunnuste_hinded <- wide_xseisma[,3:7]
xtunnuste_hinded<-xtunnuste_hinded %>% 
  rename(
    vertikaalsus = vert,
    tasakaal = tasa,
    kontakt = kont,
    v?imelisus = voim,
    liikumisv?ime = liik
  )

# ?petus siit: http://www.sthda.com/english/wiki/ggplot2-quick-correlation-matrix-heatmap-r-software-and-data-visualization

data_wide = xtunnuste_hinded

cormat <- round(cor(data_wide),2)
head(cormat)
melted_cormat <- melt(cormat)
head(melted_cormat)
# Get lower triangle of the correlation matrix
get_lower_tri<-function(cormat){
  cormat[upper.tri(cormat)] <- NA
  return(cormat)
}
# Reorder the correlation matrix
cormat <- reorder_cormat(cormat)
upper_tri <- get_upper_tri(cormat)
# Melt the correlation matrix
melted_cormat <- melt(upper_tri, na.rm = TRUE)
# Create a ggheatmap
(ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
    geom_tile(color = "white")+
    scale_fill_gradient2(low = "white", high = "gray50", mid = "gray95", 
                         midpoint = 0, limit = c(-1,1), space = "Lab", 
                         name="Korrelatsiooni-\nkordaja") +
    theme_minimal()+ # minimal theme
    theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                     size = 12, hjust = 1), 
          axis.text.y = element_text(vjust = 1, size = 12, hjust = 1))+
    coord_fixed()) + 
  geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank())