if (!require("pacman")) install.packages("pacman")
## Loading required package: pacman
## Warning: package 'pacman' was built under R version 4.0.5
if (!require("devtools")) install.packages("devtools")
## Loading required package: devtools
## Warning: package 'devtools' was built under R version 4.0.5
## Loading required package: usethis
## Warning: package 'usethis' was built under R version 4.0.5
pacman::p_load(DAAG, # für Datensatz Mireault.sav
pastecs,# für z.B. stat.desc
dplyr) # praktische Befehle wie pipes (%>%)
library(haven)
## Warning: package 'haven' was built under R version 4.0.5
Mireault <- read_sav("C:/Users/Mina/Downloads/Mireault.sav")
View(Mireault)
tinytex::install_tinytex()
pastecs::stat.desc(Mireault[,c("AgeAtLos", "SomT", "ObsessT", "DepressT", "PVTotal", "PVLoss", "SuppTotl")]) # nur für numerische Variablen
## AgeAtLos SomT ObsessT DepressT PVTotal
## nbr.val 140.0000000 3.750000e+02 3.750000e+02 3.750000e+02 3.810000e+02
## nbr.null 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## nbr.na 241.0000000 6.000000e+00 6.000000e+00 6.000000e+00 0.000000e+00
## min 1.0000000 4.100000e+01 3.800000e+01 4.200000e+01 3.200000e+01
## max 20.0000000 8.000000e+01 8.000000e+01 8.000000e+01 1.480000e+02
## range 19.0000000 3.900000e+01 4.200000e+01 3.800000e+01 1.160000e+02
## sum 1667.0000000 2.048300e+04 2.293900e+04 2.285200e+04 3.203500e+04
## median 13.0000000 5.300000e+01 6.300000e+01 6.000000e+01 8.100000e+01
## mean 11.9071429 5.462133e+01 6.117067e+01 6.093867e+01 8.408136e+01
## SE.mean 0.3782222 5.200866e-01 5.008448e-01 4.908110e-01 1.159871e+00
## CI.mean.0.95 0.7478125 1.022660e+00 9.848247e-01 9.650951e-01 2.280568e+00
## var 20.0272867 1.014338e+02 9.406705e+01 9.033580e+01 5.125592e+02
## std.dev 4.4751857 1.007143e+01 9.698817e+00 9.504515e+00 2.263977e+01
## coef.var 0.3758404 1.843864e-01 1.585534e-01 1.559685e-01 2.692602e-01
## PVLoss SuppTotl
## nbr.val 381.0000000 3.810000e+02
## nbr.null 0.0000000 0.000000e+00
## nbr.na 0.0000000 0.000000e+00
## min 6.0000000 1.200000e+01
## max 42.0000000 1.080000e+02
## range 36.0000000 9.600000e+01
## sum 7333.0000000 2.550600e+04
## median 18.0000000 6.800000e+01
## mean 19.2467192 6.694488e+01
## SE.mean 0.3285134 6.350159e-01
## CI.mean.0.95 0.6459317 1.248585e+00
## var 41.1179168 1.536364e+02
## std.dev 6.4123254 1.239502e+01
## coef.var 0.3331646 1.851526e-01
Lassen Sie sich für die sieben oben aufgeführten Variablen jeweils deskriptive Statistiken ausgeben (mindestens: Minimum, Maximum, Mittelwert, Standardabweichung) und fügen Sie den Output in Ihre Lösung ein. Begründen Sie, ob alle Werte in den Variablen AgeAtLos und DepressT (auf Basis der aus der Aufgabenstellung verfügbaren Informationen) plausibel sind.
Nur Studierende, deshalb so junger Datensatz mit Todesfall in engerer und jüngerer Altersspanne. DepressT mit BSI gemessen, höhere Werte, weil Studierende allgemein psychologische Anpassung bei Eintritt in akademisches Umfeld gefordert, Erfolgreiches Studium, vielleicht auch stichprobenspezifisch
library(DAAG)
attach(Mireault) # alle nachfolgenden Befehle beziehen sich auf diesen Datensatz
Mireault$PVLoss_z <- scale(Mireault$PVLoss) # Berechnet z-standardisierte werte
histogram(Mireault$PVLoss_z) # -3.29 als Cutoffwert für Ausreißer
Führen Sie für die Variable PVLoss eine z-Standardisierung durch und plotten Sie die Verteilung dieser z-standardisierten Variable als Histogramm. Fügen Sie das Histogramm in Ihre Lösung ein. Finden Sie visuelle Hinweise auf univariate Ausreißer? Hinweis: Für die Erstellung des Histogramms muss das Paket DAAG geladen werden.
Ab plus/minus 3 Standardabweichungen Ausreißer, scheint einen zu geben mit Abweichung zwischen 3-4 Standardabweichungen
Mireault$PVLoss_ar <- ifelse(Mireault$PVLoss_z > 3.29 | Mireault$PVLoss_z < -3.29, 1, 0) # neue Variable für Ausreißer: 1 = Ausreißer, 2 = kein Ausreißer
table(Mireault$PVLoss_ar) # Tabelle mit Häufigkeit von Ausreißern
##
## 0 1
## 380 1
#-> weiteres Vorgehen: neuen Datensatz erstellen ohne Ausreißer (ausschließen anhand der neuen Variable, z.B. mit select) & dann Analysen sowohl mit also auch ohne rechnen
Prüfen Sie für die Variable PVLoss, ob in dieser univariate Ausreißer vorliegen (gemäß der im Tutorial gezeigten Konvention). Falls Ausreißer vorliegen: Wie viele solcher Ausreißer existieren? 1 Ausreißer ist vorhanden
Mireault.multiAR <- Mireault %>% # neuer Datensatz nur mit für Analyse relevante Variablen
select(SomT, ObsessT, DepressT, SuppTotl)
Mireault.multiAR$mahal <- mahalanobis(Mireault.multiAR, # Mahalanobisdistanz berechnen
colMeans(Mireault.multiAR, na.rm = T), # nur vollständige Mittelwerte nutzen
cov(Mireault.multiAR, use = "pairwise.complete.obs")) # bei Berechnung der Covarianzen nur vollständige Werte nutzen
hist(Mireault.multiAR$mahal) # Histogramm
qchisq( # Spuckt Cutoff Wert aus
p = 0.999, # p = 1 - alpha (alpha = 0.001)
df = 4) # df = Anzahl an Variablen im Mahalanobis Datansatz
## [1] 18.46683
Mireault.multiAR[Mireault.multiAR$mahal > 18.46683, ] # Welche Werte sind Ausreißer (liegen über dem im Chi2 Test errechneten Wert)
## # A tibble: 10 x 5
## SomT ObsessT DepressT SuppTotl mahal
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NA NA NA NA NA
## 2 NA NA NA NA NA
## 3 NA NA NA NA NA
## 4 68 65 57 12 24.7
## 5 NA NA NA NA NA
## 6 52 61 54 13 22.9
## 7 NA NA NA NA NA
## 8 NA NA NA NA NA
## 9 80 74 80 18 20.2
## 10 79 38 42 12 50.9
Prüfen Sie die vier Variablen SomT, ObsessT, DepressT und SuppTotl auf (gemeinsame) multivariate Ausreißer, indem Sie deren Mahalonobis-Distanz berechnen. Fügen Sie ein Histogramm der Verteilung dieser Mahalonobis-Distanz in Ihre Lösung ein. Berechnen Sie den kritischen Wert der Mahalanobis-Distanz für ein Set von 4 Variablen (α = .001) und bestimmen Sie die Fälle (mit Angabe der Fallnummer), deren Distanz über diesem kritischen Wert liegt.
pastecs::stat.desc(Mireault[,"PVTotal"], norm = T)
## PVTotal
## nbr.val 381.0000000
## nbr.null 0.0000000
## nbr.na 0.0000000
## min 32.0000000
## max 148.0000000
## range 116.0000000
## sum 32035.0000000
## median 81.0000000
## mean 84.0813648
## SE.mean 1.1598705
## CI.mean.0.95 2.2805680
## var 512.5591518
## std.dev 22.6397693
## coef.var 0.2692602
## skewness 0.3027643
## skew.2SE 1.2110488
## kurtosis -0.2671276
## kurt.2SE -0.5356263
## normtest.W 0.9901085
## normtest.p 0.0114556
# normtest.p wichtig: ist es signifikant auf 0.05?
# -> Test allerdings fast immer signifikant, also Schiefe + Kurtosis angucken
# je größer Schiefe/skewness (Grenzwert: +/- 2) & Kortosis (Grenzwert: +/-7), desto stärker Abweichung von Normalverteilung
attach(Mireault)
## The following objects are masked from Mireault (pos = 3):
##
## AgeAtLos, AnxT, College, DepressT, filter_$, Gender, GPA, Group,
## GSIT, HostT, ID, LostPGen, ObsessT, ParT, PhobT, PsyT, PVLoss,
## PVTotal, SensitT, SomT, SuppTotl, YearColl
hist(PVTotal) # histogram für variable emotional
hist(rnorm(95, mean = 19.24, sd = 3.72)) # histogram einer normalverteilten variable mit dem selben MW und SD wie emotional als Vergleich
pastecs::stat.desc(Mireault[,"PVLoss"], norm = T)
## PVLoss
## nbr.val 3.810000e+02
## nbr.null 0.000000e+00
## nbr.na 0.000000e+00
## min 6.000000e+00
## max 4.200000e+01
## range 3.600000e+01
## sum 7.333000e+03
## median 1.800000e+01
## mean 1.924672e+01
## SE.mean 3.285134e-01
## CI.mean.0.95 6.459317e-01
## var 4.111792e+01
## std.dev 6.412325e+00
## coef.var 3.331646e-01
## skewness 5.277731e-01
## skew.2SE 2.111078e+00
## kurtosis -3.473164e-02
## kurt.2SE -6.964154e-02
## normtest.W 9.753901e-01
## normtest.p 4.522550e-06
# normtest.p wichtig: ist es signifikant auf 0.05?
# -> Test allerdings fast immer signifikant, also Schiefe + Kurtosis angucken
# je größer Schiefe/skewness (Grenzwert: +/- 2) & Kortosis (Grenzwert: +/-7), desto stärker Abweichung von Normalverteilung
attach(Mireault)
## The following objects are masked from Mireault (pos = 3):
##
## AgeAtLos, AnxT, College, DepressT, filter_$, Gender, GPA, Group,
## GSIT, HostT, ID, LostPGen, ObsessT, ParT, PhobT, PsyT, PVLoss,
## PVLoss_ar, PVLoss_z, PVTotal, SensitT, SomT, SuppTotl, YearColl
## The following objects are masked from Mireault (pos = 4):
##
## AgeAtLos, AnxT, College, DepressT, filter_$, Gender, GPA, Group,
## GSIT, HostT, ID, LostPGen, ObsessT, ParT, PhobT, PsyT, PVLoss,
## PVTotal, SensitT, SomT, SuppTotl, YearColl
hist(PVLoss) # histogram für variable emotional
hist(rnorm(95, mean = 19.24, sd = 3.72)) # histogram einer normalverteilten variable mit dem selben MW und SD wie emotional als Vergleich
Prüfen Sie die Variablen PVTotal und PVLoss jeweils auf Normalverteilung: Geben Sie den p-Wert für den Shapiro-Wilk-Test, die Schiefe und den Kurtosis-Exzess an und interpretieren Sie diese Werte. Fügen Sie die Histogramme für die Verteilungen der Variablen PVTotal und PVLoss in Ihre Lösung ein. Wie lassen sich die Verteilungen beschreiben?
Mireault$PVLoss_z <-sqrt(Mireault$PVLoss)
# Quadratwurzel
hist(Mireault$PVLoss_z)
# erneut Deskriptiv als Vergleich
pastecs::stat.desc(Mireault[,"PVLoss"], norm = T)
## PVLoss
## nbr.val 3.810000e+02
## nbr.null 0.000000e+00
## nbr.na 0.000000e+00
## min 6.000000e+00
## max 4.200000e+01
## range 3.600000e+01
## sum 7.333000e+03
## median 1.800000e+01
## mean 1.924672e+01
## SE.mean 3.285134e-01
## CI.mean.0.95 6.459317e-01
## var 4.111792e+01
## std.dev 6.412325e+00
## coef.var 3.331646e-01
## skewness 5.277731e-01
## skew.2SE 2.111078e+00
## kurtosis -3.473164e-02
## kurt.2SE -6.964154e-02
## normtest.W 9.753901e-01
## normtest.p 4.522550e-06
Führen Sie für eine selbstgewählte Variable (PVTotal oder PVLoss) eine sinnvolle Transformation durch (vgl. Tabachnik & Fidel, 2014). Vergleichen Sie Schiefe und KurtosisExzess der transformierten Variablen mit Schiefe und Kurtosis-Exzess der Ursprungsvariable. Ist die Annahme der Normalverteilung nun erfüllt?
Was müssen Sie nach der von Ihnen in 6) durchgeführten Transformation bezüglich der Interpretation der transformierten Variablen und der Rangreihenfolge der Werte beachten? (Hinweis: keine Berechungen nötig)
colSums(is.na(Mireault))
## ID Group Gender YearColl College GPA LostPGen AgeAtLos
## 0 0 0 83 81 12 241 241
## SomT ObsessT SensitT DepressT AnxT HostT PhobT ParT
## 6 6 6 6 6 6 6 6
## PsyT GSIT PVTotal PVLoss SuppTotl filter_$ PVLoss_z
## 6 6 0 0 0 0 0 0
Gibt es in den sieben oben aufgeführten Variablen fehlende Werte (missings)? Wenn ja: Wie viele fehlende Werte gibt es jeweils in welchen Variablen?