Regression mit R

Mit R markdown lassen sich R-Code und Plots in einem HTML Dokument anzeigen

wahlen <- read.csv2(file.choose(), header= T, sep= ";", dec= ",")
attach(wahlen)
names(wahlen)
## [1] "Stadtbezirk"      "Wahlbeteiligung"  "stark"           
## [4] "schwach"          "Arbeitslosigkeit" "Matura"          
## [7] "ohneAbschluss"
plot(schwach, Wahlbeteiligung)

cor(schwach, Wahlbeteiligung)
## [1] -0.8777211
plot(stark, Wahlbeteiligung)

cor(stark, Wahlbeteiligung)
## [1] 0.749697
# Anteil sozial starker Bewohner und Wahlbeteiligung
plot(Wahlbeteiligung ~ stark,  pch=19,
     xlab="Anteil stärkere Milieus", 
     ylab="Wahlbeteiligung", 
     col= "grey50")
mod <- lm(Wahlbeteiligung ~ stark)
abline(mod, lwd=2)

coef(mod)
## (Intercept)       stark 
##  67.3228377   0.2046428
fit <- data.frame(Wahlbeteiligung, yhat=fitted(mod), e=resid(mod))
round(fit, 4)
##    Wahlbeteiligung    yhat       e
## 1             65.1 70.4334 -5.3334
## 2             67.4 69.6762 -2.2762
## 3             68.5 69.2465 -0.7465
## 4             69.4 70.8836 -1.4836
## 5             69.9 71.6817 -1.7817
## 6             73.0 73.1552 -0.1552
## 7             73.6 71.2111  2.3889
## 8             74.2 69.5125  4.6875
## 9             74.4 70.7813  3.6187
## 10            75.3 73.7486  1.5514
## 11            76.6 78.2508 -1.6508
## 12            79.8 78.6191  1.1809
# ausgefeiltere Form der Darstellung
x <- stark
y <- Wahlbeteiligung
mittelw.x <- mean(x)    # arithmetisches Mittel x
mittelw.y <- mean(y)    # arithmetisches Mittel y
regression <- lm(y~x)   # Regression von y auf x
reg.param <- regression$coefficients # Vektor mit den Koeffizienten
reg.param    # Zuerst der Ordinatenabschnitt (englisch: Intercept), dann die Steigung
## (Intercept)           x 
##  67.3228377   0.2046428
#
#
# Grafikausgabe
titel <- bquote(hat(y)==.(reg.param[1])+.(reg.param[2])*x) # bquote() Setzt hier die Regressionsgleichung in die �berschrift, 
# hat(y) setzt ein Dach (circumflex) auf das y
x.achse <- "Anteil starker Milieus"
y.achse <- "Wahlbeteiligung"
plot(x, y, xlab=x.achse, ylab=y.achse, main=titel)
abline(regression, col=2)
abline(h=mittelw.y, v=mittelw.x, lty=2)     # Schwerpunkt der Punktwolke    

#Textausgabe
regression2 <- summary(regression)
print(regression2)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3334 -1.6835 -0.4508  1.7608  4.6875 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 67.32284    1.62128  41.524 1.57e-12 ***
## x            0.20464    0.05712   3.582  0.00499 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.948 on 10 degrees of freedom
## Multiple R-squared:  0.562,  Adjusted R-squared:  0.5183 
## F-statistic: 12.83 on 1 and 10 DF,  p-value: 0.004992
# Alle Berechnungen Schritt für Schritt
# Hilfstabelle
i <- 1:length(x)
x2 <- x*x
xy <- x*y
y2 <- y*y
cbind(i, x, y, x2, xy, y2)
##        i    x    y      x2      xy      y2
##  [1,]  1 15.2 65.1  231.04  989.52 4238.01
##  [2,]  2 11.5 67.4  132.25  775.10 4542.76
##  [3,]  3  9.4 68.5   88.36  643.90 4692.25
##  [4,]  4 17.4 69.4  302.76 1207.56 4816.36
##  [5,]  5 21.3 69.9  453.69 1488.87 4886.01
##  [6,]  6 28.5 73.0  812.25 2080.50 5329.00
##  [7,]  7 19.0 73.6  361.00 1398.40 5416.96
##  [8,]  8 10.7 74.2  114.49  793.94 5505.64
##  [9,]  9 16.9 74.4  285.61 1257.36 5535.36
## [10,] 10 31.4 75.3  985.96 2364.42 5670.09
## [11,] 11 53.4 76.6 2851.56 4090.44 5867.56
## [12,] 12 55.2 79.8 3047.04 4404.96 6368.04
sumx <- sum(x)
sumy <- sum(y)
# Mittelwerte
xquer <- sumx/length(x)
yquer <- sumy/length(y)
# Varianzen
s2x <- 1/length(x) * sum(x2) - xquer**2
s2y <- 1/length(y) * sum(y2) - yquer**2
# Kovarianz
cxy <- 1/length(xy) * sum(xy) - xquer*yquer
# Steigung 
b <- cxy/s2x
a <- yquer - b*xquer
a
## [1] 67.32284
b
## [1] 0.2046428