1. úloha

Načtení dat

Vážený pane zkoušející, takto jsem načetla data (penkavy_vse) a zároveň si již vytvořila soubor vyčištěný od chybějících hodnot (penkavy).

setwd("/Users/monika/Desktop/Data")
penkavy_vse<-read.table("penkavy89.txt", header=T)
penkavy<-na.omit(penkavy_vse)
attach(penkavy_vse)

1) Jejich základní charakteristiky

Tady jsou jejich základní charakteristiky:

summary(penkavy_vse)
##      penkav         popsano      
##  Min.   : 3.00   Min.   : 8.854  
##  1st Qu.:20.75   1st Qu.:13.114  
##  Median :25.00   Median :14.678  
##  Mean   :25.65   Mean   :15.129  
##  3rd Qu.:30.25   3rd Qu.:16.534  
##  Max.   :48.00   Max.   :25.465  
##                  NA's   :2
nrow(penkavy_vse)
## [1] 80
shapiro.test(penkav)
## 
##  Shapiro-Wilk normality test
## 
## data:  penkav
## W = 0.98706, p-value = 0.6027
shapiro.test(popsano)
## 
##  Shapiro-Wilk normality test
## 
## data:  popsano
## W = 0.94198, p-value = 0.00145
shapiro.test(sqrt(popsano))
## 
##  Shapiro-Wilk normality test
## 
## data:  sqrt(popsano)
## W = 0.96991, p-value = 0.06162
shapiro.test(log(popsano))
## 
##  Shapiro-Wilk normality test
## 
## data:  log(popsano)
## W = 0.98633, p-value = 0.5736
par(mfrow=c(2,2))
hist(penkav)
hist(popsano)
hist(sqrt(popsano))
hist(log(popsano))

Datový soubor obsahuje 80 pozorování, z toho dvě chybějící. Obě dvě proměnné jsou kontinuální spojité. Proměnná penkav má rozdělení blízké normálnímu (potvzeno nezamítnutím normality u Shapiro-Willkova testu). Oproti tomu proměnná popsano je pozitivně šikmá - toto dokazuje mimo jinéi fakt, že průměr je větší než medián. Shapiro-Willkův test zamítá normalitu. Při odmocnění dojde k většímu přiblížení normalitě, zlogaritmování přiblíží rozdělení tomu normálnímu ještě více.

2) Vztahy mezi proměnnými

plot(popsano~penkav)

cor(penkavy$popsano,penkavy$penkav)
## [1] 0.6892152
cor(penkavy$popsano,penkavy$penkav)^2
## [1] 0.4750176
mod<-lm(log(popsano)~penkav)
summary(mod)
## 
## Call:
## lm(formula = log(popsano) ~ penkav)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31272 -0.09977  0.02147  0.08374  0.35044 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.288921   0.052189  43.858  < 2e-16 ***
## penkav      0.015747   0.001912   8.234 3.91e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1461 on 76 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.4715, Adjusted R-squared:  0.4646 
## F-statistic: 67.81 on 1 and 76 DF,  p-value: 3.915e-12
par(mfrow=c(2,2))
plot(mod)
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(mod)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod
## BP = 0.21885, df = 1, p-value = 0.6399
detach(penkavy_vse)

Bodový graf ukazuje na pozitivní vztah mezi proměnnými. Korelační koeficient je okolo 0.69, koeficient determinace 0.48 - tj. proměnná penkav vysvětluje asi 50 % variability v datech popsano. Signifikanci tohoto vztahu ukázala lineární regrese, podrobněji viz níže.

3) Typ problému

Jde o regresní (=kauzální) problém - přilétnutí pěnkavy zapříčinilo to, že Darwin začal zapisovat. Teoreticky by to šlo i obráceně (pěnkavy jsou zvědavé a přilétají se podívat, co Darwin píše), ale toto určitě není tento případ, protože pěnkavy jsou autory experimentu a ovlivňují, kolik jich přilétne (jejich počet je tedy nezávislá proměnná).

4) Interpretace exploratorní analýzy

Nyní se vracím k předchozímu lineárnímu modelu. Závislost pospasných stran vyšla signifikantní, s p-hodnotou 3.91*10^-12. Koeficient determinace vyšel 0.4715 - znovu potvrzení, že model vysvětluje asi polovinu varibility v datech. Následně jsem zkontrolovala diagnostické grafy - linearita a absence přemíry vlivu odlehlých pozorování vyšla v pořádku (grafy 1 a 4). Problém nastával u normality reziduálů a obzvláště homoskedasticity. Tyto hodnoty nejlépe vycházely u zlogaritmovaných dat popsano, oproti ostatním vyzkoušeným: bez úpravy a po odmocnění (pro přehlednost jsem zde tyto modely neuvedla). Ani u zlogaritmovaných dat nebyla homoskedasticita ideální (drobné prohnutí uprostřed), ověřila jsem siji tedy i pomocí Breush-Paganova testu, který mi s p-hodnotou 0.6399 homoskedasticitu nezamítl.

5) Diskuse metodiky

Nevyhnutelným problémem přístupu pěnkav je to, že nejsou v izolovaném prostředí - Darwin tak může popisovat např. místní diverzitu hmyzu či rostlin. Darwin asi zkoumá pěnkavy a bude mít tedy tendence psát více o nich než např. o chrobácích, ale pěnkavy nejsou v pozici to vědět. Problém by také mohly činit denní doby: večer bude třeba Darwin unavený a psát méně než ráno na stejný počet pěnkav. (Nutno ale přiznat, že pěnkavy si počínají více vědecky než celé generace lidských vědců a třeba takový Aristoteles by si z nich mohl vzít příklad).

2. úloha

Načtení dat

Takto jsem načetla data a rovnou již převedla jména lokalit na faktory a vytvořila nové proměnné log_oxo (resp -log(oxo)) a sqrt_voda

setwd("/Users/monika/Desktop/Data")
soos<-read.table("soos89.txt", header=T, stringsAsFactors =T)
soos$log_oxo<-(-log(soos$oxonium))
soos$sqrt_voda<-sqrt(soos$voda)
attach(soos)

1) Základní charakteristiky dat

summary(soos)
##       lokalita     oxonium               voda           log_oxo      
##  okoli    :20   Min.   :2.040e-08   Min.   :0.0350   Min.   : 6.124  
##  reservace:40   1st Qu.:1.955e-06   1st Qu.:0.1422   1st Qu.: 9.300  
##                 Median :6.779e-06   Median :0.2400   Median :11.919  
##                 Mean   :1.226e-04   Mean   :0.2671   Mean   :11.494  
##                 3rd Qu.:9.145e-05   3rd Qu.:0.3812   3rd Qu.:13.147  
##                 Max.   :2.189e-03   Max.   :0.6480   Max.   :17.709  
##    sqrt_voda     
##  Min.   :0.1871  
##  1st Qu.:0.3772  
##  Median :0.4897  
##  Mean   :0.4956  
##  3rd Qu.:0.6175  
##  Max.   :0.8050
par(mfrow=c(2,2))
hist(oxonium, breaks = 40)
hist(voda, breaks = 40)
hist(voda)
hist(log_oxo, breaks = 40)

Byly naměřeny tři proměnné. Proměnná lokalita je faktorem o dvou hladinách: rezervace (40 pozorování) a okolí (20). Zbylé proměnné jsou spojité. Proměnná oxonium je extrémně pozitivně šikmá (rozdíl mediánu a průměru dva řády), proměnná voda je také pozitivně šikmá, což je lépe vidět na histogramu s méně sloupečky.

2) Jednotky

par(mfrow=c(2,2))
hist(oxonium, breaks = 40)
hist(log_oxo, breaks = 40)
qqnorm(oxonium)
qqline(oxonium)
qqnorm(log_oxo)
qqline(log_oxo)

shapiro.test(log_oxo)
## 
##  Shapiro-Wilk normality test
## 
## data:  log_oxo
## W = 0.97634, p-value = 0.2936
hist(sqrt(voda))
hist(log(voda))
qqnorm(sqrt(voda))
qqline(sqrt(voda))
qqnorm(log(voda))
qqline(log(voda))

shapiro.test(sqrt(voda))
## 
##  Shapiro-Wilk normality test
## 
## data:  sqrt(voda)
## W = 0.97696, p-value = 0.3139
shapiro.test(log(voda))
## 
##  Shapiro-Wilk normality test
## 
## data:  log(voda)
## W = 0.94954, p-value = 0.01482

Výše vytvořené jednotky jsem zvolila tak, aby maximálně odpovídaly normálnímu rozdělení. U oxonia byla volba log_oxo díky rozšířenosti pH (které je ale dekadický, ne logaritmus naturalis) této míry nasnadě. Rozdělení log_oxo se jak na histogramu, tak na QQ plotu příliš neodlišovalo od toho normálního. Shapiro-Willkův test normalitu nezamítl. U vody nebyla otázka, co s prouměnnou, takto jednoduchá. Proto jsem porovnala transformaci odmocněním a zlogaritmováním. Při zlogaritmování se pozitivně šikmá proměnná voda stala negativně šikmou, zatímco odmocněná lépe odpovídala normálnímu rozdělení (viz histogramy a QQ ploty). Proto jsem dále pracovala s sqrt_voda. Tanto výběr jsem podpořila i Shapiro-Willkovým testem, kdy normalitu zamítalmu log_voda a nezamítal u sqrt_voda.

3) Návrh pokusu

par(mfrow=c(1,1))
plot(log_oxo~sqrt_voda)

cor(log_oxo, sqrt_voda)
## [1] 0.2156651
mod1<-lm(log_oxo~sqrt_voda)
summary(mod1)
## 
## Call:
## lm(formula = log_oxo ~ sqrt_voda)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.2721 -1.7680  0.0868  1.7198  6.4573 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    9.589      1.181   8.118 3.88e-11 ***
## sqrt_voda      3.844      2.286   1.682   0.0979 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.594 on 58 degrees of freedom
## Multiple R-squared:  0.04651,    Adjusted R-squared:  0.03007 
## F-statistic: 2.829 on 1 and 58 DF,  p-value: 0.09794

Pokus není navržen dobře, protože hladina rezervace byla zastoupena 40 pozorováními, zatímco hladina okolí jen 20. Nebudeme tak moci odhadnout efekt lokality se stejnou přesností. Ale alespoň na sobě nejsou prediktory sqrt_voda a log_oxo závislé, takže by neměl být problém s korelací prediktorů (k určení závislosti prediktorů na sobě jsem použila lm, který mlčky předpokládá kauzalitu, což je asi v tomto případě statistické faux pas, ale nevím, jak by to šlo jinak).

4) Odhady statistik základního souboru

prum <- tapply(log_oxo,lokalita,mean)
se<-function(x){sd(x)/sqrt(length(x))}
chyb <- tapply(log_oxo,lokalita,se)
plot(prum~c(1:2), pch=16, xlim = c(0.5,2.5), ylim=c(5,15),col= "black", ylab="Oxonium (mol/l, -log(naturalis))", xlab="Lokalita", xaxt="n")
arrows(x0=c(1:2), x1=c(1:2), y0=prum-1.96*chyb, y1=prum+1.96*chyb,angle=90,length=0.05, code=3)
axis(1,at=1:2, label=levels(lokalita))

prum2 <- tapply(sqrt_voda,lokalita,mean)
chyb2 <- tapply(sqrt_voda,lokalita,se)
plot(prum2~c(1:2), pch=16, xlim = c(0.5,2.5), ylim=c(0,1),col= "black", ylab="Obsah vody (%, sqrt)", xlab="Lokalita", xaxt="n")
arrows(x0=c(1:2), x1=c(1:2), y0=prum2-1.96*chyb2, y1=prum2+1.96*chyb2,angle=90,length=0.05, code=3)
axis(1,at=1:2, label=levels(lokalita))

prum_oxo <- tapply(oxonium,lokalita,mean)
-log10(prum_oxo)
##     okoli reservace 
##  3.451710  5.142551

Znázorněné body jsou průměry naměřených hodnot v daných kategoriích, rozsahy 95% CI. U log_oxo jsou rozdíly mezi lokalitami výrazné - jejich intervaly se nepřekrývají. (Poznámka: důvod proč jsou hodnoty “pH” v lokalitě plné síry tak vysoké - nejde o pH (dekadický log), ale log se základem e; log_oxo průměr okolí = 8.502768 ~ 3.451710 pH, zatímco log_oxo průměr rezervace = 12.990143 ~ 5.142551). Naopak zastoupení vody v půdě se tolik nelišilo, CI se překrývají a můžeme očekkávat, že lm neuvidím mezi nimi rozdíl.

5) Interpretace výsledků

t.test(sqrt_voda~lokalita)
## 
##  Welch Two Sample t-test
## 
## data:  sqrt_voda by lokalita
## t = -1.7299, df = 41.059, p-value = 0.09115
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.14529777  0.01121952
## sample estimates:
##     mean in group okoli mean in group reservace 
##               0.4508891               0.5179282
var.test(sqrt_voda~lokalita)
## 
##  F test to compare two variances
## 
## data:  sqrt_voda by lokalita
## F = 0.84792, num df = 19, denom df = 39, p-value = 0.7163
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4045464 1.9821728
## sample estimates:
## ratio of variances 
##            0.84792
plot(sqrt_voda~lokalita)

t.test(log_oxo~lokalita)
## 
##  Welch Two Sample t-test
## 
## data:  log_oxo by lokalita
## t = -12.193, df = 54.788, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.225000 -3.749748
## sample estimates:
##     mean in group okoli mean in group reservace 
##                8.502768               12.990143
var.test(log_oxo~lokalita)
## 
##  F test to compare two variances
## 
## data:  log_oxo by lokalita
## F = 0.39482, num df = 19, denom df = 39, p-value = 0.03272
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1883696 0.9229622
## sample estimates:
## ratio of variances 
##          0.3948183
plot(log_oxo~lokalita)

Jak již jsem předeslala, analýza dat našla rozdíl mezi množství oxoniových kationtů na stanovištích (p-hodnota: < 2.2e-16), zatímco u množství vody v půdě rozdíl neviděla (p-hodnota: 0.09115). Problémem byl nestejný rozptyl v datech u log_oxo - důsledek právě jiných počtů nasbíraných dat dle lokalit.

6) Korelace a kauzality

Nijak - rezervace byla s největší pravděpodobností vyhlášena právě díky specifickým podmínkám. Dalo by se uvažovat o zmenšené míře disturbance zemědělstvím a stavbou, ale předpokládám, že na to se otázka neptá.