### STATISTICNA ANALIZA PODATKOV MERITEV JECMENA ###
### AG3, 2026-2027 ###

## POLJSKI POSKUS VKLJUCUJE: 3 OBRAVNAVANJA X 3 PONOVITVE (oz. BLOKOV) ##
## 10 MERITEV DOLZINE JECMENOVEGA LISTA ZNOTRAJ POSAMEZNEGA BLOKA ##

## 1_KORAK
# NASTAVITE POT DO MAPE (set working directory) S PODATKE V KATERO SE VAM TEKOM ANALIZE SHRANJUJEJO GRAFI, PREGLEDNICE, IPD
setwd("E:/Rvaje_jecmen26")
getwd() 
# nalozi knjiznice s paketi ki jih potrebuete pri nalaganju, analizi in prikazu podatkov
library(readr)
library(tidyverse)
library(dplyr)
library(stringr)
library(ggplot2)
library(lme4)
library(lmtest)
library(emmeans)
library(multcomp)
library(multcompView)

## 2_KORAK
# NALOZI PODATKE MERITEV JECMENA
jecmen <- read.csv(
  file = "data_jecmen.csv", 
  header = TRUE,
  na.strings = "",
  stringsAsFactors = TRUE
)

## 3_KORAK
# ustvarite povzetek po obravnavanju po posameznem bloku (9 obravnavanj)
povzetek_blok <- jecmen %>%
  group_by(obravnava, blok) %>%
  summarise(
    n = dplyr::n(),
    mean_cm = mean(dolz_lista_cm, na.rm = TRUE),
    median_cm = median(dolz_lista_cm, na.rm = TRUE),
    sd_cm = sd(dolz_lista_cm, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(obravnava, blok)
povzetek_blok

## 4_KORAK 
# ustvarite povzetek po obravnavanju (3 obravnavanja, v vsakem 3 ponovitve)
jecmen <- jecmen %>%
  mutate(
    treatment = factor(
      obravnava,
      levels = c("KONTROLA", "M60", "S600")),
    Rep = factor(blok)
  )
povz_obravn <- jecmen %>%
  group_by(treatment) %>%
  summarise(
    n = dplyr::n(),
    mean_cm = mean (dolz_lista_cm, na.rm = TRUE),
    median_cm = median(dolz_lista_cm, na.rm = TRUE),
    sd_cm = sd(dolz_lista_cm, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(treatment)
povz_obravn

blok_means <- jecmen %>%
  group_by(treatment, Rep) %>%
  summarise(mean_cm = mean(dolz_lista_cm, na.rm = TRUE), .groups = "drop")
jecmen %>% count(treatment, Rep)  #da preverimo da je res okoli 10 v vsaki celici
blok_means

## 5_KORAK
# narisi graficni prikaz s povprecji treh obravnavanj
osnoven_graf <- ggplot(
  jecmen,
  aes(x = treatment, y = dolz_lista_cm)) +
  geom_boxplot(fill = "lightgray", color = "black", width = 0.6) +
  labs(
    title = "Dolzina lista po obravnavanjih",
    x = "Obravnavanje",
    y = "Dolzina lista [cm]"
  ) +
  theme_minimal()
osnoven_graf
ggsave("osnoven_graf.tiff", osnoven_graf, device = "tiff", width = 15, height = 10, units = "cm", dpi = 600)

# oziroma ce zelim polepsati graficni prikaz uporabim spodnje ukaze
lepsi_graf <- ggplot(
  jecmen,
  aes(x = treatment, y = dolz_lista_cm)) +
  geom_boxplot(
    fill = "gray90", 
    color = "black", 
    width = 0.6,
    lwd = 0.7,           #thicker box and whisker lines
    outlier.shape = 21,
    outlier.size = 2,
    outlier.fill = "red",
    outlier.alpha = 0.5
  ) +
  stat_summary(
    fun = median, 
    geom = "crossbar",
    width = 0.6,
    color = "red",
    fatten = 0.8
  ) +
  stat_summary (
    fun.min = min, 
    fun.max = max,
    geom = "errorbar",
    color = "black",
    width = 0.25,
    linewidth = 0.8
  ) +
  labs(
    title = "Dolzina lista po obravnavanjih",
    x = "Obravnavanje",
    y = "Dolzina lista [cm]"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 13, hjust = 0.5),
    axis.title.x = element_text(face = "bold", size = 12),
    axis.title.y = element_text(face = "bold", size = 12),
    panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.8),
    axis.line = element_line(color = "black"),
    axis.ticks = element_line(color = "black"),
    axis.ticks.length = unit(2, "mm")
  )
lepsi_graf
ggsave("lepsi_graf.tiff", lepsi_graf, device = "tiff", width = 15, height = 10, units = "cm", dpi = 600)

## 6_KORAK
# STATISTICNA ANALIZA MED OBRAVNAVANJI - ANOVA
AnovaModel_1 <- aov(mean_cm ~ treatment, data = blok_means)
summary(AnovaModel_1)

##7_KORAK
#PREVERJANJE PREDPOSTAVK ANOVE
#"normalnost" ostankov
shapiro.test(residuals(AnovaModel_1))
#diagnosticni grafi
par(mfrow = c(2, 2))
plot(AnovaModel_1)
par(mfrow = c(1, 1))
#homogenost varianc med obravnavanji
fligner.test(mean_cm ~ treatment, data = blok_means)