In ToonVanDaele/metapop:
knitr::opts_chunk$set(echo = TRUE)
library(popbio)
library(Rramas)
Situering
In de paper A METAPOPULATION MODEL OF THE PEREGRINE FALCON IN CALIFORNIA: VIABILITY AND MANAGEMENT STRATEGIES beschrijven Wootton en Bell een metapopulatie model dat gebruik maakt van projectiematrices voor de modellering van een de populatie slechtvalken in California. In wat volgt zal deze modellering beschreven worden en een eigen simulatie van het model zal uitgevoerd worden.
Modellering van een enkele populatie
De auteurs geven eerst een model zonder ruimtelijke component, van de volgende vorm:
$$ \left( \begin{array}{c}
F(t+1) \
N(t+1)
\end{array}\right) = \left( \begin{array}{cc}
0 & yb \
r & s
\end{array}\right) \times \left( \begin{array}{c}
F(t) \
N(t)
\end{array}\right) + \left( \begin{array}{c}
yI/2 \
0
\end{array}\right)$$
Hierbij is $F(t)$ de populatie jongen op tijdstip $t$, $N(t)$ de populatie (broedende) volwassen vogels op tijdstip $t$ en $I$ is een vast aantal in gevangenschap gekweekte jongen die jaarlijks in de populatie worden losgelaten als middel om de populatie te ondersteunen. De parameters in dit model zijn $y$, de kans dat net uitgebroedde jongen een jaar overleven, $b$, het aantal vrouwelijke jongen dat 1 volwassen vogel per jaar voortbrengt, $r$, de kans dat de jongen nog een jaar overleven en een broedplaats bemachtigen en $s$, de jaarlijkse overlevingskans van de volwassen vogels.
We kunnen dit met de volgende code simuleren (de waarden voor de parameters en de beginwaarden zijn afkomstig uit de paper):
init_1 <- c(55,90,0)
init_2 <- c(55,90,50)
classes <- c("F","N","I")
b <- 0.69
y <- 0.36
r <- 0.72
s <- 0.77
A <- matrix(data = c(0,y*b,y/2,r,s,0,0,0,1), nrow = 3, ncol = 3, byrow = TRUE)
dimnames(A)<- list(classes,classes)
A
p.p_1<-pop.projection(A,init_1,iterations=500)
p.p_2<-pop.projection(A,init_2,iterations=500)
p.p_2$pop.sizes <- p.p_2$pop.sizes-50
Merk op dat we nu een $3\times 3$ matrix $A$ gebruikten. Vermenigvuldiging van een vector $\left( F(t) N(t) I\right)^T$ met deze matrix geeft echter hetzelfde resultaat voor $F(t+1)$ en $N(t+1)$ als de vergelijking onder de vorm waarmee we het model hebben uitgedrukt.
De reproductive rate is r p.p_1$lambda en de stable stage verdeling is r p.p_1$stable.stage.
Een sensitiviteits- en elasticiteits-analyse van de $2\times 2$ submatrix geeft de volgende resultaten:
sensitivity(A[1:2,1:2])
elasticity(A[1:2,1:2])
We kunnen dezelfde resultaten bekomen met behulp van het Rramas pakket:
v0 <- init_1[1:2]
mat <- A[1:2,1:2]
matT <- as.tmatrix(mat)
summary(matT)
We kunnen ook de evolutie van de populatie plotten (p.p_1 geeft een situatie weer waarbij er geen gekweekte jongen aan de populatie worden toegevoegd en p.p_2 geeft een situatie weer waarbij er jaarlijks 50 gekweekte jongen worden toegevoegd aan de populatie):
plot(p.p_1$pop.sizes,col=3, ylim=c(0,250))
points(p.p_2$pop.sizes,col=4)
legend("topright",legend=c("Geen toevoeging", "50 extra jongen per jaar"), col =3:4, lty=5)
Ook deze resulataten kunnen verkregen worden via Rramas.
p.p_ramas1 <- projectn(v0,mat,time=500)
plot(p.p_ramas1)
I <- 50
management <- c(y*I/2,0)
p.p_ramas2 <- projectn(v0,mat,time=500, management = management)
plot(p.p_ramas2)
Modellering van een metapopulatie
Vervolgens ontwikkelden de auteurs ook een metapopulatie model waarin de slechtvalken populatie in California in twee subpopulaties werd opgedeeld, een noordelijke en een zuidelijke populatie. Deze zullen vanaf nu aangeduid worden met subscripten $n$, resp. $s$.
Het model ziet er als volgt uit:
$$ \left( \begin{array}{c}
F_n(t+1) \
N_n(t+1) \
F_s(t+1) \
N_s(t+1)
\end{array}\right) = \left( \begin{array}{cccc}
0 & y_nb_n & 0 & 0 \
r_nh & s_n & r_nmc & 0\
0 & 0 & 0 & y_sb_s\
r_smc & 0 & r_sh & s_s
\end{array}\right) \times \left( \begin{array}{c}
F_n(t) \
N_n(t) \
F_s(t) \
N_s(t)
\end{array}\right) + \left( \begin{array}{c}
fy_nI/2 \
0 \
(1-f)y_sI/2\
0
\end{array}\right)$$
Nieuwe parameters zijn: $m$, de migratiekans van een vrouwtje dat nog geen broedgebied heeft, $c$, de kans om de migratie te overleven, $h=1-m$ en $f$, de proportie van de in gevangenschap gekweekte jongen die in de noordelijke populatie terecht komt. De andere parameters zijn dezelfde als voordien, maar zijn nu voorzien van een subscript omdat hun waarden kunnen verschillen in beide subpopulaties.
Deze matrix bouwen we als volgt op:
#Parameters
meta_classes <- c("Fn","Nn","Fs","Ns","I")
bn<-0.71
bs<-0.53
yn<-ys<-0.36
rn<-rs<-0.72
#sn<-ss<-0.77
sn<-0.91
ss<-0.77
m<-0.27
f<-1
c<-1
#Opbouw Matrix
An <- matrix(data = c(0,yn*bn,rn,sn), nrow = 2, ncol = 2, byrow = TRUE)
As <- matrix(data = c(0,ys*bs,rs,ss), nrow = 2, ncol = 2, byrow = TRUE)
An[2,1]<-An[2,1]*(1-m)
As[2,1]<-As[2,1]*(1-m)
Dsn <- matrix(data = c(0,0,rn*m*c,0), nrow = 2, ncol = 2, byrow = TRUE)
Dns <- matrix(data = c(0,0,rs*m*c,0), nrow = 2, ncol = 2, byrow = TRUE)
Fledg <- c(f*yn/2,0,(1-f)*ys/2,0)
G <- rbind(cbind(An,Dsn),cbind(Dns,As))
G <- rbind(cbind(G,Fledg),c(0,0,0,0,1))
dimnames(G)<- list(meta_classes,meta_classes)
G
We kunnen opnieuw een sensitiviteitsanalyse doen:
sensitivity(G[1:4,1:4])
elasticity(G[1:4,1:4])
# Sensitiviteit en Elasticiteit tov de parameters ipv matrix-elementen:
VitalRates <- list(bn=0.71,bs=0.53, yn=0.36, ys=0.36, rn=0.72, rs=0.72, sn = 0.77, ss=0.77, m=0.27)
matrix <- expression(
0, yn*bn, 0, 0,
rn*(1-m), sn, rn*m, 0,
0,0,0,ys*bs,
rs*m,0,rs*(1-m), ss)
vitalsens(matrix, VitalRates)
Gt <- as.tmatrix(G[1:4,1:4])
summary(Gt)
plot(Gt)
Deze analyse geeft andere resultaten dan de sensitiviteiten en elasticiteiten die gegeven worden in de paper. De auteurs bepaalden echter de sensitiviteit en elasticiteit echter via simulatie en niet via de matrixeigenschappen.
Dit kunnen we ook narekenen (De berekening van de benodigde verhogingen/verlagingen van de parameters om aan een groei-ratio van 1 te komen laten we achterwege en we nemen gewoon de waarden uit Tabel 2 uit de paper over):
dx = c(bn=1.1-0.71, bs= 1.02-0.53, yn=0.56-0.36, ys=0.69-0.36, rn=1.11-0.72,
rs=1.38-0.72, sn=0.85-0.77, ss=0.88-0.77, m=-0.18-0.27)
dR = 1-0.9438
sens<-dR/dx
sens
x<-c(bn,bs,yn,ys,rn,rs,sn,ss,m)
R=0.9438
elast<-(x/R)*sens
elast
Enkele simulaties:
init_1_met <- c(19,60,20,29,0)
init_2_met <- c(19,60,20,29,50)
p.p_1<-pop.projection(G,init_1_met,iterations=100)
p.p_2<-pop.projection(G,init_2_met,iterations=100)
p.p_2$pop.sizes <- p.p_2$pop.sizes-50
plot(p.p_1$pop.sizes,col=3)
points(p.p_2$pop.sizes,col=4)
stage.vector.plot(p.p_2$stage.vectors[1:4,], proportions = TRUE)
Opnieuw kunnen de simulaties ook uitgevoerd worden met Rramas
meta_rramas1 <- projectn(init_2_met[1:4], Gt, time=100)
meta_rramas2 <- projectn(init_1_met[1:4], Gt, time=100, management=Fledg)
plot(meta_rramas1, ylim=c(0,10000))
plot(meta_rramas2, ylim=c(0,10000))
plot(meta_rramas2, sum = FALSE, ylim=c(0,10000))
Merk op dat dit de situatie is wanneer alle gekweekte jongen in de noordelijke populatie terecht komen ($f=1$ in het begin)
Densiteitsafhankelijkheid
Wanneer de bovenstaande figuur vergelijken met de figuren in de paper, vinden we deze niet terug. De oorzaak hiervan is dat, hoewel het voorgaande model wel vermeld werd, de auteurs ook rekening hielden met een densiteitsafhankelijke populatie. Dit deden ze als volgt:
Zodra de noordelijke subpopulatie een totale grootte bereikte van $T_n=100$ individuen werd een ander matrix-model gebruikt, namelijk:
$$ \left( \begin{array}{c}
F_n(t+1) \
N_n(t+1) \
F_s(t+1) \
N_s(t+1)
\end{array}\right) = \left( \begin{array}{cccc}
s_nh & y_nb_n & s_nmc & 0 \
0 & s_n & 0 & 0\
0 & 0 & 0 & y_sb_s\
r_smc & 0 & r_sh & s_s
\end{array}\right) \times \left( \begin{array}{c}
F_n(t) \
N_n(t) \
F_s(t) \
N_s(t)
\end{array}\right) + \left( \begin{array}{c}
-(1-s_n)T_n \
(1-s_n)T_n \
0\
0
\end{array}\right)$$
Dit model stelt een situatie voor waarbij de $F_n(t)$-klasse niet enkel jongen bevat, maar ook zogenaamde 'Floaters', volwassen vogels die nog geen broedgebied konden bemachtigen. Jaarlijks zal een deel van deze floaters de broedgebieden innemen die vrijgekomen zijn door het overlijden van broedende volwassen vogels. Merk ook op dat, zodra de maximale densiteit eenmaal bereikt wordt, er geen in gekweekte jongen meer worden toegevoegd aan de populatie.
We implementeren het densiteitsafhankelijke model als volgt:
cap<-100
An2 <- matrix(data=c(sn*(1-m), yn*bn,0,sn),nrow=2,ncol=2, byrow = TRUE)
Dsn2 <- matrix(data = c(sn*m*c,0,0,0), nrow = 2, ncol = 2, byrow = TRUE)
float <- c(-(1-sn), (1-sn),0,0)
G2 <- rbind(cbind(An2,Dsn2),cbind(Dns,As))
G2 <- rbind(cbind(G2,float),c(0,0,0,0,1))
dimnames(G2)<- list(meta_classes,meta_classes)
G2
Om dit model te simuleren, kunnen we niet langer gebruik maken van het commando pop.projection uit het popbio-pakket, daarom passen we de onderliggende code van dit commando aan als volgt:
pop.projection.WB <- function(A1,A2,cap,n0,iterations=20){
x <- length(n0)
t <- iterations
stage <- matrix(numeric(x * t), nrow = x)
pop <- numeric(t)
change <- numeric(t - 1)
TN<-numeric(t)
TS<-numeric(t)
Iter <- numeric(t)
n<-n0
for (i in 1:t) {
Tn <- sum(n[2])
TN[i]<-Tn
TS[i]<-sum(n[4])
if (Tn<cap) {
stage[, i] <- n
pop[i] <- sum(n[1:x-1])
if (i > 1) {
change[i - 1] <- pop[i]/pop[i - 1]
}
Iter[i]<-1
n <- A1 %*% n
}else{
n[5]<-cap
stage[, i] <- n
pop[i] <- sum(n[1:x-1])
if (i > 1) {
change[i - 1] <- pop[i]/pop[i - 1]
}
Iter[i]<-2
n <- A2 %*% n
}
}
rownames(stage) <- rownames(A1)
colnames(stage) <- 0:(t - 1)
w <- stage[, t]
pop.proj <- list(lambda = pop[t]/pop[t - 1], stable.stage = w/sum(w),
stage.vectors = stage, pop.sizes = pop, pop.changes = change,
North = TN, South = TS, MUsage = Iter)
pop.proj}
Simulatie geeft ons de volgende resultaten:
pp_dens <- pop.projection.WB(G,G2,100,init_1_met,iterations=100)
plot(pp_dens$pop.sizes,col=3)
stage.vector.plot(pp_dens$stage.vectors[1:4,], proportions = FALSE)
NS_division <- rbind(colSums(pp_dens$stage.vectors[1:2,]),colSums(pp_dens$stage.vectors[3:4,]))
stage.vector.plot(NS_division,proportions = FALSE)
De auteurs testen ook wat er zou gebeuren moest een catastrofale gebeurtenis de gehele noordelijke subpopulatie doen uitsterven:
Local_ext <- pp_dens$stage.vectors[,100]
Local_ext[5]<-0
Local_ext[1]<-0
Local_ext[2]<-0
pp_dens2 <- pop.projection.WB(G,G2,100,Local_ext,iterations = 100)
toPlot <-cbind(pp_dens$stage.vectors[1:4,], pp_dens2$stage.vectors[1:4,])
dimnames(toPlot) <- list(meta_classes[1:4], c(0:199))
stage.vector.plot(toPlot, proportions = FALSE)
#NS_division2 <- rbind(colSums(toPlot[1:2,]),colSums(toPlot[3:4,]))
NS_division2 <- rbind(toPlot[2,], toPlot[4,])
stage.vector.plot(NS_division2,proportions = FALSE)
De rode curve stelt de noordelijke subpopulatie voor en de gele de zuidelijke. We zien dat de noordelijke populatie na de catastrofe kan gered worden door de zuidelijke populatie.
Tot slot simuleren we nog de verschillen in evolutie van de populatie wanneer er wel of geen migratie is en wat er zou gebeuren wanneer er geen source-sink interactie was tussen de noordelijke en zuidelijke subpopulaties.
### Figure 7(a)
An <- matrix(data = c(0,yn*bn,rn,sn), nrow = 2, ncol = 2, byrow = TRUE)
As <- matrix(data = c(0,ys*bs,rs,ss), nrow = 2, ncol = 2, byrow = TRUE)
Zero <- matrix(0L, nrow=2, ncol=2)
G3 <- rbind(cbind(An,Zero),cbind(Zero,As))
G3 <- rbind(cbind(G3,Fledg),c(0,0,0,0,1))
dimnames(G3)<- list(meta_classes,meta_classes)
G3
An2 <- matrix(data=c(sn, yn*bn,0,sn),nrow=2,ncol=2, byrow = TRUE)
G4 <- rbind(cbind(An2,Zero),cbind(Zero,As))
G4 <- rbind(cbind(G4,float),c(0,0,0,0,1))
dimnames(G4)<- list(meta_classes,meta_classes)
G4
pp_dens3 <- pop.projection.WB(G3,G4,100,init_1_met,iterations=100)
plot(pp_dens$North, col=1, ylim=c(0,200))
points(pp_dens$South, col=2)
points(pp_dens3$North, col=3)
points(pp_dens3$South, col=4)
legend("top",legend=c("North, migration", "South, migration",
"North, no migration", "South, no migration"), col =1:4, lty=5)
### Figure 7(b)
plot(pp_dens3$North, col=2, ylim=c(0,300))
points(pp_dens$North+pp_dens$South, col=3)
points(pp_dens3$North+pp_dens3$South, col=4)
legend("top",legend=c("No sink", "Migration", "No Migration"), col =2:4, lty=5)
ToonVanDaele/metapop documentation built on May 9, 2019, 5:11 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE) library(popbio) library(Rramas)
Situering
In de paper A METAPOPULATION MODEL OF THE PEREGRINE FALCON IN CALIFORNIA: VIABILITY AND MANAGEMENT STRATEGIES beschrijven Wootton en Bell een metapopulatie model dat gebruik maakt van projectiematrices voor de modellering van een de populatie slechtvalken in California. In wat volgt zal deze modellering beschreven worden en een eigen simulatie van het model zal uitgevoerd worden.
Modellering van een enkele populatie
De auteurs geven eerst een model zonder ruimtelijke component, van de volgende vorm: $$ \left( \begin{array}{c} F(t+1) \ N(t+1) \end{array}\right) = \left( \begin{array}{cc} 0 & yb \ r & s \end{array}\right) \times \left( \begin{array}{c} F(t) \ N(t) \end{array}\right) + \left( \begin{array}{c} yI/2 \ 0 \end{array}\right)$$
Hierbij is $F(t)$ de populatie jongen op tijdstip $t$, $N(t)$ de populatie (broedende) volwassen vogels op tijdstip $t$ en $I$ is een vast aantal in gevangenschap gekweekte jongen die jaarlijks in de populatie worden losgelaten als middel om de populatie te ondersteunen. De parameters in dit model zijn $y$, de kans dat net uitgebroedde jongen een jaar overleven, $b$, het aantal vrouwelijke jongen dat 1 volwassen vogel per jaar voortbrengt, $r$, de kans dat de jongen nog een jaar overleven en een broedplaats bemachtigen en $s$, de jaarlijkse overlevingskans van de volwassen vogels.
We kunnen dit met de volgende code simuleren (de waarden voor de parameters en de beginwaarden zijn afkomstig uit de paper):
init_1 <- c(55,90,0) init_2 <- c(55,90,50) classes <- c("F","N","I") b <- 0.69 y <- 0.36 r <- 0.72 s <- 0.77 A <- matrix(data = c(0,y*b,y/2,r,s,0,0,0,1), nrow = 3, ncol = 3, byrow = TRUE) dimnames(A)<- list(classes,classes) A p.p_1<-pop.projection(A,init_1,iterations=500) p.p_2<-pop.projection(A,init_2,iterations=500) p.p_2$pop.sizes <- p.p_2$pop.sizes-50
Merk op dat we nu een $3\times 3$ matrix $A$ gebruikten. Vermenigvuldiging van een vector $\left( F(t) N(t) I\right)^T$ met deze matrix geeft echter hetzelfde resultaat voor $F(t+1)$ en $N(t+1)$ als de vergelijking onder de vorm waarmee we het model hebben uitgedrukt.
De reproductive rate is r p.p_1$lambda en de stable stage verdeling is r p.p_1$stable.stage.
Een sensitiviteits- en elasticiteits-analyse van de $2\times 2$ submatrix geeft de volgende resultaten:
sensitivity(A[1:2,1:2]) elasticity(A[1:2,1:2])
We kunnen dezelfde resultaten bekomen met behulp van het Rramas pakket:
v0 <- init_1[1:2] mat <- A[1:2,1:2] matT <- as.tmatrix(mat) summary(matT)
We kunnen ook de evolutie van de populatie plotten (p.p_1 geeft een situatie weer waarbij er geen gekweekte jongen aan de populatie worden toegevoegd en p.p_2 geeft een situatie weer waarbij er jaarlijks 50 gekweekte jongen worden toegevoegd aan de populatie):
plot(p.p_1$pop.sizes,col=3, ylim=c(0,250)) points(p.p_2$pop.sizes,col=4) legend("topright",legend=c("Geen toevoeging", "50 extra jongen per jaar"), col =3:4, lty=5)
Ook deze resulataten kunnen verkregen worden via Rramas.
p.p_ramas1 <- projectn(v0,mat,time=500) plot(p.p_ramas1) I <- 50 management <- c(y*I/2,0) p.p_ramas2 <- projectn(v0,mat,time=500, management = management) plot(p.p_ramas2)
Modellering van een metapopulatie
Vervolgens ontwikkelden de auteurs ook een metapopulatie model waarin de slechtvalken populatie in California in twee subpopulaties werd opgedeeld, een noordelijke en een zuidelijke populatie. Deze zullen vanaf nu aangeduid worden met subscripten $n$, resp. $s$.
Het model ziet er als volgt uit: $$ \left( \begin{array}{c} F_n(t+1) \ N_n(t+1) \ F_s(t+1) \ N_s(t+1) \end{array}\right) = \left( \begin{array}{cccc} 0 & y_nb_n & 0 & 0 \ r_nh & s_n & r_nmc & 0\ 0 & 0 & 0 & y_sb_s\ r_smc & 0 & r_sh & s_s \end{array}\right) \times \left( \begin{array}{c} F_n(t) \ N_n(t) \ F_s(t) \ N_s(t) \end{array}\right) + \left( \begin{array}{c} fy_nI/2 \ 0 \ (1-f)y_sI/2\ 0 \end{array}\right)$$
Nieuwe parameters zijn: $m$, de migratiekans van een vrouwtje dat nog geen broedgebied heeft, $c$, de kans om de migratie te overleven, $h=1-m$ en $f$, de proportie van de in gevangenschap gekweekte jongen die in de noordelijke populatie terecht komt. De andere parameters zijn dezelfde als voordien, maar zijn nu voorzien van een subscript omdat hun waarden kunnen verschillen in beide subpopulaties.
Deze matrix bouwen we als volgt op:
#Parameters meta_classes <- c("Fn","Nn","Fs","Ns","I") bn<-0.71 bs<-0.53 yn<-ys<-0.36 rn<-rs<-0.72 #sn<-ss<-0.77 sn<-0.91 ss<-0.77 m<-0.27 f<-1 c<-1 #Opbouw Matrix An <- matrix(data = c(0,yn*bn,rn,sn), nrow = 2, ncol = 2, byrow = TRUE) As <- matrix(data = c(0,ys*bs,rs,ss), nrow = 2, ncol = 2, byrow = TRUE) An[2,1]<-An[2,1]*(1-m) As[2,1]<-As[2,1]*(1-m) Dsn <- matrix(data = c(0,0,rn*m*c,0), nrow = 2, ncol = 2, byrow = TRUE) Dns <- matrix(data = c(0,0,rs*m*c,0), nrow = 2, ncol = 2, byrow = TRUE) Fledg <- c(f*yn/2,0,(1-f)*ys/2,0) G <- rbind(cbind(An,Dsn),cbind(Dns,As)) G <- rbind(cbind(G,Fledg),c(0,0,0,0,1)) dimnames(G)<- list(meta_classes,meta_classes) G
We kunnen opnieuw een sensitiviteitsanalyse doen:
sensitivity(G[1:4,1:4]) elasticity(G[1:4,1:4]) # Sensitiviteit en Elasticiteit tov de parameters ipv matrix-elementen: VitalRates <- list(bn=0.71,bs=0.53, yn=0.36, ys=0.36, rn=0.72, rs=0.72, sn = 0.77, ss=0.77, m=0.27) matrix <- expression( 0, yn*bn, 0, 0, rn*(1-m), sn, rn*m, 0, 0,0,0,ys*bs, rs*m,0,rs*(1-m), ss) vitalsens(matrix, VitalRates)
Gt <- as.tmatrix(G[1:4,1:4]) summary(Gt) plot(Gt)
Deze analyse geeft andere resultaten dan de sensitiviteiten en elasticiteiten die gegeven worden in de paper. De auteurs bepaalden echter de sensitiviteit en elasticiteit echter via simulatie en niet via de matrixeigenschappen. Dit kunnen we ook narekenen (De berekening van de benodigde verhogingen/verlagingen van de parameters om aan een groei-ratio van 1 te komen laten we achterwege en we nemen gewoon de waarden uit Tabel 2 uit de paper over):
dx = c(bn=1.1-0.71, bs= 1.02-0.53, yn=0.56-0.36, ys=0.69-0.36, rn=1.11-0.72, rs=1.38-0.72, sn=0.85-0.77, ss=0.88-0.77, m=-0.18-0.27) dR = 1-0.9438 sens<-dR/dx sens x<-c(bn,bs,yn,ys,rn,rs,sn,ss,m) R=0.9438 elast<-(x/R)*sens elast
Enkele simulaties:
init_1_met <- c(19,60,20,29,0) init_2_met <- c(19,60,20,29,50) p.p_1<-pop.projection(G,init_1_met,iterations=100) p.p_2<-pop.projection(G,init_2_met,iterations=100) p.p_2$pop.sizes <- p.p_2$pop.sizes-50 plot(p.p_1$pop.sizes,col=3) points(p.p_2$pop.sizes,col=4) stage.vector.plot(p.p_2$stage.vectors[1:4,], proportions = TRUE)
Opnieuw kunnen de simulaties ook uitgevoerd worden met Rramas
meta_rramas1 <- projectn(init_2_met[1:4], Gt, time=100) meta_rramas2 <- projectn(init_1_met[1:4], Gt, time=100, management=Fledg) plot(meta_rramas1, ylim=c(0,10000)) plot(meta_rramas2, ylim=c(0,10000)) plot(meta_rramas2, sum = FALSE, ylim=c(0,10000))
Merk op dat dit de situatie is wanneer alle gekweekte jongen in de noordelijke populatie terecht komen ($f=1$ in het begin)
Densiteitsafhankelijkheid
Wanneer de bovenstaande figuur vergelijken met de figuren in de paper, vinden we deze niet terug. De oorzaak hiervan is dat, hoewel het voorgaande model wel vermeld werd, de auteurs ook rekening hielden met een densiteitsafhankelijke populatie. Dit deden ze als volgt:
Zodra de noordelijke subpopulatie een totale grootte bereikte van $T_n=100$ individuen werd een ander matrix-model gebruikt, namelijk: $$ \left( \begin{array}{c} F_n(t+1) \ N_n(t+1) \ F_s(t+1) \ N_s(t+1) \end{array}\right) = \left( \begin{array}{cccc} s_nh & y_nb_n & s_nmc & 0 \ 0 & s_n & 0 & 0\ 0 & 0 & 0 & y_sb_s\ r_smc & 0 & r_sh & s_s \end{array}\right) \times \left( \begin{array}{c} F_n(t) \ N_n(t) \ F_s(t) \ N_s(t) \end{array}\right) + \left( \begin{array}{c} -(1-s_n)T_n \ (1-s_n)T_n \ 0\ 0 \end{array}\right)$$
Dit model stelt een situatie voor waarbij de $F_n(t)$-klasse niet enkel jongen bevat, maar ook zogenaamde 'Floaters', volwassen vogels die nog geen broedgebied konden bemachtigen. Jaarlijks zal een deel van deze floaters de broedgebieden innemen die vrijgekomen zijn door het overlijden van broedende volwassen vogels. Merk ook op dat, zodra de maximale densiteit eenmaal bereikt wordt, er geen in gekweekte jongen meer worden toegevoegd aan de populatie.
We implementeren het densiteitsafhankelijke model als volgt:
cap<-100 An2 <- matrix(data=c(sn*(1-m), yn*bn,0,sn),nrow=2,ncol=2, byrow = TRUE) Dsn2 <- matrix(data = c(sn*m*c,0,0,0), nrow = 2, ncol = 2, byrow = TRUE) float <- c(-(1-sn), (1-sn),0,0) G2 <- rbind(cbind(An2,Dsn2),cbind(Dns,As)) G2 <- rbind(cbind(G2,float),c(0,0,0,0,1)) dimnames(G2)<- list(meta_classes,meta_classes) G2
Om dit model te simuleren, kunnen we niet langer gebruik maken van het commando pop.projection uit het popbio-pakket, daarom passen we de onderliggende code van dit commando aan als volgt:
pop.projection.WB <- function(A1,A2,cap,n0,iterations=20){ x <- length(n0) t <- iterations stage <- matrix(numeric(x * t), nrow = x) pop <- numeric(t) change <- numeric(t - 1) TN<-numeric(t) TS<-numeric(t) Iter <- numeric(t) n<-n0 for (i in 1:t) { Tn <- sum(n[2]) TN[i]<-Tn TS[i]<-sum(n[4]) if (Tn<cap) { stage[, i] <- n pop[i] <- sum(n[1:x-1]) if (i > 1) { change[i - 1] <- pop[i]/pop[i - 1] } Iter[i]<-1 n <- A1 %*% n }else{ n[5]<-cap stage[, i] <- n pop[i] <- sum(n[1:x-1]) if (i > 1) { change[i - 1] <- pop[i]/pop[i - 1] } Iter[i]<-2 n <- A2 %*% n } } rownames(stage) <- rownames(A1) colnames(stage) <- 0:(t - 1) w <- stage[, t] pop.proj <- list(lambda = pop[t]/pop[t - 1], stable.stage = w/sum(w), stage.vectors = stage, pop.sizes = pop, pop.changes = change, North = TN, South = TS, MUsage = Iter) pop.proj}
Simulatie geeft ons de volgende resultaten:
pp_dens <- pop.projection.WB(G,G2,100,init_1_met,iterations=100) plot(pp_dens$pop.sizes,col=3) stage.vector.plot(pp_dens$stage.vectors[1:4,], proportions = FALSE) NS_division <- rbind(colSums(pp_dens$stage.vectors[1:2,]),colSums(pp_dens$stage.vectors[3:4,])) stage.vector.plot(NS_division,proportions = FALSE)
De auteurs testen ook wat er zou gebeuren moest een catastrofale gebeurtenis de gehele noordelijke subpopulatie doen uitsterven:
Local_ext <- pp_dens$stage.vectors[,100] Local_ext[5]<-0 Local_ext[1]<-0 Local_ext[2]<-0 pp_dens2 <- pop.projection.WB(G,G2,100,Local_ext,iterations = 100) toPlot <-cbind(pp_dens$stage.vectors[1:4,], pp_dens2$stage.vectors[1:4,]) dimnames(toPlot) <- list(meta_classes[1:4], c(0:199)) stage.vector.plot(toPlot, proportions = FALSE) #NS_division2 <- rbind(colSums(toPlot[1:2,]),colSums(toPlot[3:4,])) NS_division2 <- rbind(toPlot[2,], toPlot[4,]) stage.vector.plot(NS_division2,proportions = FALSE)
De rode curve stelt de noordelijke subpopulatie voor en de gele de zuidelijke. We zien dat de noordelijke populatie na de catastrofe kan gered worden door de zuidelijke populatie.
Tot slot simuleren we nog de verschillen in evolutie van de populatie wanneer er wel of geen migratie is en wat er zou gebeuren wanneer er geen source-sink interactie was tussen de noordelijke en zuidelijke subpopulaties.
### Figure 7(a) An <- matrix(data = c(0,yn*bn,rn,sn), nrow = 2, ncol = 2, byrow = TRUE) As <- matrix(data = c(0,ys*bs,rs,ss), nrow = 2, ncol = 2, byrow = TRUE) Zero <- matrix(0L, nrow=2, ncol=2) G3 <- rbind(cbind(An,Zero),cbind(Zero,As)) G3 <- rbind(cbind(G3,Fledg),c(0,0,0,0,1)) dimnames(G3)<- list(meta_classes,meta_classes) G3 An2 <- matrix(data=c(sn, yn*bn,0,sn),nrow=2,ncol=2, byrow = TRUE) G4 <- rbind(cbind(An2,Zero),cbind(Zero,As)) G4 <- rbind(cbind(G4,float),c(0,0,0,0,1)) dimnames(G4)<- list(meta_classes,meta_classes) G4 pp_dens3 <- pop.projection.WB(G3,G4,100,init_1_met,iterations=100) plot(pp_dens$North, col=1, ylim=c(0,200)) points(pp_dens$South, col=2) points(pp_dens3$North, col=3) points(pp_dens3$South, col=4) legend("top",legend=c("North, migration", "South, migration", "North, no migration", "South, no migration"), col =1:4, lty=5) ### Figure 7(b) plot(pp_dens3$North, col=2, ylim=c(0,300)) points(pp_dens$North+pp_dens$South, col=3) points(pp_dens3$North+pp_dens3$South, col=4) legend("top",legend=c("No sink", "Migration", "No Migration"), col =2:4, lty=5)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.