Auswertung biologischer Daten

Allgemeine Statistik mit R, die Test-Methode ist noch nicht bekannt, ich habe noch keinen Plan!

Moderatoren: EDi, jogo

Antworten
FledermausR
Beiträge: 21
Registriert: Fr Mär 23, 2018 11:24 am

Auswertung biologischer Daten

Beitrag von FledermausR »

Hallo,
ich bin Masterstudentin und sitze gerade an meiner Statistik und ich brauche wirklich ganz dringend ein Feedback :cry: Ich habe u.A. untersucht, ob die Fledermäuse, die sich im Winterschlaf befinden, durch meine Arbeit im Winterquartier (Monitoring) aufwachen und es dadurch eine vermehrte Aktivität in den Lichtschranken, die am Eingang des Quartiers angebracht wurden, zu erkennen ist.
Außerdem habe ich noch andere Umweltparameter erfasst um zu Überprüfen ob die Flugaktivität auch durch andere Faktoren beeinflusst wird wie z.B. die Außentemperatur, Luftfeuchte, temperatur im Bunker, die Jahreszeit oder die Anzahl der Individuen die in einem Bunker gesessen haben.
Ich habe nun ein Model erstellt und bräuchte ganz dringend ein Feedback ob das so in Ordnung ist was ich gemacht habe.
In der Dropbox habe ich meine Dateien (einmal als txt und einmal als excel)

https://www.dropbox.com/s/tjkmjm6j20b4o ... .xlsx?dl=0
https://www.dropbox.com/s/y14zoweccru8r80/LiBa.txt?dl=0

Dazu Erklärug:
Bunker: (1,2,3,4,5,6,7,8) sind die acht verschiedenen Bunker die ich regelmäßig auf den Fledermausbesatz kontrolliert habe
Date: Jeweils das Datum des Kontrolltages (Tag vor dem Monitoring) und der Monitoringtag.
NoF: Number of flight movements .. Die Summe der Flugbewegungen in den Lichtschranken (innerhalb meines bereits definierten Zeitraumes)
Treatment: Kontrolltag oder Monitoringtag
NoI : Number of Individuals, wie Viele Tiere insgesamt an dem jeweiligen Tag in den jeweiligen Bunker gehangen haben
Temp_Bunker: Temperatur im Bunker
Hum_Bunker: Feuchtigkeit im Bunker
Temp_outside: Aussentemperatur
Season: eine von mir definierte Unterteilung in Anfangssaison des Winterschlafs und Saisonende.



Hier nun mein Skript:

Code: Alles auswählen

#Hauptfrage: Fliegen mehr Tiere am Monitoringtag durch die Lichtschranke oder am Kontrolltag (Tag davor) -
#--------------------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------------------
Dataset <- 
read.table("C:/Users/Bianca/Documents/Uni/Master/Masterarbeit/Statistik/Tabellen/LiBa.txt", 
header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)


#Data exploration
#-----------------

normalityTest(~NoF, test="shapiro.test", data=Dataset)

#	Shapiro-Wilk normality test

#data: NoF
#W = 0.77188, p-value < 2.2e-16

#Nicht normal verteilt.

scatterplotMatrix(~Hum_Bunker+NoF+NoI+Season+Temp_Bunker+Temp_outside, 
reg.line=FALSE, smooth=FALSE, spread=FALSE, span=0.5, ellipse=FALSE, 
levels=c(.5, .9), id.n=0, diagonal = 'histogram', data=Dataset)

#die Anzahl der Flugbewegungen NoF sieht poisson verteilt aus. Muss man bei solch einem Model eigentlich auch die unabhängigen Variablen auf ihre Verteilung messen ?
#Bei number of Individuals (NoI) und humidity Bunker scheint es einen Zusammenhang zu geben, das untersuche ich aber in einem zweiten model
#Temp_Bunker und Temp_ outside korrelieren (natürlich) aber um zu Testen ob es hier eine Interaktion gibt benutze ich nachher den VIF


#Model selection
------------------
LinearModel.1 <- lm(NoF ~ NoI + Season + Temp_Bunker + Temp_outside + Treatment + 
Hum_Bunker, data=Dataset)
summary(LinearModel.1)
AIC(LinearModel.1)

#Call:
#lm(formula = NoF ~ NoI + Season + Temp_Bunker + Temp_outside + 
Treatment + Hum_Bunker, data = Dataset)

#Residuals:
# Min 1Q Median 3Q Max 
#-26.103 -8.226 -2.291 5.619 65.373 

#Coefficients:
Estimate Std. Error t value Pr(>|t|) 
#(Intercept) -15.26569 11.42002 -1.337 0.182903 
#NoI 0.20971 0.04015 5.223 0.00000046 ***
#Season 7.84640 3.33686 2.351 0.019726 * 
#Temp_Bunker 3.46660 0.87401 3.966 0.000103 ***
#Temp_outside 1.00189 0.27125 3.694 0.000289 ***
#Treatment[T.Monitoring] 5.02466 1.97038 2.550 0.011558 * 
#Hum_Bunker -0.13040 0.07570 -1.723 0.086570 . 
#---
#Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Residual standard error: 13.72 on 190 degrees of freedom
# (11 observations deleted due to missingness)
#Multiple R-squared: 0.2474,	Adjusted R-squared: 0.2236 
#F-statistic: 10.41 on 6 and 190 DF, p-value: 5.714e-10


> AIC(LinearModel.1)
#[1] 1599.823

#-->Schrittweise modelselektion nach AIC hat auch keine verbesserungsvorschläge. Interaktion von Temp_outside:Temp_Bunker nicht signifikant, deshalb habe ich die Interaktion aus dem model raus gelassen.

#Suchen nach Interaktion/ Collinearität
#------------------------------------------

vif(LinearModel.1)
# NoI Season Temp_Bunker Temp_outside Treatment Hum_Bunker 
# 1.365622 2.512431 2.914203 1.233297 1.015443 1.110728

#keine Kollinearität (keines >10)

#Model validierung
#-----------------------


Dataset<- within(Dataset, {
fitted.LinearModel.1 <- fitted(LinearModel.1)
residuals.LinearModel.1 <- residuals(LinearModel.1)
rstudent.LinearModel.1 <- rstudent(LinearModel.1)
hatvalues.LinearModel.1 <- hatvalues(LinearModel.1)
cooks.distance.LinearModel.1 <- cooks.distance(LinearModel.1)
obsNumber <- 1:nrow(Dataset) 
})


scatterplot(fitted.LinearModel.1~residuals.LinearModel.1, reg.line=FALSE, 
smooth=FALSE, spread=FALSE, boxplots=FALSE, span=0.5, ellipse=FALSE, 
levels=c(.5, .9), data=Dataset)
oldpar <- par(oma=c(0,0,3,0), mfrow=c(2,2))
plot(LinearModel.1)
par(oldpar)

#Betrachtung der diagnostischen Graphen 

#Homoskedastizität: Residuals vs. fitted, ist ok oder?
#Normalverteilung: Normal Q-Q-plot, sieht nicht nach normalverteilung aus. Rechtsschief. Vielleicht wird das besser wenn ich ein glm mit einer poisson verteilung wähle?
#Cook distance auch ok, oder?

#Prüfen auf unabhängigkeit. Dafür Residuals gegen alle erklärenden Variablen geplottet. Kein Muster zu sehen.

scatterplotMatrix(~Hum_Bunker+NoF+NoI+residuals.LinearModel.1+Season+Temp_Bunker+Temp_outside, 
reg.line=FALSE, smooth=FALSE, spread=FALSE, span=0.5, ellipse=FALSE, levels=c(.5, .9), id.n=0, 
diagonal = 'density', data=Dataset)

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#Problem ist nun, dass meine Residuen nicht normal verteilt sind und ich somit die Voraussetzung eines linearen modelles nicht mehr erfülle.
#Deshalb versuche ich es jetzt mit einem glm

GLM.5 <- glm(NoF ~ NoI + Season + Temp_Bunker + Temp_outside + Treatment + 
+ Hum_Bunker, family=poisson(log), data=Dataset)

summary(GLM.5)

#Call:
#glm(formula = NoF ~ NoI + Season + Temp_Bunker + Temp_outside + 
# Treatment + Hum_Bunker, family = poisson(log), data = Dataset)

#Deviance Residuals: 
# Min 1Q Median 3Q Max 
#-6.9503 -2.4242 -0.8772 1.2364 10.8613 

#Coefficients:
# Estimate Std. Error z value Pr(>|z|) 
#(Intercept) 0.5483058 0.2491235 2.201 0.0277 * 
#NoI 0.0146728 0.0007625 19.243 < 2e-16 ***
#Season 0.5622197 0.0737185 7.627 2.41e-14 ***
#Temp_Bunker 0.2412360 0.0172221 14.007 < 2e-16 ***
#Temp_outside 0.0900575 0.0059392 15.163 < 2e-16 ***
#Treatment[T.Monitoring] 0.3418947 0.0388474 8.801 < 2e-16 ***
#Hum_Bunker -0.0108332 0.0016904 -6.409 1.47e-10 ***
#---
#Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#(Dispersion parameter for poisson family taken to be 1)

# Null deviance: 2651.5 on 196 degrees of freedom
#Residual deviance: 1785.8 on 190 degrees of freedom
# (11 observations deleted due to missingness)
#AIC: 2554.8

#Number of Fisher Scoring iterations: 5


# > exp(coef(GLM.5)) # Exponentiated coefficients
# (Intercept) NoI Season 
# 1.7303190 1.0147810 1.7545627 
# Temp_Bunker Temp_outside Treatment[T.Monitoring] 
# 1.2728213 1.0942372 1.4076120 
# Hum_Bunker 
# 0.9892252 


#Anschließend noch AIC backward/forward modellselektion gemacht aber es gab keinen besseren Vorschlag. 
#Dann Interaktionen eingebaut die ich aus einer Streuungsmatrix entnommen habe. Das Model hatte einen niedrigereren AIC als das 
#ohne Interaktionen

GLM.7 <- glm(NoF ~ NoI + Season + Temp_Bunker * Temp_outside * Hum_Bunker + Treatment, family=poisson(log), 
data=Dataset)
summary(GLM.7)

#Call:
#glm(formula = NoF ~ NoI + Season + Temp_Bunker * Temp_outside * 
# Hum_Bunker + Treatment, family = poisson(log), data = Dataset)

#Deviance Residuals: 
# Min 1Q Median 3Q Max 
#-7.308 -2.390 -0.839 1.305 11.896 

#Coefficients:
# Estimate Std. Error z value Pr(>|z|) 
#(Intercept) 2.5507609 0.4679785 5.451 5.02e-08 ***
#NoI 0.0154151 0.0007835 19.675 < 2e-16 ***
#Season 0.5836995 0.0757559 7.705 1.31e-14 ***
#Temp_Bunker -0.1409994 0.0815377 -1.729 0.0838 . 
#Temp_outside -0.1275101 0.0709328 -1.798 0.0722 . 
#Hum_Bunker -0.0374406 0.0053334 -7.020 2.22e-12 ***
#Treatment[T.Monitoring] 0.3751392 0.0396056 9.472 < 2e-16 ***
#Temp_Bunker:Temp_outside 0.0347641 0.0169447 2.052 0.0402 * 
#Temp_Bunker:Hum_Bunker 0.0048477 0.0010007 4.844 1.27e-06 ***
#Temp_outside:Hum_Bunker 0.0028211 0.0008712 3.238 0.0012 ** 
#Temp_Bunker:Temp_outside:Hum_Bunker -0.0004497 0.0002028 -2.218 0.0266 * 
#---
#Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#(Dispersion parameter for poisson family taken to be 1)

# Null deviance: 2651.5 on 196 degrees of freedom
#Residual deviance: 1742.0 on 186 degrees of freedom
# (11 observations deleted due to missingness)
#AIC: 2519

#Number of Fisher Scoring iterations: 5

vif(GLM.7)
# NoI Season 
# 1.623205 2.870362 
# Temp_Bunker Temp_outside 
# 69.868462 171.922953 
# Hum_Bunker Treatment 
# 12.750599 1.054920 
# Temp_Bunker:Temp_outside Temp_Bunker:Hum_Bunker 
# 313.505445 92.461581 
# Temp_outside:Hum_Bunker Temp_Bunker:Temp_outside:Hum_Bunker 
# 173.216313 328.841561
#Modelvalidierung
-----------------------
#Diagnostische Graphen anschauen

oldpar <- par(oma=c(0,0,3,0), mfrow=c(2,2))
plot(GLM.7)
par(oldpar)


#Homostedastizität ist erfüllt
#Normalverteilung ist nun auch erfüllt (außer die paar Ausreißer aber um diese Ausreißer ging es mir ja)
#Bei der cook distance fällt ein einziger punkt ganz leicht aus der Grenze von 1. 

scatterplotMatrix(~Hum_Bunker+NoF+NoI+residuals.LinearModel.1+Season+Temp_Bunker+Temp_outside, reg.line=FALSE, 
smooth=FALSE, spread=FALSE, span=0.5, ellipse=FALSE, levels=c(.5, .9), id.n=0, diagonal = 'histogram', 
data=Dataset)

#Prüfen auf unabhängigkeit. Dafür Residuals gegen alle erklärenden Variablen geplottet. Kein Muster zu sehen.

#Post Hoc
#----------------
#Die faktoriellen Variablen waten treatment und season. Da sich es hier aber nur um 2 Gruppen handelt die verglichen werden (Monitoring und control bzw 1 und 2)
#ist hier kein post hoc von nöten sondern der signifikante Unterschied dieser zwei Gruppen wird in einem t Test (oder Wilcox, je nach Verteilung) getestet.
#Die anderen Variablen wie Temp_Bunker, Temp_outside usw kann ich soweit ich das verstanden habe nicht mit einem post hoc untersuchen weil es sind ja keine Gruppen 
#sondern metrische Daten. Die würde ich dann eher plotten und optisch interpretieren.
Vielen Dank, ich würde mich wirklich sehr über Antwort freuen bzw wenn jemand, der Ahnung hat, mal drüber schauen könnte. :D
LG
Bianca
jogo
Beiträge: 2085
Registriert: Fr Okt 07, 2016 8:25 am

Re: Auswertung biologischer Daten

Beitrag von jogo »

Hallo Bianca,

willkommen in diesem Forum!
Für die anderen Leser noch die Anmerkung, dass ich empfohlen hatte, diese Frage hier zu posten, weil es doch ein paar mehr Spezialisten hier gibt.
In dem anderen Forum gibt es eine Antwort von Hufeisen:
http://www.r-forum.de/deskriptive-stati ... html#p2453

Gruß, Jörg
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: Auswertung biologischer Daten

Beitrag von EDi »

Man könnte NoI auch als 'offset' ins model nehmen. Das würde dann auf die NoI skalieren, also modelliert man NoF pro NoI.

Bei deinem glm hast du ziemliche "overdispersion", d.h mehr Varianz als vom Poission Model angenommen. Das ist normal mit ökologischen Daten, probiert mal ein negativ-binomial Model.
Bitte immer ein reproduzierbares Minimalbeispiel angeben. Meinungen gehören mir und geben nicht die meines Brötchengebers wieder.

Dieser Beitrag ist lizensiert unter einer CC BY 4.0 Lizenz
Bild.
FledermausR
Beiträge: 21
Registriert: Fr Mär 23, 2018 11:24 am

Re: Auswertung biologischer Daten

Beitrag von FledermausR »

ok vielen Dank :) Woran hast du die Overdispersion abgelesen?
Athomas
Beiträge: 768
Registriert: Mo Feb 26, 2018 8:19 pm

Re: Auswertung biologischer Daten

Beitrag von Athomas »

Vier Bemerkungen:

1) In einer "poisson regression" sollte der Logarithmus der Anzahl der Individuen (aus nachvollziehbaren Gründen :) ) als "offset" eingehen. Also ähnlich wie von EDi angeregt.

2) Bei einer Poisson-Verteilung liegt mit dem Mittelwert auch die Varianz fest - sie hat nämlich denselben Wert. In der Praxis kommt es häufig vor, dass die Varianz größer als der Mittelwert ist, dann spricht man von "overdispersion". Das kann viele Gründe haben, z.B. können im Modell wesentliche "erklärende" Variablen fehlen - oder aber das Poisson-Modell stimmt für jedes einzelne Individuum, aber die Größe des Verteilungsparameters variiert.
Wenn dieser Mittelwert selbst wieder eine gammaverteilte Zufallsgröße ist, entsteht aus der "Mischung" dieser Verteilungen eine negative Binomialverteilung.

Eine andere Methode, die "overdispersion" in den Griff zu bekommen, ist das "quasi-Poisson-Model", was auch im Rahmen der bekannten glm-Funktion zur Verfügung steht. "Im Prinzip" werden dadurch nicht die Schätzungen selbst, sondern die geschätzten Fehler verändert. Bei hoher Overdispersion werden die geschätzten Fehler dadurch deutlich größer - die ursprüngliche Poisson-Regression gaukelt also viel höhere Signifikanzen vor!

3) Was läßt Dich glauben, die Bunker alle in einen Topf schmeißen zu können? Dass mit Bunker 6 irgendwas anders ist, sieht man doch schnell!?
Wenn man in einem quasi-Poisson-Modell den Bunker als zusätzlichen Faktor mit aufnimmt, sinkt die Overdispersion um fast 60%...

4) Was soll der "Normalitätstest" am Anfang? Ob die Zielgröße normal (oder wie auch immer) verteilt ist, ist völlig irrelevant - die Residuen zählen!
Antworten