Monthly Archives: February 2014

VP-PA: Konfidenzintervall von Mittelwerten mit R



Konfidenzintervall am Beispiel der mittleren Jahresniederschlagssummen

Der vorliegende Datensatz ist eine Teilmenge des Lageparameterdatensatzes und
beinhaltet die mittlere Jahressummen der Niederschläge zwischen
1961 und 1990 für einige der im Gesamtdatensatz vorhandenen Messtationen des DWD.
Die Leitfrage dieses Beispiels ist:

//Wie verhält sich die Stichprobe des Datensatzes zur Grundgesamtheit in Form des
Lageparameter-Datensatzes? //

Zur Beantwortung dieser Frage sind mehrere Schritte notwendig, die im folgenden
kurz erläutert werden.

Einlesen des Datensatztes

Bevor wir die Frage mit Hilfe einiger Befehle und Rechnunge beantworten, müssen
die Daten zunächst wieder eingelesen werden. Hierfür definieren wir zunächst das
Arbeitsverzeichnis, in dem die Datei liegt (alternativ kann natürlich auch der
Pfad + Dateinamen beim öffnen angegeben werden). Beim anschließenden Lesen der
Daten ist auf die korrekte Angabe des Trennzeichens zu achten (hier: “;”).

setwd("D:/active/vp-pa/2013_winter/seminarsitzungen/statistik")
dataset <- read.table("konfidenzintervall.csv", header = TRUE, sep = ";")
summary(dataset)
##                 Station.des.DWD  Hoehe.NN..m. Geogr..Breite..Grad.
##  ABENSBERG-SANDHARLANDEN: 1     Min.   :  1   Min.   :47.7        
##  ALBSTADT-LAUTLINGEN    : 1     1st Qu.:106   1st Qu.:49.1        
##  ASCHAU                 : 1     Median :298   Median :50.1        
##  BERLIN-FRIEDRICHSFELDE : 1     Mean   :311   Mean   :50.4        
##  BORKEN/HESSEN-GOMBETH  : 1     3rd Qu.:429   3rd Qu.:51.5        
##  BOXBERG                : 1     Max.   :805   Max.   :54.0        
##  (Other)                :24                                       
##  Geogr..Laenge..Grad. Jahresniederschlagssumme Erwartungswert
##  Min.   : 7.20        Min.   : 545             Min.   :798   
##  1st Qu.: 8.85        1st Qu.: 643             1st Qu.:798   
##  Median :11.40        Median : 722             Median :798   
##  Mean   :10.88        Mean   : 822             Mean   :798   
##  3rd Qu.:12.07        3rd Qu.: 897             3rd Qu.:798   
##  Max.   :14.60        Max.   :1994             Max.   :798   
## 

T-Test gegen Grundgesamtheit

Da es sich bei dem Datensatz um eine Teilmenge des aus dem vorherigen Beispiel
bekannten Lageparameterdatensatzes handelt, können wir letzteren als Grund-
gesamtheit betrachten und den vorliegenden Datensatz als Stichprobe.

Um die Frage zu beantworten, ob sich die Stichprobe von der Grundgesamtheit
statistisch signifikant unterscheidet, führen wir einen T-Test durch, bei dem
der Mittelwert der Niederschläge der Stichprobe gegen den Mittelwert der
Grundgesamtheit getestet wird. Letzterer ist in der Spalte “Erwartungswert”
im Datensatz enthalten.

t <- t.test(dataset$Jahresniederschlagssumme, dataset$Erwartungswert)
print(t)
## 
##  Welch Two Sample t-test
## 
## data:  dataset$Jahresniederschlagssumme and dataset$Erwartungswert 
## t = 0.4173, df = 29, p-value = 0.6795
## alternative hypothesis: true difference in means is not equal to 0 
## 95 percent confidence interval:
##  -94.13 142.39 
## sample estimates:
## mean of x mean of y 
##     821.6     797.5

Wie man sieht, sind die Mittelwerte der beiden Gruppen lediglich mit einer
Irrtumswahrscheinlichkeit von 0.6795 verschieden. Anders gesagt: die Mittelwerte
der beiden Datensätze sind statistisch nicht signifikant unterschiedlich, also
gleich.

Das ganze kann man dahingehend vereinfachen, als dass der Mittelwert der
Erwartungswerte gleich dem einzelnen Erwartungswert ist, weil in der entsprechenden
Spalte lauter gleiche Werte stehen.

t2 <- t.test(dataset$Jahresniederschlagssumme, mu = 797.5)
print(t2)
## 
##  One Sample t-test
## 
## data:  dataset$Jahresniederschlagssumme 
## t = 0.4173, df = 29, p-value = 0.6795
## alternative hypothesis: true mean is not equal to 797.5 
## 95 percent confidence interval:
##  703.4 939.9 
## sample estimates:
## mean of x 
##     821.6

Bei der Berechnung eines solchen Test mit nur einer Stichprobe erhält man
zudem das Konfidenzintervall (hier 95%) für den Mittelwert. Im aktuellen Beispiel
liegt es zwischen 703,37 und 939,89 (der Wert von mu spielt keine Rolle, solange
man nur das Konfidenzintervall des Mittelwerts berechnen möchte und nicht die
alternative Hypothese testet.)

VP-PA: Lageparameter mit R



Lageparameter am Beispiel der mittleren jährlichen Niederschläge

Der vorliegende Datensatz enthält mittlere Jahressummen der Niederschläge zwischen
1961 und 1990 für diverse Messtationen des DWD. Die Leitfrage dieses
Beispiels ist:

//Wie kann man den Datensatz hinsichtlich allgemeiner Lage- und Streuparameter
beschreiben? //

Zur Beantwortung dieser Frage sind mehrere Schritte notwendig, die im folgenden
kurz erläutert werden.

Einlesen des Datensatztes

Bevor wir die Frage mit Hilfe einiger Befehle und Rechnunge beantworten, müssen
die Daten zunächst eingelesen werden. Hierfür definieren wir zunächst das
Arbeitsverzeichnis, in dem die Datei liegt (alternativ kann natürlich auch der
Pfad + Dateinamen beim öffnen angegeben werden). Beim anschließenden Lesen der
Daten ist auf die korrekte Angabe des Trennzeichens zu achten (hier: “;”).

setwd("D:/active/vp-pa/2013_winter/seminarsitzungen/statistik")
dataset <- read.table("lageparameter.csv", header = TRUE, sep = ";")

Lage- und Streuparameter

Als Lage- und Streuparameter sollen die folgenden Kenngrößen berechnet werden:
mean: Mittelwert
mod: Modus
med: Median
qxx: Quartile
qa: Quartilsabstand
rg: Range
s: Standardabweichung
s2: Varianz
d: Mittlere Abweichung vom Mittelwert
v: Variationskoeffizient
g: Schiefe
ku: Kurtosis

Als Datengrundlage dienen die mittleren Jahressummen des Niederschlags, die in
der Spalte “Jahr” abgelegt sind.

Viele der genannten Informationen sind bereits über den Befehl summary verfügbar.

summary(dataset)
##                  STATION      Hoehe...NN       B.Grad     B.Minuten   
##  AACH,KR.KONSTANZ    :  1   Min.   :   3   Min.   :47   Min.   : 0.0  
##  AACHEN-BRAND        :  1   1st Qu.:  90   1st Qu.:48   1st Qu.:17.0  
##  AACHEN (KLAERANLAGE):  1   Median : 269   Median :50   Median :35.0  
##  AACHEN (WST)        :  1   Mean   : 308   Mean   :50   Mean   :32.2  
##  AALEN-UNTERKOCHEN   :  1   3rd Qu.: 473   3rd Qu.:51   3rd Qu.:47.0  
##  AARBERGEN-PANROD    :  1   Max.   :1045   Max.   :54   Max.   :59.0  
##  (Other)             :187                                             
##      L.Grad     L.Minuten       Jahr     
##  Min.   : 6   Min.   : 0   Min.   : 399  
##  1st Qu.: 9   1st Qu.:13   1st Qu.: 637  
##  Median :10   Median :29   Median : 747  
##  Mean   :10   Mean   :29   Mean   : 812  
##  3rd Qu.:11   3rd Qu.:45   3rd Qu.: 913  
##  Max.   :14   Max.   :59   Max.   :2187  
## 

Um diese Informationen in einzelnen Variablen zu speichern, muss direkt auf die
entsprechenden Einträge zugegriffen werden. Hierzu empfiehlt es sich, die
Zusammenfassung zum einen als Variable zu speichern, zum anderen sollte eine
Begrenzung auf die relevante Spalte erfolgen. Anschließend kann über den
Befehl str die Struktur der Variablen abgefragt und die einzelnen Parameter
mittels dieser Informationen ausgelesen werden. Offenbar sind die Informationen
in einer Tabelle mit 6 Spalten gespeichert, wobei die erste Spalte das Minimum,
die zweite das 1. Quartil usw. enthält.

summary.as.variable <- summary(dataset$Jahr)
str(summary.as.variable)
## Classes 'summaryDefault', 'table'  Named num [1:6] 399 637 747 812 913 2190
##   ..- attr(*, "names")= chr [1:6] "Min." "1st Qu." "Median" "Mean" ...
min <- summary.as.variable[1]  # equals 0. quartile
q01 <- summary.as.variable[2]
med <- summary.as.variable[3]
mean <- summary.as.variable[4]
q03 <- summary.as.variable[5]
max <- summary.as.variable[6]  # equals 4. quartile
print(min)
## Min. 
##  399
print(q01)
## 1st Qu. 
##     637
print(med)
## Median 
##    747
print(mean)
## Mean 
##  812
print(q03)
## 3rd Qu. 
##     913
print(max)
## Max. 
## 2190

Um den Quartilsabstand zwischen Q 25% (q01) und Q 75% (q03) zu berechnen, reicht
eine einfache Differenz. Gleiches gilt für die Spannweite (Range), bei dem die
Differenz zwischen Maximum und Minimum berechnet wird. Da die für die Differenz
verwendeten Variablen über ein Attributwert verfügen (den sie aus dem
Summary-Objekt geerbt haben), überschreiben wir den Attributwert anschließend
mittels attr und setzten ihn auf die Bedeutung des Werts (z. b. IQR für
Inter Quartile Range).

qa <- q03 - q01
attr(qa, "names") <- list("IQR")
print(qa)
## IQR 
## 276

rg <- max - min
attr(rg, "names") <- list("Range")
print(rg)
## Range 
##  1791

Ob man es glaubt oder nicht, im Basis-Paket von R ist keine Funktion zur
Berechnung des Modalwerts. Daher lesen wir zunächst jeden im Datensatz
vorkommenden Wert aus und bestimmen anschließend, wie häufig dieser Wert vorkommt.
Umständlich, geht aber.

unique.values <- unique(dataset$Jahr)
mod <- unique.values[which.max(tabulate(match(dataset$Jahr, unique.values)))]
attr(mod, "names") <- list("Modus")
print(mod)
## Modus 
## 740.1

Standardabweichung und Varianz sind wieder einfach, da es hier eine fertige
Funktion gibt. Zuvor muss noch die entsprechende Bibliothek (stats) geladen
werden.

library(stats)
s <- sd(dataset$Jahr)
s2 <- var(dataset$Jahr)
print(s)
## [1] 290.3
print(s2)
## [1] 84268

Für die Berechnung des Variantionskoeffizienten würde z. B. das raster-Paket
eine Funktion zur Verfügung stellen. Da sich ein umwandeln des Datensatzes dafür
aber kaum lohnt, berechnen wir diese von Hand. Ebenfalls von Hand berechnen
wir die mittlere Abweichung vom Mittelwert.

v <- s/mean * 100
attr(v, "names") <- list("Var.Coef")
print(v)
## Var.Coef 
##    35.75

d <- sum(abs(dataset$Jahr - mean))/length(dataset$Jahr)
attr(d, "names") <- list("Mean.Dif.Mean")
print(d)
## Mean.Dif.Mean 
##         201.5

Für die Schiefe und Kurtosis stehen mit dem fBasics-Paket wieder Funktionen zur
Verfügung, die nach dem Laden der Bibliothek angewendet werden können.

library(fBasics)
## Warning: package 'fBasics' was built under R version 2.15.3
## Loading required package: MASS
## Loading required package: timeDate
## Warning: package 'timeDate' was built under R version 2.15.3
## Loading required package: timeSeries
## Warning: package 'timeSeries' was built under R version 2.15.3
## 
## Attaching package: 'fBasics'
## 
## The following object(s) are masked from 'package:base':
## 
##     norm
g <- skewness(dataset$Jahr)
attr(g, "names") <- list("Skweness")
print(g)
## Skweness 
##     1.95 
## attr(,"method")
## [1] "moment"

ku <- kurtosis(dataset$Jahr, method = "excess")
attr(ku, "names") <- list("Kurtosis")
print(ku)
## Kurtosis 
##    5.168 
## attr(,"method")
## [1] "excess"

Abschließend visualisieren wir noch viele der oben berechneten Variablen in einem
Box-Whisker-Plot.

boxplot(dataset$Jahr, main = "Mittlerer Jahresniederschlag 1961 - 1990")

plot of chunk unnamed-chunk-9