Popis
V jednom mysliveckém revíru se po dobu 7 let lovili zajíci. Známe počty honů a počty
ulovených zajíců pro každý rok. Nově se začala lovit i jelení zvěř. Chceme zjistit, jestli je lov
zajíců a jelení zvěře trvale udržitelný.
Data
Pro lov zajíců máme pro každý rok záznamy o počtu honů a ulovených zajíců:
Rok Počet honů
Počet ulovených
zajíců
1 23 95
2 18 85
3 20 79
4 11 58
5 37 73
6 27 86
7 5 29
Jelikož lov jelení zvěře započal před rokem, známe pouze počet ulovených jedinců v loňském
roce. To bylo 72 jedinců. Poznáme však následující parametry: nosnou kapacitu prostředí, K =
450, konečnou míru populačního růstu, λ = 1.63, a délku života jelenů, L = 17 let.
Úkoly
1/ Zjistěte hodnotu trvale udržitelného lovu pro zajíce přepočtenou na jeden hon.
2/ Odhadněte hodnotu maximálního trvale udržitelného odlovu (MSY) pro jelení zvěř a
porovnejte s hodnotou skutečného lovu. Spočtěte také hodnotu MSH.
14 Trvale udržitelný lov
Řešení
1/ Do vektoru s názvem hunt vložte počet honů a do vektoru s názvem hare počty
ulovených zajíců. Vykreslete závislost hare na hunt do bodového grafu.
> hunt <- c(23,18,20,11,37,27,5)
> hare <- c(95,85,79,58,73,86,29)
> plot(hunt,hare)
Překreslete graf do sešitu.
5 10 15 20 25 30 35
30405060708090
Závislost ulovených zajíců na počtu honů
hunt
hare
Roste počet ulovených zajíců lineárně s počtem honů? …………………………………
Daty proložte Schaeferův model (kvadratický model bez průsečíků) ve tvaru:
2
cxbxy += ,
kde y je počet ulovených zajíců a x je počet honů. Použijte příkaz lm. Hodnoty parametrů
zjistěte příkazem coef.
> m <- lm(hare~hunt+I(hunt^2)-1)
> coef(m)
Výsledný model má tvar: …………………………………………
14. TRVALE UDRŽITELNÝ LOV
Odhadnutý model překreslete do grafu uvedeného výše s použitím kombinace příkazů lines
a predict.
> x <- seq(0,40,1)
> lines(x,predict(m,list(hunt=x)))
Hodnotu trvale udržitelného lovu z odhadnutého modelu zjistíte jako lokální maximum
paraboly. První derivaci výsledného modelu proveďte analyticky a vyjádřete x:
x = ……………
> 6.8857977/(2*0.1327181)
Hodnota trvale udržitelného lovu přepočtena na jeden hon je ………………
Zjistěte, kolik je maximální počet zajíců pro trvale udržitelný lov.
> 6.8857977*25.94144-0.1327181*25.94144^2
Je to ………………… zajíců.
2/ Hodnotu MSY pro jelení zvěř spočtěte podle vzorce Robinsona & Redforda (1991):
⎟
⎠
⎞
⎜
⎝
⎛ −
=
2
KK
aMSY
λ
,
kde hodnota a je závislá na délce života organismů následovně: a = 0.6 pro L < 5, a = 0.4 pro
L = (5,10), a = 0.2 pro L > 10. V našem případě je a = 0.2.
> 0.2*(1.63*450-450)/2
Odhad zaznamenejte do tabulky:
MSY MSH
Hodnotu MSH spočteme podle vzorce
4
rK
MSH = ,
kde r = ln(λ).
> log(1.63)*450/4
14. TRVALE UDRŽITELNÝ LOV
Odhad zapište do tabulky výše.
Je současný lov zajíců a jelení zvěře udržitelný?
…………………………………………………………………………………………………
Poznámka
Kromě uvedeného modelu existují další praxí ověřené postupy pro výpočet trvale
udržitelného lovu, jako je model Foxův, Gullandův nebo model Cadimy atd. (Sparre &
Venema 1998).
14. TRVALE UDRŽITELNÝ LOV
Popis
V porostu mladé slunečnice byla sledována početnost síťových pavouků. Pavouci si staví sítě
na spodní straně listů. Bylo prostudováno 5 rostlin, každá měla 5 listů. Početnost pavouků je
uvedena jako součet jedinců na jeden list.
Data
Byly zjištěny tyto početnosti (N) na list, seřazené od spodního listu k vrchnímu:
Rostlina N
1 0, 0, 1, 5, 7
2 0, 1, 1, 4, 1
3 0, 7, 2, 3, 0
4 0, 1, 6, 1, 1
5 1, 2, 0, 3, 2
Úkoly
1/ Pomocí χ2
-testu a koeficientu disperze zjistěte, jestli je rozmístění pavouků náhodné,
pravidelné nebo agregované vzhledem k listům i k celým rostlinám.
2/ Pokud je rozmístění agregované, odhadněte stupeň agregace.
Řešení
1/ Do vektoru s názvem spider vložte početnost pavouků. Rozložení pavouků na listech
zobrazte histogramem, příkaz hist.
> spider <- c(0,0,1,5,7,0,1,1,4,1,0,7,2,3,0,0,1,6,1,1,1,2,0,3,2)
> hist(spider)
Prostorové rozmístění
populace15
Překreslete histogram do sešitu.
Nyní otestujte hypotézu o pravidelném rozmístění pavouků na listech použitím χ2
-testu
homogenity:
∑=
−
=
I
i
i
x
xx
1
2
2 )(
χ , kde I = 25 a x je aritmetický průměr.
Nulová hypotéza testu předpokládá, že pravděpodobnost nalezení pavouka na nějakém listu je
stejná pro všechny listy v souboru. Tento test je dostupný v příkazu chisq.test.
> chisq.test(spider)
Je rozmístění pavouků na listech náhodné?
…………………………………………………………………………………………………
2/ Koeficient disperze (CD) a stupeň agregace (k) spočítejte podle vzorců:
x
s
CD
2
= ,
xs
x
k
−
= 2
2
,
kde s2
je výběrový rozptyl.
Histogram pavouků
spider
Frequency
0 1 2 3 4 5 6 7
0
5
10
15
15. PROSTOROVÉ ROZMÍSTĚNÍ POPULACE
> var(spider)/mean(spider)
> mean(spider)^2/(var(spider)-mean(spider))
CD = ……………
k = ……………
Popište, jaké mají pavouci rozmístění na listech? ………………………………………………
Nyní otestujte rozložení pavouků vzhledem k celým rostlinám. Do vektoru s názvem plant
vložte sumy jedinců z jedné rostliny kombinací příkazů c a sum. Součty jedinců na
slunečnicích otestujte opět χ2
-testem homogenity a spočtěte koeficient disperze (CD) na
rostlinu:
> plant <- c(sum(spider[1:5]),sum(spider[6:10]),sum(spider[11:15]),
+ sum(spider[16:20]),sum(spider[21:25]))
> chisq.test(plant)
> var(plant)/mean(plant)
Jaká je hodnota CD vzhledem k rostlinám? ……………………………………………
Jak tedy klasifikujete rozmístění pavouků na rostlinách?
……………………………………………………………………………………………
Poznámka
Příkaz pro výpočet CD je dostupný v balíčku lawstat. Kromě koeficientu disperze bylo
navrženo několik dalších a přesnějších metod pro stanovení typu rozmístění populace, jako je
Lloyd mean crowding, Iwao patchiness regression nebo metody založené na měření
vzdálenosti a použité především mezi sedentérními jedinci (viz Southwood & Henderson
2000).
15. PROSTOROVÉ ROZMÍSTĚNÍ POPULACE
Popis
Snůšky pavouka jednoho druhu byly zavlečeny na nové území, kde se mladí jedinci začali
šířit do okolí. Hned po zavlečení bylo sledováno jejich šíření, a to v průběhu 3 měsíců na
podzim a 3 měsíců na jaře (v zimě jsou pavouci neaktivní, proto nebylo šíření studováno),
celkem tedy 6 měsíců před prvním rozmnožováním.
Data
Zjištěná data o ploše, na které byli pavouci nalezeni v průběhu 6 měsíců, jsou v tabulce:
Měsíc Plocha
[km2
]
0 0
1 0.02
2 0.05
3 0.11
4 0.17
5 0.26
6 0.30
Úkoly
1/ Odhadněte koeficient šíření a míru expanze, pokud znáte konečnou míru růstu, λ = 1.4.
2/ Namodelujte šíření pavouků s použitím Skellamova modelu a odhadněte, do jaké
vzdálenosti se rozšíří za 5 let, pokud by na počátku z kokonů vylezlo 50 jedinců.
Řešení
1/ Do vektoru s názvem area vložte plochu a do vektoru month měsíce. Ty však
musíte přepočíst na rok (vydělit 12), protože později budeme pracovat s časovou jednotkou
Šíření populace v prostoru16
rok. Vyneste závislost area na čase t do bodového grafu.
> plot(month,sqrt(area) ,xlab="year")
Překreslete graf do sešitu.
0.0 0.1 0.2 0.3 0.4 0.5
0.00.10.20.30.40.5
Šíření pavouků za 6 měsíců
year
sqrt(area)
Koeficienty šíření pavouků (D) odhadněte pomocí lineárního modelu závislosti odmocniny
plochy na čase:
tarea β=
Lineární model, příkaz lm, aplikujte bez průsečíku (-1), protože na počátku nebyl na novém
území přítomen žádný pavouk. Odhad β zjistíte příkazem coef. Odhadnutý model přidejte
příkazem abline.
> month <- 0:6/12
> area <- c(0,0.02,0.05,0.11,0.17,0.26,0.3)
> m <- lm(sqrt(area)~month-1)
> coef(m)
> abline(m)
Odhad koeficientu šíření D spočítejte podle:
4
ˆβ
=D
> 1.195874/4
16. ŠÍŘENÍ POPULACE V PROSTORU
Odhadnutý koeficient šíření D (bez příspěvku množení) je …… km2
/rok.
Odhad míry expanze, c, která zohledňuje i množení, spočítejte podle
rDc 2= ,
kde )ln(λ=r .
> log(1.4)
> 2*sqrt(0.34*0.3)
Odhadnutá míra expanze je ………… km2
/rok.
2/ Pro simulaci šíření použijte Skellamův model. Ten je rozšířením Gaussovy funkce o
vnitřní míru růstu populace (r) a kruhovitý dvourozměrný prostor:
Dt
rt
e
Dt
N
tN 4
2
4
)0,0(
),(
ρ
π
ρ
−
= ,
kde ρ je poloměr plochy rozšíření. Výsledkem Skellamova modelu je početnost populace
v daném čase a vzdálenosti od počátku šíření.
Vytvořte matici početností N příkazem matrix tak, že každý řádek bude jeden časový krok a
sloupce budou odpovídat lineární vzdálenosti od počátku šíření. Proto nejprve vytvořte vektor
rho s rozsahem 0 až 5 km a vektor months s rozsahem 0 až 5 let použitím příkazu seq.
Simulaci proveďte příkazem for. Použijte výše odhadnuté hodnoty parametrů r a D a
počáteční počet jedinců 50. Graf, který vytvoříte příkazem matplot, bude poněkud
netypicky modelovat početnost pavouků vzhledem ke vzdálenosti od počátku šíření. Čas
nebude na žádné z os, bude prezentován jednotlivými křivkami.
> rho <- seq(0,5,by=0.2)
> months <- seq(0.083,5,by=1)
> N <- matrix(0,ncol=length(rho),nrow=length(months))
> r=0.34; D=0.3; i=1
> for(t in months){N[i,] <- 50*exp(r*t-rho^2/(4*D*t))/(4*pi*D*t);i=i+1}
> matplot(rho,t(N),type="l")
> legend("topright",c("1","2","3","4","5"),lty=1:5,col=1:5)
Překreslete graf do sešitu.
16. ŠÍŘENÍ POPULACE V PROSTORU16. ŠÍŘENÍ POPULACE V PROSTORU
0 1 2 3 4 5
050100150
Simulace šíření pavouků
rho
N
Do jaké vzdálenosti od epicentra se pavouci dostanou za 5 let?
…………………………………………………………………………………………………
Poznámka
Modelování šíření druhů je detailně popsáno v Shigesada et al. (1995) a Shigesada &
Kawasaki (1997), kde lze nalézt model i pro kombinaci šíření na krátkou a delší vzdálenost,
tzv. stratified diffusion model. Lubina & Levin (1988) ukazují, jak pořídit odhad D z hodnot
rozptylu šířící se populace. V prostředí R lze vyzkoušet balíček MigClim pro predikci
rozšíření druhů třeba v závislosti na změně klimatu.
16. ŠÍŘENÍ POPULACE V PROSTORU
Popis
Populace skokanů byla výstavbou dálnice rozdělena na dvě nerovnoměrně veliké části.
Demografickým výzkumem obou populací se zjistilo, že větší populace vyčerpala potravní
zdroje, proto její početnost klesá. Menší populace má k dispozici nadbytek potravy, proto její
početnost stoupá.
Data
Aktuální počet jedinců (N) a konečná míra růstu (λ) pro obě populace jsou uvedeny v tabulce:
Populace 1 Populace 2
N 100 10
λ 0.8 1.2
Úkoly
1/ Použijte hustotně nezávislý model pro diskrétně se množící organismy a nasimulujte
populační dynamiku obou populací po dobu 20 generací. Zjistěte, jestli jsou izolované
populace dlouhodobě životaschopné.
2/ Pokud je alespoň jedna populace v ohrožení, simulací zjistěte, jestli by pomohlo
vytvoření koridoru. Jak velký by měl být ve smyslu propustnosti jedinců?
Řešení
1/ K simulaci použijte diskrétní model exponenciálního růstu pro větší (N1) i menší (N2)
populaci:
))1(( ,2,111,1 ttt dNNdN +−=+ λ
))1(( ,1,221,2 ttt dNNdN +−=+ λ ,
Metapopulační dynamika17
kde d je koeficient propustnosti neboli výměny jedinců mezi populacemi.
Vytvořte matici N12 složenou z 2 vektorů délky 20. Na první pozici jednoho vektoru umístěte
hodnotu 100, u druhého 10, jako počáteční početnosti. Do objektu N12 ukládejte
nasimulované početnosti obou populací. Pro izolované populace bude d = 0, tj. 0%
propustnost. Simulaci proveďte příkazem for a početnosti vyneste do grafu příkazem
matplot.
> N12 <- data.frame(N1=rep(0,20),N2=rep(0,20))
> N12[1,1] <- 100
> N12[1,2] <- 10
> d=0
> for(t in 1:20) {
+ N12[t+1,1] <- ((1-d)*N12[t,1]+d*N12[t,2])*0.8
+ N12[t+1,2] <- ((1-d)*N12[t,2]+d*N12[t,1])*1.2}
> matplot(N12,type="l",xlab="generation",ylab="N",lty=1:2)
> legend("topleft",legend=c("N1","N2"),lty=1:2,col=1:2)
Překreslete graf do sešitu.
5 10 15 20
0100200300
Dynamika izolovaných populací
generation
N
Popište, jaký je osud izolovaných populací?
…………………………………………………………………………………………………
…………………………………………………………………………………………………
2/ Nyní proveďte simulaci dynamiky propojených populací. Použijte stejný systém
17. METAPOPULAČNÍ DYNAMIKA
rovnic a stejnou syntaxi, ale s malou hodnotou koeficientu propustnosti, například d = 0.2, tj.
s 20% propustností.
> d=0.2
> for(t in 1:20) {
+ N12[t+1,1] <- ((1-d)*N12[t,1]+d*N12[t,2])*0.8
+ N12[t+1,2] <- ((1-d)*N12[t,2]+d*N12[t,1])*1.2}
> matplot(N12,type="l", xlab="generation",ylab="N",lty=1:2)
> legend("topleft",c("N1","N2"),lty=1:2,col=1:2)
Překreslete graf do sešitu.
5 10 15 20
50100150
Dynamika málo propojených populací
generation
N
Pomohlo to k udržení obou populací? …………………………………………………………
Nyní nasimulujte vysokou propojenost populací, například d = 0.8, tedy 80% propustnost.
> d=0.8
> for(t in 1:20) {
+ N12[t+1,1] <- ((1-d)*N12[t,1]+d*N12[t,2])*0.8
+ N12[t+1,2] <- ((1-d)*N12[t,2]+d*N12[t,1])*1.2}
> matplot(N12,type="l",xlab="generation",ylab="N",lty=1:2)
> legend("topright",c("N1","N2"),lty=1:2,col=1:2)
Překreslete graf do sešitu.
17. METAPOPULAČNÍ DYNAMIKA
5 10 15 20
20406080100
Dynamika hodně propojených populací
generation
N
Jaká by měla být propustnost koridoru, aby obě populace byly životaschopné?
…………………………………………………………………………………………………
17. METAPOPULAČNÍ DYNAMIKA
Popis
V epigeonu pole se vyskytují dva druhy pavouků, jeden z rodu Pardosa a druhý z rodu
Pachygnatha. Oba druhy jsou pravými predátory a konzumují pouze bezobratlé živočichy.
Detailním pozorováním lovu kořisti v terénu a složením dostupné kořisti bylo zjištěno,
kterými druhy bezobratlých živočichů se oba pavouci živí.
Data
Relativní četnosti kořisti z 5 skupin bezobratlých pro oba druhy pavouků, relativní dostupnost
kořisti a počet pozorování (N) jsou uvedeny v tabulce:
Collembola Hemiptera Ensifera Diptera Isopoda N
Pardosa 0.61 0.15 0.12 0.07 0.05 100
Pachygnatha 0.93 0.05 0.01 0.00 0.01 100
dostupnost 0.60 0.10 0.10 0.10 0.10 50
Úkoly
1/ Odhadněte šířku trofické niky pro oba druhy pavouků a otestujte, jestli jsou si podobné.
2/ Zjistěte, jak moc se jejich niky překrývají.
Řešení
1/ Do vektoru s názvem Par a Pach umístěte relativní četnosti zastoupení jednotlivých
řádů kořistí. Do vektoru s názvem av pak umístěte relativní četnost dostupné kořisti. Pak
vektory Par a Pach spojte příkazem rbind do objektu both a relativní četnosti obou
druhů pavouků vyneste do sloupcového grafu příkazem barplot.
Šířka niky18
> Par <- c(0.61,0.15,0.12,0.07,0.05)
> Pach <- c(0.93,0.05,0.01,0,0.01)
> av <- c(0.6,0.1,0.1,0.1,0.1)
> both <- rbind(Par,Pach)
> barplot(both,beside=T,ylim=c(0,1),legend.text=c("Pardosa","Pachygnatha"),
+ names.arg=c("Collembola","Hemiptera","Ensifera","Diptera","Isopoda"))
Překreslete graf do sešitu.
Collembola Hemiptera Ensifera Diptera Isopoda
0.00.20.40.60.81.0
Relativní četnosti ulovené kořisti
Pro odhad šířky trofické niky použijte Smithův index (Smith 1982):
∑=
=
5
1k
kk qpFT ,
kde pk je relativní četnost uloveného taxonu k a qk je relativní dostupnost tohoto taxonu.
Smithův index se interpretuje vzhledem k teoretické minimální (0) a maximální (1) šířce niky.
> F1 <- sum(sqrt(Par*av)); F1
> F2 <- sum(sqrt(Pach*av)); F2
Odhady indexů šířky niky jsou:
Pardosa Pachygnatha
18. ŠÍŘKA NIKY
Který druh má širší niku a proč?
…………………………………………………………………………………………………
Je rozdíl mezi druhy statisticky významný? K zodpovědění otázky je potřeba nejprve spočíst
rozptyly obou odhadů podle vzorce:
N
FT
FTVar
4
1
)(
2
−
= ,
kde N je celkový počet zaznamenaných jedinců kořisti.
> V1 <- (1-F1^2)/400
> V2 <- (1-F2^2)/400
Nyní proveďte oboustranný test nulové hypotézy o rovnosti šířky nik. Testovou statistiku
Z spočtěte podle vzorce:
)()( 21
21
FTVarFTVar
FTFT
Z
+
−
=
Tato statistika má standardní normální rozdělení Z ~ N(0,1). Pravděpodobnost chyby prvního
druhu spočtěte z distribuční funkce pro normální rozdělení (pnorm):
> Z <- (F1-F2)/sqrt(V1+V2)
> 2*(1-pnorm(abs(Z)))
P = ………
Je tento rozdíl mezi druhy i statisticky významný? …………………………………
2/ Překryv nik pro oba druhy spočtěte podle vzorce Hulberta (1978):
∑=
k
kk
q
pp
L 21
> sum(Par*Pach/av)
Odhadnutý překryv nik obou druhů je: .……………………
Poznámka
Existuje celá řada indexů pro měření šířky niky a překryvu nik. Liší se v několika aspektech:
především v typu dat, ze kterých se počítají (kategorické či kontinuální dimenze niky), ale liší
se také v tom, jestli zohledňují dostupnou kořist (MacArthur & Levins 1967, Colwell &
Futuyma 1971, Hurlbert 1978, Feinsinger et al. 1981).
18. ŠÍŘKA NIKY
Popis
Dva druhy sekáčů, Opilio a Phalangium, se vyskytovaly společně na jednom stanovišti. Na
vybrané ploše byly sledovány jejich početnosti jednou ročně po dobu 12 let. Chceme zjistit,
zda mezi nimi dochází ke kompetici.
Data
Byly zaznamenány tyto početnosti (N) sekáčů v průběhu 12 let:
Rok N Opilio N Phalangium
1 50 120
2 66 81
3 59 76
4 74 71
5 83 63
6 84 65
7 96 55
8 98 37
9 101 28
10 126 22
11 135 16
12 136 10
Úkoly
1/ Odhadněte parametry míry populačního růstu, nosné kapacity prostředí a koeficientu
kompetice pro oba druhy.
2/ Nasimulujte populační dynamiku obou druhů po dobu 20 generací s použitím
odhadnutých parametrů a počátečními početnostmi 20 jedinců pro každý druh.
Odhady parametrů modelu
kompetice19
Řešení
1/ Do vektoru s názvem opi umístěte početnosti druhu Opilio a do vektoru s názvem
phal početnosti druhu Phalangium. Oba vektory spojte do objektu dat. Do grafu vyneste
početnosti pavouků v závislosti na čase příkazem ts.plot.
> opi <- c(50,66,59,74,83,84,96,98,101,126,135,136)
> phal <- c(120,81,76,71,63,65,55,37,28,22,16,10)
> dat <- data.frame(opi,phal)
> ts.plot(dat,ylab="N")
Překreslete graf do sešitu.
Početnosti sekáčů v průběhu 12 let
Time
N
2 4 6 8 10 12
20406080100120140
Jsou dynamiky obou druhů podobné? Popište, jaké jsou.
…………………………………………………………………………………………………
…………………………………………………………………………………………………
Předpokládáme, že sekáči se množí diskrétně. Proto použijte např. Rickerův model pro
diskrétní růst závislý na hustotě pro dva konkurenční druhy:
⎟⎟
⎠
⎞
⎜⎜
⎝
⎛ −−
+ = 1
,212,11
1
,11,1
K
NNK
r
tt
tt
eNN
α
a
⎟⎟
⎠
⎞
⎜⎜
⎝
⎛ −−
+ = 2
,121,22
2
,21,2
K
NNK
r
tt
tt
eNN
α
,
kde ri je vnitřní míra populačního růstu, Ki je nosná kapacita a ijα je koeficient kompetice.
19. ODHADY PARAMETRŮ MODELU KOMPETICE
Odhady parametrů ri, Ki, a ijα získáme tzv. dynamickou regresí. Ta zahrnuje nejprve
linearizaci těchto rovnic (zvlášť pro každý druh):
tj
i
iji
ti
i
i
i
ti
ti
N
K
r
N
K
r
r
N
N
,,
,
1,
ln
α
−−=⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛ +
⇒ tjijtiii
ti
ti
NcNba
N
N
,,
,
1,
ln ++=⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛ +
pro t = 1, …, 11.
Je to rovnice roviny v prostoru s osami Ni,t, Nj,t a ⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛ +
ti
ti
N
N
,
1,
ln . Pro grafické znázornění bychom
potřebovali trojrozměrný graf. Pomocí regresního modelu se hledá nejlepší proložení roviny
naměřenými body zvlášť pro druh Opilio a Phalangium. Z regresních odhadů parametrů ai, bi
a cij odvoďte odhady parametrů ri, Ki a αij podle vzorců:
ii ar ˆ=
i
i
i
i
i
b
a
b
r
K
ˆ
ˆ
ˆ
−=−=
i
ij
i
iji
ij
b
c
r
cK
ˆ
ˆˆ
=−=α pro i = 1, 2 a j = 2, 1.
Hodnoty
ti
ti
N
N
,
1, +
si musíte nejprve dopočíst a uložit do vektorů opi1 a pahl1. Použijte
lineární model, příkaz lm, a odhady parametrů zjistěte příkazem coef.
> opi1 <- opi[-1]/opi[-12]
> phal1 <- phal[-1]/phal[-12]
> m1 <- lm(log(opi1)~opi[-12]+phal[-12])
> coef(m1)
> -0.2827058214/-0.0019409899
> -145.6503*-0.0003509306/0.2827058214
> m2 <- lm(log(phal1)~phal[-12]+opi[-12])
> coef(m2)
> -1.253529913/-0.007836019
> -159.9702*-0.011631378/1.253529913
Odhady parametrů zapište do tabulky y:
Opilio Phalangium
r
K
α
Porovnejte odhady všech parametrů mezi druhy. Mají podobnou vnitřní míru růstu a nosnou
kapacitu?
…………………………………………………………………………………………………
Jsou oba druhy podobně silní konkurenti?
…………………………………………………………………………………………………
19. ODHADY PARAMETRŮ MODELU KOMPETICE
2/ K simulaci dynamiky použijte systém rovnic uvedených výše s počátečními
hodnotami početností 20 jedinců pro oba druhy. Výsledky simulací ukládejte do matice
s názvem N12. Simulaci proveďte příkazem for a její výsledek vyneste příkazem matplot.
> N12 <- matrix(ncol=2,nrow=21)
> N12[1,] <- c(20,20)
> for(t in 1:20) {
+ N12[t+1,1] <- N12[t,1]*exp(0.28*(145.7-N12[t,1]-0.18*N12[t,2])/145.7)
+ N12[t+1,2] <- N12[t,2]*exp(1.25*(160-N12[t,2]-1.48*N12[t,1])/160)}
> matplot(N12,type="l",xlab="generation",lty=1:2)
> legend("right",c("Opilio","Phalangium"),col=1:2,lty=1:2)
Překreslete graf do sešitu.
5 10 15 20
020406080100120140
Simulovaná populační dynamika kompetitorů
generation
N12
Jaké jsou osudy obou druhů? ………………………………………………………………
……………………………………………………………………………………………….
Poznámka
Dalšími metodami používanými k odhadu parametrů kompetice je tzv. statická regrese
(Schoener 1974, Pfister 1995), která používá data o abundanci konkurentů z několika lokalit,
nebo manipulativní experiment, ve kterém se sleduje početnost konkurentů a jejich míra růstu
po narušení jednoho z druhů (Bender et al. 1984). V dynamické regresi, kterou jsme si
ukázali, by měl být odhad parametrů (především jejich chyb) proveden regresní metodou,
která dokáže modelovat i časovou závislost v datech (třeba GLS nebo GEE, viz Pekár &
19. ODHADY PARAMETRŮ MODELU KOMPETICE
Brabec 2012). Odhad parametrů z dynamické regrese je lepší než ze statické, ve které je nutné
abundance standardizovat a zahrnout do modelu i faktory prostředí (Fox & Luo 1996).
19. ODHADY PARAMETRŮ MODELU KOMPETICE
Popis
Na jedné lokalitě se začal šířit invazní druh mravence a zdá se, že vytlačí jeden z původních
druhů mravenců. Předchozí studie ukázala, že oba druhy mají podobnou niku. Zjistěte, jestli
je původní druh v ohrožení.
Data
V předchozí studii se zjistily tyto hodnoty parametrů původního (1) a invazního (2) druhu:
Původní Invazní
r 0.8 0.7
K 200 300
α 1.6 2.1
Úkoly
1/ Odhadněte výsledek interakce analyticky.
2/ Nasimulujte vývoj početnosti mravenců obou druhů, pokud má populace původního
druhu počáteční početnost 20 jedinců a populace invazního druhu početnost 10 jedinců.
Simulaci proveďte pro dobu 30 let.
3/ Jaký bude výsledek simulace, obrátí-li se počáteční početnosti ve prospěch invazního
druhu (tj. N1,t0 = 10, N2,t0 = 20), např. na lokalitách, kde je původní druh málo hojný?
Řešení
1/ Pro analytické řešení využijte metodu nulových izoklin. Hodnoty K1,
12
1
α
K
, K2 a
21
2
α
K
určují 4 body, kterými jsou izokliny definovány. Spočtěte jejich hodnoty.
Mezidruhová kompetice20
> 200/1.6
> 300/2.1
Vyneste izokliny do grafu:
0 50 100 150 200 250
050100150200250300350
původní mravenec
invaznímravenec
Nulové izokliny konkurenčních druhů
Jaký bude výsledek interakce?
…………………………………………………………………………………………………
2/ K simulaci použijte systém diferenciálních rovnic kompetičního modelu Lotky-
Volterra:
⎟⎟
⎠
⎞
⎜⎜
⎝
⎛ +
−=
1
2121
11
1
1
d
d
K
NN
rN
t
N α
⎟⎟
⎠
⎞
⎜⎜
⎝
⎛ +
−=
2
1212
22
2
1
d
d
K
NN
rN
t
N α
Diferenciální rovnice vyčísluje změnu početnosti Ni za jednotku času. Výsledkem simulace
má být ovšem očekávaná početnost v daném čase. Bylo by potřeba integrovat uvedené
diferenciální rovnice, abyste získali funkci popisující přímo početnost (analytické řešení). To
je často těžko řešitelný úkol. Proto zvolte cestu numerického řešení obyčejných
diferenciálních rovnic (ordinary differential equations, ve zkratce ode). Numerický postup
aproximuje výslednou početnost bez nutnosti integrování rovnic. Nejjednodušším typem
numerického řešení je tzv. Eulerova metoda.
20. MEZIDRUHOVÁ KOMPETICE
Výpočet postupuje po časových krocích typicky menších, než je základní časová jednotka
(např. Δ = 0.1 časové jednotky), aby nepřesnost způsobená aproximací nebyla zbytečně velká.
V každém kroku spočítejte změny početností
t
N
d
d 1
a
t
N
d
d 2
podle posledních známých
početností a k těmto je přičtěte. Tím dostanete odhady početností v čase t+Δ. Délku kroku
musíte zohlednit i v připočtených změnách (za desetinu časové jednotky přibude/ubude jen
desetina spočítané změny početnosti):
t
N
NNN i
tititi
d
d
1,1,2, ⋅Δ+== Δ+
Do objektu delta vložte délku kroku. Ve vektoru time specifikujeme čas pro simulaci.
Definujeme hodnoty všech parametrů a vytvoříme matici N se dvěma sloupci. Na jejich
prvních pozicích budou počáteční početnosti. Simulaci proveďte příkazem for a výsledek
vyneste příkazem matplot.
> delta=0.1
> time <- seq(0,30,by=delta)
> r1=0.8; r2=0.7; a12=1.6; a21=2.1; K1=200; K2=300
> N=matrix(ncol=2,nrow=length(time))
> N[1,]=c(20,10)
> for(i in 1:(length(time)-1)){
+ N[i+1,1] <- N[i,1]+delta*r1*N[i,1]*(1-(N[i,1]+a12*N[i,2])/K1)
+ N[i+1,2]<-N[i,2]+delta*r2*N[i,2]*(1-(N[i,2]+a21*N[i,1])/K2)}
> matplot(time,N,type="l",xlab="years",lty=1:2,col=1)
> legend("topleft",c("N1","N2"),lty=1:2)
Překreslete graf do sešitu.
0 5 10 15 20 25 30
050100150
Simulace početnosti mravenců
years
N
20. MEZIDRUHOVÁ KOMPETICE
Vytlačí invazní druh původního mravence? Pokud ano, jak dlouho to bude trvat?
…………………………………………………………………………………………………
3/ Jak dopadne simulace, bude-li mít na počátku převahu invazní druh mravence?
Použijte stejný systém rovnic, ale řešení tentokrát nalezněte pomocí složitějších a přesnějších
metod. Přístupné jsou v balíku deSolve. Příkaz ode je nejobecnějším řešitelem obyčejných
diferenciálních rovnic. Výpočet probíhá ve dvou krocích. V prvním kroku specifikujte model,
nazvěte jej competition, tj. vytvořte funkci, do které zapište potřebné diferenciální
rovnice. Do kulatých závorek v definici funkce zapište trojici objektů, se kterými bude funkce
počítat: ti bude čas; ini bude vektor s počátečními hodnotami stavových proměnných a
par bude vektor se jmény a hodnotami parametrů. Příkaz with zpřístupní „mateřské“ funkci
všechny parametry a proměnné pod jejich jmény ve formě seznamu (as.list). Zápis
funkce zakončete příkazem return, který předá výsledky rovnic k dalšímu zpracování. Zde
je důležité pořadí: výsledky dN1, dN2 musí odpovídat pořadí N1, N2 ve vektoru počátečních
hodnot.
> competition <- function(ti,ini,par){
+ with(as.list(c(ini,par)),{
+ dN1 <- r1*N1*(1-(N1+a12*N2)/K1)
+ dN2 <- r2*N2*(1-(N2+a21*N1)/K2)
+ return(list(c(dN1,dN2)))})}
Zaveďte parametry a jejich hodnoty (do vektoru param) a zadejte jejich počáteční početnosti
do vektoru initial. Načtěte balík deSolve. Definujte časy, ve kterých má být proveden
výpočet (vektor time). Samotný výpočet obstará příkaz ode, který pojmenuje vstupní
objekty jako (y, times, func, parms). Výslednou matici uložte do objektu sim.
Výsledek simulace vykreslete příkazem matplot.
> param <- c(r1=0.8,r2=0.7,a12=1.6,a21=2.1,K1=200,K2=300)
> initial <- c(N1=10,N2=20)
> library(deSolve)
> time <- seq(0,30,by=0.1)
> sim <- ode(y=initial,times=time,func=competition,parms=param)
> matplot(time,sim[,-1],type="l", xlab="years", ylab="N",lty=1:2)
> legend("topleft",c("N1","N2"),lty=1:2,col=1:2)
Překreslete graf do sešitu.
20. MEZIDRUHOVÁ KOMPETICE
0 5 10 15 20 25 30
050100150200250300
Simulace početnosti mravenců
years
N
Vytlačí nyní invazní druh původního mravence? Pokud ano, jak dlouho to bude trvat?
…………………………………………………………………………………………………
Poznámka
Jak přesné je použití jednoduché Eulerovy numerické metody proti jiným složitějším
metodám z příkazu ode (defaultně je nastavena metoda lsoda), lze zjistit třeba grafickým
porovnáním výsledků simulace. Pokud to uděláte, uvidíte, že rozdíl je v podstatě
zanedbatelný. O použití příkazu ode v simulačním modelování se lze dozvědět více např. v
Petzold (2003).
20. MEZIDRUHOVÁ KOMPETICE
Popis
Funkční odpověď střevlíků konzumujících semena byla zkoumána v miskách o ploše 10 cm2
.
Pro každou hustotu semen bylo uskutečněno několik opakování. Hustota semen byla
konstantní po celou dobu experimentu, tj. 6 hodin, protože zkonzumovaná semena byla
neustále nahrazována novými. Na konci experimentu byl stanoven počet sežraných a
nesežraných semen.
Data
Průměrné počty zkonzumovaných semen (Ha) pro jednotlivé hustoty semen (H) jsou uvedeny
v následující tabulce:
H Ha
1 2.5
5 6.1
10 7.9
20 10.5
40 12.3
50 11.8
Úkoly
1/ Jaký typ funkční odpovědi mají střevlíci?
2/ Odhadněte parametry vyhledávací účinnosti, času zpracování a asymptoty nasycení.
Řešení
1/ Do vektoru s názvem H vložte hustoty semen a do vektoru s názvem Ha průměrné
počty sežraných semen. Do bodového grafu vyneste počty zkonzumovaných semen
v závislosti na jejich hustotě.
Funkční odpověď21
> H <- c(1,5,10,20,40,50)
> Ha <- c(2.5,6.1,7.9,10.5,12.3,11.8)
> plot(H,Ha)
Překreslete graf do sešitu.
0 10 20 30 40 50
4681012
Funkční odpověď střevlíků
H
Ha
Který typ funkční odpovědi střevlíci vykazují? ………………………………………
2/ Hodnoty parametrů odhadněte pomocí příslušné Hollingovy rovnice:
h
a
aHT
aHT
H
+
=
1
Tato rovnice je inherentně nelineární. Abychom k odhadu parametrů mohli použít lineární
regresi, musí se linearizovat. Lineární tvar rovnice získáte převrácením:
T
T
aHTH
h
a
+=
11
a úpravou:
αβ +=
HHa
11
Zjistěte hodnoty parametrů a, Th i asymptotu nasycení (max) podle vzorců:
21. FUNKČNÍ ODPOVĚĎ
TTh αˆ=
βˆ
1
T
a =
αˆ
1
==
hT
T
max
Připravte převrácené hodnoty H a Ha a vyneste je do bodového grafu. Pak použijte lineární
regresní model, příkaz lm. Výslednou přímku vynesete příkazem abline.
> y <- 1/Ha
> x <- 1/H
> plot(x,y,xlab="1/H",ylab="1/Ha")
> m <- lm(y~x)
> abline(m)
Překreslete graf do sešitu.
0.0 0.2 0.4 0.6 0.8 1.0
0.100.150.200.250.300.350.40
Závislost převrácených hodnot
1/H
1/Ha
Odhady parametrů modelu zjistíte příkazem coef. Vypočítejte z nich a, Th a max podle
vzorců uvedených výše.
> coef(m)
> 0.08445655*6
> 1/(0.31904090*6)
> 1/0.08445655
21. FUNKČNÍ ODPOVĚĎ
Odhady parametrů zaznamenejte do tabulky:
a cm2
/h
Th h
max
Poznámka
Přesnější odhady parametrů pro funkční odpověď Typu II lze získat s pomocí GLM (viz
Pekár & Brabec 2009). Dále pro popis funkčních odpovědí bylo navrženo několik
nelineárních modelů (viz Juliano 2001). Jedním z obecných modelů je )1(
c
bx
eay −
−= . Zde
parametr a definuje asymptotu (nasycení), kladné b určuje rychlost růstu a c prohnutí křivky.
Touto funkcí můžeme modelovat funkční odpověď Typu II i III. Pokud c ≤ 1, jde o typ II,
pokud c > 1, jde o typ III.
21. FUNKČNÍ ODPOVĚĎ
Popis
Saranče zkonzumují za svůj život určité množství potravy. V laboratorním experimentu byly
chovány saranče od vylíhnutí až po dospění a následné vykladení při různém konstantním
množství biomasy. Chceme zjistit, jak množství potravy ovlivní jejich fitness.
Data
Ze zjištěných dat o přežívání a plodnosti byla pro každé množství biomasy sestavena Leslieho
matice a z ní spočtena míra populačního růstu (r). Hodnoty vnitřní míry růstu v závislosti na
množství potravy (V) jsou uvedeny v tabulce:
V(g) r
50 -1.0
100 -0.6
200 -0.1
500 0.3
1000 0.5
2000 0.7
4000 1.0
Úkoly
1/ Zobrazte závislost vnitřní míry růstu (r) na množství zkonzumované potravy (V).
Proložte daty model přeměny potravy na potomstvo a odhadněte jeho parametry.
Stanovte minimální množství potravy na jedince potřebné pro vitální populaci.
Řešení
1/ Do vektoru s názvem r vložte hodnoty vnitřní míry růstu a do vektoru s názvem v
množství zkonzumované potravy. Do bodového grafu vyneste závislost míry růstu na
množství potravy:
Numerická odpověď22
> v <- c(50,100,200,500,1000,2000,4000)
> r <- c(-1,-0.6,-0.1,0.3,0.5,0.7,1)
> plot(v,r)
Překreslete graf do sešitu.
0 1000 2000 3000 4000
-1.0-0.50.00.51.0
Závislost míry růstu na množství potravy
v
r
Pro kolik hodnot je míra růstu kladná a pro kolik záporná?
kladné záporné
Použijte Ivlevův model v modifikaci pro numerickou odpověď, který má tři parametry:
dear fV
−−= −
)1(
Parametr a popisuje vyhledávací účinnost konzumenta (search rate), f udává účinnost
přeměny potravy na potomstvo (konverzní účinnost) a d je maximální mortalita konzumenta.
Jelikož je tvar modelu nelineární, v parametrech k odhadu použijte nelineární regresi. Před
tím musíte odhadnout startovací hodnoty všech parametrů. Odvoďte je analyticky pro
podmínky V = 0, V = ∞ a pro r = 0.
22. NUMERICKÁ ODPOVĚĎ
Startovací hodnoty budou: a = ………, f = ………, d = ………
Pro fit nelineární regrese použijte příkaz nls. Hodnoty parametrů zjistíte příkazem coef.
> m <- nls(r~a*(1-exp(-f*v))-d,start=list(a=2,f=0.001,d=1))
> coef(m)
Odhadnuté hodnoty parametrů zapište do tabulky:
a
f
d
Modelovou křivku zakreslete do grafu výše použitím příkazů lines a predict.
> x <- seq(0,4000,1)
> lines(x,predict(m,list(v=x)))
Jaké je minimální množství potravy potřebné k nezápornému populačnímu růstu? Tuto
hodnotu udává kořen Ivlevovy rovnice, tedy:
dea fV
−−= −
)1(0
Ten zjistíte příkazem uniroot z balíčku rootSolve.
> library(rootSolve)
> uniroot(function(x) 1.94*(1-exp(-0.003*x))-1.17,lower=0,upper=1000)
Minimální množství potravy je ………………
Poznámka
Další modely používané pro modelování numerické odpovědi, jmenovitě lineární a
asymptotický, odvozený z Hollingovy rovnice funkční odpovědi Typu II, lze nalézt
v Beddington et al. (1976).
22. NUMERICKÁ ODPOVĚĎ