rm(list=ls())
sys <- Sys.info()["sysname"]
if (sys == "Windows") {
Sys.setenv(LANG = "en")
} else {
Sys.setlocale("LC_MESSAGES", "en_US.utf8")
}
if (!require(car)) install.packages("car"); library(car) # Companion to Applied Regression: Package by John Fox accompanying his book; provides additional functions and plots for lm's and glm's
## Loading required package: car
## Warning: package 'car' was built under R version 4.1.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.1.3
if (!require(effects)) install.packages("effects"); library(effects) # graphical and tabular effect displays for fitted models
## Loading required package: effects
## Warning: package 'effects' was built under R version 4.1.3
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
if (!require(corrplot)) install.packages("corrplot"); library(corrplot) # used to plot a correlation matrix
## Loading required package: corrplot
## Warning: package 'corrplot' was built under R version 4.1.3
## corrplot 0.92 loaded
if (!require(ggplot2)) install.packages("ggplot2"); library(ggplot2) # graphics package
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.1.3
if(!require(MVA)) install.packages("MVA"); library(MVA) # bivariate boxplot
## Loading required package: MVA
## Warning: package 'MVA' was built under R version 4.1.3
## Loading required package: HSAUR2
## Warning: package 'HSAUR2' was built under R version 4.1.3
## Loading required package: tools
if(!require(KernSmooth)) install.packages("KernSmooth"); library(KernSmooth) # kernel density estimates
## Loading required package: KernSmooth
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
if (!require(performance)) install.packages("performance"); library(performance) # posterior predictive checks
## Loading required package: performance
## Warning: package 'performance' was built under R version 4.1.3
if (!require(see)) install.packages("see"); library(see) # posterior predictive checks
## Loading required package: see
## Warning: package 'see' was built under R version 4.1.3
# library(foreign) # Read Data Stored by Minitab, S, SAS, SPSS, Stata, Systat, Weka, dBase, ...
library(MASS)
library(lattice) # graphics package
#----------------------------------1. Simple Linear Regression---------------------------------------------
# CEO salary, sales, and profit data:
dataceo <- read.table("http://statmath.wu.ac.at/~wurzer/ADAR/Data/dataceo.txt", header = TRUE)
# salary = Gehalt + Boni
# totcomp = CEO-Gesamtvergütung
# tenure = # Anzahl der Jahre als CEO (=0, wenn weniger als 6 Monate)
# age = Alter des CEO
# sales = Gesamtverkaufserlös von Unternehmen i
# profits = Gewinn für Unternehmen i
# assets = Gesamtvermögen des Unternehmens i
# 53 der Fortune-500-Firmen wurden ausgeschlossen, weil Daten zu einer oder mehreren Variablen fehlten oder weil die CEOs keine Vergütung erhielten.
#---1.1 Preliminary Descriptive Analysis:---
names(dataceo) # get variable names
## [1] "salary" "totcomp" "tenure" "age" "sales" "profits" "assets"
summary(dataceo) # Informationen über die Variablen im Datenrahmen -> Ergebnisse des Befehls 'summary' hängen von dem in Klammern angegebenen Objekt ab
## salary totcomp tenure age
## Min. : 100 Min. : 100 Min. : 0.000 Min. :34.00
## 1st Qu.: 1084 1st Qu.: 1576 1st Qu.: 2.000 1st Qu.:52.00
## Median : 1600 Median : 2951 Median : 5.000 Median :57.00
## Mean : 2028 Mean : 8340 Mean : 7.834 Mean :56.47
## 3rd Qu.: 2348 3rd Qu.: 6043 3rd Qu.:10.000 3rd Qu.:61.00
## Max. :15250 Max. :589101 Max. :60.000 Max. :84.00
## sales profits assets
## Min. : 2896 Min. :-2669.0 Min. : 717.8
## 1st Qu.: 4184 1st Qu.: 108.5 1st Qu.: 3856.9
## Median : 6704 Median : 333.1 Median : 7810.8
## Mean : 11558 Mean : 700.5 Mean : 27054.3
## 3rd Qu.: 12977 3rd Qu.: 738.0 3rd Qu.: 21105.6
## Max. :161315 Max. :22071.0 Max. :668641.0
dim(dataceo) # Anzahl der Zeilen und Spalten im Datensatz
## [1] 447 7
complete.cases(dataceo) # Sind die Einzelfälle (= Beobachtungseinheiten = Zeilen) vollständig (keine NA's)?
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [136] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [151] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [166] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [181] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [196] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [211] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [226] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [241] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [256] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [271] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [286] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [301] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [316] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [331] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [346] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [361] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [376] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [391] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [406] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [421] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [436] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
table(complete.cases(dataceo)) # Tabelle mit TRUE (vollständig) und FALSE (NA's vorhanden) Fällen
##
## TRUE
## 447
#####Punktplot----
attach(dataceo) # Objekt 'dataceo' an den R-Suchpfad anhängen -> direkter Zugriff auf die Variablen (salary, totcomp etc.) wird möglich
plot(assets, salary) # scatterplot ,Streudiagramm
options(scipen = 10) # Strafe für die wissenschaftliche Notation, d.h. die feste Notation wird bis zu einem gewissen Grad bevorzugt (siehe ?Optionen)
plot(assets, salary)
abline(lm(salary ~ assets)) # a-b-line': fügt eine Regressionslinie durch das aktuelle Diagramm hinzu. Funktioniert nur, wenn bereits ein Diagramm erstellt wurde.
#####Balkendiagramm----
hist(salary, freq = F) # histogram
########Boxplot----
Boxplot(~ salary) # Boxplot wird erstellt; Ausreißer (d. h. die Zeilennummern der entsprechenden Beobachtungen) werden in der Konsole zurückgegeben
## [1] "192" "27" "185" "5" "7" "6" "23" "96" "390" "207"
#####Punktediagramm mit 25 und 75 Quantil----
scatterplot(salary ~ assets) # komplexere Scatterplot-Version (einschließlich marginaler Boxplots,
# eine nicht-parametrische Regressionsgerade, nicht-parametrische Schätzung der Streuung), die im Paket "car" enthalten ist
scatterplot(salary ~ assets, ellipse = TRUE, smooth = list(style="lines")) # eine Konzentrationsellipse enthalten (d. h. ein "2-dim-Konfidenzintervall"; die Ellipse zeigt
# der Bereich, der 95 % aller Stichproben enthält, die aus der zugrunde liegenden Normalverteilung gezogen werden könnten)
smoothScatter(assets, salary) # geglättetes Streudiagramm, um übereinander liegende Beobachtungen zu berücksichtigen
bvbox(dataceo[, c("assets", "salary")], xlab = "assets", ylab = "salary") # bivariater Boxplot zur Identifizierung von 2-dim Ausreißern
# Alternative zum bivariaten Boxplot - konvexe Hülle:
(hull <- with(dataceo, chull(assets, salary)))
## [1] 11 100 338 423 358 376 157 27 192 7
with(dataceo,
plot(assets, salary, pch = 1, xlab = "assets", ylab = "salary"))
with(dataceo,
polygon(assets[hull], salary[hull], density = 15, angle = 30))
# Korrelation unter Verwendung aller Punkte gegenüber dem Entfernen der Punkte, die auf der konvexen Hülle liegen:
with(dataceo, cor(assets, salary))
## [1] 0.431036
with(dataceo, cor(assets[-hull], salary[-hull]))
## [1] 0.4977579
# Plots auf der Grundlage von bivariaten Kernel-Dichte-Schätzungen:
image(kde2d(dataceo$assets, dataceo$salary, n = 250)) # colored image
persp(kde2d(dataceo$assets, dataceo$salary, n = 50)) # perspective plot
contour(kde2d(dataceo$assets, dataceo$salary)) # contour plot
xyplot(salary ~ assets) # Streudiagramm unter Verwendung des Lattice-Pakets (könnte auch erweitert werden)
##lineare Dartesllung----
ggplot(dataceo, aes(x = assets, y = salary)) + # Streudiagramm mit ggplot2-Paket
geom_point(shape = 1) + # use circles die Zahl verändert die Zeichen der einzelnen Datenpunkte
geom_smooth(method = lm) # lineare Regressionslinie hinzufügen
## `geom_smooth()` using formula 'y ~ x'
## Kurvendarstellung----
ggplot(dataceo, aes(x = assets, y = salary)) +
geom_point(shape = 1) + # Use circles
geom_smooth() # add loess smooth
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
# # # # # Stichwort: Parametrische, nichtparametrische und semiparametrische Verfahren -> Begleitfolien# # # # #
# # # # # Schlüsselwort: Lokale polynomiale Regression -> Begleitfolien# # # # #
#---1.2 Parameter-Schätzung:---
#dataceo[,7]
regfit1.1 <- lm(dataceo[, 1] ~ dataceo[, 7]) # linke Seite der "~": abhängige Variable; rechte Seite: unabhängige Variable(n)
regfit1.1 <- lm(dataceo$salary ~ dataceo$assets) # the same
regfit1.1 <- lm(salary ~ assets, data = dataceo) # the same
regfit1.1 # geschätzte Regressionskoeffizienten
##
## Call:
## lm(formula = salary ~ assets, data = dataceo)
##
## Coefficients:
## (Intercept) assets
## 1716.84888 0.01148
names(regfit1.1) # Namen der Listenelemente des Objekts regfit1.1
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
regfit1.1$coefficients
## (Intercept) assets
## 1716.84887504 0.01148313
summary(regfit1.1) # detailliertere Informationen über das Modellobjekt Gesamte Auswertung der Regressionsrechnung
##
## Call:
## lm(formula = salary ~ assets, data = dataceo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5309.7 -812.4 -290.1 315.7 13032.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1716.84887 79.79613 21.52 <2e-16 ***
## assets 0.01148 0.00114 10.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1556 on 445 degrees of freedom
## Multiple R-squared: 0.1858, Adjusted R-squared: 0.184
## F-statistic: 101.5 on 1 and 445 DF, p-value: < 2.2e-16
#---1.3 Data Transformation:----
par(mfrow = c(1, 2))
plot(density(salary), xlab = "salary", main = " (a) ") #Dichtefunktion von Gehalt
plot(density(log(salary)), xlab = "log of salary", main = " (b) ") # Log-Transformation zur Normalität #Verwendung des logarithmus zur ANnäherung an Normalfunktion
par(mfrow = c(1, 1))
scatterplot(salary ~ assets) # beide Variablen sind verzerrt
scatterplot(salary ~ assets, log = "x") # Beachten Sie die logarithmische Skala auf der x-Achse log scale verwendet z.B. Faktor 100
scatterplot(salary ~ assets, log = "y") # logarithmische Skala auf der y-Achse
# durch die Logarithmierung beider Achsen kommt es zu einer linearen Darstellung
scatterplot(salary ~ assets, log = "xy") # logarithmische Skala auf beiden Achsen -> Umwandlung in Linearität
tenure.1 <- tenure + 1 #andernfalls würde log von 0 Probleme verursachen
scatterplot(salary ~ tenure.1, log = "xy")
scatterplot(salary ~ age, log = "y")
regfit1.2 <- lm(log(salary) ~ log(assets), data = dataceo)
regfit1.2 # ein Anstieg des Vermögens ist nun mit einem proportionalen (und nicht einem absoluten wie in regfit1.1) Anstieg des Gehalts verbunden
##
## Call:
## lm(formula = log(salary) ~ log(assets), data = dataceo)
##
## Coefficients:
## (Intercept) log(assets)
## 5.3133 0.2267
1.01 ^ coef(regfit1.2)[2] # Beim Vergleich von CEO-Paaren, die sich bei den Vermögenswerten ihres Unternehmens um 1 % unterscheiden, sind die CEOs der Unternehmen mit den um 1 % höheren
## log(assets)
## 1.002258
# Vermögenswerte haben ein um 0,2 % höheres Gehalt # -> von Ökonomen "Elastizität" genannt
#andere Transformationen: Quadratwurzel, quadratisch, kubisch, ... -> abhängig von der Struktur der Daten und dem Ziel der Transformation
# (Beziehungen linearisieren, Streuung ausgleichen usw.)
#----------------------------------2. Multiple Linear Regression-------------------------------------------
#---2.1 Preliminary Descriptive Analysis:---
#Streudiagramm mehrerer Variablen----
pairs(dataceo[c(1, 3, 4, 7)]) # Streudiagramm-Matrix der beteiligten Variablen (paarweise)
pairs(dataceo[c(1, 3, 4, 7)], panel = panel.smooth) # jedem Panel des Bildschirms eine glatte Oberfläche hinzufügen
pairs(dataceo[c(1, 3, 4, 7)], upper.panel = panel.smooth, lower.panel = NULL) # keine Paneele oberhalb der Diagonale
# um zu sehen, wie man die Funktionen des Panels weiter anpassen kann, siehe ?pairs für Beispiele
scatterplotMatrix(~ salary + tenure + age + assets) # 'car' package version
scatterplotMatrix(~ salary + tenure + age + assets, diagonal = list(method = "boxplot")) # boxplots instead of
# adaptive (= mit variabler Bandbreite) Kernel-Dichte-Schätzungen in den diagonalen Feldern des Diagramms
#Korrelatinsmatrix erstellen----
cor.dataceo <- cor(dataceo[c(1, 3, 4, 7)]) # eine Korrelationsmatrix zu berechnen
corrplot(cor.dataceo, method = "ellipse") # Korrelationsmatrix visualisieren (Paket 'corrplot')
corrplot(cor.dataceo, method = "number", type = "lower") # siehe ?corrplot für mögliche Anpassungen
#---2.2 "Naive" Schätzung eines multiplen Regressionsproblems (d. h. es werden keine Annahmen über Linearität usw. im Voraus getroffen)
# -> Regression des Gehalts auf Betriebszugehörigkeit, Alter und Vermögen:---
regfit2.1 <- lm(salary ~ tenure + age + assets)
summary(regfit2.1)
##
## Call:
## lm(formula = salary ~ tenure + age + assets)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5429.5 -770.5 -275.5 338.6 12718.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1071.28459 638.36341 1.678 0.09402 .
## tenure 28.39738 9.66493 2.938 0.00347 **
## age 7.62357 11.72896 0.650 0.51604
## assets 0.01121 0.00113 9.917 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1538 on 443 degrees of freedom
## Multiple R-squared: 0.2084, Adjusted R-squared: 0.203
## F-statistic: 38.88 on 3 and 443 DF, p-value: < 2.2e-16
round(confint(regfit2.1), 3) # Berechnung der Konfidenzintervalle für alle Koeffizienten die 3 gibt die Anzahl der Nachkommastellen an
## 2.5 % 97.5 %
## (Intercept) -183.312 2325.882
## tenure 9.403 47.392
## age -15.428 30.675
## assets 0.009 0.013
#---2.3 Diagnostics:---
qqnorm(residuals(regfit2.1)) # QQ-Plot zur Identifizierung nicht normaler Fehler
qqline(residuals(regfit2.1)) # dem QQ-Plot eine Linie hinzufügen, die eine Normalverteilung darstellt
shapiro.test(residuals(regfit2.1)) # hochsignifikant -> Hinweis auf nicht-normale Fehler
##
## Shapiro-Wilk normality test
##
## data: residuals(regfit2.1)
## W = 0.71725, p-value < 2.2e-16
ks.test(residuals(regfit2.1), "pnorm") # hochsignifikant -> auch Hinweis auf nicht-normale Fehler
##
## One-sample Kolmogorov-Smirnov test
##
## data: residuals(regfit2.1)
## D = 0.63087, p-value < 2.2e-16
## alternative hypothesis: two-sided
# 'car' package diagnostics:
qqPlot(regfit2.1, line = "none") # eine annähernde punktweise 95%ige Konfidenzhüllkurve auf der Grundlage eines parametrischen Bootstraps verwendet
## [1] 27 192
# (berechnet unter der Annahme, dass die Fehler normalverteilt sind) wenn sich die Punkte innerhalb des Konfidenzintervalls befinden wäre das Modell in Ordnung
# Identify non-linear relationships:
residualPlots(regfit2.1) # Diagramm unten rechts: gekrümmter Trend (sollte linear sein) -> Anpassung nicht angemessen; alle Tests auf mangelnde Anpassung (für Linearität)
## Test stat Pr(>|Test stat|)
## tenure -2.1628 0.0310896 *
## age -2.0328 0.0426727 *
## assets -3.4680 0.0005757 ***
## Tukey test -3.3912 0.0006959 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# der einzelnen Variablen zeigen das gleiche Ergebnis (p-Werte < 0,05) die Grafik zeigt eine Kurve darum ist es ein Problem
plot(fitted.values(regfit2.1), rstandard(regfit2.1)) # Nicht konstante Fehlervarianz identifizieren -> hier schwer zu erkennen
#wenn man das naive Modell zuerst schätzt sieht man wo die Probleme liegen
# -> Modellverfeinerung 1: Nehmen Sie den Logarithmus der Antwort:
regfit2.2 <- lm(log(salary) ~ tenure + age + assets) #die Krümmung wird versucht mit der Logarithmusfunktion wegzubekommen
summary(regfit2.2)
##
## Call:
## lm(formula = log(salary) ~ tenure + age + assets)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3610 -0.3797 0.0023 0.3646 2.1555
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.7842961490 0.2436236184 27.847 < 2e-16 ***
## tenure 0.0068813935 0.0036885029 1.866 0.0628 .
## age 0.0081186011 0.0044762139 1.814 0.0704 .
## assets 0.0000035201 0.0000004314 8.160 3.49e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5869 on 443 degrees of freedom
## Multiple R-squared: 0.1579, Adjusted R-squared: 0.1522
## F-statistic: 27.69 on 3 and 443 DF, p-value: < 2.2e-16
qqPlot(regfit2.2, line = "none") #die Punkte liegen schon näher an dem gemeinsam geschätzeten Konfidenzintervall, es gibt nicht mehr viele problematische Beobachtungen
## [1] 27 100
residualPlots(regfit2.2) # Residueenplot ist nach wie vor ein Problem, die Signifikanz hat sich im Vergleich zu 2.1 verbessert (siehe Statistik)
## Test stat Pr(>|Test stat|)
## tenure -2.2732 0.0234946 *
## age -2.5724 0.0104253 *
## assets -3.8652 0.0001276 ***
## Tukey test -4.2438 0.00002198 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted.values(regfit2.2), rstandard(regfit2.2))
# -> model refinement 2: transform predictors: es wird der logarithmus von den Prädiktoren assets und tenure gerechnet
regfit2.3 <- lm(log(salary) ~ log(tenure + 1) + age + log(assets)) # 1 hinzufügen, um -Inf-Werte aufgrund von log(0) zu vermeiden
regfit2.3 <- lm(log(salary) ~ log1p(tenure) + age + log(assets)) # the same
# andere mögliche Transformationen bei Vorhandensein von Nullen: Box-Cox-Transformationen (BC), inverse hyperbolische Sinustransformation (IHS), Mischungsmodelle
summary(regfit2.3)
##
## Call:
## lm(formula = log(salary) ~ log1p(tenure) + age + log(assets))
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6089 -0.2961 -0.0168 0.3218 1.9555
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.922576 0.275706 17.854 < 2e-16 ***
## log1p(tenure) 0.135574 0.032170 4.214 0.0000304 ***
## age 0.002826 0.004234 0.667 0.505
## log(assets) 0.225227 0.020302 11.094 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5524 on 443 degrees of freedom
## Multiple R-squared: 0.2539, Adjusted R-squared: 0.2488
## F-statistic: 50.25 on 3 and 443 DF, p-value: < 2.2e-16
qqPlot(regfit2.3) #sieht schon viel besser aus
## [1] 27 100
residualPlots(regfit2.3) # alle Diagramme und Tests zeigen jetzt lineare Muster; im Allgemeinen können die erforderlichen Transformationen komplizierter sein
## Test stat Pr(>|Test stat|)
## log1p(tenure) -0.6861 0.4930
## age -1.5827 0.1142
## log(assets) 0.2438 0.8075
## Tukey test -0.5616 0.5744
#die Krümmung von tenure und assets hat sich stark verbessert
plot(fitted.values(regfit2.3), rstandard(regfit2.3)) # die Varuianz hat sich auch stark verbessert, konstant über alle fitted values, haben die gleiche Variation
# Ausreißerdiagnose (ein einflussreicher Ausreißer scheint vorhanden zu sein - Beobachtung Nr. 100):
residualPlots(regfit2.3, id = TRUE) # id=TRUE, die Ausreißer werden gelabelt values outside -2,+2
## Test stat Pr(>|Test stat|)
## log1p(tenure) -0.6861 0.4930
## age -1.5827 0.1142
## log(assets) 0.2438 0.8075
## Tukey test -0.5616 0.5744
qqPlot(regfit2.3, id = TRUE, line = "none")
## [1] 27 100
#Test auf Ausreißer----
outlierTest(regfit2.3) # zeigt die outlier an, kann sich jedoch aufgrund der Normalverteilungskurve andern!!! Signifikanzlevel
## rstudent unadjusted p-value Bonferroni p
## 100 -6.940743 0.000000000013968 0.0000000062436
outlierTest(regfit2.3)$p * 447 # == "Bonferroni p" -> adjustment for multiple testing 'bonferronie behebt das Problem das in der Vorzeile angesprochen wurde
## 100
## 0.000000006243589
#447 hängt mit der EInstellung des Konfidenzintervalls an
#Plot zur ANsicht von Ausreißer----
infIndexPlot(regfit2.3, id = TRUE) # residuen sind eingezeichnet, jeder Punkt ist getestet und Bonferroni korrigiert 1Zeile obere AUsreißer
# 2 Zeile Mittele werte, 3 Zeile untere Ausreißer, 4. Zeile identifizerit Werte mit großem Einfluss
# (Regressions-)Ausreißer: y-Werte, die in Abhängigkeit von den Werten der Prädiktoren ungewöhnlich sind;
# große student Restwerte (Panels 2 und 3) deuten auf Ausreißer hin
# ---
# leverage values: Beobachtungen, die relativ weit vom Regressionsraum entfernt sind;
# sie haben potenziell einen großen Einfluss auf die Regressionskoeffizienten
# large hat-values (panel 4) indicate leverage values
# ---
# einflussreiche Punkte sind Beobachtungen, die sowohl abweichend sind als auch einen hohen leverage
# große Cook's distances (panel 1) einflussreiche Punkte angeben
influencePlot(regfit2.3, id = TRUE) # all in one: x-axis: hat values; y-axis: studentized residuals;
## StudRes Hat CookD
## 27 3.5972258 0.005499006 0.017418228
## 67 1.0147136 0.042694920 0.011479521
## 97 -2.5894799 0.024753055 0.042007006
## 100 -6.9407430 0.019699244 0.218723617
## 184 -0.4029077 0.041768248 0.001772347
#1. wenn es ein outlier ist sieht man es , man hat die Werte der hat values und die gestrichelten Linien zeigt an welche Werte als problematisch anzusehen sind
#Cooks distance ist anhand des Durchmessers der Blasen dargestellt, es zeigt das Pkt 100 ein großes Problem ist
# Flächen von Kreisen: proportional zum Wert von Cook's distance,
# Punkte mit größter hat values (Linien mit dem Doppelten und Dreifachen des Durchschnitts hat value), Cook's distances,
# or studentized residuals are flagged
# Faustregeln: Hat werte, die größer sind als das Dreifache des Durchschnitts, deuten auf Hebelwirkungen hin, Cook's Distanzen
# > 1 indicate influential points
# model refinement 4: Neueinschätzung des Modells ohne den Ausreißer
regfit2.4 <- update(regfit2.3, subset = -c(100)) # take previous model regift2.3, but without observation no. 100
# hier wird die Reihe 100 ausgeschlossen
compareCoefs(regfit2.3, regfit2.4) # Vergleichen Sie die Koeffizienten der beiden Modelle, bei ALter erhöht sich der Koeffizient am meisten
## Calls:
## 1: lm(formula = log(salary) ~ log1p(tenure) + age + log(assets))
## 2: lm(formula = log(salary) ~ log1p(tenure) + age + log(assets), subset =
## -c(100))
##
## Model 1 Model 2
## (Intercept) 4.923 4.734
## SE 0.276 0.264
##
## log1p(tenure) 0.1356 0.1486
## SE 0.0322 0.0306
##
## age 0.00283 0.00396
## SE 0.00423 0.00403
##
## log(assets) 0.2252 0.2372
## SE 0.0203 0.0194
##
cC <- as.data.frame(compareCoefs(regfit2.3, regfit2.4, print = FALSE))
round(cbind(coef = cC$'Model 2' / cC$'Model 1', se = cC$'SE 2' / cC$'SE 1') * 100, 2) # heier sieht man die Koeffizietne sind größer als zuvor was sehr gut ist
## coef se
## [1,] 96.16 95.58
## [2,] 109.61 95.25
## [3,] 140.27 95.15
## [4,] 105.30 95.44
#das R² hat sich ebenfalls von 25 auf 29% gesteigert. summary()
# Eliminiert man nur diesen einen von 447 CEOs, erhöht sich der log(Assets)-Koeffizient um 5%
# log1p(tenure) Koeffizient um fast 10 % und der Alterskoeffizient um über 40 %. Standardfehler
# sind nicht so stark betroffen, gehen aber alle um etwa 5% zurück -> höhere Genauigkeit der Schätzung
# die Fehlernormalität der vier Modelle zu vergleichen, um die Auswirkungen der Verfeinerungen zu erkennen:
par(mfrow = c(2, 2)) #den Parameter mfrow ändern
qqPlot(regfit2.1)
## [1] 27 192
qqPlot(regfit2.2)
## [1] 27 100
qqPlot(regfit2.3)
## [1] 27 100
qqPlot(regfit2.4)
## 27 157
## 27 156
par(mfrow = c(1, 1)) # restore the initial settings for mfrow
#jetzt schaut die Normalverteilung gut aus
residualPlots(regfit2.4)
## Test stat Pr(>|Test stat|)
## log1p(tenure) -0.0560 0.9554
## age -1.2634 0.2071
## log(assets) 0.6021 0.5474
## Tukey test 0.3434 0.7313
# linearität ist weitgehenst gegeben
# die Anpassung der vier Modelle zu vergleichen:
par(mfrow = c(2, 2))
plot(fitted.values(regfit2.1), rstandard(regfit2.1))
plot(fitted.values(regfit2.2), rstandard(regfit2.2))
plot(fitted.values(regfit2.3), rstandard(regfit2.3))
plot(fitted.values(regfit2.4), rstandard(regfit2.4))
par(mfrow = c(1, 1))
vif(regfit2.4) # Erkennen von Multikollinearität---- - Varianzinflationsfaktoren ################Multikollinearitätscheck
## log1p(tenure) age log(assets)
## 1.198365 1.208011 1.010450
#Multikolinearität wie ist die Vorhersagekraft des Modells durch einen einzelnen Prädiktor auf alle vorhergesagt werden
#in der Regressionsrechneung verwendet man 4-5 spezifische Untersuchungen die zur Abschätzung der Modellierung angewandt werden.
#Wie verändert sich der Effekt beim hinzufügen anderer Variablen
# -> Faustregel: vifs sollten 4 nicht überschreiten (-> Breite der Konfidenzintervalle)
# wird verdoppelt / die Genauigkeit der Schätzung wird halbiert)
sqrt(vif(regfit2.4)) # -> Multikollinearität ist hier kein Problem, die prozentualen Veränderungen sind gering
## log1p(tenure) age log(assets)
## 1.094699 1.099095 1.005212
# Übliche Methoden zur Behebung von Multikollinearität: Modellwiederholung, Variable
# Es gibt verschiedene Wege um mit Multikollinearität umzugehen, z.B. vorgelagerte PCA
# Auswahl, voreingenommene Schätzung, Einbeziehung zusätzlicher Vorabinformationen, vorgelagerte Analyse mit PCA
# Vorsicht! Keine dieser Maßnahmen sollte mechanisch durchgeführt werden, da sie möglicherweise nicht
# vertretbar. Eine schnelle Lösung gibt es nicht!
#---2.4 Interpretation:---
# Konfidenzintervalle für das Modell regfit2.4:
confint(regfit2.4)
## 2.5 % 97.5 %
## (Intercept) 4.215747669 5.25153754
## log1p(tenure) 0.088381151 0.20881873
## age -0.003953078 0.01188138
## log(assets) 0.199076642 0.27524121
cbind(Estimate = coef(regfit2.4), confint(regfit2.4))
## Estimate 2.5 % 97.5 %
## (Intercept) 4.73364261 4.215747669 5.25153754
## log1p(tenure) 0.14859994 0.088381151 0.20881873
## age 0.00396415 -0.003953078 0.01188138
## log(assets) 0.23715893 0.199076642 0.27524121
# Effect Plots erleichtern die Interpretation von Koeffizienten in komplexen Modellen (d. h. Modelle, die Folgendes enthalten
# transformierte Variablen, Interaktionen usw.):
#effect plots age ist nicht logarithmiert, es zeigt die geschätzte REgressionslinie von Age auf salary, es zeigt keinen wirklich starken effect, age war nicht signifikant
plot(effect("age", regfit2.4))# Wie wirkt sich das Alter nach dem oben beschriebenen Modell auf das log(Gehalt) aus?
# die übrigen Prädiktoren auf typische Werte (z. B. Mittelwerte) setzen?
eff.regfit2.4 <- allEffects(regfit2.4) # zeigt alle effekte der REgressionsrechnung an
eff.regfit2.4
## model: log(salary) ~ log1p(tenure) + age + log(assets)
##
## tenure effect
## tenure
## 0 20 30 40 60
## 7.375866 7.828282 7.886156 7.927703 7.986742
##
## age effect
## age
## 30 50 60 70 80
## 7.593986 7.673269 7.712911 7.752552 7.792194
##
## assets effect
## assets
## 700 200000 300000 500000 700000
## 6.833990 8.175122 8.271282 8.392429 8.472226
class(eff.regfit2.4) # object of type 'efflist'
## [1] "efflist"
class(eff.regfit2.4$age) # object of type 'eff'
## [1] "eff"
plot(allEffects(regfit2.4)) # x-axis is on the original scale in every plot,
#es zeigt eine Kurve an weil die Prädiktoren logarithmiert sind. tenure hat enenstarken effekt konfidenzintervall ist nicht so weit, stärkster effect bei assets
#wenn assets steigen, steigt salary stark an;;
# aber zwei der Prädiktoren wurden im Modell transformiert -> scheinbar nicht-lineare Beziehungen!
plot(effect("tenure", regfit2.4)) # Note - effect is not called tenure, but log1p(tenure)
## NOTE: tenure does not appear in the model
plot(effect("log1p(tenure)", regfit2.4))
plot(allEffects(regfit2.4)[1]) # the same
# uDie logarithmische Skala auf der x-Achse führt zu einer grafischen Darstellung der im Modell verwendeten linearen Beziehung
plot(effect("log(assets)", regfit2.4), transform.x = list(assets = list(trans = log, inverse = function(x) exp(x))),
ticks.x = list(assets = list(at = c(1e2, 1e3, 1e4, 1e5))))
#hier wird wieder zurücklogarithmiert für die Grafik, linearer Zusammenhang
#aber log und exp verursachen Probleme mit dem Plot für Tenure:
plot(effect("log1p(tenure)", regfit2.4), transform.x = list(tenure = list(trans = log, inverse = function(x) exp(x))),
ticks.x = list(tenure = list(at = c(10, 20, 30, 40, 50, 60))))
# Lösung: log1p anstelle von log und expm1 anstelle von exp(x) verwenden
plot(effect("log1p(tenure)", regfit2.4), transform.x = list(tenure = list(trans = log1p, inverse = function(x) expm1(x))),
ticks.x = list(tenure = list(at = c(10, 20, 30, 40, 50, 60))))
# all effects, same y limits
plot(allEffects(regfit2.4), transform.x = list(tenure = list(trans = log1p, inverse = function(x) expm1(x)),
assets = list(trans = log, inverse = function(x) exp(x))),
ylim = c(6.5, 9),
ticks.x = list(tenure = list(at = c(0, 10, 20, 30, 40, 60)), assets = list(at = c(1e2, 1e3, 1e4, 1e5))))
# die STeigung zeigt die Stärke des Effekts an
#---2.5 Prediction:---
predict(regfit2.4) # prediction for all observations in the data set
## 1 2 3 4 5 6 7 8
## 8.239068 7.871394 7.915165 7.972891 8.452412 7.943661 8.564511 7.850641
## 9 10 11 12 13 14 15 16
## 7.658088 7.742206 8.566611 7.770588 7.744820 7.705575 7.408011 8.125529
## 17 18 19 20 21 22 23 24
## 7.477561 8.434612 8.287656 7.594916 7.768171 7.811271 8.385632 7.512667
## 25 26 27 28 29 30 31 32
## 7.587097 7.373760 7.520433 7.525725 7.843131 7.399033 7.688896 7.827345
## 33 34 35 36 37 38 39 40
## 7.526742 7.549279 7.642985 7.554458 8.322039 7.810138 7.556566 7.432221
## 41 42 43 44 45 46 47 48
## 7.425861 7.316710 7.740658 7.628451 7.892725 7.532024 8.307004 8.083665
## 49 50 51 52 53 54 55 56
## 7.374324 7.126370 7.564849 7.884972 8.138529 7.392757 7.664425 7.588028
## 57 58 59 60 61 62 63 64
## 8.065699 7.401105 7.367759 7.598069 7.417328 8.047607 7.441371 7.600908
## 65 66 67 68 69 70 71 72
## 8.136913 7.557353 7.375952 8.348046 8.040059 7.525159 7.580606 7.477470
## 73 74 75 76 77 78 79 80
## 7.439437 7.412046 7.381494 7.649282 7.839062 7.585681 7.153812 7.376199
## 81 82 83 84 85 86 87 88
## 7.362312 7.408813 7.624837 7.369778 7.681059 7.172187 7.077835 7.109611
## 89 90 91 92 93 94 95 96
## 7.603022 6.842794 7.926618 7.572617 7.466052 7.724806 7.633587 7.737979
## 97 98 99 101 102 103 104 105
## 7.716766 7.362253 7.393001 7.667408 7.002067 7.277403 7.099295 7.717770
## 106 107 108 109 110 111 112 113
## 7.713393 7.460036 7.227967 7.506033 7.470968 7.625392 8.009696 7.339863
## 114 115 116 117 118 119 120 121
## 7.147152 8.119914 7.584925 7.375389 7.462492 7.377277 7.158513 7.404762
## 122 123 124 125 126 127 128 129
## 7.488136 7.462492 7.234089 7.963833 7.358634 7.437166 7.513418 7.302118
## 130 131 132 133 134 135 136 137
## 7.253721 7.304196 7.803433 7.692087 7.272140 7.690717 7.167197 7.302125
## 138 139 140 141 142 143 144 145
## 7.562712 7.470210 7.911772 7.109070 7.692574 7.201065 7.273892 8.137601
## 146 147 148 149 150 151 152 153
## 7.677565 7.484303 7.371536 7.037107 7.546033 7.226718 7.509584 7.166938
## 154 155 156 157 158 159 160 161
## 7.824511 7.302014 7.386841 7.151727 7.488681 7.286614 7.505371 7.719379
## 162 163 164 165 166 167 168 169
## 7.491362 6.977822 7.248367 7.331342 7.781415 7.043502 7.516930 6.845556
## 170 171 172 173 174 175 176 177
## 7.121546 7.182849 7.334245 7.324170 7.358155 7.389458 7.296786 6.990053
## 178 179 180 181 182 183 184 185
## 7.309083 7.556536 7.412942 7.292205 7.883642 7.236198 7.814120 8.114031
## 186 187 188 189 190 191 192 193
## 7.561205 7.328476 7.201811 8.050881 7.167953 7.284719 7.929682 7.980326
## 194 195 196 197 198 199 200 201
## 7.557947 7.182578 7.852215 7.403394 7.174774 7.229544 7.143303 7.168723
## 202 203 204 205 206 207 208 209
## 7.015776 7.062048 7.760359 7.658560 7.300132 8.017544 7.176640 7.525308
## 210 211 212 213 214 215 216 217
## 7.454980 7.469544 7.420358 7.236562 7.716157 7.868365 7.572249 7.314619
## 218 219 220 221 222 223 224 225
## 6.712702 7.228025 7.298597 7.051842 7.331689 6.934388 6.730781 7.567837
## 226 227 228 229 230 231 232 233
## 7.316990 7.272089 7.279092 7.255703 7.156636 7.094387 7.598281 7.315938
## 234 235 236 237 238 239 240 241
## 7.253819 7.173870 7.753901 7.720411 7.885395 7.682278 7.591253 7.267535
## 242 243 244 245 246 247 248 249
## 7.250748 7.305012 7.236270 7.381814 7.734301 7.642811 7.141550 7.675437
## 250 251 252 253 254 255 256 257
## 6.895731 7.850861 7.069005 7.333114 7.264327 7.245901 7.656981 7.525299
## 258 259 260 261 262 263 264 265
## 6.980997 7.729983 7.471852 7.043927 7.331932 7.336529 7.266556 7.143375
## 266 267 268 269 270 271 272 273
## 7.860962 7.143819 7.339614 7.318245 7.281275 7.039983 7.068146 7.561698
## 274 275 276 277 278 279 280 281
## 7.102904 7.553592 7.217998 7.375642 7.664377 7.123439 6.863714 6.895893
## 282 283 284 285 286 287 288 289
## 7.435198 7.730933 7.525595 7.428011 7.666426 7.563632 7.462686 7.612323
## 290 291 292 293 294 295 296 297
## 7.156261 7.062889 7.321024 7.117594 6.907713 7.410624 7.408148 7.165703
## 298 299 300 301 302 303 304 305
## 6.979370 7.426612 7.189744 7.338120 6.733712 7.065277 7.190917 7.149525
## 306 307 308 309 310 311 312 313
## 6.794055 7.509515 7.209560 7.620135 7.148963 7.129611 7.301932 6.786533
## 314 315 316 317 318 319 320 321
## 7.285719 7.389998 7.344067 6.879372 7.348492 7.166498 7.598641 7.329307
## 322 323 324 325 326 327 328 329
## 7.036507 7.396825 7.278637 7.009430 7.145603 7.258979 7.146036 7.097420
## 330 331 332 333 334 335 336 337
## 7.030880 7.426829 7.186451 7.828278 7.260174 7.093515 7.283360 7.446808
## 338 339 340 341 342 343 344 345
## 6.607673 7.061825 6.896061 6.620364 7.140483 7.175217 7.160072 7.125302
## 346 347 348 349 350 351 352 353
## 7.725400 7.453996 7.419908 7.242592 7.263169 7.530294 7.112906 7.135712
## 354 355 356 357 358 359 360 361
## 7.164497 7.078723 7.603765 7.657317 7.103129 7.178735 6.951316 7.510510
## 362 363 364 365 366 367 368 369
## 7.163889 7.419603 7.300447 7.296361 7.694986 6.988951 7.072079 7.297803
## 370 371 372 373 374 375 376 377
## 7.386860 6.924590 7.196563 6.972699 6.722807 7.192293 7.281539 7.435121
## 378 379 380 381 382 383 384 385
## 6.934721 7.429700 6.718275 7.371539 6.973920 7.507789 7.740742 7.235060
## 386 387 388 389 390 391 392 393
## 6.885521 7.737115 7.808918 7.418168 7.754334 7.152038 7.425065 7.255515
## 394 395 396 397 398 399 400 401
## 6.686934 6.954325 7.165889 7.393019 7.242633 7.178372 7.120853 7.003526
## 402 403 404 405 406 407 408 409
## 7.107373 7.081541 7.081883 7.249014 6.983906 6.971627 6.970274 6.977050
## 410 411 412 413 414 415 416 417
## 7.235711 7.173681 7.117257 7.823807 6.694830 7.202549 7.318970 6.971498
## 418 419 420 421 422 423 424 425
## 7.211300 6.831632 6.911005 8.043616 7.184507 6.935172 7.122111 7.560082
## 426 427 428 429 430 431 432 433
## 7.331249 7.602936 7.361365 6.777894 6.761430 6.813436 7.638508 6.750751
## 434 435 436 437 438 439 440 441
## 7.159163 7.765682 7.270522 7.011500 7.275226 7.381188 7.918719 6.821923
## 442 443 444 445 446 447
## 7.032031 7.810634 7.275977 7.178933 7.261196 7.132691
comp.salary <- data.frame(obs = log(dataceo$salary)[-100], est = predict(regfit2.4)) # compare observed and fitted/predicted values
head(comp.salary)
## obs est
## 1 8.016318 8.239068
## 2 8.707814 7.871394
## 3 8.180601 7.915165
## 4 8.101678 7.972891
## 5 9.210340 8.452412
## 6 9.145802 7.943661
head(resid(regfit2.4)) # difference between observed and predicted values
## 1 2 3 4 5 6
## -0.2227502 0.8364194 0.2654361 0.1287872 0.7579281 1.2021409
comp.salary2 <- data.frame(obs = dataceo$salary[-100], est = exp(predict(regfit2.4))) # comparison on the original scale
head(comp.salary2) # zeigt an ob man über oder unterschätzt
## obs est
## 1 3030 3786.010
## 2 6050 2621.217
## 3 3571 2738.498
## 4 3300 2901.232
## 5 10000 4686.364
## 6 9375 2817.657
head(predict(regfit2.4, interval = "confidence")) # predictions including confidence intervals
## fit lwr upr
## 1 8.239068 8.103023 8.375113
## 2 7.871394 7.707141 8.035648
## 3 7.915165 7.822757 8.007573
## 4 7.972891 7.872251 8.073530
## 5 8.452412 8.293110 8.611714
## 6 7.943661 7.846362 8.040960
#ceo kriegt in wirklichkeit 9,21 millionen sollte eigentlich 8,45 Millionen bekommen
newceo <- data.frame(tenure = 1, age = 60, assets = 500000) # suppose there is a new CEO with the given characteristics
# -> what is the predicted salary?
predict(regfit2.4, newdata = newceo)
## 1
## 8.186579
round(exp(predict(regfit2.4, newdata = newceo))) # 3.59 Millionen sollten dem neuen CEO bezahlt werden
## 1
## 3592
# Posterior predictive checks: Posterior predictive checks means âsimulating replicated data under the fitted model
# and then comparing these to the observed dataâ. Posterior predictive checks can be used to âlook
# for systematic discrepancies between real and simulated dataâ.
# direct usage of the function when using original (non-transformed) data: ####Multiple SImulationen auf Daten angewandt
#?ppcheck
#pp_check(regfit2.1)
# doesn't work after transformations took place:
#pp_check(regfit2.4)
# Transform data outside the model formula, then build and check it:
dat.trans <- regfit2.4$model
names(dat.trans) <- c("salary", "tenure", "age", "assets")
regfit.trans <- lm(salary ~ tenure + age + assets, dat.trans)
#pp_check(regfit.trans)
#---2.6 Checking model fit:---
# Posterior predictive checks: Posterior predictive checks means â???osimulating replicated data under the fitted model
# and then comparing these to the observed dataâ???. Posterior predictive checks can be used to â???olook
# for systematic discrepancies between real and simulated dataâ???.
# Model using the original (non-transformed) data:
check_posterior_predictions(regfit2.1)
## Warning: Maximum value of original data is not included in the replicated data.
## Model may not capture the variation of the data.
check_posterior_predictions(regfit2.1) #yrep = Y=immer und immer wieder geschätztes Modell Y = der Wirklichkeit # kein gute Modell, negative Gehälter anzugeben macht keinen Sinn
## Warning: Maximum value of original data is not included in the replicated data.
## Model may not capture the variation of the data.
# after model refinements:
check_posterior_predictions(regfit2.2) # response transformation # über und unterschätzung von Werten in Modell 2
## Warning: Minimum value of original data is not included in the replicated data.
## Model may not capture the variation of the data.
check_posterior_predictions(regfit2.3) # Log TRansformation der Prädiktoren die Kurve hinten ist besser geschätzt, trotzdem gibt es wieder über und unterschätzung
## Warning: Minimum value of original data is not included in the replicated data.
## Model may not capture the variation of the data.
check_posterior_predictions(regfit2.4) # outlier #outlier wird entfernt darum gibt es keine Warnmeldung mehr
#R² kann nicht unbedingt weiterhelfen