-------------------------
#  Regression-Examples  #
-------------------------

x = seq(from=-5,to=5,by=0.5)
x
n = length(x)
n
noise = rnorm(n)
y = 3 - 0.5 * x + noise           # Gleichung (1)
plot(x,y)

# Regression: Nehmen wir an, die Daten x und y sind gegeben, 
# wir kennen aber nicht den genauen Zusammenhang zwischen 
# x und y. Man schaut sich den plot an und vermutet einen 
# linearen Zusammenhang:
#
#   y = beta0 + beta1 * x           Gleichung (2)
#
# Die lm()-Funktion bestimmt diejenigen beta0 und beta1, 
# die am besten, sum_i (y_i^observed - y_i^model)^2 -> min, 
# zu den gegebenen Daten passt. Dabei ist hier y^observed 
# durch (1) gegeben und y^model durch Gleichung (2).

lm( y ~ x )

# der Rueckgabetyp der lm()-Funktion ist eine Liste mit 
# 12 Elementen, die default-maessig nicht alle angezeigt 
# werden.

res = lm( y ~ x )
res

summary(res)

names(res)

# structural information:
mode(res)
str(res)
class(res)

# list elements:
res[[1]]       # Listen-Elemente werden mit doppelten 
               # eckigen Klammern aufgerufen
res[1]         # ok, scheint hier auch zu gehen..
names(res)
res$coefficients
res$coef
res$c          # es gibt noch ein res$call
res$co
res$ca
res$call

# Regression-Fit:
res$fit
plot(x,y)
lines(x,res$fit,col="red")


# Ueberpruefen wir die Formeln aus der Vorlesung:
# Die beta's  beta=(beta0,beta1)  waren gegeben 
# durch die Loesung des Gleichungssystems 
#
# A beta = b
#
# mit 
#
# a11 = e^2,  a12 = e*x
# a21 = x*e,  a22 = x^2
#
# und
#
# b1 = e*y,  b2 = x*y
#
# checken wir das:

n     # 21
e = rep(1,n)
e

ee = sum(e*e)
ex = sum(e*x)
xx = sum(x*x)
ey = sum(e*y)
xy = sum(x*y)

a1 = c(ee,ex)  # erste  Zeile der Matrix A
a2 = c(ex,xx)  # zweite Zeile der Matrix A 

A = rbind(a1,a2)      # rbind fuer rowbind, zeilenweises 
                      # Zusammensetzen von A
A

b = c(ey,xy)

# Loesen von A beta = b :

beta = solve(A,b)

beta
res$coeff      # ok, das passt


# Berechnen wir noch den Regression-Fit, die 
# Regressionsgerade, von Hand und checken 
# gegen die res$fit-Zahlen: 

beta0 = beta[1]
beta1 = beta[2]
beta0
beta1

yfit = beta0 + beta1*x

plot(x,y)
lines(x,res$fit,col="red")
lines(x,yfit, col="green")    # ok, das passt auch



----------------
#  Beispiel 2  #
----------------
rm(list=ls())                 # loescht alles

x = seq(from=-5,to=5,by=0.5)
n = length(x)
n
noise = rnorm(n)
y = 3 - 0.5 * x + 0.25 * x^2 + noise         
plot(x,y)

# durch Betrachten des plots vermutet man 
# einen quadratischen Zusammenhang:
#
# y = beta0 + beta1 * x + beta2 * x^2         (3)
#
# Das Modell (3) ist nicht mehr linear in x, 
# aber immer noch linear in den gesuchten 
# Modell-Parametern beta0, beta1, beta2 -> 
# wir koennen die betas mit linearer Regression, 
# mit der lm()-Funktion, bestimmen. Wir haben 
# jetzt 3 Regressionskoeffizienten und damit 
# auch 3 Regressoren: den konstanten Vektor 
# e=(1,1,...,1) fuer das beta0 und dann die 
# Vektoren x und x^2 fuer beta1 und beta2. 
# Man muss alle Regressoren bis auf das e 
# in eine Regressionsmatrix, sagen wir, X 
# packen:

X = cbind(x,x^2)     # cbind fuer "columnbind", 
                     # spaltenweises Zusammensetzen 
                     # von X
X

# wir koennen die 2.Spalte noch "x^2" nennen (das muss 
# man aber nicht machen):
colnames(X)
colnames(X)[2] = "x^2"
colnames(X)
X

# die lineare Regression kann man dann ganz einfach  
# so machen:

res = lm( y ~ X )
res


# die lm()-Funktion nimmt automatisch immer einen
# konstanten Term mit, sie produziert automatisch 
# immer ein beta0 (auch "intercept" genannt); moechte 
# man das unterdruecken, also ein Modell 
#
# y = beta1 * x + beta2 * x^2                 (4)
#
# fitten, dann muss man folgende Syntax benutzen:

res2 = lm( y ~ -1 + X )
res2

# In diesem Fall kommt das natuerlich nicht so gut 
# hin:

plot( x , y )
lines( x , res$fit , col="red" )          # mit  beta0
lines( x , res2$fit , col="green" )       # ohne beta0

# vielleicht noch die y-Achse etwas rescalieren:

plot( x , y , ylim = c(0,14) )
lines( x , res$fit , col="red" )          # mit  beta0
lines( x , res2$fit , col="green" )       # ohne beta0

# fuegen wir noch eine x-Achse hinzu:

lines( x , rep(0,length(x)) )

# also res2$fit an der Stelle x=0 ist 0, wie es sein 
# muss.