library(ggplot2)
library(dplyr)
library(ez)

###########################################################################
# 1. 'Uebungen'
###########################################################################

# 1.
# Fuer den in R vorhandenen Data-Frame 'trees' pruefen Sie
# inwiefern 'Height' aus 'Volume'  vorhersagbar ist. 
# Schaetzen Sie 'Height' ein bei einem Volumen von 110.

head(trees)

plot(Height ~ Volume, data = trees)

# Regression
trees.lm = lm(Height ~ Volume, data = trees)
summary(trees.lm)
# Es gibt eine signifikante lineare Beziehung zwischen
# Volume und Height (R^2 = 0.36, F[1,29] = 16.2, p < 0.001)

# Bei lm() koennen wir die Voraussetzungen erst nach
# erfolgter Modellierung pruefen:
# Weichen die Residuals von einer Normalverteilung ab?
shapiro.test(resid(trees.lm))
# Nein
# W = 0.9822, p-value = 0.8707

# Gleichmaessige Verteilung um die 0-Linie?
# Mehr oder weniger, ja.
plot(resid(trees.lm))

# Autokorrelation?
# Nein - die meisten Werte - und vor allem den zweiten Wert -
# liegen innherhalb den blauen Linien
acf(resid(trees.lm))

# Daher koennen wir das Problem loesen
# Abbildung mit Regressionslinie
plot(Height ~ Volume, data = trees)
abline(trees.lm)

# Alternativ mit ggplot()
coef = coef(trees.lm)
p1 = ggplot(trees.lm) + aes(y = Height, x = Volume) + geom_point()
p2 = geom_abline(intercept = coef[1], slope = coef[2])
p1 + p2

#Das Folgende ist eine Wiederholung (da wir jetzt wissen, dass wir das machen duerfen)
# Statistik
summary(trees.lm)

# Es gibt eine signifikante lineare Beziehung zwischen Height und Volume
# (R^2 = 0.36, F[1,29] = 16.2, p < 0.001).

# Die eingeschaetzte Hoehe bei einem Volumen von 110:
predict(trees.lm, data.frame(Volume = 110))

#oder die Formel SLOPEWERT * x + INTERCEPTWERT (wobei hier x=110), also hier:
0.2319 * 110 + 69.0034

# ggf. Bild neu malen mit diesem Wert
xlim = c(10, 110)
ylim = c(60, 100)
plot(Height ~ Volume, data = trees, xlim=xlim, ylim = ylim)
abline(trees.lm)
points(100, 92.19335, col = 2)
abline(v = 100, h = 92.19335, lty=2, col=2)



# 2. 
# Fuer diese Daten wurde F2 - F1 (Hz) in einem Vokal
# zwischen 1910 und 1997 gemessen. aendert sich F2-F1 mit der Zeit?
# Wenn ja, schaetzen Sie den Wert von F2-F1 ein im Jahr 2012.

jahr = c(1910, 1920, 1930, 1940, 1950, 1959, 1969, 1978, 1987, 1997)
fdat = c(139, 149, 157, 175, 216, 303, 390, 449, 462, 487)

plot(fdat~jahr)
abline(lm(fdat~jahr))
shapiro.test(resid(lm(fdat~jahr)))

plot(resid(lm(fdat~jahr)))
abline(h = 0)

acf(resid(lm(fdat~jahr)))

#Okay, daher:
summary(lm(fdat~jahr))
# Der Wert F2-F1 aendert sich in Abhaengigkeit von
# der Zeit signifikant (R^2=0.93, F[1,8]=107.2, p < 0.001)

predict(lm(fdat~jahr),data.frame(jahr=2012))
#oder m*x+k
coef(lm(fdat~jahr))[2]*2012+coef(lm(fdat~jahr))[1]
# 2012 soll laut Modell der Wert bei 566 Hz liegen. 
# (Hinterfragen Sie jedoch die Sinnhaftigkeit solcher Vorhersagen!)

# 3. 
# Die Grundfrequenz wurde in der selben Person
# in einem Zeitraum von 10 Jahren gemessen (d.h. eine Messung pro Jahr).
# Der erste Werte ist aus den Jahr 1987, der letzte aus dem Jahr 1996:
# 137.0   131.2   127.1   123.4   119.2   114.6   109.6   104.5   99.4    95.3
# Aendert sich die Grundfrequenz mit der Zeit?
# Wenn ja, welchen Wert muesste f0 im Jahr 2000 gehabt haben?

grund = c(137, 131.2, 127.1, 123.4, 119.2, 114.6, 109.6, 104.5, 99.4, 95.3)
jahr = 1987:1996
plot(grund~jahr)
abline(lm(grund~jahr))
shapiro.test(resid(lm(grund~jahr)))
plot(resid(lm(grund~jahr)))
acf(resid(lm(grund~jahr)))
summary(lm(grund~jahr))
# Die Grundfrequenz sank signifikant mit den Jahren
# (R^2=0.99 (sic!, da die Daten fast auf der Linie liegen), 
# F[1,8]=4276, p <0.001)

# Gemaess der Vorhersage des Modells haette die Grundfrequenz 
# im Jahre 2000 bei einem Werte von
predict(lm(grund~jahr),data.frame(jahr=2000))
# ca. 77 Hz liegen muessen. 

#Waere auch durch m*x+k zu ermitteln gewesen:
coef(lm(grund~jahr))
# -->
-4.587273 * 2000 + 9251.683636

# Man beachte auch nochmal das Intercept: da x zwischen 1987 und 1996 schwankt,
# ist das Intercept der Wert, denn das Modell für das Jahr 0 (!!) vorhersagen würde.
coef(lm(grund~jahr))[1]
predict(lm(grund~jahr),data.frame(jahr=0))

# Die Versuchsperson lebte damals aber noch nicht, noch hat Sie damals eine Grundfrequenz von 9252 Hz produziert...