---
title: "Chi2 Test Phosphatkonzentrationen fuer Hangklassen"
output: html_notebook
---


```{r warning=FALSE, message= FALSE, eval=FALSE, include=FALSE}
library("rgdal")
download.url <- "https://zenodo.org/record/3469889/files/supplement-10.5281zenodo.3469889.zip?download=1"
download.file(paste0(download.url), destfile = file.path(as.character(getwd()), "supplement.zip"), method = "auto", mode = "wb")

unzip(file.path(getwd(), "supplement.zip"), exdir = getwd())

settlement.phos <- readOGR(file.path(getwd(), "supplement", "spot_test_cleaned_32635.shp"))
slope <- readOGR(dsn = file.path(getwd(), "supplement", "slope_cat_32635.shp"))
settlement.phos.slope <- settlement.phos; settlement.phos.slope@data["cat"] <-  over(settlement.phos.slope, slope["cat"])
writeOGR(settlement.phos.slope, file.path(getwd(), "supplement", "spot_test_slope.shp"), driver = "ESRI Shapefile", layer = "Slope", overwrite_layer = TRUE)
```


```{r}
phos.slope.chi2 <- chisq.test(settlement.phos.slope@data$phsEstm, settlement.phos.slope@data$cat)
phos.slope.chi2
```

Theoretisch kann die Nullhypothese also verworfen werden, da der p-Wert < 0.05 ist. Allerdings sollte hier etwas genauer hingeschaut werden. 

Welche Zellen haben den groeßten Einfluss auf den P-Wert der Chi2 Statistik?

```{r}
round(phos.slope.chi2$residuals, 3)
```


Beobachtete Phosphateinschaetzungen nach Hangklasse:

```{r}
phos.slope.chi2$observed
paste("Ingsgesamt fliessen", sum(phos.slope.chi2$observed), "Beobachtungen  in die chi2 Berechnung ein.")
```
Insgesamt werden also nur 53 Proben mit einer Hangklasse assoziiert.

Da die Menge der Pav Proben im Bereich des klassifizierten Hanges sind hier natuerlich problematisch, da > 1/3 der Daten herausfaellt. 
```{r}
paste("Gesamtmenge der Proben: ", length(settlement.phos.slope@data[,1]))
paste("Proben ohne Hangkategoriezuweisung:", length(which(is.na(settlement.phos.slope@data["cat"]) == TRUE)))
```

Ausserdem sollte die erwarteten Haeufigkeiten pro Zelle > 5 sein, um die Aussagekraft einzuschaetzen 

```{r}
phos.slope.chi2$expected
```
nicht in allen Zellen erwartete Wert hoeher als 5, daher ist der chi2 Test wenig aussagekraeftig.

Versucht Korrektur:
```{r}
phos.slope.fisher <- fisher.test(settlement.phos.slope@data$phsEstm, settlement.phos.slope@data$cat)
paste0("p-Wert nach ",phos.slope.fisher$method, ": ", phos.slope.fisher$p.value)
```

Die Null-Hypothese kann also nicht verworfen werden, fuer eine genauere Untersuchung des Verhaeltnisses der Phpsphatkonzentrationen zur Hangneigung sind die Informationen nicht aussagekrae?ftg genug und es muessen weitere Faktoren fuer die Hangkategoriesierung herangezogen werden. Dann kann moeglicherweise ein Zusammenhang zwischen hoeheren Konzentrationen am Hang festgestellt und die H0 verworfen werden.

##Test mit interpolierten Werten

Krieging
```{r}
kr.res.path <- file.path(getwd(), "supplement", "Kriging_results_32635.shp")
kr.res <- readOGR(kr.res.path)

#classify
kr.res$var1_pred[which(kr.res$var1_pred >= 1 & kr.res$var1_pred < 2)] <- 1
kr.res$var1_pred[which(kr.res$var1_pred >= 2 & kr.res$var1_pred < 3)] <- 2
kr.res$var1_pred[which(kr.res$var1_pred >= 3 & kr.res$var1_pred < 4)] <- 3
kr.res$var1_pred[which(kr.res$var1_pred >= 4)] <- 4

settlement.phos.kr.slope <- kr.res; settlement.phos.kr.slope@data["cat"] <-  over(settlement.phos.kr.slope, slope["cat"])

phos.slope.kr.chi2 <- chisq.test(settlement.phos.kr.slope@data$var1_pred, settlement.phos.kr.slope@data$cat)
phos.slope.kr.chi2

round(phos.slope.kr.chi2$residuals, 3)
phos.slope.kr.chi2$observed

paste("Gesamtmenge der Proben: ", length(settlement.phos.kr.slope@data[,1]))
paste("Proben ohne Hangkategoriezuweisung:", length(which(is.na(settlement.phos.kr.slope@data["cat"]) == TRUE)))

phos.slope.kr.chi2$expected
#p-value = 0.7024
```

IDW
```{r}
#test mit interpo?
idw.res.path <- file.path(getwd(), "supplement", "IDW_results_32635.shp")
idw.res<- readOGR(dsn = idw.res.path)
#classify

idw.res$var1_pred[which(idw.res$var1_pred >= 1 & idw.res$var1_pred < 2)] <- 1
idw.res$var1_pred[which(idw.res$var1_pred >= 2 & idw.res$var1_pred < 3)] <- 2
idw.res$var1_pred[which(idw.res$var1_pred >= 3 & idw.res$var1_pred < 4)] <- 3
idw.res$var1_pred[which(idw.res$var1_pred >= 4)] <- 4

settlement.phos.idw.slope <- idw.res;settlement.phos.idw.slope@data["cat"] <-  over(settlement.phos.idw.slope, slope["cat"])

settlement.phos.idw.chi2 <- chisq.test(settlement.phos.idw.slope@data$var1_pred, settlement.phos.idw.slope@data$cat)
settlement.phos.idw.chi2

round(settlement.phos.idw.chi2$residuals, 3)
settlement.phos.idw.chi2$observed

paste("Gesamtmenge der Proben: ", length(settlement.phos.idw.slope@data[,1]))
paste("Proben ohne Hangkategoriezuweisung:", length(which(is.na(settlement.phos.idw.slope@data["cat"]) == TRUE)))

settlement.phos.idw.chi2$expected

#p-value = 0.2576
```

