ich bin Masterstudentin und sitze gerade an meiner Statistik und ich brauche wirklich ganz dringend ein Feedback
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:
#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.
LG
Bianca
