#' ---
#' lesson: Statistik3
#' knitr:
#'   opts_chunk: 
#'     collapse: false 
#' ---
#' 
#' 
#' # Statistik 3: Demo
#' 
#' [Demoscript herunterladen (.R)](../purl/Statistik3_Demo.R){.dld}
#' 
#' [Demoscript herunterladen (.qmd)](../purl/Statistik3_Demo.qmd){.dld}
#' 
#' 
#' ## Korrelation vs. Regression
#' 
## -----------------------------------------------------------------------------
library(readr)
library(ggplot2)

#' 
#' 
## -----------------------------------------------------------------------------
# Datensatz zum Einfluss von Stickstoffdepositionen auf den Pflanzenartenreichtum
df <- read_delim("datasets/stat/Nitrogen.csv", delim = ";")
summary(df)

# Plotten der Beziehung
ggplot(df, aes(N_deposition, Species_richness)) +
  geom_point()

#' 
## -----------------------------------------------------------------------------
# Pearson Korrelation
# zuerst Species_richness dann N_deposition
cor.test(df$Species_richness, df$N_deposition, method = "pearson")

# Pearson Korrelation
# zuerst N_deposition dann Species_richness  
cor.test(df$N_deposition, df$Species_richness, method = "pearson") 

# Rang-Korrelation Spearman
cor.test(df$Species_richness, df$N_deposition, method = "spearman")

# Rang-Korrelation Kendall
cor.test(df$Species_richness, df$N_deposition, method = "kendall")

# Jetzt als Regression

# zuerst Species_richness dann N_deposition
lm_1 <- lm(Species_richness ~ N_deposition, data = df) 
lm_1

# zuerst N_deposition dann Species_richness  
lm(N_deposition ~ Species_richness, data = df) 

summary(lm_1) 

# Signifikantes Ergebnis visualisieren
ggplot(df, aes(N_deposition, Species_richness)) +
  geom_point() +
  geom_abline(intercept = lm_1$coefficients["(Intercept)"], slope = lm_1$coefficients["N_deposition"], 
              color = "blue")

#' 
#' 
#' ## Einfache und Polynomische Regression
#' 
## -----------------------------------------------------------------------------
# Daten generieren 

# "pred" sei unsere unabhängige Variable
pred <- c(20, 19, 25, 10, 8, 15, 13, 18, 11, 14, 25, 39, 38, 28, 24) 

# "resp" sei unsere abhängige Variable
resp <- c(12, 15, 10, 7, 2, 10, 12, 11, 13, 10, 9, 2, 4, 7, 13) 

# Dataframe erstellen
df_2 <- data.frame(pred, resp) 

# Daten anschauen
ggplot(df_2, aes(pred, resp)) +
  geom_point()

#' 
#' 
## -----------------------------------------------------------------------------
# Modell definieren und anschauen

# Einfaches lineares Modell
lm_1 <- lm(resp ~ pred, data = df_2) 

# Modell anschauen
summary(lm_1) 

#' 
#' -> kein signifikanter Zusammenhang im einfachen linearen Modell und entsprechend 
#' kleines Bestimmtheitsmass (adj. R^2^ = 0.07)
#' 
#' 
## -----------------------------------------------------------------------------
# Polynomische Regression Modell definieren und anschauen 

# lineares Modell mit quadratischem Termsummary
lm_quad <- lm(resp ~ pred + I(pred^2), data = df_2) 

# lineares Modell mit quadratischem Term anschauen
summary(lm_quad) 

#' 
#' -> Signifikanter Zusammenhang und viel besseres Bestimmtheitsmass 
#' (adj. R^2^ = 0.60)
#' 
#' 
## -----------------------------------------------------------------------------
# Modelle darstellen

# Vorhersagen der Modelle generieren

# 100 x-Werte, mit denen man die Modelle "füttern" kann 
xv <- seq( min(df_2$pred), max(df_2$pred), length = 100) 

# Vorhersagen des linearen Modells für die y-Werte
y_lm_1 <- predict(lm_1, data.frame(pred = xv) ) 

# Vorhersagen des quadratischen Modells für die y-Werte
y_lm_quad <- predict(lm_quad, data.frame(pred = xv)) 

ModPred <- data.frame(xv, y_lm_1, y_lm_quad)

# Modellvorhersagen plotten
ggplot(df_2, aes(pred, resp)) +
  geom_point() +
  geom_line(data = ModPred, aes(xv, y_lm_1), 
            color = "red", linetype = "dashed") + 
  geom_line(data = ModPred, aes(xv, y_lm_quad), color = "blue")

# Mit der funktion geom_smooth welche predict im Hintergrund verwendet
ggplot(df_2, aes(pred, resp)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x, color = "red", linetype = "dashed",
              se = FALSE) +
  geom_smooth(method = "lm", formula = y ~ x + I(x^2), 
              color = "blue", se = FALSE)

#' 
## -----------------------------------------------------------------------------
# Residualplots
par(mfrow = c(2, 2))
plot(lm_1, main = "Lineares Modell")
plot(lm_quad, main = "Quadratisches  Modell")

#' 
#' -> Die Plots sehen beim Modell mit quadratischem Term besser aus
#' 
#' ## ANCOVA
#' 
#' Experiment zur Fruchtproduktion ("Fruit") von Ipomopsis sp. in Abhängigkeit von 
#' der Beweidung ("Grazing" mit 2 Levels: "Grazed", "Ungrazed") und korrigiert für 
#' die Pflanzengrösse vor der Beweidung (hier ausgedrückt als Durchmesser an der 
#' Spitze des Wurzelstock: "Root")
#' 
## -----------------------------------------------------------------------------
# Daten einlesen und anschauen
compensation <- read_delim("datasets/stat/ipomopsis.csv")
compensation$Grazing <- as.factor(compensation$Grazing)

str(compensation)
summary(compensation)

# Plotten der vollständigen Daten/Information
ggplot(compensation, aes(Root, Fruit, color = Grazing)) +
  geom_point()

#' 
#' -> Je grösser die Pflanze, desto grösser ihre Fruchtproduktion.
#' -> Die grössere Fruchtproduktion innerhalb der beweideten Gruppe scheint ein 
#' Resultat von unterschiedlichen Pflanzengrössen zwischen den Gruppen zu sein.
#' 
## -----------------------------------------------------------------------------
# Lineare Modelle definieren und anschauen

# Volles Modell mit Interaktion
aoc_1 <- lm(Fruit ~ Root * Grazing, data = compensation) 
summary.aov(aoc_1)

# Finales Modell ohne die (nicht signifikante) Interaktion
aoc_2 <- lm(Fruit ~ Grazing + Root, data = compensation) 
summary.aov(aoc_2) # ANOVA-Tabelle
summary(aoc_2) # Parameter-Tabelle

# Residualplots anschauen
par(mfrow = c(2, 2) ) # 2x2 Plots pro Grafik
plot(aoc_2)

#' 
#' -> Das ANCOVA-Modell widerspiegelt die Zusammenhänge wie sie aufgrund der 
#' grafisch dargestellten Daten zu vermuten sind gut. Die Residual-Plots zeigen 3 Ausreisser (Beobachtungen 27, 34 und 37), welche "aus der Reihe tanzen".
#' 
## -----------------------------------------------------------------------------
# Visualisierung
ggplot(compensation, aes(Root, Fruit, color = Grazing)) +
  geom_point() +
  stat_smooth(formula = y ~ x, method = "lm" ) +
  ylim(0, 125)