Category Archives: Project Work

VP-PA: Korrelation zweier Parameter mit R






Korrelation

Der vorliegende Datensatz enthält Daten zur Bodenfeuchte und zum Chlorophyll-
Gehalt. Die Daten wurden auf Flächen von 1 Quadratmeter in Tibet aufgenommen.
Die Leitfrage dieses Beispiels ist:
//Korreliert die Bodenfeuchte mit dem Chlorophyllgehalt auf den Flächen?//

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

Einlesen des Datensatztes

Bevor wir die Frage mit Hilfe eines statistischen Tests 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")

data <- read.table("korrelation.csv", header = TRUE, sep = ";", stringsAsFactor = FALSE)

Jetzt da die Daten in den Arbeitsspeicher von R geladen sind, kann sich eine
z. b. die Struktur des Datensatzes oder eine kurze Zusammenfassung anzeigen
lassen. Der Zugriff auf den Datensatz erfolgt über den Variablennamen
(hier: data).

str(data)
## 'data.frame':    31 obs. of  3 variables:
##  $ Plot        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ SoilMoisture: num  28.37 8.75 14.8 36 12.84 ...
##  $ Chlorophyll : num  8.5 12.9 12.9 13.9 14.2 14.2 15.2 15.9 16.9 17.4 ...
summary(data)
##       Plot       SoilMoisture    Chlorophyll  
##  Min.   : 1.0   Min.   : 0.18   Min.   : 8.5  
##  1st Qu.: 8.5   1st Qu.: 4.61   1st Qu.:16.4  
##  Median :16.0   Median : 9.30   Median :20.8  
##  Mean   :16.0   Mean   :13.65   Mean   :20.5  
##  3rd Qu.:23.5   3rd Qu.:17.81   3rd Qu.:23.8  
##  Max.   :31.0   Max.   :49.85   Max.   :32.6

Schnelle Visualisierung des Datensatzes

Es gibt viele Möglichkeiten, Daten in R zu visualisiern. In diesem Beispiel
verwenden wir das lattice-Paket, dass zunächst importiert werden muss.
Anschließend lassen wir uns die Variablen in einem Streudiagramm anzeigen,
um einen Eindruck über die Abhägnigkeit zu bekommen. Die Angabe von
Titel und Axenbeschriftungen ist optional (und sollte bei Unsicherheiten darüber,
was auf der X- und was auf der Y-Achse dargestellt wird, unterbleiben).

library(lattice)
## Warning: package 'lattice' was built under R version 2.15.3

my.xyplot <- xyplot(SoilMoisture ~ Chlorophyll, data = data, main = "Bodenfeuchte als Funktion des Chlorophyllgehalts", 
    xlab = "Chlorophyllgehalt", ylab = "Bodenfeuchte", col = "red")
print(my.xyplot)

plot of chunk unnamed-chunk-3

Auch mit viel Phantasie sieht man hier keinen sehr guten, linearen Zusammenhang.
Aber jetzt rechnen wir es einfach trotzdem einmal.

Korrelation und Regression

Für die Berechnung der Korrelation kann man verschiedene Koeffizienten
verwenden. Das erste, was hier häufig genannt wird, ist der
Korrelationskoeffizient nach Pearson.

cor.pearson <- cor.test(data$Chlorophyll, data$SoilMoisture, method = "pearson")
cor.pearson
## 
##  Pearson's product-moment correlation
## 
## data:  data$Chlorophyll and data$SoilMoisture 
## t = -1.47, df = 29, p-value = 0.1523
## alternative hypothesis: true correlation is not equal to 0 
## 95 percent confidence interval:
##  -0.5650  0.1003 
## sample estimates:
##     cor 
## -0.2634

Das Problem dabei ist, dass der Korrelationskoeffizient nach Pearson nur lineare
Zusammenhänge untersucht. Ggf. ist ein verteilungsfreier Koeffizient, z. B.
der von Spearmann oder von Kendall, hier die bessere Wahl. Der grundsätzliche
Befehl hierfür bleibt gleich, lediglich die Methode wird verändert.

cor.spearman <- cor(data$Chlorophyll, data$SoilMoisture, method = "spearman")
cor.spearman
## [1] -0.3061

Wie man sieht, ist der Spearman-Koeffizient etwas größer, was bei der Verteilung
der Daten (wie oben visualisiert) nicht wundert.

Kommen wir jetzt zur Regression.

data.lm <- lm(SoilMoisture ~ Chlorophyll, data = data)
summary(data.lm)
## 
## Call:
## lm(formula = SoilMoisture ~ Chlorophyll, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15.26  -8.84  -3.39   3.86  37.30 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   25.931      8.639    3.00   0.0055 **
## Chlorophyll   -0.600      0.408   -1.47   0.1523   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Residual standard error: 12.3 on 29 degrees of freedom
## Multiple R-squared: 0.0694,  Adjusted R-squared: 0.0373 
## F-statistic: 2.16 on 1 and 29 DF,  p-value: 0.152

Wie bei der Pearson-Korrelation zu erwarten (ja mathematisch zwingend), ist die
Regression mi teinem R-Quadrat von 0.069 nur sehr schwach. Den Verlauf der
Regressionsgerade kann man sich trotzdem noch einmal anschauen. Hierzu wird der
obige Plot um ein weiteres “Panel” ergänzt. Das abline-Panel erzeugt eine
Gerade, die dem Ergebnis der Regerssion (gespeichert in der Variablen data.lm)
folgt.

my.xyplot <- xyplot(SoilMoisture ~ Chlorophyll, data = data, main = "Bodenfeuchte als Funktion des Chlorophyllgehalts", 
    xlab = "Chlorophyllgehalt", ylab = "Bodenfeuchte", col = "red", panel = function(x, 
        y, ...) {
        panel.xyplot(x, y, ...)
        panel.abline(data.lm, col = "blue")
    })
print(my.xyplot)

plot of chunk unnamed-chunk-7

Bleibt noch ein Blick in die Residuenverteilung, insbesondere deren
Normalverteilung. Hierzu bietet das lm-Objekt 4 vordefinierte Visualisierungen.
Wir schauen uns nur die ersten beiden an.

plot(data.lm, which = c(1:2))

plot of chunk unnamed-chunk-8 plot of chunk unnamed-chunk-8

Die Normal-Q-Q-Visualisierung zeigt, dass insbesondere die Datenpunkte mit der
ID 14, 18 und 22 von zur Abweichung der Residuen von der Normalverteilung führen.
Bei nicht-normalverteilten Residuen kann man sich auf die Angaben zum
Korrelationskoeffizienten nicht verlassen.

VP-PA: Sonnenstunden mit R






Sonnenscheindauer 1900-1910 vs. 1990-1999

Der vorliegende Datensatz enthält tägliche Messungen der Sonnenscheindauer in
Potsdam zwischen 1900 und 1909 sowie 1990 und 1999. Die Leitfrage dieses
Beispiels ist:

//Unterscheidet sich die mittelere Sonnenscheindauer der Dekade 1900 bis 1909
statistisch signifikant von der Dekate 1900 bis 1999.//

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

Einlesen des Datensatztes

Bevor wir die Frage mit Hilfe eines statistischen Tests 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")

ssd <- read.table("sonnenscheindauer.csv", header = TRUE, sep = ";", stringsAsFactor = FALSE)

Jetzt da die Daten in den Arbeitsspeicher von R geladen sind, kann sich eine
z. b. die Struktur des Datensatzes oder eine kurze Zusammenfassung anzeigen
lassen. Der Zugriff auf den Datensatz erfolgt über den Variablennamen
(hier: ssd).

str(ssd)
## 'data.frame':    7304 obs. of  8 variables:
##  $ Jahr             : int  1900 1900 1900 1900 1900 1900 1900 1900 1900 1900 ...
##  $ Monat            : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Tag              : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Sonnenscheindauer: num  1.9 0 0 4.1 0 0 1.3 0 0.2 2.1 ...
##  $ X                : logi  NA NA NA NA NA NA ...
##  $ X.1              : logi  NA NA NA NA NA NA ...
##  $ X.2              : logi  NA NA NA NA NA NA ...
##  $ X.3              : logi  NA NA NA NA NA NA ...
summary(ssd)
##       Jahr          Monat            Tag       Sonnenscheindauer
##  Min.   :1900   Min.   : 1.00   Min.   : 1.0   Min.   : 0.00    
##  1st Qu.:1905   1st Qu.: 4.00   1st Qu.: 8.0   1st Qu.: 0.30    
##  Median :1950   Median : 7.00   Median :16.0   Median : 3.60    
##  Mean   :1950   Mean   : 6.52   Mean   :15.7   Mean   : 4.51    
##  3rd Qu.:1994   3rd Qu.:10.00   3rd Qu.:23.0   3rd Qu.: 7.60    
##  Max.   :1999   Max.   :12.00   Max.   :31.0   Max.   :16.10    
##     X             X.1            X.2            X.3         
##  Mode:logical   Mode:logical   Mode:logical   Mode:logical  
##  NA's:7304      NA's:7304      NA's:7304      NA's:7304     
##                                                             
##                                                             
##                                                             
## 

Vorbereiten des Datensatztes

Im Unterschied zur ebenfalls vorliegenden Excel-Tabelle ist der Datensatz nicht
in Spaltengruppen untergliedert. Für das weitere Vorgehen wäre aber eine
Gruppeninformation (Dekade x, Dekade y) hilfreich. Eine Möglichkeit, dem
Datensatz eine Gruppeninformation hinzuzufügen ist, eine neue Spalte zu
definieren, in die der Wert 0 geschrieben wird, wenn das Datum zwischen 1900
und 1909 liegt und der Wert 1 geschrieben wird, wenn das Datum zwischen 1990
und 1999 liegt. Vereinfachend kann man auch einfach das Jahr 1910 als Trennjahr
für die Gruppen angeben.

ssd$Gruppe <- ifelse(ssd$Jahr < 1910, "1900-1909", "1990-1999")
str(ssd)
## 'data.frame':    7304 obs. of  9 variables:
##  $ Jahr             : int  1900 1900 1900 1900 1900 1900 1900 1900 1900 1900 ...
##  $ Monat            : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Tag              : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Sonnenscheindauer: num  1.9 0 0 4.1 0 0 1.3 0 0.2 2.1 ...
##  $ X                : logi  NA NA NA NA NA NA ...
##  $ X.1              : logi  NA NA NA NA NA NA ...
##  $ X.2              : logi  NA NA NA NA NA NA ...
##  $ X.3              : logi  NA NA NA NA NA NA ...
##  $ Gruppe           : chr  "1900-1909" "1900-1909" "1900-1909" "1900-1909" ...

Schnelle Visualisierung des Datensatzes

Es gibt viele Möglichkeiten, Daten in R zu visualisiern. In diesem Beispiel
verwenden wir das lattice-Paket, dass zunächst importiert werden muss.
Anschließend lassen wir uns die Verteilung der Sonnenscheindauer getrennt nach
den beiden Dekaden (also getrennt nach 0 und 1 in der Spalte “Gruppe”) als
Grafik anzeigen.

library(lattice)
## Warning: package 'lattice' was built under R version 2.15.3

my.densityplot <- densityplot(~Sonnenscheindauer | Gruppe, data = ssd, plot.points = FALSE, 
    main = "Sonnenscheindauer in Potsdam")
print(my.densityplot)

plot of chunk unnamed-chunk-4

Sieht ja doch sehr ähnlich aus. Aber ist das tatsächlich so?

T-Test

Zur Überprüfung, ob die mittlere Sonnenscheindauer der beiden Dekaden
statistisch signifikant unterschiedlich ist, bietet sich ein T-Test an.

t <- t.test(Sonnenscheindauer ~ Gruppe, data = ssd, conf.level = 0.99)
print(t)
## 
##  Welch Two Sample t-test
## 
## data:  Sonnenscheindauer by Gruppe 
## t = -5.414, df = 7285, p-value = 6.371e-08
## alternative hypothesis: true difference in means is not equal to 0 
## 99 percent confidence interval:
##  -0.7974 -0.2831 
## sample estimates:
## mean in group 1900-1909 mean in group 1990-1999 
##                   4.245                   4.785

Wie man sieht, sind die Mittelwerte der beiden Gruppen mit einer
Irrtumswahrscheinlichkeit von lediglich 0.000000063 verschieden.

Visualisierung der T-Test-Ergebnisse in der Grafik

Zu guter letzt noch einmal die Grafik von oben, aber diesmal mit der Information
aus dem T-Test.

library(lattice)

my.densityplot <- densityplot(~Sonnenscheindauer | Gruppe, data = ssd, plot.points = FALSE, 
    main = paste("Sonnenscheindauer der Dekaden ist ", "unterschiedlich mit p = ", 
        t$p.value), panel = function(x, ...) {
        panel.densityplot(x, ...)
        panel.text(labels = paste(attributes(t$estimate[panel.number()]), round(t$estimate[panel.number()], 
            2), sep = ": "), x = 7, y = 0)
    })
print(my.densityplot)

plot of chunk unnamed-chunk-6