Analyse zum Datensatz “Hirsch”

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':

    date

R> 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: nlme

This 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 = 7618

R> 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 = 7211

R> 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 = 14829

R> 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 = 14829

R> 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.65

R> AIC(gam.fm)

[1] 27203.76

Der 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 ' ' 1

Laut 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 = 14829

R> plot(gam.fm2, pages=1)

R> AIC(gam.b)

[1] 27338.65

R> AIC(gam.fm)

[1] 27203.76

R> AIC(gam.fm2)

[1] 27248.75

R> 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 ' ' 1

Mit 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!