Statistik mit R The R Bootcamp |
Am Ende des Practicals wirst du wissen…
Öffne dein TheRBootcamp
R project.
Öffne ein neues R Skript. Schreibe deinen Namen, das Datum und “Robuste Statistik Practical” als Kommentare an den Anfang des Skripts.
## NAME
## DATUM
## Robuste Statistik Practical
Speichere das neue Skript unter dem Namen
robuste_statistik_practical.R
im 2_Code
Ordner.
Lade die nötigen Pakete. Siehe unten.
# Lade die nötigen Pakete
library(tidyverse)
library(lubridate)
Attaching package: 'lubridate'
The following objects are masked from 'package:base':
date, intersect, setdiff, union
library(lme4)
Loading required package: Matrix
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
library(boot)
library(quantreg)
Loading required package: SparseM
Attaching package: 'SparseM'
The following object is masked from 'package:base':
backsolve
library(Rfit)
read_csv()
Funktion um
avocado
und avocado_cali
einzulesen.# Lade die Daten
avocado <- read_csv(file = "1_Data/avocado.csv")
avocado_cali <- read_csv(file = "1_Data/avocado_cali.csv")
Printe die Datensätze.
Verwende names(XX)
, summary(XX)
, und
View(XX)
um einen weiteren Überblick über die Daten zu
bekommen.
Wiederum, führe den Code unten aus um sicherzustellen, dass alle
character
Variablen als Faktoren vorliegen, was den
statistischen Modellen hilft kategoriale Variablen richtig zu
interpretieren.
# Konvertiere alle character zu factor
avocado <- avocado %>% mutate_if(is.character, factor)
avocado_cali <- avocado_cali %>% mutate_if(is.character, factor)
avocado_cali
Datensatzes, welcher
Verkaufsdaten für Avocados in Kalifornien über die letzten drei Jahre
enthält. Das Ziel ist es zunächst Verkaufsvolumen
über
Preis
vorherzusagen mit einer normalen Regression
(lm
) vorherzusagen.# Regression Verkausvolumen auf Preis
cali_lm <- lm(formula = YY ~ XX,
data = ZZ)
# Regression Verkausvolumen auf Preis
cali_lm <- lm(formula = Verkaufsvolumen ~ Preis,
data = avocado_cali)
summary()
des Modells. Ist
Preis
ein guter Prädiktor für Verkaufsvolumen
?
Achte auch auf R-squared
.summary(cali_lm)
Call:
lm(formula = Verkaufsvolumen ~ Preis, data = avocado_cali)
Residuals:
Min 1Q Median 3Q Max
-4692333 -1486900 351468 1342770 4276676
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11577703 371529 31.2 <2e-16 ***
Preis -6115691 256440 -23.9 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1840000 on 336 degrees of freedom
Multiple R-squared: 0.629, Adjusted R-squared: 0.628
F-statistic: 569 on 1 and 336 DF, p-value: <2e-16
Preis
ist ein signifikanter Prädiktor mit einem sehr
extremen t-Wert. Zudem weist R-squared
einen hohen
Wert von 62.8%
erklärte Varianz auf. Bedeutet dies nun,
dass das Modell, welches einen linearen Zusammenhang annimmt, gut den
Zusammenhang von Preis
und Verkaufsvolumen
beschreibt? Verwende den Code unten um die Daten und die Regression zu
plotten.# Plotte Modell und Daten
ggplot(avocado_cali, aes(Preis, Verkaufsvolumen)) +
geom_point() + geom_smooth(method = 'lm')
# Residualplot
ggplot(mapping = aes(fitted(cali_lm), resid(cali_lm))) + geom_point()
R-squared
Werts liegen offensichtliche
Verletzungen der Annahmen vor: Linearität, Normalität, und
Homoskedastizität scheinen nicht zu halten. Ein Grund für diese
Verletzungen sind fehlende Variablen. Ergänze das Modell um einen
weiteren Prädiktor, nämlich Typ
, welcher die Art der
Avocado kodiert und plotte erneut Residuen gegen gefittete Werte.
Speichere das Modell als cali_lm
.# Regression Verkausvolumen auf Preis & Typ
cali_lm <- lm(formula = XX ~ XX + XX, data = XX)
ggplot(mapping = aes(fitted(XX), resid(XX))) + geom_point()
# Regression Verkausvolumen auf Preis & Typ
cali_lm <- lm(Verkaufsvolumen ~ Preis + Typ, data = avocado_cali)
ggplot(mapping = aes(fitted(cali_lm), resid(cali_lm))) + geom_point()
Verkaufsvolumen
. Es gibt nun mehrere Wege damit umzugehen.
1) Transformation von Verkaufsvolumen
z.B. mit
log()
, 2) Verwendung der alternativen Variable
Verkaufsvolumen_index
, oder 3) die Modellierung mit einer
passenden Verteilung. Probiere alle drei aus. Beginne mit 1), der
Transformation. Überprüfe zunächst mit dem Template wie sich die Schiefe
darstellt und wie sie sich durch log()
verändert.# Histogramm Verkaufsvolumen
ggplot(avocado_cali, aes(Verkaufsvolumen)) + geom_histogram()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(avocado_cali, aes(log(Verkaufsvolumen))) + geom_histogram()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
log(Verkaufsvolumen)
als Kriterium und
kreiere den typischen Residualplot.# Regression Verkausvolumen auf Preis & Typ
cali_lm <- lm(formula = log(XX) ~ XX + XX, data = XX)
ggplot(mapping = aes(fitted(XX), resid(XX))) + geom_point()
# Regression Verkausvolumen auf Preis & Typ
cali_lm <- lm(log(Verkaufsvolumen) ~ Preis + Typ, data = avocado_cali)
ggplot(mapping = aes(fitted(cali_lm), resid(cali_lm))) + geom_point()
Verkaufsvolumen_index
(ohne log()
) und
überprüfe, ob das vielleicht die verbleibenden Probleme löst.# Regression Verkausvolumen auf Preis & Typ
cali_lm <- lm(formula = XX ~ XX + XX, data = XX)
ggplot(mapping = aes(fitted(XX), resid(XX))) + geom_point()
# Regression Verkausvolumen auf Preis & Typ
cali_lm <- lm(Verkaufsvolumen_index ~ Preis + Typ, data = avocado_cali)
ggplot(mapping = aes(fitted(cali_lm), resid(cali_lm))) + geom_point()
# Regression Verkausvolumen auf Preis & Typ
tapply(resid(cali_lm), fitted(cali_lm) > 7, sd)
FALSE TRUE
0.229 0.178
Verkaufszahlen
(einer Häufigkeitsvariable) als Kriterium.
Achtung: Die Warnungen könnt ihr ignorieren.# Regression Verkausvolumen auf Preis & Typ
cali_glm <- glm(formula = XX ~ XX + XX, data = XX, family = 'XX')
ggplot(mapping = aes(fitted(XX), resid(XX))) + geom_point()
# Regression Verkausvolumen auf Preis & Typ
cali_glm <- glm(Verkaufsvolumen ~ Preis + Typ, data = avocado_cali, family = 'poisson')
ggplot(mapping = aes(fitted(cali_glm), resid(cali_glm))) + geom_point()
Verkausvolumen
über die Jahre 2016 und 2017. Führe zunächst
den Code unten aus um zwei vergleichbare Datensätze für die beiden Jahre
zu kreieren.# Jahr 2016
avocado_cali_2016 <- avocado_cali %>%
filter(year(Datum) == 2016)
# Jajr 2017 mit Wochen gematched
avocado_cali_2017 <- avocado_cali %>%
filter(year(Datum) == 2017,
week(Datum) %in% week(avocado_cali_2016$Datum))
t.test()
das
Verkaufsvolumen
und den Verkaufsvolumen_index
zwischen den beiden Jahren. Setze dabei paired = TRUE
um zu
berücksichtigen, dass es sich hier um abhängige Daten handelt. Siehe
Template.# Verkaufsvolumen 2016 vs 2017
t.test(avocado_cali_2016$XX,
avocado_cali_2017$XX, paired = TRUE)
# Verkaufsvolumen-Index 2016 vs 2017
t.test(avocado_cali_2016$XX,
avocado_cali_2017$XX, paired = TRUE)
# Verkaufsvolumen 2016 vs 2017
t.test(avocado_cali_2016$Verkaufsvolumen,
avocado_cali_2017$Verkaufsvolumen, paired = TRUE)
Paired t-test
data: avocado_cali_2016$Verkaufsvolumen and avocado_cali_2017$Verkaufsvolumen
t = 2, df = 103, p-value = 0.05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3290 307840
sample estimates:
mean of the differences
152275
# Verkaufsvolumen-Index 2016 vs 2017
t.test(avocado_cali_2016$Verkaufsvolumen_index,
avocado_cali_2017$Verkaufsvolumen_index, paired = TRUE)
Paired t-test
data: avocado_cali_2016$Verkaufsvolumen_index and avocado_cali_2017$Verkaufsvolumen_index
t = 2, df = 103, p-value = 0.03
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.0056 0.0965
sample estimates:
mean of the differences
0.051
Verkaufsvolumen
zeigte keinen signifikanten
Unterschied, Verkaufsvolumen_index
dagegen schon. Dies ist
wiederum auf die Schiefe in Verkaufsvolumen
zurückzuführen.
Vergleiche nun wie sich ein robuster wilcox.test()
in den
beiden Fällen schlägt.# Verkaufsvolumen 2016 vs 2017
wilcox.test(avocado_cali_2016$XX,
avocado_cali_2017$XX, paired = TRUE)
# Verkaufsvolumen-Index 2016 vs 2017
wilcox.test(avocado_cali_2016$XX,
avocado_cali_2017$XX, paired = TRUE)
# Verkaufsvolumen 2016 vs 2017
wilcox.test(avocado_cali_2016$Verkaufsvolumen,
avocado_cali_2017$Verkaufsvolumen, paired = TRUE)
Wilcoxon signed rank test with continuity correction
data: avocado_cali_2016$Verkaufsvolumen and avocado_cali_2017$Verkaufsvolumen
V = 3491, p-value = 0.01
alternative hypothesis: true location shift is not equal to 0
# Verkaufsvolumen-Index 2016 vs 2017
wilcox.test(avocado_cali_2016$Verkaufsvolumen_index,
avocado_cali_2017$Verkaufsvolumen_index, paired = TRUE)
Wilcoxon signed rank test with continuity correction
data: avocado_cali_2016$Verkaufsvolumen_index and avocado_cali_2017$Verkaufsvolumen_index
V = 3474, p-value = 0.02
alternative hypothesis: true location shift is not equal to 0
# Verkaufsvolumen 2016 vs 2017
# Vorzeichentest Vorbereitung
signs <- avocado_cali_2017$XX > avocado_cali_2016$XX
n_plus <- sum(signs)
n <- length(signs)
# p-Wert für Vorzeichentest
pbinom(q = n_plus, size = n, prob = .5)
# Verkaufsvolumen-Index 2016 vs 2017
# Vorzeichentest Vorbereitung
signs <- avocado_cali_2017$XX > avocado_cali_2016$XX
n_plus <- sum(signs)
n <- length(signs)
# p-Wert für Vorzeichentest
pbinom(q = n_plus, size = n, prob = .5)
# Verkaufsvolumen 2016 vs 2017
# Vorzeichentest Vorbereitung
signs <- avocado_cali_2017$Verkaufsvolumen > avocado_cali_2016$Verkaufsvolumen
n_plus <- sum(signs)
n <- length(signs)
# p-Wert für Vorzeichentest
pbinom(q = n_plus, size = n, prob = .5)
[1] 0.0475
# Verkaufsvolumen-Index 2016 vs 2017
# Vorzeichentest Vorbereitung
signs <- avocado_cali_2017$Verkaufsvolumen_index > avocado_cali_2016$Verkaufsvolumen_index
n_plus <- sum(signs)
n <- length(signs)
# p-Wert für Vorzeichentest
pbinom(q = n_plus, size = n, prob = .5)
[1] 0.0475
Verkaufsvolumen
durch Preis
und
Typ
. Ein weiterer Weg gegenüber den bisher eingeschlagenen
um mit potentiellen Verletzungen der Regression umzugehen sind robuste
Regressionen. Rechne zunächst noch einmal das ursprüngliche Modell mit
einer normalen Regression.# Regression Verkausvolumen auf Preis
cali_lm <- lm(formula = YY ~ XX + XX,
data = ZZ)
# Regression Verkausvolumen auf Preis
cali_lm <- lm(formula = Verkaufsvolumen ~ Preis + Typ,
data = avocado_cali)
rq()
(quantreg
Paket) und eine Rang-basierte Regression mit
rfit()
(Rfit
Paket).# Quantil-basierte Regression
cali_rq <- rq(formula = YY ~ XX + XX,
data = ZZ)
# Rang-basierte Regression
cali_rfit <- rfit(formula = YY ~ XX + XX,
data = ZZ)
# Quantil-basierte Regression
cali_rq <- rq(formula = Verkaufsvolumen ~ Preis + Typ,
data = avocado_cali)
# Rang-basierte Regression
cali_rfit <- rfit(formula = Verkaufsvolumen ~ Preis + Typ,
data = avocado_cali)
summary()
der drei Modelle. Welche
Unterschiede und Gemeinsamkeiten fallen dir auf?# regular regression
summary(cali_lm)
Call:
lm(formula = Verkaufsvolumen ~ Preis + Typ, data = avocado_cali)
Residuals:
Min 1Q Median 3Q Max
-2029523 -338620 -65747 347099 4692280
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7418161 181045 40.97 < 2e-16 ***
Preis -1338574 155354 -8.62 2.8e-16 ***
Typorganic -5012181 121165 -41.37 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 745000 on 335 degrees of freedom
Multiple R-squared: 0.939, Adjusted R-squared: 0.939
F-statistic: 2.59e+03 on 2 and 335 DF, p-value: <2e-16
# quantile regression
summary(cali_rq)
Call: rq(formula = Verkaufsvolumen ~ Preis + Typ, data = avocado_cali)
tau: [1] 0.5
Coefficients:
coefficients lower bd upper bd
(Intercept) 6.06e+06 5.68e+06 1.80e+308
Preis -1.54e+05 -3.16e+05 -4.18e+04
Typorganic -5.64e+06 -6.39e+06 -5.41e+06
# rank-based regression
summary(cali_rfit)
Call:
rfit.default(formula = Verkaufsvolumen ~ Preis + Typ, data = avocado_cali)
Coefficients:
Estimate Std. Error t.value p.value
(Intercept) 6364261 79959 79.59 < 2e-16 ***
Preis -455900 68628 -6.64 1.2e-10 ***
Typorganic -5453944 53525 -101.89 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Multiple R-squared (Robust): 0.93
Reduction in Dispersion Test: 2231 p-value: 0
Verkaufsvolumen
auf Preis
und Typ
zu bestimmen. Dazu haben wir
für dich schon die bootstrap Funktion geschrieben die wir im nächsten
Schritt in die boot
-Funktion stecken können. Versuche
zunächst jedoch nachzuvollziehen was unsere boot_lm
Funktion tut.# Bootstrap Funktion für lm
boot_lm <- function(data, indices, formula){
data <- data[indices,] # Bootstrap sample
data_lm <- lm(formula = formula,
data = data)
coefficients(data_lm)
}
# Lasse Bootstrap laufen
cali_boot<- boot(data = XX, # Der Datensatz
statistic = XX, # Die Funktion
formula = XX ~ XX + XX, # Die Modellgleichung
R = 1000) # Die Anzahl der Ziehungen
# Lasse Bootstrap laufen
cali_boot<- boot(data = avocado_cali,
statistic = boot_lm,
formula = Verkaufsvolumen ~ Preis + Typ,
R = 1000)
cali_boot
und evaluiere das Ergebnis. Die
Spalte original
enthält die Gewichte für den gesamten
Datensatz. Die Spalte bias
den Unterschied zwischen dem
Mittel der Gewichte über die Bootstrap Samples und der Spalte
original
. Diese sollte möglichst klein sein. Die Spalte
std. error
enthält das Ergebnis von Interesse: den über
Bootstrap geschätzten Standardfehler der Gewichte.# Zeige Bootstrap Ergebnisse
cali_boot
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = avocado_cali, statistic = boot_lm, R = 1000, formula = Verkaufsvolumen ~
Preis + Typ)
Bootstrap Statistics :
original bias std. error
t1* 7418161 -10142 267313
t2* -1338574 7983 211733
t3* -5012181 -3101 121897
z = b / se
, d.h., das Gewicht
geteilt durch den Standardfehler folgt der z-Verteilung. Der Code unten
sieht die relevanten Werte aus dem Objekt und berechnet die
z-Werte für Intercept und die zwei Gewichte.# Berechne z-Werte
cali_boot$t0[1] / sd(cali_boot$t[,1])
(Intercept)
27.8
cali_boot$t0[2] / sd(cali_boot$t[,2])
Preis
-6.32
cali_boot$t0[3] / sd(cali_boot$t[,3])
Typorganic
-41.1
# Normale Regression
summary(lm(Verkaufsvolumen ~ Preis + Typ, avocado_cali))
Call:
lm(formula = Verkaufsvolumen ~ Preis + Typ, data = avocado_cali)
Residuals:
Min 1Q Median 3Q Max
-2029523 -338620 -65747 347099 4692280
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7418161 181045 40.97 < 2e-16 ***
Preis -1338574 155354 -8.62 2.8e-16 ***
Typorganic -5012181 121165 -41.37 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 745000 on 335 degrees of freedom
Multiple R-squared: 0.939, Adjusted R-squared: 0.939
F-statistic: 2.59e+03 on 2 and 335 DF, p-value: <2e-16
library(tidyverse)
# Normale regression ----------------
# Modell
m <- lm(formula = hwy ~ displ,
data = mpg)
# Evaluiere Modell
summary(m) # Summary
fitted(m) # gefittete Werte
resid(m) # Residuen
# Residualplot
plot(fitted(m), resid(m))
# non-parametric tests ----------------
# Teile Daten
mpg_suv <- mpg %>% filter(class == 'suv')
mpg_compact <- mpg %>% filter(class == 'compact')
# Wilcoxon
wilcoxon.test(mpg_suv$hwy, mpg_compact$hwy)
# Robuste Regression ----------------
library(quantreg)
library(Rfit)
# Quantil Regression
m_q <- rq(hwy ~ displ, data = mpg)
# Rang Regression
m_rb <- rfit(hwy ~ displ, data = mpg)
File | Zeilen | Spalten |
---|---|---|
avocado.csv | 17573 | 16 |
avocado_cali.csv | 338 | 16 |
Die Datensätze entstammen der Hass Avocado Board Webseite und enthalten wochenweise Avocado Verkaufsvolumina von April 2015 bis Mai 2018 für verschiedene Regionen der USA. Der kleinere Datensatz enthält nur die Daten für Kalifornien.
Name | Description |
---|---|
Verkaufsvolumen |
Das totale Verkaufsvolumen für Avocados in einer Woche. |
Verkaufsvolumen_index |
Ein Index von 0 bis 10 der das Verkaufsvolumen repräsentiert. |
Preis |
Der durchschnittliche Preis einer Avocado. |
Typ |
Der Typ der Avocado: “conventional” oder “organic”. |
Temparatur |
Die regionale Temperatur. |
Temperatur_USA |
Die landesweite Temperatur |
Luftfeuchtigkeit |
Die regionale Luftfeuchtigkeit. |
Luftfeuchtigkeit_USA |
Die landesweite Luftfeuchtigkeit_USA. |
Niederschlag |
Der regionale Niederschlag. |
Niederschlag_USA |
Der landesweite Niederschlag. |
Jahr |
Jahr. |
Jahreszeit |
Jahreszeit. |
Datum |
Datum. |
Region |
Region. |
Längengrad |
Längengrad. |
Breitengrad |
Breitengrad. |
Package | Installation |
---|---|
tidyverse |
install.packages("tidyverse") |
lubridate |
install.packages("lubridate") |
quantreg |
install.packages("quantreg") |
boot |
install.packages("boot") |
Rfit |
install.packages("Rfit") |
lme4 |
install.packages("lme4") |
rsq |
install.packages("rsq") |
Function | Package | Description |
---|---|---|
lm , glm |
stats |
Regressionsmodelle. |
summary |
base |
Modellergebnisse |
fitted , resid |
stats |
Extrahiere gefittete Werte und Residuen. |
rq |
quantreg |
Quantil-Regression. |
rfit |
Rfit |
Rangregression. |
wilcox.test |
stats |
Wilcoxon-Test. |
pbinom |
stats |
Kumulative Binomialverteilung (für den Vorzeichentest). |