Statistik in der Archaeologie (2020): Beispiel-Skripte fertig fuer das Arbeiten mit dem R-Commander (Teil 3/4: Fallstudien)
Frank Siegmund
2020-08-05
Vorab:
Wegen des Wechsels von Datensaetzen, Modellen und dem Editieren der Daten etc. laeuft dieses Skript nicht “in einem Satz” ganz durch, sondern muss Stueck fuer Stueck ausgefuehrt werden.
Kap. 25.6: Regression am Beispiel Koerperehoehenschaetzung
Einlesen des Datensatzes Rollet 1888:
> setwd("D:/RData")
> #
> Rollet1888 <- readSPSS("D:/RData/5_Rollet1888.sav", rownames=FALSE,
+ stringsAsFactors=TRUE, tolower=FALSE)Aufteilen der Serie in zwei Stichproben: Referenzserie und Testserie
> Rollet1888$dummy <- with(Rollet1888, runif(100, min=0, max=10))
> #
> Rollet1888 <- with(Rollet1888, Rollet1888[order(Sex, dummy,
+ decreasing=FALSE), ])! Nun haendische Nachbearbeitung des Datensatzes, wie S. 357 beschrieben, um je 5 per Zufallszahl bestimmte Maenner und je 5 Frauen zu markieren, indem “dummy” auf “20” gesetzt wird. Anschliessend rausschreibend der beiden Teildatensaetze “Referenserie” mit 90 Faellen und “Testserie” mit 10 Faellen:
> Referenzserie <- subset(Rollet1888, subset=dummy < 15)
> #
> Testserie <- subset(Rollet1888, subset=dummy > 15)
> #
> # exemplarische Kontrolle fuer Referenzserie:
> numSummary(Referenzserie[,"F1", drop=FALSE], groups=Referenzserie$Sex,
+ statistics=c("mean", "sd", "quantiles"), quantiles=c(0,.25,.5,.75,1)) mean sd 0% 25% 50% 75% 100% F1:n
female 413.46 22.55081 372.5 398.125 412.5 425.375 478 50
male 450.43 23.73235 394.0 434.750 446.0 470.500 500 50Lineare Regression auf Basis Referenzserie anhand des Oberschenkelknochens (F1)
> RegF1 <- lm(Khoehe~F1, data=Referenzserie)
> summary(RegF1)
Call:
lm(formula = Khoehe ~ F1, data = Referenzserie)
Residuals:
Min 1Q Median 3Q Max
-10.6781 -2.4242 -0.0569 2.1813 10.6220
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 52.8077 5.6715 9.311 3.79e-15 ***
F1 0.2487 0.0131 18.986 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.857 on 98 degrees of freedom
Multiple R-squared: 0.7863, Adjusted R-squared: 0.7841
F-statistic: 360.5 on 1 and 98 DF, p-value: < 2.2e-16Ergebnis: geschaetzte Koerperhoehe = F1 mal 0.2487 plus 52.8077.
Achtung: da eingangs die Individuen, die der Referenzserie resp. der Testserie zugeordnet wurden, zufaellig ausgewaehlt wurden, werden beim Nacharbeiten die Ergebnisse jeweils leicht unterschiedlich ausfallen!
Nun ein Blick auf die Residuen = Schaetzfehler:
> Referenzserie<- within(Referenzserie, {
+ fitted.RegF1 <- fitted(RegF1)
+ residuals.RegF1 <- residuals(RegF1)
+ rstudent.RegF1 <- rstudent(RegF1)
+ hatvalues.RegF1 <- hatvalues(RegF1)
+ cooks.distance.RegF1 <- cooks.distance(RegF1)
+ })
> #
> numSummary(Referenzserie[,"residuals.RegF1", drop=FALSE],
+ statistics=c("mean", "sd", "quantiles"), quantiles=c(0,.25,.5,.75,1)) mean sd 0% 25% 50% 75% 100% n
-4.676814e-16 3.837362 -10.67805 -2.424205 -0.056882 2.181264 10.62196 100Zu den Zahlen die "grundlegenden diagnostischen Grafiken:
> oldpar <- par(oma=c(0,0,3,0), mfrow=c(2,2))
> plot(RegF1)
> par(oldpar)Das Bild “Residuals vs. Fitted” ist (zumindest bei meinen Beispielsberechnungen) unauffaellig, im Q-Q-Plot sind einige wenige Faelle markiert, die deutlicher abnweichen. Gemaess “Scale-Location” mit den standardisierten Residuen sowie anhand der Cook-Distanz sind ihre Abweichungen jedoch nicht extrem. Anschliessend wendet man die Formeln auf die Testserie an und ermittelt an ihr, wie gross die Fehler ausfallen, wenn sie jenseits der Referenzserie angewendet werden.
Das Vorgehen ist hier sehr vereinfacht umrissen. Tatsaechlich wuerde man neben dem Oberschenkelknochen auch Schaetzformeln fuer die anderen Langknochen entwickeln, also fuer H1, R1, F1 und T1, und diese nutzen. Fuer Individuen, von denen mehr als ein Langknochen verfuegbar ist, verringert sich der Schaetzfehler dadurch.
Kap. 26.6.6: Weitere Diagnose-Optionen
Wiederherstellen eines Ausgangspunktes fuer die Uebung:
> Rollet1888 <- readSPSS("D:/RData/5_Rollet1888.sav", rownames=FALSE, stringsAsFactors=TRUE, tolower=FALSE)
> #
> RegF1 <- lm(Khoehe~F1, data=Rollet1888)
> summary(RegF1)
Call:
lm(formula = Khoehe ~ F1, data = Rollet1888)
Residuals:
Min 1Q Median 3Q Max
-10.6781 -2.4242 -0.0569 2.1813 10.6220
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 52.8077 5.6715 9.311 3.79e-15 ***
F1 0.2487 0.0131 18.986 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.857 on 98 degrees of freedom
Multiple R-squared: 0.7863, Adjusted R-squared: 0.7841
F-statistic: 360.5 on 1 and 98 DF, p-value: < 2.2e-16> #
> Rollet1888<- within(Rollet1888, {
+ fitted.RegF1 <- fitted(RegF1)
+ residuals.RegF1 <- residuals(RegF1)
+ rstudent.RegF1 <- rstudent(RegF1)
+ hatvalues.RegF1 <- hatvalues(RegF1)
+ cooks.distance.RegF1 <- cooks.distance(RegF1)
+ obsNumber <- 1:nrow(Rollet1888)
+ })Bonferroni-Test auf Ausreisser:
> outlierTest(RegF1)No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferroni p
51 -2.911925 0.0044572 0.44572> #
> # generiert Abb. 115:
> influenceIndexPlot(RegF1, id=list(method="y", n=3), vars=c("Cook", "Studentized", "Bonf", "hat"))
> #
> # generiert Abb. 116:
> influencePlot(RegF1, id=list(method="noteworthy", n=3))
StudRes Hat CookD
15 -1.88556783 0.02804904 0.04999755643
33 -1.51993070 0.04888039 0.05857999452
40 2.86813449 0.01000355 0.03870721209
50 -0.04437313 0.06342830 0.00006735933
51 -2.91192464 0.02705026 0.10951412429
52 -0.65139559 0.05076449 0.01141311919
85 -2.47239202 0.01463977 0.04315764862Lineare Regression mit Resampling
Start wieder mit “frischem” Datensatz:
> Rollet1888 <- readSPSS("D:/RData/5_Rollet1888.sav", rownames=FALSE, stringsAsFactors=TRUE, tolower=FALSE)
> #
> RegF1 <- lm(Khoehe~F1, data=Rollet1888)
> summary(RegF1)
Call:
lm(formula = Khoehe ~ F1, data = Rollet1888)
Residuals:
Min 1Q Median 3Q Max
-10.6781 -2.4242 -0.0569 2.1813 10.6220
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 52.8077 5.6715 9.311 3.79e-15 ***
F1 0.2487 0.0131 18.986 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.857 on 98 degrees of freedom
Multiple R-squared: 0.7863, Adjusted R-squared: 0.7841
F-statistic: 360.5 on 1 and 98 DF, p-value: < 2.2e-16> #
> # konventionelle Konfidenzintervalle:
> library(MASS, pos=17)
> Confint(RegF1, level=0.95) Estimate 2.5 % 97.5 %
(Intercept) 52.8076803 41.5527726 64.0625881
F1 0.2487176 0.2227216 0.2747136Nun wie p. 367 beschrieben Bootstrap-Konfidenzintervalle berechnen lassen. Wir bestellen das Ziehen von 2.000 Zufallsstichproben aus unserem Datensatz, weshalb das Rechnen eine kurze Weile dauern kann:
> .bs.samples <- Boot(RegF1, R=2000, method="case")
> plotBoot(.bs.samples)
> confint(.bs.samples, level=0.95, type="bca")Bootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) 43.7351544 62.2162861
F1 0.2271214 0.2691816> remove(.bs.samples)Der akribische Vergleich der Schaetzwerte der konventionellen Konfidenzintervalle und der Bootstrapping-Intervalle zeigt: Die Bootstrap-Konfidenzintervalle fuer den Schnittpunkt mit der y-Achse (“Intercept”) und die Steigung der Regressionsgeraden sind - im Vergleich zu den Konfidenzintervallen des klassischen Vorgehens - minim enger, d.h. sie liegen naeher am eigentlichen Schaetzwert. Somit unterstreicht das Bootstrapping die (interne) Validitaet und Staibilitaet des Ergebnisses.