library(ggplot2)
library(gridExtra)

epg = read.table(file.path(pfadu, "epg.txt"))
amp = read.table(file.path(pfadu, "dbdauer.txt"))

source(file.path(pfadu, "ggplotRegression.R"))

#############################################################

head(epg)
# theme_set(theme_grey(base_size = 18))

#############################################################
# Abbildungen: entweder mit plot()
############################################################# entweder mit plot()

# 1 Spalte, 3 Reihen
par(mfrow = c(1, 3))

# F2: abhaengige Variable, COG = unabhaengige Variable
plot(F2 ~ COG, data = epg)

# F1: abhaengige Variable, COG = unabhaengige Variable
plot(F1 ~ COG, data = epg)

# 4 Vokale auswaehlen
temp = with(epg, V %in% c("i", "I", "E", "a"))

# F1 abhaengige Variable, SUM1278 unabhaengige Variable
plot(F1 ~ SUM1278, data = epg[temp, ])


# oder mit ggplot()


p1 = ggplot(epg) + 
  aes(y = F2, x = COG) + 
  geom_point()

p2 = ggplot(epg) + 
  aes(y = F1, x = COG) + 
  geom_point()

temp = with(epg, V %in% c("i", "I", "E", "a"))
p3 = ggplot(epg[temp, ]) + 
  aes(y = F1, x = SUM1278) + 
  geom_point()

# um die Bilder in 3 Spalten und einer Reihe zu arrangieren
grid.arrange(p1, p2, p3,  ncol = 3, nrow = 1)




#############################################################
# Kovarianz
#############################################################

# Kovarianz

# Abhaengige Variable
y = epg$F2

# Unabhaengige Variable
x = epg$COG

# Anzahl der Stichproben
n = length(y)
#oder
n = nrow(epg)

# Mittelwert von x
mx = mean(x)

# Mittelwert von y
my = mean(y)

# Stichproben-Abweichung vom Mittelwert
dx = x - mean(x)
dy = y - mean(y)

# Kovarianz = Produkt dieser Abweichungen
# dividiert durch den Stichprobenanzahl minus 1
covxy = sum(dx * dy) / (n - 1)
covxy
# das gleiche
cov(x, y)

#############################################################
# Korrelation
#############################################################

xgross = x * 1000

cov(x, y)
cov(xgross, y)

# Korrelation: die Kovarianz dividiert durch
# den Produkt der Standardabweichungen der beiden Variablen
# normalisiert (r zwischen -1 und 1)!
r = cov(x, y) / (sd(x) * sd(y))

# Mit der cor() Funktion:
cor(x, y)
cor(xgross, y)

#############################################################
# Regression
#############################################################

########################## Die Steigung: entweder
b = r * sd(y) / sd(x)
b
# oder
b = cov(x, y) / var(x)
b
########################## Das Intercept
k = my - b * mx
k
########################## Die eingeschaetzen Werte
yhut = b * x + k
yhut
########################## Abbildung mit ueberlagerter Regressionsline
# entweder:
par(mfrow = c(1, 1))
plot(y ~ x)
# Regressionslinie ueberlagern
abline(k, b)
# die eingeschaetzen Werte ueberlagern
points(x, yhut, col = 2)

# oder mit ggplot() (um Einiges komplizierter!)
# Data-frame
df = data.frame(y, x, yhut)
# Die Abbildung
p1 = ggplot(df) + aes(y = y, x = x) + geom_point()
# Eingeschaetze Werte
p2 = geom_point(aes(y = yhut), color = "red")
# Regressionsline
p3 = geom_abline(intercept = k,
                 slope = b,
                 col = "red")
# Alle zusammen
p1 + p2 + p3



############################ SSE
error =  y - yhut
SSE = sum(error ^ 2)

#############################################################
# Regression und die lm() Funktion
#############################################################

########################### Regression berechnen
reg = lm(y ~ x)
########################### Regression ueberlagern
plot(y ~ x)
abline(reg)


# oder mit der Funktion ggplotRegression()
# (Quelle: https://rdrr.io/github/G-Thomson/gthor/man/ggplotRegression.html)
ggplotRegression(reg)

########################### Steigung, Intercept
coef(reg)
########################### die eingeschaetzten Werte
yhut = predict(reg)

# Der Error (Abstand zwischen den tatsaechlichen und eingeschaetzen Werten)
residuals(reg)
# oder auch
resid(reg)
# SSE
deviance(reg)
# das gleiche (siehe oben)
sum(error ^ 2)

#############################################################
# SSY, SSR, SSE
#############################################################
# SSY, SSR, SSE
SSY = sum((y - my) ^ 2)
SSR = sum((yhut - my) ^ 2)
# bestaetigen, dass SSR + SSE = SSY
SSR + SSE

# R-squared
SSR / SSY
# das gleiche
cor(x, y) ^ 2

#############################################################
# Die Pruefstatistik
#############################################################
# Der Standard-Error von 'r'
rsb = sqrt((1 - r ^ 2) / (n - 2))
# Die t-Statistik
tstat = r / rsb
# Die Wahrscheinlichkeit, dass die Werte durch eine Regressionslinie
# modelliert werden koennen
2 * (1 - pt(tstat, n - 2))
# Das gleiche mit der Die F-Statistik
fstat = tstat ^ 2
1 - pf(fstat, 1, n - 2)

# Diese Informationen sind auch in summary() enthalten:
summary(reg)
# Residual standard error: 300 on 43 degrees of freedom
# Multiple R-squared:  0.7952,	Adjusted R-squared:  0.7905
# F-statistic:   167 on 1 and 43 DF,  p-value: < 2.2e-16
#
# Es gibt eine signifikante lineare Beziehung zwischen
# COG und F2 (R^2 = 0.80, F[1, 43] = 167.0, p < 0.001).


#############################################################
# Kriterien fuer die Durchfuehrung einer Regression
#############################################################
# Die Residuals:
# sollen einer Normalverteilung folgen:
shapiro.test(resid(reg))

#	(a). Shapiro-Wilk normality test
# data:  resid(reg)
# Wenn p > 0.05, ist es OK also die Residuals sind mit einer
# Normalverteilung konsistent
# W = 0.9704, p-value = 0.2987

#   (b). Residuals beobachten
# Die Residuals sollen auf eine randomisierte Weise
# um die Null-Linie verteilt sein - dies ist eventuell nicht der Fall
plot(resid(reg))
abline(h = 0, lty = 2)

#   (c). Keine Autokorrelation
# Vor allem sollen Werte bei lag 1 und lag 2 innerhalb des
# Konfidenzintervalls liegen. Dies ist nicht der Fall...
acf(resid(reg))