1. Przygotowanie zbioru danych
Pobranie dwóch zbiorów danych:
- Pierwszego - ftp://h05-ftp.jrc.it/panagpa/LUCAS.SOIL_corr.Rdata.zip (0.5 GB) - zawierającego 19036 obserwacji i 88 zmiennych. W części zmiennych były braki wartości (mniejsza liczba N w poniższej tabeli). Uzyskanie zmienne numerycznie obejmują takie cechy jak:
| Statistic | N | Mean | St. Dev. | Min | Max |
| coarse | 19,036 | 13.159 | 10.629 | 0 | 83 |
| clay | 17,939 | 18.884 | 13.005 | 0 | 79 |
| silt | 17,939 | 38.226 | 18.300 | 0 | 92 |
| sand | 17,939 | 42.883 | 26.107 | 1 | 99 |
| pH.in.CaCl2 | 19,036 | 5.594 | 1.426 | 2.570 | 9.250 |
| pH.in.H2O | 19,036 | 6.198 | 1.353 | 3.210 | 10.080 |
| OC | 19,036 | 50.003 | 91.307 | 0.000 | 586.800 |
| CaCO3 | 19,036 | 51.601 | 125.318 | 0 | 944 |
| N | 19,036 | 2.925 | 3.755 | 0.000 | 38.600 |
| P | 19,036 | 30.076 | 32.855 | 0.000 | 1,366.400 |
| K | 19,035 | 197.052 | 229.299 | 0.000 | 7,342.000 |
| CEC | 19,036 | 15.757 | 14.485 | 0.000 | 234.000 |
| POINT_ID | 19,036 | 42,757,642.000 | 7,506,485.000 | 26,521,974 | 57,981,484 |
| OBSERVED | 19,036 | 1.002 | 0.048 | 1 | 2 |
| OBS_TYPE | 19,036 | 1.034 | 0.194 | 1 | 3 |
| GPS_PREC | 19,036 | 32.310 | 489.697 | 1 | 8,888 |
| GPS_LAT | 19,036 | 50.547 | 7.953 | 34.984 | 88.889 |
| GPS_LONG | 19,036 | 10.682 | 11.406 | -9.746 | 88.889 |
| OBS_DIST | 19,036 | 283.467 | 4,898.888 | 0 | 88,888 |
| OBS_DIRECT | 19,036 | 1.031 | 0.208 | 1 | 3 |
| OBS_RADIUS | 19,036 | 1.587 | 0.492 | 1 | 2 |
| LC1_PERCENT | 19,036 | 4.269 | 1.261 | 1 | 5 |
| LC2_PERCENT | 18,062 | 0.281 | 1.068 | 0 | 5 |
| AREA_SIZE | 19,036 | 3.023 | 0.834 | 1 | 4 |
| TREES_HEIGHT | 19,001 | 0.605 | 0.906 | 0 | 2 |
| FEATURES_WIDTH | 19,000 | 0.572 | 0.899 | 0 | 2 |
| LAND_MNGT | 19,036 | 1.926 | 0.352 | 1 | 3 |
| WM_WATER_MNGT | 19,036 | 5.702 | 1.880 | 1 | 8 |
| WM_SRC_IRRIGATION | 19,036 | 0.313 | 1.881 | 0 | 16 |
| WM_TYP_IRRIGATION | 19,036 | 0.297 | 1.777 | 0 | 16 |
| WM_DELIVERY_SYST | 19,036 | 0.226 | 1.058 | 0 | 8 |
| SOIL_SURVEY | 19,036 | 1.006 | 0.107 | 1 | 3 |
| SOIL_PLOUGH | 19,036 | 1.701 | 0.467 | 0 | 2 |
| SOIL_CROP | 19,036 | 1.598 | 1.042 | 0 | 4 |
| SOIL_STONES | 19,036 | 1.395 | 0.766 | 0 | 4 |
| X_LAEA | 19,036 | 4,275,455.000 | 750,602.800 | 2,652,000 | 5,798,000 |
| Y_LAEA | 19,036 | 3,089,891.000 | 842,917.900 | 1,464,000 | 5,214,000 |
| STRATA | 19,036 | 2.333 | 1.354 | 1 | 5 |
| STRATA2 | 18,332 | 0.173 | 0.779 | 0 | 7 |
| area0 | 19,036 | 352,862.900 | 185,687.900 | 2,596 | 638,341 |
| area1 | 19,036 | 115,072.300 | 103,440.500 | 387 | 335,767 |
| area2 | 19,036 | 40,000.290 | 36,056.480 | 387 | 165,075 |
| area3 | 19,036 | 11,095.510 | 15,762.700 | 43 | 105,862 |
| PHOTO_POINT | 19,036 | 1.020 | 0.345 | 1 | 8 |
| PHOTO_CROP | 19,036 | 1.226 | 1.160 | 1 | 8 |
| PHOTO_N | 19,036 | 1.001 | 0.035 | 1 | 2 |
| PHOTO_E | 19,036 | 1.001 | 0.037 | 1 | 2 |
| PHOTO_S | 19,036 | 1.001 | 0.035 | 1 | 2 |
| PHOTO_W | 19,036 | 1.003 | 0.050 | 1 | 2 |
| PHOTO_I | 19,036 | 7.631 | 1.562 | 1 | 8 |
| PHOTO_T | 19,036 | 1.061 | 0.363 | 1 | 8 |
| PHOTO_SOIL | 19,036 | 1.036 | 0.450 | 1 | 8 |
| PHOTO_CONFLICT | 19,036 | 1.109 | 0.845 | 0 | 2 |
| PCAREA | 19,018 | 54.129 | 15.346 | 0 | 100 |
| PARMADO1 | 19,018 | 4.595 | 1.889 | 0 | 8 |
| PARMADO2 | 19,018 | 47.761 | 18.861 | 0 | 80 |
| PARMADO3 | 19,018 | 478.577 | 188.260 | 0 | 800 |
| PARMASE1 | 19,018 | 2.622 | 2.637 | 0 | 8 |
| PARMASE2 | 19,018 | 27.474 | 27.231 | 0 | 80 |
| PARMASE3 | 19,018 | 274.339 | 272.532 | 0 | 800 |
Zbiór posiadał także zmienne nienumeryczne, między innymi takie jak:
- ID
- Date
- Notes
- country
- GPS_PROJ
- LC1
- LC2
- LU1
- LU2
- LC1_SPECIES
- LC2_SPECIES
- Transect
- CLIMANR
- MAIN_CLIMA
- WRBFU
- WRBLV1
- WRBADJ1
- MIN_TOP
- mineral
Dodatkowo, w zbiorze znajdowały się pomiary spektralne 4200 długości fali dla każdej obserwacji. Nie zawierały one jednak powtórzeń.
- Drugiego - ftp://h05-ftp.jrc.it/panagpa/LUCAS_TOPSOIL_v1_spectral.csv (3 GB) - drugi zbiór danych znajdował się w częściowo uszkodzonym pliku tekstowym. Udało się z niego odzyskać 17814 obserwacji, każda zawierająca po dwa powtórzenia dla 4200 długości fali. Niestety nie każda z tych obserwacji miała ID zgodny z tym w Pierwszej bazie danych.
2. Obliczenie RD oraz stworzenie dwóch zbiorów do analizy
RD wyliczono poprzez:
- odjęcie drugiego pomiaru od pierwszego (diff)
- wyliczenie średniej wartości (średnia z pomiar pierwszy i pomiar drugi) (rowmean)
- podzielenie wartość różnicy poprzez wartość średnią (diff/rowmean)
- wyliczenie wartości absolutnej z uzyskanego powyżej wyniku
Wydzielono dwa zbiory do analizy:
ZBIÓR A - Zawierający wartość RD dla każdej długości fali w każdej obserwacji. Przykład dla jednego z pomiarów:
| obserwacja | wavelength | powtorzenie_1 | powtorzenie_2 | diff | mean | RD |
|---|---|---|---|---|---|---|
| 1 | 400.0 | 0.7564347 | 0.7557260 | 0.0007087 | 0.7560803 | 0.0009373 |
| 1 | 400.5 | 0.7655868 | 0.7648741 | 0.0007128 | 0.7652305 | 0.0009315 |
| 1 | 401.0 | 0.7746894 | 0.7739715 | 0.0007179 | 0.7743305 | 0.0009271 |
| 1 | 401.5 | 0.7837104 | 0.7829874 | 0.0007230 | 0.7833489 | 0.0009230 |
| 1 | 402.0 | 0.7926247 | 0.7918956 | 0.0007291 | 0.7922602 | 0.0009203 |
| 1 | 402.5 | 0.8014045 | 0.8006683 | 0.0007362 | 0.8010364 | 0.0009191 |
ZBIÓR B Zawierający medianę, średnią, odchylenie standardowe, oraz 95% przedział ufności średniej dla wartości RD w każdej obserwacji (brak długości fali w zbiorze). Do tego zbioru dołączono metadane z Pierwszej bazy danych. W efekcie połączenia uzyskano 16934 obserwacje. Przykład:
| obserwacja | RD_mean | RD_median | RD_sd |
|---|---|---|---|
| 1 | 0.0012313 | 0.0012444 | 0.0004211 |
| 10 | 0.0062178 | 0.0062494 | 0.0011757 |
| 100 | 0.0016846 | 0.0019665 | 0.0005974 |
| 10002 | 0.0010547 | 0.0010784 | 0.0005268 |
| 10004 | 0.0054591 | 0.0056815 | 0.0005605 |
| 10005 | 0.0082551 | 0.0085443 | 0.0012659 |
3. Analiza względnych różnic pomiędzy pomiarami
Do analizy różnic pomiędzy parami wykorzystano ZBIÓR B.
Rozkład wartości RD
Podstwowe statystyki:
| RD_mean | RD_median | RD_sd | |
|---|---|---|---|
| Min. :0.0000712 | Min. :0.00004446 | Min. :0.00006993 | |
| 1st Qu.:0.0013685 | 1st Qu.:0.00128928 | 1st Qu.:0.00046620 | |
| Median :0.0026944 | Median :0.00270712 | Median :0.00070979 | |
| Mean :0.0037372 | Mean :0.00382124 | Mean :0.00088616 | |
| 3rd Qu.:0.0049934 | 3rd Qu.:0.00513453 | 3rd Qu.:0.00108702 | |
| Max. :0.0293527 | Max. :0.03173210 | Max. :0.01031784 |
Powyższe statystyki i wykresy pokazują, że istnieją sytuacje w których średnia wartość RD przekracza 0,01. Nie występują one jednak bardzo często. W 94.2% przypadków średnia wartość RD jest niższa niż 0,01, a więc odwrotna sytuacja występuje w 5.8% przypadków.
Sprawdzono korelację pomiędzy RD_mean a kilkoma wybranymi zmiennymi, aby określić która z nich najbardziej wpływa na wartości RD.
| Korelacja ze zmienną RD_mean | |
|---|---|
| RD_mean | 1.0000000 |
| coarse | 0.1287916 |
| clay | 0.1882575 |
| silt | 0.0106097 |
| sand | -0.1018697 |
| pH.in.CaCl2 | 0.0570795 |
| pH.in.H2O | 0.0625079 |
| OC | 0.1539108 |
| CaCO3 | 0.0952680 |
| N | 0.1834410 |
| P | -0.0494880 |
| K | 0.0811864 |
| CEC | 0.2189456 |
Wyniki korelacji przedstawiono także na wykresie.
Wsród powyższych zmiennych, żadna nie wykazała znaczacego wpływu na wartość RD_mean. Najwyższe korelacje dodanie miały zmienne CEC, clay, N i OC - wynosiły one od ok. 0,15 do 0,22. W przypadku korelacji ujemnych najniższą wartość miała zmienna sand - ale znów była ona dość niewielka (-0,10).
Następnie nastąpiło bliższe przyjrzenie się relacjom RD do zmiennej clay oraz OC.
Relacja RD do zmiennej clay
##
## Spearman's rank correlation rho
##
## data: danee$RD_mean and danee$clay
## S = 554360000000, p-value < 0.00000000000000022
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1782827Test korelacji rangowej Spearmana dał podobny wynik do powyższej korelacji liniowej Pearsona - ok. 0,18. Mimo dość niewielkiej wartości, jest ona istotna statystycznie.
Na pierwszym (“surowym”) wykresie relacji RD do clay nie można zbyt wiele powiedzieć o ich zależności. Zmienna clay ma dość skośny rozkład - większość jej wartości jest niższa niż ok. 30. Wyliczono więc logarytm ze zmiennej clay. Jego porówanenie do wartości RD daje trochę wyraźniejsze wyniki. Na ostatnim wykresie powyżej widać, że wartości RD poniżej 0,01 występują przy każdej możliwej wartości logarytmu clay. Z drugiej jednak strony - obserwacje z dużą wartością RD pojawiają się częściej wraz z wyższą wartością logarytmu clay.
Wartości zmiennej clay pogrupowano do ośmiu kategorii i porównano z wartościami RD. Tutaj również widoczny jest najniższy błąd (najniższa wartość RD) w kategoriach najmniejszej zawartości gliny, następnie błąd wzrasta do wartości gliny około 30-40, od tego momentu się stabilizujac.
Wykonano również odwrotną operację - pogrupowano wartości błędu na dwie kategorie - poniżej i powyżej 0,01. Na poniższej rycinie widać, że zawartość gliny jest wyraźnie niższa dla niższej grupy wartości RD. Wykres zawiera nacięcia - jeżeli one się na siebie nie nakładają to sugeruje to, że mediany pomiędzy kategoriami są istotnie statystycznie różne.
Relacja RD do zmiennej OC
Podobą analizę wykonano w stosunku do zmiennej OC.
##
## Spearman's rank correlation rho
##
## data: danee$RD_mean and danee$OC
## S = 661080000000, p-value < 0.00000000000000022
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1831785Tutaj również korelacja była dość niewielka (ok. 0,18), ale istotna statystycznie.
Wartości OC były bardzo skośne, więc również zostały zlogaritmizowane.
Uzyskane w ten sposób wykresy również ukazują, że wartości OC są wyższe w przypadkach większej wartości RD, jednak tutaj ta relacja nie jest tak widoczna jak dla zmiennej clay.
4. Analiza względnych różnic pomiędzy zmiennością wartości RD a długością fali
Do analizy różnic pomiędzy pomiędzy zmiennością wartości RD a długością fali wykorzystano ZBIÓR A.
Dla każdej obserwacji wyliczono zmienność RD w zależności od długości fali. Poniżej można zobaczyć przykład dla obserwacji numer 10:
Następnie takie przejścia wartości RD zostały uśrednione oraz wyliczono dla nich odchylenie standardowe. Na poniższym wykresie, czarna linia oznacza wartość średnią, a kolejne szare pasy - jedno, dwa i trzy odchylenia standardowe.
5. Modele relacji zmiennych coarse, clay, silt, sand, pH.in.CaCl2, pH.in.H2O, OC, CaCO3, N, P, K, CEC a średniej wartości RD
Zbudowano model próbujący określać średnią wartość RD na podstawie wybranych cech numerycznych. Poniżej widać wyniki prostego modelu liniowego:
##
## Call:
## lm(formula = log10(RD_mean) ~ ., data = ddd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.42060 -0.25923 0.02627 0.27423 1.01463
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.735e+00 5.983e-01 -4.571 4.90e-06 ***
## coarse 2.730e-03 2.871e-04 9.509 < 2e-16 ***
## clay 2.677e-03 5.986e-03 0.447 0.654702
## silt -2.180e-03 5.986e-03 -0.364 0.715677
## sand -6.812e-04 5.981e-03 -0.114 0.909324
## pH.in.CaCl2 -1.136e-01 1.520e-02 -7.471 8.40e-14 ***
## pH.in.H2O 1.107e-01 1.590e-02 6.966 3.40e-12 ***
## OC 9.709e-04 2.463e-04 3.941 8.13e-05 ***
## CaCO3 2.018e-04 2.951e-05 6.840 8.19e-12 ***
## N 1.572e-02 4.549e-03 3.456 0.000549 ***
## P -5.025e-04 1.012e-04 -4.964 6.98e-07 ***
## K 3.036e-05 1.548e-05 1.961 0.049921 *
## CEC 3.106e-03 5.045e-04 6.156 7.62e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3756 on 15924 degrees of freedom
## Multiple R-squared: 0.07752, Adjusted R-squared: 0.07683
## F-statistic: 111.5 on 12 and 15924 DF, p-value: < 2.2e-16Uzyskane wyniki są bardzo niskie, z R2 na poziomie około 0.08.
6. Modele relacji zmiennych coarse, clay, silt, sand, pH.in.CaCl2, pH.in.H2O, OC, CaCO3, N, P, K, CEC a grup wartości RD
Druga próba modelowania miała na celu określenie grup wielkości RD w zależności od cech numerycznych. Użyto w niej metodę lasów losowych.
| TrainAccuracy | TrainKappa | method |
|---|---|---|
| 0.3799605 | 0.0412619 | ranger |
Wyniki, podobnie jak powyżej, okazały się niezadowalające. Dokładność predykcji klas wynosiła ok. 0.38, wartości Kappy równa była 0.04.