Datensatz "Rothirsch"
Analyse zum Datensatz “Hirsch”
Martin Hofbauer, Friedrich Leisch
Dieser Datensatz wurde im Zuge der Lehrveranstaltung “Vertiefung in statistische Methoden der Wildtierforschung” an der Universität für Bodenkultur Wien im Wintersemester 2018/19 von den Studierenden Lukas Graf, Jörg Fabian Knufinke und Merit Finia Pokriefke zur Analyse eingereicht.
Zu Beginn lesen wir den vorliegenden Datensatz als csv-file ein. Mit dem hier angeführten Befehl funktioniert das nur, wenn sich die Datei direkt im aktuellen Arbeitsverzeichnis (“Working Directory”) befindet. Zum Festlegen des Arbeitsverzeichnisses kann man in R den Befehl setwd() verwenden.
R> hirsch <- read.csv("hirsch.csv")
R> summary(hirsch) collar sex sl t2
c1828:7211 female:7618 Min. : 0.00 2017-06-07T14:00:42Z: 2
c3555:7618 male :7211 1st Qu.: 15.65 2017-06-07T20:00:16Z: 2
Median : 56.29 2017-06-09T17:00:42Z: 2
Mean : 124.20 2017-06-16T18:00:36Z: 2
3rd Qu.: 155.77 2017-07-04T11:00:48Z: 2
Max. :2952.90 2017-07-09T04:00:42Z: 2
(Other) :14817 Aufbereitung der Daten
Wie aus der Zusammenfassung ersichtlich, handelt es sich um einen Datensatz mit 7618 Beobachtungen des weiblichen Individuums mit dem Halsband “c3555” und 7211 Beobachtungen des männlichen Individuums mit dem Halsband “c1828” (in Summe also 14829 Beobachtungen).
Die Spalte sl gibt dabei die zurückgelegte Distanz im Zeitblock t2 an. Zum leichteren Verständnis, benennen wir sl in distance um (optional).
R> colnames(hirsch)[3] <- "distance"Die Spalte t2 ist hier im Format year-month-dayThour:minute:secondZ angegeben. Um dieses unhandliche Format zur weiteren Analyse in ein R-Datumsformat umzuwandeln, verwenden wir das Paket lubridate und die zum Ausgangsformat passende Funktion ymd_hms(). Um den originalen Datensatz nicht zu überschreiben, speichern wir das passende Datumsformat in eine neue Spalte date.
R> # install.packages("lubridate") (falls das package noch nicht installiert wurde)
R> library("lubridate")
Attaching package: 'lubridate'The following object is masked from 'package:base':
dateR> hirsch$date <- ymd_hms(hirsch$t2)Zur Analyse des Aktivitätsverlaufs innerhalb eines Tages und im Laufe eines Jahres, erstellen wir nun noch zwei weitere Variablen. Eine Variable hour, die ungeachtet des Tages die aktuelle Stunde von 0 bis 23 angibt und eine Variable yday, die jeden Tag des Jahres von 1 bis 365 durchnummeriert, ausgibt. Dabei ist die Verwendung der gleichnamigen Funktionen yday() bzw. hour(), beide wieder aus dem Package lubridate, sinnvoll.
R> hirsch$yday <- yday(hirsch$date)
R> hirsch$hour <- hour(hirsch$date)
R> summary(hirsch) collar sex distance t2
c1828:7211 female:7618 Min. : 0.00 2017-06-07T14:00:42Z: 2
c3555:7618 male :7211 1st Qu.: 15.65 2017-06-07T20:00:16Z: 2
Median : 56.29 2017-06-09T17:00:42Z: 2
Mean : 124.20 2017-06-16T18:00:36Z: 2
3rd Qu.: 155.77 2017-07-04T11:00:48Z: 2
Max. :2952.90 2017-07-09T04:00:42Z: 2
(Other) :14817
date yday hour
Min. :2017-05-30 22:00:26 Min. : 1.0 Min. : 0.00
1st Qu.:2017-08-29 23:00:21 1st Qu.: 96.0 1st Qu.: 5.00
Median :2017-11-20 14:00:41 Median :205.0 Median :12.00
Mean :2017-11-23 10:50:39 Mean :193.4 Mean :11.53
3rd Qu.:2018-02-22 14:00:54 3rd Qu.:283.0 3rd Qu.:18.00
Max. :2018-05-18 04:01:12 Max. :365.0 Max. :23.00
Analyse des Datensatzes hinsichtlich des Vergleichs der beiden Individuen
Um einen ersten Überblick über die jeweiligen Verteilungen der zurückgelegten Distanz zu bekommen, betrachten wir die Histogramme (benötigtes Package: lattice) der Variable distance separat für die beiden Individuen bzw. die Gruppen male/female.
R> library("lattice")
R> histogram(~distance|sex, data=hirsch)
Aufgrund der extremen Rechtsschiefe der Verteilungen, empfiehlt es sich die Beobachtungen der zurückgelegten Distanzen zu logarithmieren und die Betrachtung dann mit den transformierten Daten zu wiederholen. Dabei ist zu beachten, das sich nur Beobachtungen \(\gneq 0\) logarithmieren lassen. Daher addieren wir zu jeder Beobachtung den Wert 1 und können somit auch Beobachtungen darstellen, die vor der Transformation den Wert \(0\) hatten. Durch diese Addition erhalten alle Werte, die vor der Transformation den Wert 0 hatten, auch nach der Transformation den Wert 0, da \(\log(1)=0\) gilt.
R> hirsch$logdistance <- log10(hirsch$distance+1)
R> histogram(~logdistance|sex, data=hirsch)
Die vorliegenden Histogramme sind nun leichter analysierbar und es ist bereits auf einen Blick eine gewissen Ähnlichkeit erkennbar. Zur weiteren Analyse betrachten wir die Boxplots, wieder getrennt nach Geschlecht.
R> boxplot(logdistance~sex, data=hirsch)
Auch im Boxplot ist kein wirklich großer Unterschied zwischen den beiden Gruppen ersichtlich. Um Unterschiede zwischen den beiden Gruppen mittels eines parametrischen Verfahrens zu identifizieren, würde man auf einen Zwei-Stichproben t-Test zurückgreifen. Die Ergebnisse sind hier aufgrund des großen Stichprobenumfangs aber mit Vorsicht zu genießen, da der Test selbst minimalste Unterschiede als signifikant kennzeichnet:
R> t.test(logdistance~sex, data=hirsch)
Welch Two Sample t-test
data: logdistance by sex
t = 2.8085, df = 14798, p-value = 0.004984
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.008828845 0.049628010
sample estimates:
mean in group female mean in group male
1.722298 1.693069 Analyse des Aktivitätmusters mittels additiver Modelle
Additive Modelle werden in R mit der Funktion gam() aus dem Package mgcv erstellt, dieses muss zunächst geladen werden. Anschließend erstellen wir jeweils ein additives Modell für das weibliche und das männliche Individuum, das die Bewegungsaktivität in Abhängigkeit des Tages und der Tageszeit modelliert.
R> library("mgcv")Loading required package: nlmeThis is mgcv 1.8-27. For overview type 'help("mgcv-package")'.R> gam.f <- gam(logdistance~s(yday)+s(hour), data=hirsch, subset=(sex=="female"))
R> summary(gam.f)
Family: gaussian
Link function: identity
Formula:
logdistance ~ s(yday) + s(hour)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.722298 0.006959 247.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(yday) 8.554 8.944 14.71 <2e-16 ***
s(hour) 8.790 8.987 70.88 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.0906 Deviance explained = 9.26%
GCV = 0.36984 Scale est. = 0.36895 n = 7618R> gam.m <- gam(logdistance~s(yday)+s(hour), data=hirsch, subset=(sex=="male"))
R> summary(gam.m)
Family: gaussian
Link function: identity
Formula:
logdistance ~ s(yday) + s(hour)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.693069 0.007087 238.9 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(yday) 8.814 8.990 22.99 <2e-16 ***
s(hour) 8.822 8.991 54.07 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.0877 Deviance explained = 9%
GCV = 0.36313 Scale est. = 0.36219 n = 7211R> plot(gam.f, page=1, main="female")
R> plot(gam.m, page=1, main="male")
Insgesamt ähneln sich die Grafiken der jeweiligen Effekte im Großen und Ganzen. Wir erstellen nun ein Modell, in das zwar alle Beobachtungen einfließen, das aber trotzdem noch zwischen male und female unterscheidet, wodurch auch eine bessere Vergleichbarkeit aufgrund von gleicher Skalierung der Grafiken garantiert ist.
R> gam.fm <- gam(logdistance~sex+s(yday, by=sex) + s(hour, by=sex), data=hirsch)
R> summary(gam.fm)
Family: gaussian
Link function: identity
Formula:
logdistance ~ sex + s(yday, by = sex) + s(hour, by = sex)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.723556 0.006938 248.434 < 2e-16 ***
sexmale -0.031344 0.009953 -3.149 0.00164 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(yday):sexfemale 8.558 8.945 14.84 <2e-16 ***
s(yday):sexmale 8.810 8.989 22.70 <2e-16 ***
s(hour):sexfemale 8.793 8.987 71.52 <2e-16 ***
s(hour):sexmale 8.821 8.990 53.54 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.0896 Deviance explained = 9.18%
GCV = 0.36659 Scale est. = 0.36568 n = 14829R> plot(gam.fm, pages=1)
Betrachten wir die Grafiken der Effekte der bisherigen Modelle, so fällt auf, dass die Aktivitätsmuster des weiblichen und männlichen Individuums im Verlauf eines Tages (untere Grafiken) sehr ähnlich aussehen und maximal ein geringer Unterschied besteht. Im Bezug auf den Jahresverlauf fällt auf, dass ein Anstieg des Aktivitätsmusters im Spätsommer/Herbst (Brunftzeit) beim männlichen Individuum deutlich später und steiler ausfällt als beim weiblichen. Daher macht es durchaus Sinn, den Verlauf des Aktivitätsmusters unabhängig vom Geschlecht zu modellieren und dann mit dem getrennten Modell zu vergleichen.
R> gam.b <- gam(logdistance~s(yday) + s(hour), data=hirsch)
R> summary(gam.b)
Family: gaussian
Link function: identity
Formula:
logdistance ~ s(yday) + s(hour)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.708084 0.004992 342.2 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(yday) 8.807 8.989 26.47 <2e-16 ***
s(hour) 8.899 8.997 117.28 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.0801 Deviance explained = 8.12%
GCV = 0.36994 Scale est. = 0.36947 n = 14829R> plot(gam.b, pages=1)
Aber welches Modell ist nun “besser”? Wir versuchen durch das AIC abzuschätzen welches Modell “besser” ist:
R> AIC(gam.b)[1] 27338.65R> AIC(gam.fm)[1] 27203.76Der Wert des AIC ist für das getrennte Modell besser (niedrigerer Wert). Die Anpassung durch das Modell, das zwischen den Geschlechtern unterscheidet ist also besser. Es ist aber zu beachten, dass die Verbesserung des AIC nur sehr gering ausfällt und daher fraglich ist, ob dieses komplexere Modell nun tatsächlich vorzuziehen wäre. Um einen Unterschied im Verlauf der Effektgrafiken zu bestätigen, führen wir eine Varianzanalyse mit den beiden Modellen durch:
R> anova(gam.b, gam.fm, test="F")Analysis of Deviance Table
Model 1: logdistance ~ s(yday) + s(hour)
Model 2: logdistance ~ sex + s(yday, by = sex) + s(hour, by = sex)
Resid. Df Resid. Dev Df Deviance F Pr(>F)
1 14810 5472.0
2 14791 5409.1 18.926 62.9 9.0886 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Laut der durchgeführten Analyse besteht ein signifikanter Unterschied zwischen den Abweichungen der Modelle.
Da sich der Unterschied der Individuen auf den Jahresverlauf zu beschränken scheint, erstellen wir noch ein Modell, das nur im Jahresverlauf zwischen den Geschlechtern unterscheidet, im Tagesverlauf aber von einem gemeinsamen Effekt ausgeht. Dieses neue Modell analysieren wir dann in gleicher Weise wie vorher.
R> gam.fm2 <- gam(logdistance~ sex + s(yday, by=sex) + s(hour), data=hirsch)
R> summary(gam.fm2)
Family: gaussian
Link function: identity
Formula:
logdistance ~ sex + s(yday, by = sex) + s(hour)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.723528 0.006950 247.984 < 2e-16 ***
sexmale -0.031339 0.009971 -3.143 0.00167 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(yday):sexfemale 8.558 8.945 14.71 <2e-16 ***
s(yday):sexmale 8.808 8.989 22.68 <2e-16 ***
s(hour) 8.900 8.997 118.27 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.0863 Deviance explained = 8.8%
GCV = 0.3677 Scale est. = 0.367 n = 14829R> plot(gam.fm2, pages=1)
R> AIC(gam.b)[1] 27338.65R> AIC(gam.fm)[1] 27203.76R> AIC(gam.fm2)[1] 27248.75R> anova(gam.b, gam.fm2, gam.fm, test="F")Analysis of Deviance Table
Model 1: logdistance ~ s(yday) + s(hour)
Model 2: logdistance ~ sex + s(yday, by = sex) + s(hour)
Model 3: logdistance ~ sex + s(yday, by = sex) + s(hour, by = sex)
Resid. Df Resid. Dev Df Deviance F Pr(>F)
1 14810 5472.0
2 14800 5431.9 9.9451 40.082 11.022 < 2.2e-16 ***
3 14791 5409.1 8.9810 22.818 6.948 4.76e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mit dem neuen Modell zeichnet sich ein ähnliches Bild, wie bei der analogen Analyse vorher. Der Wert des AIC liegt zwischen den Werten der vorherigen Modelle, was nicht allzu verwunderlich ist, da es sich um ein Mischmodell der beiden vorherigen Modelle handelt. Die Varianzanalyse ergibt wieder einen signifikanten Unterschied zwischen den Modellen.
Insgesamt können wir somit festhalten, dass durch das getrennte Modell bessere Ergebnisse erzielt werden, wodurch die Unterscheidung der beiden Individuen eine Rolle spielt. Da sich die Stichprobe aber nur auf jeweils ein Individuum der Gruppen male/female beschränkt, kann keine Aussage darüber getroffen werden, ob es sich um Unterschiede der Individuen oder des Geschlechts handelt!