Popis
Mandelinky kladou vajíčka v různě velikých snůškách. V laboratoři bylo sledováno kladení
jednoho druhu. Deset samic bylo umístěno jednotlivě do nádob s dostatečným prostorem a
množstvím kořisti. Po dobu dvou týdnů byly každý den zaznamenávány počty nakladených
vajíček u každé samice.
Data
V průběhu 15 dní byly zaznamenány tyto průměrné počty (N) nakladených vajíček:
Den N
1 0
2 0
3 3.2
4 3.3
5 3.1
6 2.1
7 1.5
8 1
9 0.8
10 0.6
11 0.3
12 0.1
13 0
14 0
15 0
Úkoly
1/ Vytvořte graf kladení v čase. Popište průběh kladení modelem a odhadněte parametry
modelu.
2/ Zjistěte, který den nastala maximální plodnost.
7Model kladení
Řešení
1/ Průměrný počet vajíček na den umístěte do vektoru egg, dny do vektoru day. Do
bodového grafu vyneste závislost egg na day.
Překreslete graf do sešitu.
0 5 10 15
0123456
Změna plodnosti v čase
day
egg
K namodelování kladení se používá taková funkce z exponenciální třídy, která nejdříve
prudce roste a pak pomalu klesá a tím popisuje přirozenou reprodukční senescenci. V
základním tvaru je to funkce: x
xey −
= . Funkce je prosta parametrů a tudíž ji nelze přizpůsobit
datům. Toho lze docílit použitím obecného tvaru
decxay cxb
+−= −− )(
)( ,
se 4 neznámými parametry: a, b, c a d. Kombinace a a b udává výšku vrcholu, c definuje
průsečík s vodorovnou asymptotou a d průsečík s osou y. Jaké ale budou startovací hodnoty
všech parametrů? Navrhněte a zdůvodněte analyticky pro podmínky x = 0 a x = ∞.
R
7. MODEL KLADENÍ
Startovací hodnoty budou: a = ………, b = ………, c = ………, d = ………
Před prokládáním křivky je třeba zajistit, aby se parametr posunutí c vyskytoval s každým
výskytem x. Zahrnete-li do modelu první den, posunete průsečík funkce s osou x k nule více,
než odpovídá našim datům. Proto data pro odhad modelu upravte: první záznam o snůškách i
první den vypusťte. K proložení použijte funkci pro nelineární regresní model, příkaz nls, ve
které startovací hodnoty definujete argumentem start. Výsledky odhadů zjistíte příkazem
coef. Křivku modelu přidejte do grafu příkazem curve.
Model s odhadnutými parametry má tvar ..................................................................................
Překreslete model do grafu výše.
2/ Odhadovaný počátek kladení zjistíte nalezením kořene rovnice
decxa cxb
+−= −− )(
)(0
Použijte příkaz uniroot z balíčku rootSolve.
Který den začaly samice klást? …………………………………
R
R
7. MODEL KLADENÍ
Hodnota maximální plodnosti odpovídá lokálnímu maximu funkce. K tomu je potřeba
spočítat nejprve první derivaci funkce (pomocí příkazu D) podle proměnné day2 a potom její
nulové body.
Derivace funkce je …………………………………………………
Odhadněte lokální maximum pomocí příkazu uniroot:
Maximální plodnost nastala …… den.
Poznámka
Použitý odhad modelu ignoruje pravděpodobnou časovou závislost v datech (v důsledku
měření stejných samic), proto pro přesnější odhad je nutné použít metodu jako je Non-linear
Generalised Least Squares, která modeluje i autokorelaci (Pinheiro & Bates 2000).
Další typy exponenciálních funkcí používaných pro modelování kladení lze nalézt např.
v Bieri et al. (1983). Zajímavou triangulární funkci navrhl Kindlmann et al. (2001), ve které je
plodnost vztažena k několika fenotypovým charakteristikám, jmenovitě velikosti těla a
velikosti gonád.
Tahák
day<-1:15
egg<-c(0,0,3.2,3.3,3.1,2.1,1.5,1,0.8,0.6,0.3,0.1,0,0,0)
plot(day,egg, xlim=c(0,15),ylim=c(0,4))
day2<-2:15
egg2<-egg[-1]
m<-nls(egg2~a*(day2-c)*exp(-b*(day2-c))+d,start=list(a=5,b=1,c=2,d=0))
coef(m)
curve(5.49*(x-2)*exp(-0.57*(x-2))-0.03,xlim=c(1,15),add=T)
library(rootSolve)
start<-uniroot(function(x) 5.5*(x-2)*exp(-0.6*(x-2))-0.03,
lower=0,upper=10);start
mo<-expression(a*(day2-c)*exp(-b*(day2-c))+d)
D(mo,"day2")
max<-uniroot(function(x) 5.5*exp(-0.6*(x-2))-5.5*(x-2)
*(exp(-0.6*(x-2))*0.6),lower=0,upper=10);max
R
R
7. MODEL KLADENÍ
Popis
Na jaře byla na jedné lokalitě provedena demografická studie tchoře. U každého odchyceného
jedince byl stanoven věk. U samic byl zaznamenán i počet nově narozených mláďat
(postreprodukční census). Tchoř se rozmnožuje v pulzech jednou ročně.
Data
Ze zjištěných počtů jedinců byly stanoveny hodnoty standardizovaného přežívání (lx) a
plodnosti (mx) jako průměrný počet mláďat na jedince u každé věkové kategorie (x):
x lx mx
0 1.00 0.0
1 0.90 3.0
2 0.70 6.0
3 0.40 8.0
4 0.07 4.0
Úkoly
1/ Vykreslete do grafu standardizované přežívání v závislosti na věku.
2/ Odhadněte čistou reprodukční rychlost a generační dobu.
3/ Sestrojte Leslieho přechodovou matici.
Řešení
1/ Do vektoru x umístěte věkové kategorie a do vektoru lx hodnoty standardizovaného
přežívání. Vytvořte graf závislosti lx na x s logaritmickou škálou na ose y.
8
Demografie věkově
strukturované populace
Překreslete graf do sešitu.
0 1 2 3 4
0.10.20.51.0
Křivka přežívání
x
lx
Kterému typu křivky přežívání (I, II, III) zjištěná data odpovídají? …………………………
2/ Čistou reprodukční rychlost, R0, tj. kolikrát se populace za generaci namnoží, spočtete
podle vztahu
∑=
=
n
x
xxmlR
0
0
a generační dobu, T, podle vztahu
0
0
R
mxl
T
n
x
xx∑=
=
Do vektoru s názvem mx umístěte hodnoty plodnosti.
8. DEMOGRAFIE VĚKOVĚ STRUKTUROVANÉ POPULACE
R
R
Čistá reprodukční rychlost, R0 = ………
Generační doba, T = ………
3/ K sestavení Leslieho matice spočtěte míry přežívání (px) a fertilitu (Fx) z hodnot lx a mx
podle vztahů pro postreprodukční census:
x
x
x
l
l
p 1+
= a 1+= xxx mpF
Zakreslete odhady px a Fx do schématu životního cyklu. Připomeňme si, že F4 = 0.
Sestavte z odhadů px a Fx Leslieho přechodovou matici:
















=
...............
...............
...............
...............
...............
A
Poznámka
Demografické metody jsou podrobně popsány v Caswell (2001). V prostředí R je pár balíčků
pro detailní demografickou analýzu: BaSTA, demography a Rramas, který si představíme
v následující kapitole. Sofistikovanější modelování přežívání je dostupné v příkazu survreg
(balíček survival), která však vyžaduje individuální historii všech jedinců v populaci (viz.
např. Therneau & Grambsch 2000).
R
0 1 2 3 4
8. DEMOGRAFIE VĚKOVĚ STRUKTUROVANÉ POPULACE
Tahák
x<-c(0,1,2,3,4)
lx<-c(1,0.9,0.7,0.4,0.07)
plot(x,lx,log="y",type="b")
mx<-c(0,3,6,8,4)
R0<-sum(lx*mx);R0
T<-sum(x*lx*mx)/R0;T
lx[-1]/lx[-5]
lx[-1]/lx[-5]*mx[-1]
8. DEMOGRAFIE VĚKOVĚ STRUKTUROVANÉ POPULACE
Popis
Populace myši vykazuje přemnožení. Demografickou analýzou odhadněte osud populace.
Myš se množí kontinuálně a má několik generací v roce. Při detailním studiu životního cyklu
byly věkové kategorie rozlišeny v 3měsíčním intervalu.
Data
V tabulce jsou hodnoty standardizovaného přežívání (lx) a plodnosti (mx) pro každou věkovou
kategorii (x).
x lx mx
0 1 0
1 0.8 4
2 0.5 6
3 0.3 9
4 0 0
Úkoly
1/ Porovnejte plodnost a reprodukční hodnotu každé kategorie. Nalezněte stabilní věkové
rozdělení.
2/ Zjistěte hodnotu vnitřní míry růstu.
3/ Použitím simulačního modelu předpovězte, jak se vyvine velikost populace v průběhu
následujících 10 let, pokud nyní (v čase t = 0) je populační početnost následující: N0 =
(30, 20, 10, 5).
Řešení
1/ Do vektoru s názvem lx vložte hodnoty standardizovaného přežívání. K výpočtu
9
Analýza věkověstrukturované
populace
reprodukční hodnoty potřebujete znát přechodovou matici A. Pro její sestavení je nutné
spočítat px a Fx. Tyto se pro kontinuální rozmnožování počítají podle vztahů:
xx
xx
x
ll
ll
p
+
+
=
−
+
1
1
a
( )
2
1
1
++
= xxx
x
mpm
lF pro x > 0.
Hodnoty p0 a F0 nejsou definovány. Dále p4 = 0 a F4 = 0.
Zakreslete odhady px a Fx do schématu životního cyklu:
Sestavte z odhadů px a Fx přechodovou matici A:














=
............
............
............
............
A
Reprodukční hodnoty vx se počítají z levého charakteristického vektoru (eigenvector)
přechodové matice A, který přísluší k dominantnímu charakteristickému číslu (eigenvalue).
Načtěte balíček Rramas. Sestrojte matici A příkazem rbind a převeďte ji na přechodovou
matici (tA) příkazem as.tmatrix.
Vytvořte graf závislosti mx a vx na x. Hodnoty vx zjistíte příkazem summary a jejich grafické
zobrazení příkazem plot. Reprodukční hodnota vx nejmladší věkové třídy je rovna 1 a
hodnoty ostatních tříd jsou vyjádřeny v jednotkách reprodukční hodnoty nejmladší věkové
třídy.
9. ANALÝZA VĚKOVĚ-STRUKTUROVANÉ POPULACE
R
R
1 2 3 4
Překreslete graf do sešitu.
0 1 2 3 4
02468
Plodnosti a reprodukční hodnoty
x
mx,vx
Shodují se trendy mx a vx? Pokud se neshodují, odůvodněte proč?
…………………………………………………………………………………………………
…………………………………………………………………………………………………
Stabilní věkové rozdělení (SCD), tj. vektor poměru relativní početnosti jednotlivých stádií, se
počítá z pravého charakteristického vektor přechodové matice příslušející k dominantnímu
charakteristickému číslu. Opište je ze sumarizačního výstupu této matice.
Odhady SCD zapište do tabulky:
R
R
9. ANALÝZA VĚKOVĚ-STRUKTUROVANÉ POPULACE
x SCD
1
2
3
4
2/ Vnitřní míru populačního růstu za předpokladu stabilního věkového rozložení spočtěte
jako logaritmus dominantního charakteristického čísla (λ) přechodové matice tA. Hodnotu λ
vyčtete ze sumarizačního výstupu:
)ln(λ=r
Vnitřní míra populačního růstu, r = ……………
3/ K simulaci velikosti populace po dobu 10 let použijte maticový model pro
exponenciální růst:
0
t
t NAN = , 10.,2,1, K=t
Početnost v roce nula je N0 = (30, 20, 10, 5). K simulaci použijte příkaz projectn, ve které
nastavíte čas. Příkazem plot a legend vytvořte graf.
Překreslete graf do sešitu.
R
R
9. ANALÝZA VĚKOVĚ-STRUKTUROVANÉ POPULACE
0 2 4 6 8 10
0.0e+005.0e+071.0e+081.5e+082.0e+08
Simulace početnosti v průběhu 10 let
time
abundance
Poznámka
U organismů s reprodukčními pulsy, postreprodukčním sčítáním a s věkovou třídou 0 se levý
charakteristický vektor Leslieho matice nerovná vektoru reprodukčních hodnot (vx), ale
reziduálním reprodukčním hodnotám. Reprodukční hodnoty získáte až přičtením současné
reprodukce (mx).
Analýza strukturovaných populací je podrobně popsána v Caswell (2001) a Morris & Doak
(2002). Pro detailní analýzu populace z přechodové matice v prostředí R doporučujeme
balíček popbio, který obsahuje mnoho dalších užitečných příkazů k uvedeným pracím.
Tahák
lx<-c(1,0.8,0.5,0.3,0)
mx<-c(0,4,6,9,0)
px<-(lx[2:4]+lx[3:5])/(lx[1:3]+lx[2:4])
Fx<-sqrt(0.8)*(mx[2:4]+px*mx[3:5])/2
A<-rbind(c(Fx,0),diag(px,ncol=4))
library(Rramas)
tA<-as.tmatrix(A)
summary(tA)
x<-0:4
plot(x,mx,ylab="mx, vx")
9. ANALÝZA VĚKOVĚ-STRUKTUROVANÉ POPULACE
plot(tA)
log(4.61773)
N0<-c(30,20,10,5)
sim<-projectn(v0=N0,mat=tA,time=10)
plot(sim,sum=F)
legend("topleft",legend=c(1:4),lty=1:4,col=1:4)
9. ANALÝZA VĚKOVĚ-STRUKTUROVANÉ POPULACE
Popis
Populace vzácného druhu motýla vykazuje pokles. Byla provedena podrobná analýza
životního cyklu tohoto druhu, aby se zjistilo, který faktor způsobuje nejvýraznější pokles
početnosti. Na základě tohoto zjištění je třeba sestavit plán ochrany druhu.
Data
Pro každý instar (x) jsou v tabulce uvedeny hodnoty standardizovaného přežívání (lx) a
plodnosti (mx) a název faktoru, který způsobil největší mortalitu instaru. Odhad plodnosti byl
proveden rozborem ovarií, jde tedy o předreprodukční census.
x lx mx Faktor
Vajíčko 1.0 0 přezimování
Larva 1 0.7 0 parazitoidi
Larva 2 0.3 0 predace ptáky
Kukla 0.25 0 redukce biotopu
Dospělec 0.02 800
Úkoly
1/ Vytvořte přechodovou matici. Zjistěte hodnotu konečné míry růstu populace.
2/ Proveďte analýzu senzitivity a elasticity a určete nejdůležitější faktory. Navrhněte plán
na ochranu motýle.
Řešení
1/ Do vektoru s názvem lx umístěte hodnoty standardizovaného přežívání a do vektoru
mx plodnosti. Spočtěte px a Fx podle vzorců pro předreprodukční census:
10
Analýza stádiověstrukturované
populace
x
x
x
l
l
p 1+
= , xx mpF 0= .
Zakreslete odhady px a Fx do schématu životního cyklu:
Sestavte z odhadů px a Fx přechodovou matici A:
















=
....................
....................
....................
....................
....................
A
Hodnota konečné míry růstu populace, λ, je dominantní charakteristické číslo přechodové
matice (λ). Načtěte balíček Rramas. Sestrojte matici A příkazem rbind a převeďte ji na
přechodovou matici (tA) příkazem as.tmatrix. Hodnotu λ zjistíte příkazem summary.
Konečné míra růstu λ = ……………
2/ Sensitivita, neboli míra změny v populačním růstu λ vzhledem ke změně v přechodové
matici A, se počítá pro každý z elementů aij s využitím vektoru reprodukčních hodnot (vx) a
vektoru stabilního věkového rozdělení (SCD). Výsledkem je matice S se stejnou dimenzí jako
přechodová matice A. Ze sensitivity se pak spočte elasticita e, která měří relativní příspěvky
k populačnímu růstu λ. Hodnoty sensitivit i elasticit vyčtěte ze sumarizačního výstupu.
Přepište hodnoty elasticit pro jednotlivé parametry přechodové matice do tabulky:
R
R
V L1 L2 K D
10. ANALÝZA STÁDIOVĚ-STRUKTUROVANÉ POPULACE
eij
p1
p2
p3
p4
F4
Který faktor ovlivňuje růst populace nejvíce? ……………………………
Popište možný záchranný plán:
…………………………………………………………………………………………………
…………………………………………………………………………………………………
…………………………………………………………………………………………………
Tahák
lx<-c(1,0.7,0.3,0.25,0.02)
mx<-c(0,0,0,0,800)
px<-lx[2:5]/lx[1:4]
Fx<-px[1]*mx[2:5]
A<-rbind(c(0,Fx),diag(x=px,ncol=5,nrow=4))
library(Rramas)
tA<-as.tmatrix(tA)
summary(tA)
10. ANALÝZA STÁDIOVĚ-STRUKTUROVANÉ POPULACE
Popis
V laboratorních podmínkách byl sledován vývoj žlabatky v závislosti na teplotě. Z rozmezí
odpovídajícímu teplotám v přirozeném prostředí bylo nastaveno 6 konstantních teplot ze
škály 5 – 30 °C. Pro každou teplotu byla zaznamenána délka celkového vývoje u několika
jedinců.
Data
Průměrná délka celkového vývoje (d) ve dnech, tj. od nakladení vajíček po imágo, pro 6
teplot (t) je zaznamenaná v tabulce. Při teplotě 5 °C se žlabatky nevyvíjely.
t [°C] d
5 -
10 49
15 22
20 16
25 12
30 9
Úkoly
1/ Proložte daty lineární model. Určete spodní práh vývoje a sumu efektivních teplot.
2/ Zjistěte, za kolik dnů se objeví imága tohoto škůdce, pokud jste právě zaznamenali
kladení a znáte průměrné denní teploty v následujících 15 dnech: 17, 18, 21, 23, 24, 25,
23, 24, 21, 25, 22, 25, 26, 22 a 23 °C.
Řešení
1/ Do vektoru s názvem t umístěte teploty a do vektoru s názvem d délku vývoje. Pro
každou teplotu spočtěte rychlost vývoje, v, podle vztahu:
Lineární teplotní model11
d
v
1
= .
Rychlost vývoje pak vyneste do grafu v závislosti na teplotě. Pro teplotu 5 °C bude rychlost
vývoje 0.
Překreslete graf do sešitu.
5 10 15 20 25 30
0.000.020.040.060.080.10
Závislost vývoje na teplotě
temperature
v
Daty proložte lineární model, příkaz lm. Přímku do grafu vložíte příkazem abline.
Hodnoty parametrů zjistíte příkazem coef.
Překreslete výsledný model do grafu výše. Z odhadnutého modelu spočítejte spodní práh a
sumu efektivních teplot podle vztahů:
β
α
ˆ
ˆ
min −=t a
βˆ
1
=S ,
R
R
11. LINEÁRNÍ TEPLOTNÍ MODEL
kde αˆ a βˆ jsou odhady koeficientů proloženého modelu.
Spodní práh, tmin = ………… °C.
Suma efektivních teplot, S = ……… denních stupňů.
2/ Spočtěte hodnoty efektivních teplot pro každý den a do grafu je vyneste kumulativně
v závislosti na čase, tedy pro 15 sledovaných dnů podle vzorce:
∑=
−=
n
i
ttS
1
min
Do vektoru s názvem temp vložte zaznamenané teploty. Odhadněte za kolik dnů se objeví
imága s použitím příkazu cumsum.
Překreslete graf do sešitu.
2 4 6 8 10 12 14
50100150200250
Kumulativní suma teplot
days
cumulativetemperature
11. LINEÁRNÍ TEPLOTNÍ MODEL
R
R
Imága se objeví za ………… dnů.
Tahák
t<-c(5,10,15,20,25,30)
d<-c(0,49,22,16,12,9)
v<-1/d
v[1]<-0
plot(t,v, xlab="temperature")
m<-lm(v~t)
coef(m)
abline(m)
-(-0.022336/0.004351)
1/0.004351
temp<-c(17,18,21,23,24,25,23,24,21,25,22,25,26,22,23)
ET<-temp-5.13
plot(cumsum(ET),type="s", xlab="day", ylab="cumulative temperature")
abline(229,0,lty=2)
11. LINEÁRNÍ TEPLOTNÍ MODEL
Popis
V laboratoři byl sledován vývoj klopušky za konstantních teplotních podmínek. Bylo
nastaveno 7 teplot v rozmezí 18 až 32 °C. Pro každou teplotu bylo uskutečněno několik
opakování a ze zjištěných hodnot byla spočtena délka larválního vývoje.
Data
Průměrná délka vývoje (d) ve dnech pro 7 teplot (t) je zaznamenaná v tabulce:
t [°C] d
18 23.5
20 18.5
22 13
25 7.3
28 5.5
30 5
32 10.9
Úkoly
1/ Vykreslete data do grafu a proložte je Lactinovým modelem. Odhadněte spodní a horní
práh vývoje a optimální teplotu pro vývoj.
Řešení
1/ Do vektoru s názvem t vložte teploty a do vektoru d délku vývoje. Vykreslete rychlost
vývoje (v) v závislosti na teplotě, kterou z dat v tabulce spočtete podle vztahu:
d
v
1
=
12Nelineární teplotní model
Překreslete graf do sešitu:
18 20 22 24 26 28 30 32
0.050.100.150.20
Závislost rychlosti vývoje na teplotě
temperature
v
Rychlost vývoje modelujte jako funkci času pomocí Lactinova modelu se čtyřmi parametry
(ρ, Tmax, ∆, φ).
ϕ
ρ
ρ
+−= ∆
−
−
tT
T
t
eetv
max
max
)(
Parametr ρ ovlivňuje míru růstu, Tmax je teplota horního prahu vývoje, ∆ určuje tvar křivky a φ
posunuje křivku po ose y a tím umožňuje odhadnout teplotu spodního prahu vývoje.
Z naměřených hodnot odvoďte startovací hodnoty parametrů Tmax a φ analyticky.
12. NELINEÁRNÍ TEPLOTNÍ MODEL
R
Startovací hodnoty parametrů ρ a ∆ zjistěte iteracemi použitím příkazu curve:
Startovací hodnoty budou: ρ = ………, Tmax =………, ∆ = ………, φ = ………
K odhadu použijte nelineární regresi, příkaz nls. Odhady parametrů zjistíte příkazem coef.
Odhadnutou křivku vložíte do grafu příkazem lines s argumentem predict.
Modelovou křivku přidáme do grafu výše.
Spodní a horní práh vývoje vypočítejte z odhadnutého modelu jako průsečíky s osou x, tedy
teploty, při nichž je rychlost vývoje nulová (v(t) = 0). Využijte balíček rootSolve a příkaz
uniroot.all:
ϕ
ρ
ρ
+−= ∆
−
−
tT
T
t
ee
max
max
0
Spodní práh vývoje je ………°C.
Horní práh vývoje je ………°C.
Optimální teplotu vývoje odhadněte pomocí první derivace modelové funkce. Extrém funkce
najděte pomocí příkazu uniroot:
∆
−
−
∆
−=
∂
∂
tT
T
t
ee
t
tv max
max1)( ρ
ρ
ρ
Optimální teplota je ………°C.
R
R
R
R
12. NELINEÁRNÍ TEPLOTNÍ MODEL
Poznámka
Pro popis nelineární závislosti rychlosti vývoje na teplotě bylo navrženo mnoho modelů
s různými počty parametrů. Jedním ze základních modelů je x
xey −= , v obecném tvaru
decxay cxb
+−−= −+ )(
)( , kde c definuje průsečík s vodorovnou asymptotou. Tato teplotní
křivka je vlastně zrcadlově otočená křivka kladení (kap. 7) – všimněte si opačných znamének
u parametrů a a b. Další běžně používané nelineární modely jsou Janischův, Loganův-10,
Taylorův, nebo Briereův (Kontodimas et al. 2004, Roy et al., 2002, Walgama & Zalucki
2006).
Tahák
t<-c(18,20,22,25,28,30,32)
d<-c(23.5,18.5,13,7.3,5.5,5,10.9)
v<-1/d
plot(t,v, xlab="temperature")
curve(exp(1*x)-exp(1*32-((32-x)/1))-1, xlim=c(0,40),ylim=c(0,1))
curve(exp(0.1*x)-exp(0.1*32-((32-x)/1))-1, xlim=c(0,40),ylim=c(0,1))
curve(exp(0.01*x)-exp(0.01*32-((32-x)/1))-1, xlim=c(0,40),ylim=c(0,1))
m<-nls(v~exp(rho*t)-exp(rho*Tm-(Tm-t)/delta)+phi,
start=c(rho=0,Tm=32,delta=1,phi=0))
coef(m)
x<-seq(15,40,0.1)
plot(t,v,xlim=c(15,35),ylim=c(0,0.22))
lines(x,predict(m,list(t=x)))
abline(h=0,lty=2)
library(rootSolve)
tminmax<-uniroot.all(function(x) exp(0.011*x)-exp(0.011*33.681
-(33.681-x)/0.74)-1.195,lower=0,upper=40);tminmax
topt<-uniroot(function(x) 0.011*exp(0.011*x)-(1/0.74)
*exp(0.011*33.681-(33.681-x)/0.74),lower=0,upper=40);topt
12. NELINEÁRNÍ TEPLOTNÍ MODEL
Popis
V průběhu 15 let byla na jaře na poli pšenice zaznamenávána velikost populace mšic. Na 20
rostlinách byli vždy na začátku léta spočteny všichni jedinci. Stejný postup se opakoval každý
rok. Zajímá nás jaký typ dynamiky populace mšic vykazovala.
Data
V jednotlivých letech byly zaznamenány následující početnosti (N):
Rok N
1 95
2 134
3 180
4 531
5 277
6 296
7 528
8 329
9 397
10 572
11 625
12 318
13 567
14 681
15 386
Úkoly
1/ Vytvořte graf populační dynamiky mšic a zjistěte, je-li početnost závislá na hustotě.
2/ Odhadněte maximální konečnou míru růstu a nosnou kapacitu prostředí.
3/ Nasimulujte dynamiku populace mšic po dobu 20 generací. Použijte model pro diskrétní
růst populace závislý na hustotě s odhadnutou nosnou kapacitou prostředí. Hodnotu
konečné míry růstu simulujte náhodně z odhadnutých hodnot. Početnost mšic z
posledního roku použijte jako počáteční početnost.
Diskrétní hustotně-závislý
růst populace13
Řešení
1/ Do vektoru s názvem aphid vložte početnosti mšic. Vykreslete početnosti do
spojnicového grafu.
Překreslete graf do sešitu:
2 4 6 8 10 12 14
100200300400500600700
Populační dynamika mšic
year
aphid
Ke zjištění závislosti početnosti na hustotě je potřeba vynést logaritmy roční populační míry
růstu tλ proti početnosti populace, Nt. Tu spočtěte podle vztahu:
t
t
t
N
N 1+
=λ , kde t = 1, ..., 14.
Daty pak proložte lineární regresní model, příkaz lm. Odhady parametrů zjistíte příkazem
coef.
Překreslete graf do sešitu:
13. DISKRÉTNÍ HUSTOTNĚ-ZÁVISLÝ RŮST POPULACE
R
R
100 200 300 400 500 600 700
-0.50.00.51.0
Závislost ln(lambda) na N
aphid[-15]
log(lambda2)
Je v populační dynamice přítomná hustotní závislost? Pokud ano, vysvětlete proč.
…………………………………………………………………………………………………
…………………………………………………………………………………………………
2/ Parametry λmax a K lze odvodit z odhadů parametrů lineárního modelu podle vzorců:
α
λ ˆ
max e= a
β
α
ˆ
ˆ
−=K
kde α a β jsou parametry lineárního modelu. Body proložte odhadnutý lineární model
příkazem abline.
Hodnota maximální konečné míry růstu, λmax = ………
Hodnota nosné kapacity prostředí, K = …………
R
13. DISKRÉTNÍ HUSTOTNĚ-ZÁVISLÝ RŮST POPULACE
3/ K simulaci populační dynamiky použijte model pro diskrétní logistický růst závislý na
hustotě:





 −
+
=+
K
N
N
N
t
t
t
1
1
1
λ
λ
Míra populačního růstu λ se bude v každém kroku generovat náhodně z rovnoměrného
rozdělení na intervalu (1, λmax). Využijte příkaz runif.
Překreslete graf do sešitu:
0 5 10 15 20
Simulace početnosti v průběhu 20 generací
generation
N
Jaký typ populační dynamiky simulace vykazuje?
…………………………………………………………………………………………………
R
13. DISKRÉTNÍ HUSTOTNĚ-ZÁVISLÝ RŮST POPULACE
Poznámka
K modelování hustotní závislosti lze použít i různé typy logistických funkcí, jako například
d
e
a
y cxb
+
+
= −− )(
1
. Tato funkce docela dobře popisuje sigmoidní charakter početností
populace: v počáteční fázi je růst přibližně exponenciální, později s rostoucím nasycením se
zpomaluje a nakonec se asymptoticky zastaví. Parametr a definuje asymptotu, b udává
rychlost růstu křivky a jeho znaménko mění růst funkce na klesání, parametry c a d udávají
posunutí.
Tahák
aphid<-c(95,134,180,531,277,296,528,329,397,572,625,318,567,681,386)
plot(aphid,type="b", xlab="yeras")
lambda2<-aphid[-1]/aphid[-15]
plot(aphid[-15],log(lambda2))
m<-lm(log(lambda2)~aphid[-15])
coef(m)
abline(m,col=2)
exp(0.903429139)
-0.903429139/-0.002033643
N<-1:21
N[1]<-386
for(t in 1:20) N[t+1]<-N[t]*runif(1,1,2.47)/(1+N[t]*(runif(1,1,2.47)
-1)/444)
plot(0:20,N,type="b", xlab="generation")
13. DISKRÉTNÍ HUSTOTNĚ-ZÁVISLÝ RŮST POPULACE