From 7e5187d7a94267fe51ed9a16f56143e7778aaec3 Mon Sep 17 00:00:00 2001 From: Steffen Date: Fri, 17 Jul 2020 11:09:54 +0100 Subject: [PATCH] upload code and data to supplement manuscript manuscript to be submitted to Journal of Applied Ecology must be accompanied by reproducible code and data --- Balkan_EGVU_PVA_DATA.RData | Bin 0 -> 17192 bytes EGVU_IPM_Balkans.jags | 348 ++++++++++++++++++++ Oppel_etal_JApplEcol_PublishedCode.r | 465 +++++++++++++++++++++++++++ 3 files changed, 813 insertions(+) create mode 100644 Balkan_EGVU_PVA_DATA.RData create mode 100644 EGVU_IPM_Balkans.jags create mode 100644 Oppel_etal_JApplEcol_PublishedCode.r diff --git a/Balkan_EGVU_PVA_DATA.RData b/Balkan_EGVU_PVA_DATA.RData new file mode 100644 index 0000000000000000000000000000000000000000..af3692fd70c7845ebc0abd79a022f8ed514432f1 GIT binary patch literal 17192 zcmeIZXF!wPvNkNDMn#BgCJt z!&~;Fa@h4Z)>;DguHMy^zpyh;dOY;z>YH%)3pYxBJbJLVcQK}^<&{1U zORpA40)ZAYvcgF1qZfWWRV_1$n3bKLt%ZQf&kemvfa`tLoLDpv+H9wLe-V=YI6#Lt z!>;>D?s{rdgD;R?05vr<=JvkGX=q_-Uq_`}yCyhxVw1jg`eCzf@sd0Y60?&!UH&^u9!%Xu1*-fFhm$wiVWbud8UQpV!r7*RG z#;UXDbXYf;&%c6rd1B}G>y7P@DHjj*{1aX66mt7+NP?z5)I(WabK#YTzxb5B*<+;4 zyq`CHnE_3^tyrgDXq|x_TjR$DeEbT2ZR+fz#yM%xd2wFF7BzcEN9|H&l45~YTX2Lh zsA>2BpE{;BZF~@5<(It|udw;{me1~y+@=Nf#<5YXON8yI%j$=nDvPV-mGzrgpQ+$v zx0en-G8eQk(3z?7)}fC*_jNaLT8V0&E72#mt((>5A=W>t;lATreP1GLeOE!gG&o_a zXDhp`5|R?nRs)^i6#we{T0!t?KExRm4+Jyi&3Zm#9*1H2UJR~d!Deoa+DdQiIY-jn zjO?-YT<($H%??JISvvk`Fehw1r859sXg=|-+2Og-!I&H|Kt*e#*-nccvs`yS+*&!x z*aXEn;13@d+xnH_ILMNPoHr}~5ZF+q725h(X^*RJ<4JiqE`Mzw7v*#~pE?vK@5Yzi`6|sp1(6h=uyuo`7bH3eq`|oSTCgRW`cvaK6U~C1Wq92s@#Y zjVMz($9F&8cW)~=9_%91*)<9+^Yj89O&6I)n)Y3)=3=dO(p_{63+1HYt@UKkf3&qY zyEoLhZ{IL#{}N$-=9!idx3)%W{$j;=Uio?38HHZ|%CuEQpmRPW4?6{@C{BI6@#nBm zNc!aMgJ$EUDrT(Og{MnjkA&2I5ZUFAf>S~3e??d=MX9Y5i5 z;{d&)cK60KUd)eTu2GFTb^mM!hw`RGsnBnZ9_KyIZw1q)c^EP=`S?7BKb@!R$p3N1 z`_azOVSY70#snzH^o!5!>DeK9{6^J=(5mP}&KN$Kx^~TUoY!g-Syjw2Z5YV3;9ZnM zh}?BDiO{7%?GyZoKe(5tgR%mBHtqANbcZq3!mwR>p} zE*@Qoy|oe(+{_VaxV2<{SWCq+Voj%#-t226f!&b=KD)~0T5s^M!Mf#QQy(|f7>PVu zibM*GYS}-b)Ad}fgjx*EJ~i7DS??iUfDXHt7;^swh#hskK;pbtt`wu zST7$1x2$f!vlbbxLT5LL{Bxa|rNOwcQmj%3$HVyiLE_PV0#diMojQ-Mye@Mz|AWbC z4(HD_s`@Yog=vk1(og=7b%2v=Vgq(EGn$zJ&M0hYg>QTl>o*Gx8`(Fy=26I~S3Q#Z zbQ5B3{BCfSQS_I7`Gy*us-12o3=zmim2lAz#YCB5oVE0v*aX*9E-9mrAAv2`EtJN$ znzVAM^C8Pro?LXwhjqL(M;b0zc_i)3zebKlLyK3#rah~+z=_S~-=toMIK~~Z+nHu0 zS){8eQYtC>6AJ9?vS~Jd`sHtj9{;!%sK1`HL^G94ZT3CZaf5a?veNDz`bC|c6EZU? zsxQP-2`^&aV5c&L_7%M(SIy!}B1+F_f2= z$<9#mR#W(fsqE@XnD71A$GpNhpE$)QD4XhQ{9Jnl!gQ_z%S~|&V?IqwG4X5rRga~&VbT3_;SOIo*Z=d(&{|X z`B8=a4)2OM>dpZSmQ<_$I?feaQScC-8ws#a!JiEb!c@;zynS7EIFa`|nwnYybax54 z--+XlfvjBBU+>h-H6cv4(;9Ph&(?9k4?h-0maEt5o$6vcc+@8E=@r8Q!b+qyJ6n1w z)Og(_SdK6v__LyqZ+cQ%z^!mB|HGP*xkVb26D6)`&ocC$aN&}}y70G+`;AwansOxz z@^3fn8tmULa<-fC6F`(~+0Nvy?b#7tRuUGtXeCPYJzsuFv5@R}9LmAtrEFihF-JJ? z@!h*sazog&@5l*oN*xu6+E<(!$g)=-A>V6Z4|Pae)Y@K-9yLK-%F8)j4%m`AiDPLh z(oLS};2kzEP@^l3It786h5fzdQsB%O0MFs1^-OYI^Aw%8+$1+c9icJUqbdtIccmU(Q=iFd-TK#^0N?M8>^w)N-Mr3zKdI2!<7lOXS_#87l|*Edp-Z~`nu0wT z5iiByK=8-x*3J)b`umnE>Uv&F9?gq3srOqW)^=sAYU5zrex6>1Z-OFDMlaG>5n2Ut z=5GtsV<_8t5T;^Yu7*GOUnpx(=0Em~Rs(z7q)_!OQ#`W%Y!r$~ki8`Q;UrlX%D=WG zy4kwVh`=>h%)A-3l(iSrqx45T$nZr33na@ERv!(VW*nWeoqB&W0w3}3wv|?kz2T*C zl1hN7UC$YxZS}D<^j@GIRiP_?P$bs5a~G_gexHHp9^$^TP@xREC3Rys!UkqzWjs-@ z3%0tDLsxLvV$?TN4ggZMx{IBtDt^UU<@QZ5w>)^EwL_x?+bg zkmI6}vENh+^L-?2tDhvKkn4hd7nZEbJhGM+Vy&D%!s%Zd$s~B>K2Y>f;E8@&+O&(Y zyZVry+$`|m%x8GVUux)bb3~c#4;E)%TbSBVEB&U`o?($T&&sX3Gx-Hp3N;9?H`^V0|w5c^ngt;OB^uv#X>(Nn!v zk=a1nXWzaPUtFH%9}Ep_#WAJ~BWmX567G)9Z`)I?ZDeeFYxJj{qA$&ey~JHVsXqLY zB`}LNJFIDk?^h{bkl#vqtDU@ERy)wsI>w;sxJjF9Vy2O2yZJuXNlz`ZT#0gT!=jCd znOquKZlsS;<|%_fV6p4#5i(68+$HP7LjhhH8f5}tvov*oPvrNIrqsx4o(lGChi#OA z*I<(7j?eBuY>Q5sU25|{x#!v8$nq-16E(eS8>Xx43%#E@nYGFT9M@jLa`b{71}hW# z*36ot!B?Jo!$K}AGtnFW^ql|I#+ccXZaJ-L;^VI8uoYe?zNCYav2+drTnQVlzOM9x z)mOqWd!C-d8KMn(a1v-#@7#dFhkwOgFtfhL>&rqIFQ|9+i$%2JzuaaB<>aADuV|-1 z|2zmx{>U)kL1ZN`DGX%0t8QOQQwVEy*=>&uMFJlI(zcp!2Mu5+1N^;-mHg{c)rU^=#P z0o9W=@7(n`tTo^r%4I5rBfVjx&e)}4Rqz&+1C$_nRnH7^3!1Qu3vC}N9JDp&k@7Y* z^Nj^7?M8OQg|>Py>oup&<8bgJ$@VySD*~Wi{|?pGPBWzLXMrE}`ghR-<8Adl<6}ya zr^i21_ZsHuBZ{3MJbqfIKk3DX;(%OVz)I@%TM*>w!~BeFdV8R~{_d%9JSSKiJ?tNn zC3@7o(2f*Qv)HabdJsqgP(Od0VG~VzbJW-q(=Mdl;iB#@=Mq^nh6Gy$D!vuoZ4}z` zRnj_5b&kwlt6ytyR(E;EEmrF0dUX7{wr>TSV0mxXIj)j#$TbEa)Uqwz@zkEtZhU&QQQW1b*7od>1Z6EtNnI4NgMeJb`WYa?*I5q?ait)f%DF zH)8e%jbvK^#`%uWt7{TRS0A6O?YsFMm(hPo%Szs3|J0&dR(9aZ?+ThM)G<@3u5d`x z`k87s%Cj55G9>Q#4Sxtq=Q}*yJx}yB4NT5MhE?kdn`~Ciy&fbi9w5H50dB$t0-h2} zR=}x(>60(BZjueY08VN|OQxuD#A&34zo+sy3Df~{R4c#@SbK!--4z!_t8$ry*1uXZ zlXTH~*j+CD&5xJY+G!pa^&wQt2YCO@t9t>C`?@jBwD~kUuX*La)zhS>ft=Xp;}Nba z)8EtgSU-y!U@{P8prmrETw&Z zpv9edwIyEHM{DMtRD^ODZc!^NRS4~9?SDr$`f^kY-r(jjw2;=M=-xcWiSg&_Z7Md= zoAPxWT@qL@!#2ATU5{_uI<0fG{2~A=>aD}haw87->yg5t`Q;SID>16~~0|d@$iLb@6vc$>dsZLD1GEuMgF(XIO(Qp!{ zw{7QXX(1(VfMgl9S`7H=kiN>bbn&*ppq?8!GfQwv@naF*6BzQG-v(Tr7 zINR7Lesf4sC8veTn_7GR!{9;RPR*i`lE^*BPRZ5lik&IHBha}{KOoC$T4QnK0$yGB z=i6r=BW|$z&$>)&+~s(7XsIMcX9cmC@yR!EtJUn9oSZr>b>bo{KFdMB{w!pxf`Sj; zHkz^o)GTP*PKVaozkjr{$gs)3dcE!V>UfpchK$K<_0WBTon)o0y3^J`vn$FlG~;|5DfmIx60Jx42s5sLFqs1J(@kJ0=KybO3l^hWvfNmkrmOZ=P@G4NNG4^;F zI|qq5tp;@9%s5e4_6VwK*WzV#QloH>?P2RNtH+ig?x`fV*7sFYHQy%+Cg^0%ifxHi zS)2mytn%ckJad`)o$7sw`L;X>iZ@<}oX_Hd5FuAv07tIO5=B~pYBf0eCnmnq{p{T8 zd|!yJ$xA{(PRQe8|4Mr4(#xxr+}X-kx9@o?zPXnLav@Hg87V}RvU~>k@gTnh@Ncud z*ej8Agv_5+vZ=67@+d#{N|`CQ$=B-4p9ba+EsBkKJ_`_ZY)?zmvt%cKkJa5wsCD}x zdm#|&lhuDH!yx$~9MxJkRU3wUAN#cKh%MbwiL3JP*hMBO^@m>W{$O`*BOeS<7co__ zS;*^Bh5UjvgQRHqmV>7cvf3G~&XfooZ9%Mf+*W56JwjqqODQ#7^l8rqCIwohtfG1^ zwujAl`N)mW<)rse7rilb|5J9MxWtq49>r01 zReQz#70BhKya9hB{bm+Z+7!mRfbH9tfb$6z$o{eIL>L1AI`l)U%6Ky4EU&>4`W36J zTg}Nz;fj`#G_6HAj~3P`&dOp+KA}Q&GR>V%QbBroq5{N^`vJxsv)P3yZar&LlbMnp z9+CFBSb>hh0xRDc>dRgnnyUf@J>ksSixJb~g*=pn-xk#RuDMcvzWd~;(Nn6p5!@|; z8nxQs?oQw4WSIHR9>7aefqmY4aQ|Y`_T{dv(0#Po=euW6v)bV>sPdNYQD3bwjV!lMm6>roz)@5p{R7 zPh9us$C;P%t1BmzG;gfHQI5xlap5I2$$Po=;F!Q$ny<#lSAJ-E~%_j1n15E4`)` zapNZGGp`>8=c%*kA&HXszyNJ;?Z-_whSCKVynH238*0XcSDz5MbMr#I`9!@Y^=ybT zv4e|l`L_5LFQa#S&7=N4>D|_ALTAorTd>cQ8lmu2ox7l?cAC*Gfc^5-3CfbOY_{!E z>)D*UvbCb?AkR++*QYJ>UDqJ1rrs}Wt^#1=EeOP3(~0vHnoSQ`?Lw3xN5J=mWWClK z+CPyj6jaxX8cW(I$Z0}D$$TTD>{0dF?VotH`>oODI~m}_ZGH5C)>u;`-|4{wAzV2O z-x;F@>hKLy2Iu!$*4rZ$^dGHbw`}Df6v9@B33U|#8^|6ro&{wYO z?riTBo5XyHN{>Z29D}sNOW<@g61(?2O523zO+Rg!qKp+7GQv%!T|3&eaADOJkvzJu zZpn2eP%@;XPsWA&J~Xr!DB-2CZ|G*r-nr7vQ}g5`{`H_VuJ-U<|ON8XsU#OA%-c%al}u$S)5>6yx7c4{4?CGnsT zFmQaYofhD2-CD1COjJY^`LkuVR7u9Y+`gCB-?P&1J)f58=zGgLjMu#GtoKuq z{FBx?e9x0!j}z+Cn*zpw5BDfi4V#$hV&|1!a} zbAeMJ;hl)q^*Yv^aPkbW`u?Ty+rFMF?BYDOLSM5pR2<#`M(13T$P!}vJcfdubO?{} z(({KK5hg&LKDDc;%K=|Bf^!HFr7yB)^%d8Qqp=NYrUXgtEnjhCX%}`RT4O3 z@$BY8!Z7D<54EN8-Jcq_H0NY>I6SKl@@$^HkQL(>YE3qSeNEV`?4G33ON2Vmoz$`i zPA9cL0#>uNLj1R7c5``M(-Lug7Qe1ggx=v09-Qu|o#HfUeGvLuCr12c@umAyvCYqA z?(Y+PQl&jlwFGY%6gR1Dckp$>H~n9=IQ1)_6jT-IUBb-JHb3ic@Bi`(Fzda$Ui!0K zhgpUe;>Mq;({-_v(e_F9S2y;sSX(^qSNvnQc|!tE}MG+3}Y0wwXR2 zF&V#>ssdyfHlkqBN8(5nh)B~dL0gPn&XN}rqq|rSq3B3H{9bY6({RO3`{*~(i-EQ4sGJSmm?cZFkXO1>ZcS$v-87Xrou9(XYRS|eu5&Nu zSJIxpfD5V0uYP$0b$Zu-AFADS3Tv-Blt0S6yLh&A_mkdWrT4_+!zAzg-ZRRpj@pl0 zA{s%GTp5@k2ztCoC}F=|Gd&z-ba485M4{lvU>%&Q-ZO5bx+PByGWf&jvyz8iH_<=i zl0Wd)W|77IMrEbnh%PYEMHDIt^0z8PPryy)T<$J*epU!>S_)m@$68LmewEbf%R#g0 zk5fUlU(astsCIID25-24XWQ{Jvxr)3cmz} zyYn;>yNm6$PlDzu#oaKt)Qn2EXtK1E>68s|*EAC!B2OTPr)gZ9Rv}2+=04)s0`0Bw zwEg@=uWreljiDL7R;&(i6;rzrho+>;zlnhMd28aX#kp61pRh|GD>(jC+_$HG{Ll_= zIWphE9Hv!zbj?_IgI@qX+}Y8)!l~Bjl(XTB0gG1*XmYd6T%Ka#3fO)CY&0x!osc+hktk5hh@D&R5O4?HpmRZRX_?KQ&DSNJx=H!qJiUFNt{Zf3VNXkH7dlpOhrk3+~4O1?1Rg#69$wx2Zj6J?l z%4d-xV}?77^#*5$9W(gSfx&gJyp+w0^RQfl5qNWvdRP=Sw)1VvOK(w~Oz)WGf>MF~@#edGmqBJi9MYcQ{y_zkNN*T`0s|aomos z5GeR@$t7d{NQym*Yu)epZr4d?yyt#rv2TvL3qdI|3}hA-=GKSC^7_uH>PWWOcl-*% z7xe8nsrc6QzPSpl#K4ov8`H~SyPJ5Mb2O%#!xYntBHAD70uNvNN{|u3aIi1b<&J@~ zs?S^l8F@n)#VW{5oAAAGV7%o}iw?-^Aa_4~3S^?NeEFbx7~0(4>ZI*=%&r585Z5q5 zz1k_F`SPLBFsI|RE5%ns16VfTq%QU`yo6<0@ct6S!izi0x4Z%gI$`ODlw=ExzNH5F+~_eGG=)iagftpW5i@Q49GIKzDakLyai+3> z5*%2sy|>E|055KPnY1>E%pq&$%`%LUFxNjKtCxrW$VzpbXc|I`vgT4@=Wk?z`%5uO*Xv>4L*8G;5U8Dd+ZH zTur}UXHfga;LR4c&z#fyoONY)Y1{+RKiE4`)Z085uLQ~;XdolfVkOvUvm9HU#M~N& z9&ns|TZto5&#*Zf!?yJz5L|@?ul>>#($MlDs`Z2{?d;QMKexi~AM>&TG>RqSDPrU| zv&;>Dx+nc|OKz0d+um@AkTu6@SmzbdLr4KAdUq1!(FBN+UTL<74!h%0py_nz!}DarINjLQ*Zcnnp!L zox6iXo>afc_N%GLYe`E_7CmHe9~-WOo*et<4H_f1+!hxT9!L=IUKazbPL0#MhJKzh zrB7Cc>;-@}Ow;9zB_S<6u%jaI`vtM<<#9dfk`ET2>Ea}pe@?$GdwIC6b*YtmRgU&% z!$9~eOF3~*4z`SYNou-ge8iW2rICxTr)20Sjl4c>93vgUt@l#b28ciXN9q$^4s|lL z(q~6}T5r|m6<*cXtL$yNyG`nOKmGiq73epelC$V0UpSUh{54P!9)RLY;7frC*$@fV zlud|RDm=YL`0l{zQ>~o(l{8$-x?0R&aNP;b8Q+Hz=Ou%M(OoXUwabsZ9D#!_`!_aBYv-T`Jz$ z4|kdbeN6IoyJTJ}^sK~9n$c$u_O4oQENXcy5((}$yyIGF4+aSC557nJ_4Xfzxs@>q zkz9WkV!Z*OBnSwWWg$beaQ#|D+TIgNO?8hS+v94IV+sNBwV~ZVLUXNBpQav+jV#9F;&I|_$WqwWZlwK# z^oqkITw}<3zsJmKN=8GOk=BBnt=~-1{d^?1-r1LKq{E&Iq3g7++{iz${KN>!jr1Ux zAskga{E20_Jv1ue7jE(k{wva95`WOBb9`Eci;n<~i_;f1?8pvuV7q8~wEgJADTWvs z5|mSp8-)Z-RgqKbzPm9$lZN6kLrgJ+QJPmLq;|7zvg=OYNOl`Ppwnee(@~(C7}0dI zQ_)eNo|r2VzgHxe6P9tvRjfberNJzuqaZY(A>$US$($?HC~rjtAymBy5e%(i^Q2bk z?>{&{p`qq6KEhO!X{=ln#c$yJF$wWlSgp#)CcT6E$GH8Mrv1xaN5Bzb{5}I_qi?s= z@Yu)4$J@>)L?UvVS)^H0!3#LVV!@wQO`PY%>zO`H5*MqCej+kvb92P02+CL)mWRH* zlaF12GNa=ZA{Y|ho)`|FF{nuUqgq+i7&G>eo$O2MW-)LK%!y)t^la9ey8H2su*bgp z>9NM4^*a{on%a0zdHKTO0742tZpb?2)A6$tzZ_P0oZ5@9iv8A~d6B0xc(aj2KtDN# z4F7s`cQXu~*0VYJPX_G3X-E2jN{W}`QrkGp0`v&@_sCbs$%c1gXmb~NmLOcxp~q5f zr~}5BOcxH2!#@%NI)mQrJYORO!M;LncaZmXK%W1*@;7?m7v}G9dnPEaH&bqd_b?q& zAGB)%)P|ta-AAJU8SEddxqt;YD*Q)+D24?G>kqD#!jOfNFA?&0UVjIW*oWqz%Gw@? z4IJwMvVru=P2v*)4#ehe@=`va!LbD>Q&inF5pe3PG7#K)kQlm#ZzIWCfT~>{0S%@F zwqgH}@xcVahDl}$gfKEBMY~9tx_y7Jp#^B10iQ_r^`u-B(&4xs0O+R$n+hUM7BaO_On<&NIr z!@VRJlYyeZ>TbxL`OFhBXFlK*$&+CllJi7mh;w&P4@)F$n1uA#d%%iRB%o#C>jK2k zeHLwE;PGJyb2_siR30*j*7{og_L*{Lb>}nX_UfEx%H7qA$+|)5essf&!FlHdLoy$p zSK+f!lGeJK%s+g==KO26bAKeVuh}lx{E6Is&34}APsIE++a;Snk(;mC$ZY;ZG+(n_ zwD}X^f6Yd2^C#l^nhjv{pCo^TOIvD&r$jhog7ynaR@bQ0F2ca!N-?)(jpe`EcwENN ze|Gf#>!Z$$*|?XI??9S+=Bl$6Gpmk*?SzJjN4>qEq@JOGl%asBiGXhMM1BwdLhVE~ zK-R^ksbBw4pKsZ?8HU(cLqxf%^Xz(Cy#lMGMD}f1gUnTadX^lD7j>9<9N4`=*}7I} zGY=L8mp9HxYnu(!i zkk}5$$$ZvL%p4t1=eJ;&9A$q8FUW356?8e*7SD!TR0S!fhrtRGlwprtB&Ef% znjj{X{$ORErM6DOx!$2%NMb$c&c%jU05y_4=TX9rdT%G>z3W3DTNyynpNb+}>+>T&l#tcqSufR1mGhb9xUwTFGDw!O9Gz?c1Ig}I{ zbnN55ZY#m#+f7vj_)seH?H%eM6(-8?EjWE|Ym7V9ljl0{5+;|rE?O9B4Z5^}Bh*-@ z^Ur^=Ud#0P`okyqRmmT%u>H1jyS^NWOs`$2{HO4}>TNN6!OJ^9AJmwq@!kh!ZQyMq zr;j;*o6cXaU!#{O06nUmz4J7L_MbEUJ3Ga0!431`S1YN9{!ec#aFmsWvJ@EtK=OUU z(SI{N@O#GLHQbflTIQdj>vJxs^vP;wCa^7gsP}Ny?CS+G&@BT@(4tILIXg_9$#gQH zVDqAl?-OLOAOtwwAba35$>Y&&~FbUeg^YHod#vnohV> z6mTbnv$CJOH#nku@J9GTVhwoJmna{@P={ya9ly(fB;YMLy;HOYND(A>?Bh+0 zM8UIL>0a+|KQn#+9WCQ~TW;MFR0RG;8_tE7-pQOob{sp62QU#@&NaiW=ViB}Jyq`b zN#ZEAPo-=7#ew&L$^Mef#3$=L8fc2I$KGr_x4^av(!T*090k29{Ko&ki2UW#Z;`nF z1&XQ`Fw^1XuTE65tbHpcOP5n7y|^itsOg2sK?e zyxRWK9L~TglQfgczXzBKN%lq3nD;Bm=Yd z?KGoXkXJ~ic5H!&dc8PAnv|nVDQxIIzCipli#bynJA*7kdg!_!;oMbH_&*0c_kO)2 zy(1ksxOu+Ir56P(BY*B&Z~A(Hb=^?DClERr3c+T~DL@W?|9)e)u-Ej&yiJu zqz}bi5&yl=@;ics|4QyT?_pB@|8{oZ#)7lt2-EMp0PtV`(+1hZ-kl;r^>`9s&POy&&*v%|D{0&U|jVgK(Zy!~Yuy?!Q3X{Qm|bMxChuIF0)dgK|my zz0#Tc$^TC3K(H*e0@)l#=za~m{IR>uAd@(UK5#`=Gu#s0r(nLNAz zEFhtj8-%G#*N1bRlGbkZ{PIrbH^UCkrRw8U(w~!1Cj7ZgtcJ&ha6oJpE(?!-4}e50a-VOzDjM*T`V- z`F5R`E1?*fEEP;Eeu(dOmvI5`A7tKv%u^= zE+7gfeflq+W`P;^Z(b&jPVVsSN&jH*wP+Cqs1IIUIa}=DUhLFyQW`*PygX#N0~yu| z875Gi5WSvDnC=BwVcTd{kutk;IC(SB@=-2rXp6^hZWhaUdFbELtHjM=6rZnD6@K{g zss!uNmd9saAP+<25mBPr^ZV@38bt8j5zl+5U)pjXxd|fkj7U|gY`JOe2}+jshkgc{ zu^;V7BYAn4_LeI51tfqM#_8@Hsgv-Lm<2?H98T#kgg1kRqT#c@;r*VFVKzfL0E3l@ z7J~kEg(NDN@%}V%@a)1~@BzOVP$jGN_DCpRrj)g0&1O%J>>^-7N1f9X)7) z$6Fy2#IWDpJqz(;>L0-QzD3B_A(b`{Iu7u-c)MZVh+xLn;V5aI-&*2?_amfr+>QMX zsEPWLhoehY>kN2i9eDn@&2Rey#7V#7U#iJvW;W0p$9!q0SJnV{ONY;RYj5?Jk)%s< zcFZqS{vGI{tKd%a3rZ=Gi%63KKXu$?`^JE*!?!~TlUu+2IDCVGMPn)L4;miU zZ^$kM8!G6>FsR_e{y7?954L{<3@_Pk*pVsCrUE7(2xb(V-for(I(^<+zgO|j?q8aU zcyHQ6So3))*s2$_7=@T+mo$}hj21I}CXjx(Vc%osw(4b$-F*1mP&V(rxuJR9{TGIs zdG{?0f0ZuzZB4pDz3s)o%{8gOS;Ikq6E`-vE*{@z)OYyLVQK1ff@FUt=LK1$ssBXC z1aC-F|A~+b@<~(wfyA5>ydwQaa$fMEbj+U!nIOG%%%2FkpqO;bf09V6tfLjr7d3|^ zwXY{=D^YiKjr#4B^&K7-b6eCH)ceu$WQ6ot&+Yurk7i{JTH-HfpRN=o)7Hap#-_EN zb%$rc+^}(GT?#LrMp3taDqcc2I#xVydvlS3mP`N3^R^xAl`_Ei8<;u;MuwYoKqQ?+ z|K~$rU>xbC0L}~lyv)y5!XFXJPX~O5dBNl7Hie}non9~EGigc3!AC&5|MtKq1VrE( zxZHaBR`B*1S_}M}{cYo+>7-Sds1}kvZ@x|Y46`>T0>ZFZZ}R~$+ah(O%axCSVi9oI zBRHbrbT={-_RYYrKNxm0u%7QkI=Tv^ID@5~NB<^z}cphh$b~Z4Fehk$=v*syHuJU?3XyOKiKw!5QiQh{|=rInKVUhgi6KiIRG2nGJ5M0F^AC$jAD9K})Ft1c1!+vdsJ{a&ULXUSR!3W-29vay@| z+&9K<3UFT;`^3*pKX%@UlT?p!__pcqhrb+1u4~k8r)=OrUp(;F_1ynF*OM;o>g>}I()jU?$?*Kv zlu?HXlQ?Xqc+s|?`6}9w$TBkDtGxqWl!(dqdr4h@zO+iQ>=@LK#_%M&Sl#LavwEAH zJ!`}5Sbs0~kkPRho6!G2TUA0~YMp&wuS6MSj5Wi_+IcqTjs6e5HKB1_-RoIy#cggY z1G?tHC(32|36T21Rgc2ljYK62KEwV>9`Cfsc;gBlm~+|$PE%oqFRt6XI?>W(OPrCr z!Q5o$mv}KUH@(qBx!fz2u<%;41* z&8p{?R&7fCZ{iXza&lAX>e80vaOsN~e9^3YZuzoJ>D!w$_*+ZX7gLJ*!kY1sv6ifs zQ{VK3`Qsx$Sh8AAeKQc=d6r_PcbgZik{p`Ql#zjP%?3Dx*2CSrK&vc3j0j_+RW{EFQr6<{J;qa}E?22f;h5Qjz*1&G5{%iEu z70cYrNC~h9du+|DnE?tRoYgNFaU}C#sYNF~`fcTkpY1Wu_p;@*8+li!w(JyPDfOy# z;|qt>OE`0}*6N+&RCQGA%{;%`6qEVhKfe^0gKm7UOtVT=>wD*iS+BH9{LsUE&u0uz*HFu@SF_fBNfn)ij3SPkiSCq>A*O(2uVm|F zq!cT<`7@m8*xiKvSP&r`=~CH%a#}iClsA*gXCmw#Sm3VbY_2iMBt6@{SHQ@Kz+vif&JqpdCZeWYW=!`%}v#spHCiS)ED0MJea|r=G%Ya z*%cpd9)pZEpV1#;XdDhU>3q@zk^YT3LKDnu;JcD;1_BgaK@dMVmWx=ba=7MH!I-V3sYvtNQdIFDPpn%-2scMsSN(xOlQ_b5^h zT`LcK(h7SV?2#o!{2YiFJLWw)qltf5mw5c#Cq#gwFZ@orTHc)gU4Q?)S6@-0r0AcW zM3jM72w9tEcz#Q$K6}wodtJ*zXHsN~zwxq+;wTMW*4_NSI0hq1ZLiPE*hFcj#kFmOeg(*ZQ>Um}nhY3yBRfk*Ib&`kAU#=? zWhhH73x>;L^*im3q`!7FSVWarmvlkwR_gK^UpLy9X=c4_qGk<4I$Q;BR9kqg1tcA+ zy6FUSG(;E`K7!F1mt8oWaz=^S!CHU%3W69zxF96ut>0|0tuxHRTvD_k2BB7o*3 z<*lph>+rXMKrh4%i{?}JD)lanw7y5yG=gqX9MU%ZwaRsa>xROcv$byQ1By4KfND-y zvxdVW^EXQsL?J%&YIs!!V$<96*|*{3iWg<{7$r`@u(*f3h}Cxn`uwLz(#Fd{AIq(> zGrrS=KbT>_zb184mhEK2NWcDWr_3P+gA4K9|6-8(IIx2m(r$v(ocTuHD*%^F1N?r! zign=Pebm|t`>lyN5o`bN8s9h}%P9a35EEcf4h_2nrQ?CXe>XOkNE8~_FSke1t6p!| z;iDs+oJqTt1!vR5(2^bvQd6jMh=}MWw*%obNiB$PP$H27_~5t3p#^BSBZLHGs%Int z;7{5NPGN&#!=O-nJLF{skhDqpA5x?ajXetVkCu;bB){{KjjbQr~2*f9OrL}SSI6qZN^p@G_@|GoMe5O@_>P@gv+)JNK%3W5Hu z@$#R27=OIn)YvJheIhrcfEIn=+eqrskT&E9QLzD}F=SDY>N{!AmiL_w|J?WtgONew z1$9~z0Si0Ca{R|J(zdo_VVM!$m(Ceg0qZ#H?mjNxdS0T77;QL*4|kr>Aa!E&;6^W_ x=*B(}e-)b1s3Lg8|FDJoG`11 + + ### probabilistic formulation + nestlings.f[ncr,is,fut] <- (fut.fec[fut,ncr] * 0.5 * Nterr.f[ncr,is,fut]) ### number of local recruits calculated as fecundity times number of territorial pairs + N1.f[ncr,is,fut] ~ dbin(min((imp.surv[fut,is]*ann.phi.juv.telemetry),1),round(nestlings.f[ncr,is,fut-1])) ### number of 1-year old survivors + N2wild.f[ncr,is,fut] ~ dbin(min((imp.surv[fut,is]*ann.phi.sec.telemetry),1),round(N1.f[ncr,is,fut-1])) ### number of 2-year old wild survivors + N2released.f[ncr,is,fut] ~ dbin(min((imp.surv[fut,is]*ann.phi.capt.rel.first.year),1),round(capt.release[fut,ncr])) ### number of 2-year old captive-released survivors + N2.f[ncr,is,fut] <- N2wild.f[ncr,is,fut] + N2released.f[ncr,is,fut] ### sum of the N2 cohort derived from wild and captive-bred juveniles + N3.f[ncr,is,fut] ~ dbin(min((imp.surv[fut,is]*ann.phi.third.telemetry),1),round(N2.f[ncr,is,fut-1])) ### number of 3-year old survivors + N4.f[ncr,is,fut] ~ dbin(fut.survival[ncr,is,fut],round(N3.f[ncr,is,fut-1])) ### number of 4-year old survivors + N5.f[ncr,is,fut] ~ dbin(fut.survival[ncr,is,fut],round(N4.f[ncr,is,fut-1])) ### number of 5-year old survivors + N6.f[ncr,is,fut] ~ dbin(fut.survival[ncr,is,fut],round((N5.f[ncr,is,fut-1]+N6.f[ncr,is,fut-1]))) ### number of 6-year or older (adult) birds + Nterr.f[ncr,is,fut] <- round((N4.f[ncr,is,fut] * 0.024) + (N5.f[ncr,is,fut] * 0.124) + (N6.f[ncr,is,fut])) + + + } # fut + + for (fut2 in 1:(PROJECTION-1)){ + lambda.t.f[ncr,is,fut2] <- Nterr.f[ncr,is,fut2+1] / max(Nterr.f[ncr,is,fut2],1) + loglambda.t.f[ncr,is,fut2]<-log(lambda.t.f[ncr,is,fut2])## for calculating geometric mean of overall population growth rate + } # fut2 + + + #### FUTURE POPULATION GROWTH RATE ######### + fut.lambda[ncr,is]<-exp((1/(PROJECTION-1))*sum(loglambda.t.f[ncr,is,1:(PROJECTION-1)])) # Geometric mean + + } # end scenario of imp.surv + + } # end scenario of n capt released + + } # END of the JAGS model + diff --git a/Oppel_etal_JApplEcol_PublishedCode.r b/Oppel_etal_JApplEcol_PublishedCode.r new file mode 100644 index 0000000..07e842e --- /dev/null +++ b/Oppel_etal_JApplEcol_PublishedCode.r @@ -0,0 +1,465 @@ +########################################################################## +# +# EGYPTIAN VULTURE POPULATION VIABILITY ANALYSIS ON BALKAN PENINSULA +# +########################################################################## +# original model based on Lieury et al. 2015: Relative contribution of local demography and immigration in the recovery of a geographically-isolated population of the endangered Egyptian vulture. Biological Conservation 191: 349-356. + +library(jagsUI) +library(tidyverse) +load("https:\\github.com\\steffenoppel\\vultures\\Balkan_EGVU_PVA_DATA.RData") +# load("C:\\STEFFEN\\MANUSCRIPTS\\in_prep\\EGVU_papers\\PVA_CaptiveRelease\\EGVU_IPM2020_output_v4_FINAL.RData") +# rm(list=setdiff(ls(), c("z.telemetry","INPUT"))) +# head(INPUT) +# INPUT$nrep.terrvis<-NULL +# INPUT$rand.phi.offset<-NULL +# INPUT$J.fec.red<-NULL +# save.image("Balkan_EGVU_PVA_DATA.RData") + + + +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# EXPLANATION OF INPUT DATA +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +### the object 'INPUT' is a list with the following elements: +str(INPUT) +# y.terrvis: matrix of territory occupancy for 87 territories (rows) over 14 years (columns) with maximum number of adult birds observed per territory and year +# z.terrvis: matrix of initial values for territory occupancy for 87 territories (rows) over 14 years (columns) to initiate the estimation +# eff.terrvis: matrix of territory monitoring effort for 87 territories (rows) over 14 years (columns) scaled to the range -2 to +2 +# f.obs: vector with one value per monitored territory specifying the first year of monitoring (1-13) +# nsite.terrvis: number of observed territories for territory monitoring (87) +# nprim.terrvis: number of observed primary occasions for territory monitoring (14) +# phase: vector of indicator (1 or 2) whether survival in a particular year was 'good' or 'bad' +# y.count: matrix with the number of adult breeding birds counted in each country (column, 4) in each year (row, 14) +# T.count: number of years with count data (14) +# countries: number of countries with count data (4) +# country.prop: vector of length 'countries' with proportion of the Balkan population in each country +# R.fec: vector of the number of monitored breeding pairs in each year for calculation of fecundity +# J.fec: vector of the number of fledglings produced by monitored breeding pairs in each year for calculation of fecundity +# y.telemetry: matrix of observed states (1-5) for 33 tracked juvenile birds (rows) over 118 months (columns): +# Observation states are as follows: +# 1 Tag ok, bird moving +# 2 Tag ok, bird not moving (dead) +# 3 Tag failed, bird observed alive +# 4 Dead bird recovered +# 5 No signal (=not seen) +# nind.telemetry: number of individual juveniles tracked with satellite telemetry (33) +# n.occasions.telemetry: number of months over which individuals were tracked (118) +# f.telemetry: vector with one value per tracked individual specifying the first month of tracking +# l.telemetry: vector with one value per tracked individual specifying the last month of tracking +# age.telemetry: matrix of individual age per month for 33 tracked juvenile birds (rows) over 119 months (columns) +# capt.telemetry: vector with one value per tracked individual specifying whether the bird was captive-bred (1) or wild (0) +# mig.telemetry: matrix specifying whether an individual migrated (1) or not (0) in a given month for 33 tracked juvenile birds (rows) over 119 months (columns) +# migprog.age: vector specifying at which age first autumn migration occurs in wild juveniles +# captprog.age: vector specifying at which age first autumn migration occurs in captive-released juveniles +# PROJECTION: number of years over which population is projected into the future (30) +# scen.capt.release: number of scenarios of captive-released birds (48) +# scen.imp.surv: number of scenarios of survival improvement in the wild (18) +# capt.release: matrix of the number of captive birds released every year (rows, 50 - only 30 are used if PROJECTION=30) for each scenario of captive release strategies (columns, 48) +# imp.surv: matrix of survival improvement (multiplier for estimated survival) for every year (rows, 30 - only 30 are used if PROJECTION=30) for each scenario of survival improvement (columns, 18) + + + + +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# SPECIFY INTEGRATED POPULATION MODEL FOR JAGS +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +sink("EGVU_IPM_Balkans.jags") +cat(" + model { + #------------------------------------------------- + # integrated population model for the Balkan Egyptian Vulture population + # - age structured model with 6 age classes + # - age-specific probability to recruit at ages 4 or 5 + # - age-specific survival derived from tracking data for the first 3 years + # - adult survival based on territory occupancy + # - pre breeding census, female-based assuming equal sex ratio & survival + # - FUTURE PROJECTION COMBINES SCENARIOS: baseline scenario (=do nothing, no chick removal) + # - productivity supplemented by captive-bred juveniles (0-15 per year, for 10,20, or 30 years) + # - improvement of survival 0-10% with time lag for full improvement to occur after 10 years + #------------------------------------------------- + + + #------------------------------------------------- + # 1. PRIORS FOR ALL DATA SETS + #------------------------------------------------- + + # Priors and constraints FOR FECUNDITY + mu.fec ~ dunif(0,2) # Priors on fecundity can range from 0- 2 chicks per pair (uninformative) + + # Priors and constraints FOR POPULATION COUNTS OBSERVATION + for (s in 1:countries){ ### start loop over every country + sigma.obs.count[s] ~ dunif(0.1,100) #Prior for SD of observation process (variation in detectability) + tau.obs.count[s]<-pow(sigma.obs.count[s],-2) + } + + # Priors and constraints FOR JUVENILE SURVIVAL FROM TELEMETRY + #### MONTHLY SURVIVAL PROBABILITY + for (i in 1:nind.telemetry){ + for (t in f.telemetry[i]:(n.occasions.telemetry)){ + logit(phi.telemetry[i,t]) <- lp.mean.telemetry + ### intercept for mean survival + b.phi.capt*(capt.telemetry[i]) + ### survival dependent on captive-release (captive-raised or other) + b.phi.mig*(mig.telemetry[i,t]) + ### survival dependent on first migration across the sea + b.phi.age*(age.telemetry[i,t]) ### survival dependent on age (juvenile or other) + } #t + } #i + + #### BASELINE FOR SURVIVAL PROBABILITY (wild adult stationary from east) + mean.phi.telemetry ~ dunif(0.9, 1) # uninformative prior for all MONTHLY survival probabilities + lp.mean.telemetry <- log(mean.phi.telemetry/(1 - mean.phi.telemetry)) # logit transformed survival intercept + + #### SLOPE PARAMETERS FOR SURVIVAL PROBABILITY + b.phi.age ~ dnorm(0, 0.001) # Prior for changes in survival probability with age on logit scale + b.phi.capt ~ dnorm(0, 0.001) # Prior for captive release on survival probability on logit scale + b.phi.mig ~ dunif(-2,0) # Prior for COST OF MIGRATION on survival probability on logit scale + + + #### TAG FAILURE AND LOSS PROBABILITY + for (i in 1:nind.telemetry){ + for (t in f.telemetry[i]:(n.occasions.telemetry)){ + logit(p.obs.telemetry[i,t]) <- base.obs.telemetry + logit(tag.fail.telemetry[i,t]) <- base.fail.telemetry + logit(p.found.dead.telemetry[i,t]) <- base.recover.telemetry + } #t + } #i + + + ##### SLOPE PARAMETERS FOR OBSERVATION PROBABILITY + base.obs.telemetry ~ dnorm(0, 0.001) # Prior for intercept of observation probability on logit scale + base.fail.telemetry ~ dnorm(0, 0.001) # Prior for intercept of tag failure probability on logit scale + base.recover.telemetry ~ dnorm(0, 0.001) # Prior for intercept of recovery probability on logit scale + + + # Priors and constraints FOR ADULT SURVIVAL FROM TERRITORY MONITORING + ## Priors for detection probability + lmu.p.terrvis <- log(mean.p.terrvis/(1 - mean.p.terrvis)) # logit transformed survival intercept + mean.p.terrvis ~ dunif(0.5, 1) + + # Det prob relationship with observation effort + beta.obs.eff.terrvis ~ dnorm(0, 0.0001) + + # Priors for survival + for (nypterr in 1:2){ ## only 2 survival 'phases' - good or bad years + mean.phi.terrvis[nypterr] ~ dunif(0.5, 1) # informative prior for annual survival probability + mean.rec.terrvis[nypterr] ~ dunif(0, 1) # informative prior for annual recruitment probability after territory abandonment + } + + + ### RANDOM OBSERVATION AND SURVIVAL EFFECT FOR EACH TERRITORY AND YEAR + + for (nyRpterr in 1:nprim.terrvis){ + rand.obs.terrvis[nyRpterr] ~ dnorm(0, tau.obs.terrvis) + } + + tau.obs.terrvis <- 1 / (sd.obs.terrvis * sd.obs.terrvis) + sd.obs.terrvis ~ dunif(0, 3) + + + # Priors for population process and immigration + + # Initial population sizes for first year of monitoring + nestlings[1] ~ dunif(50, 100) ##changed from JUV + N1[1] ~ dunif(10, 35) + N2[1] ~ dunif(5, 20) + N3[1] ~ dunif(3, 15) + N4[1] ~ dunif(0, 10) + N5[1] ~ dunif(0, 5) + N6[1] ~ dunif(200, 250) + + + #------------------------------------------------- + # 2. LIKELIHOODS AND ECOLOGICAL STATE MODEL + #------------------------------------------------- + + # ------------------------------------------------- + # 2.1. System process: female based matrix model + # ------------------------------------------------- + + for (tt in 2:T.count){ + + nestlings[tt] <- mu.fec * 0.5 * Nterr[tt] ### number of local recruits + N1[tt] ~ dbin(ann.phi.juv.telemetry, round(nestlings[tt-1])) ### number of 1-year old survivors - add CAPT.ADD in here + N2[tt] ~ dbin(ann.phi.sec.telemetry, round(N1[tt-1])) ### number of 2-year old survivors + N3[tt] ~ dbin(ann.phi.third.telemetry, round(N2[tt-1])) ### number of 3-year old survivors + N4[tt] ~ dbin(mean.phi.terrvis[phase[tt]], round(N3[tt-1])) ### number of 4-year old survivors + N5[tt] ~ dbin(mean.phi.terrvis[phase[tt]], round(N4[tt-1])) ### number of 5-year old survivors + N6[tt] ~ dbin(mean.phi.terrvis[phase[tt]], round((N5[tt-1]+N6[tt-1]))) ### number of 6-year or older (adult) birds + + } # tt + + + + # ------------------------------------------------- + # 2.2. Observation process for population counts: state-space model of annual counts + # ------------------------------------------------- + + for (tlc in 1:T.count){ + Nterr[tlc] <- N4[tlc] * 0.024 + N5[tlc] * 0.124 + N6[tlc] ### number of observable territorial birds + for (s in 1:countries){ ### start loop over every country + Nterr.country[tlc,s] ~ dbin(country.prop[s],round(Nterr[tlc])) + y.count[tlc,s] ~ dnorm(Nterr.country[tlc,s], tau.obs.count[s]) # Distribution for random error in observed numbers (counts) + } + + + # ------------------------------------------------- + # 2.3. Likelihood for fecundity: Poisson regression from the number of surveyed broods + # ------------------------------------------------- + + J.fec[tlc] ~ dpois(rho.fec[tlc]) + rho.fec[tlc] <- R.fec[tlc]*mu.fec + } # close loop over every year in which we have count and fecundity data + + + + # ------------------------------------------------- + # 2.4. Likelihood for juvenile survival from telemetry + # ------------------------------------------------- + + # ------------------------------------------------- + # Define state-transition and observation matrices + # ------------------------------------------------- + + for (i in 1:nind.telemetry){ + + for (t in f.telemetry[i]:(n.occasions.telemetry-1)){ + + # Define probabilities of state S(t+1) [last dim] given S(t) [first dim] + + ps[1,i,t,1]<-1 ## dead birds stay dead + ps[1,i,t,2]<-0 + ps[1,i,t,3]<-0 + + ps[2,i,t,1]<-(1-phi.telemetry[i,t]) + ps[2,i,t,2]<-phi.telemetry[i,t] * (1-tag.fail.telemetry[i,t]) + ps[2,i,t,3]<-phi.telemetry[i,t] * tag.fail.telemetry[i,t] + + ps[3,i,t,1]<-(1-phi.telemetry[i,t]) + ps[3,i,t,2]<-0 + ps[3,i,t,3]<-phi.telemetry[i,t] + + # Define probabilities of O(t) [last dim] given S(t) [first dim] + + po[1,i,t,1]<-0 + po[1,i,t,2]<-p.obs.telemetry[i,t] * (1-tag.fail.telemetry[i,t]) * (1-p.found.dead.telemetry[i,t]) + po[1,i,t,3]<-0 + po[1,i,t,4]<-p.found.dead.telemetry[i,t] + po[1,i,t,5]<-(1-p.obs.telemetry[i,t]) * tag.fail.telemetry[i,t] * (1-p.found.dead.telemetry[i,t]) + + po[2,i,t,1]<-p.obs.telemetry[i,t] * (1-tag.fail.telemetry[i,t]) + po[2,i,t,2]<-0 + po[2,i,t,3]<-0 + po[2,i,t,4]<-0 + po[2,i,t,5]<-(1-p.obs.telemetry[i,t]) * tag.fail.telemetry[i,t] + + po[3,i,t,1]<-0 + po[3,i,t,2]<-0 + po[3,i,t,3]<-0 + po[3,i,t,4]<-0 + po[3,i,t,5]<-1 + + } #t + } #i + + # Likelihood + for (i in 1:nind.telemetry){ + # Define latent state at first capture + z.telemetry[i,f.telemetry[i]] <- 2 ## alive when first marked + for (t in (f.telemetry[i]+1):n.occasions.telemetry){ + # State process: draw S(t) given S(t-1) + z.telemetry[i,t] ~ dcat(ps[z.telemetry[i,t-1], i, t-1,]) + # Observation process: draw O(t) given S(t) + y.telemetry[i,t] ~ dcat(po[z.telemetry[i,t], i, t-1,]) + } #t + } #i + + + # ------------------------------------------------- + # 2.5. Likelihood for adult survival from territory monitoring + # ------------------------------------------------- + + ### ECOLOGICAL STATE MODEL WITH ESTIMATE OF SURVIVAL + + for (ilterr in 1:nsite.terrvis){ + x.terrvis[ilterr,f.obsvis[ilterr]] <- 2 #firstobs[ilterr] + + for (klterr in (f.obsvis[ilterr]+1):nprim.terrvis){ + z.terrvis[ilterr,klterr] ~ dbin(mean.phi.terrvis[phase[klterr]],x.terrvis[ilterr,klterr-1]) + mu.x[ilterr,klterr] ~ dbin(mean.rec.terrvis[phase[klterr]],(2-z.terrvis[ilterr,klterr])) + x.terrvis[ilterr,klterr] <- mu.x[ilterr,klterr] + z.terrvis[ilterr,klterr] + + } # close klterr loop over primary period - years + } # close ilterr loop over sites + + ### OBSERVATION MODEL WITH RANDOM EFFECT OF YEAR + + for (iobs in 1:nsite.terrvis){ + for (kobs in f.obsvis[iobs]:nprim.terrvis){ + linpred.terrvis[iobs,kobs] <- lmu.p.terrvis+beta.obs.eff.terrvis*eff.terrvis[iobs,kobs] + rand.obs.terrvis[kobs] + p.terrvis[iobs,kobs] <- 1 / (1 + exp(-linpred.terrvis[iobs,kobs])) + y.terrvis[iobs,kobs] ~ dbin(p.terrvis[iobs,kobs], x.terrvis[iobs,kobs]) + + } # close kobs loop over year + } # close iobs loop over sites + + + + + # ------------------------------------------------- + # 3. DERIVED PARAMETERS + # ------------------------------------------------- + + ### 3.1 TELEMETRY DERIVED SURVIVAL ESTIMATES + + ## for WILD BIRDS + for (ageprog in 1:36){ + logit(phi.wild.telemetry[ageprog]) <- lp.mean.telemetry + ### intercept for mean survival + b.phi.mig*(migprog.age[ageprog]) + ### survival dependent on migration (first-time crossing of sea) + b.phi.age*ageprog ### survival dependent on age (juvenile or other) + } + + ann.phi.juv.telemetry<-phi.wild.telemetry[1]*phi.wild.telemetry[2]*phi.wild.telemetry[3]*phi.wild.telemetry[4]*phi.wild.telemetry[5]*phi.wild.telemetry[6]*phi.wild.telemetry[7]*phi.wild.telemetry[8]*phi.wild.telemetry[9]*phi.wild.telemetry[10]*phi.wild.telemetry[11]*phi.wild.telemetry[12] ### multiply monthly survival from age 1-12 + ann.phi.sec.telemetry<-phi.wild.telemetry[13]*phi.wild.telemetry[14]*phi.wild.telemetry[15]*phi.wild.telemetry[16]*phi.wild.telemetry[17]*phi.wild.telemetry[18]*phi.wild.telemetry[19]*phi.wild.telemetry[20]*phi.wild.telemetry[21]*phi.wild.telemetry[22]*phi.wild.telemetry[23]*phi.wild.telemetry[24] ### multiply monthly survival from age 13-24 + ann.phi.third.telemetry<-phi.wild.telemetry[25]*phi.wild.telemetry[26]*phi.wild.telemetry[27]*phi.wild.telemetry[28]*phi.wild.telemetry[29]*phi.wild.telemetry[30]*phi.wild.telemetry[31]*phi.wild.telemetry[32]*phi.wild.telemetry[33]*phi.wild.telemetry[34]*phi.wild.telemetry[35]*phi.wild.telemetry[36] ### multiply monthly survival from age 25-36 + + ## for CAPTIVE-REARED DELAYED RELEASE BIRDS + for (captageprog in 1:36){ + logit(phi.capt.telemetry[captageprog]) <- lp.mean.telemetry + ### intercept for mean survival + b.phi.capt + ### survival dependent on captive-release (captive-raised or other) + b.phi.mig*(captprog.age[captageprog]) + ### survival dependent on migration (first-time crossing of sea) + b.phi.age*captageprog ### survival dependent on age (juvenile or other) + } + + ann.phi.capt.rel.first.year<-phi.capt.telemetry[8]*phi.capt.telemetry[9]*phi.capt.telemetry[10]*phi.capt.telemetry[11]*phi.capt.telemetry[12]*phi.capt.telemetry[13]*phi.capt.telemetry[14]*phi.capt.telemetry[15]*phi.capt.telemetry[16]*phi.capt.telemetry[17]*phi.capt.telemetry[18]*phi.capt.telemetry[19]*phi.capt.telemetry[20]*phi.capt.telemetry[21]*phi.capt.telemetry[22]*phi.capt.telemetry[23]*phi.capt.telemetry[24] ### first year for delayed-release bird is longer, from month 8 to 24 + + ### 3.2 POPULATION GROWTH DERIVED FROM COUNTS + + # Annual population growth rate + for (ipop in 1:(T.count-1)){ + lambda.t[ipop] <- Nterr[ipop+1] / max(Nterr[ipop],1) + loglambda.t[ipop]<-log(lambda.t[ipop])## for calculating geometric mean of overall population growth rate + } + + #### OVERALL POPULATION GROWTH RATE ######### + mean.lambda<-exp((1/(T.count-1))*sum(loglambda.t[1:(T.count-1)])) # Geometric mean + + + + + # ------------------------------------------------- + # 4. PREDICTION INTO THE FUTURE + # ------------------------------------------------- + ### INCLUDE SCENARIOS FOR N CAPTIVE BIRDS AND SURVIVAL IMPROVEMENT + + # CAPTIVE RELEASE OF JUVENILE BIRDS + for (ncr in 1:scen.capt.release){ + + ### FUTURE FECUNDITY IS BASED ON mu.fec - removed the chick removal scenario + + for (fut in 1:PROJECTION){ + fut.fec[fut,ncr] <-mu.fec ## no chicks taken after 5 years anymore + } + + + + # SPECIFY IMPROVEMENT OF SURVIVAL + for (is in 1:scen.imp.surv){ + + + ## POPULATION PROCESS + ### DOES NOT WORK WITH SAME ARRAY OF POP TRAJECTORY AS 3 dimensions required + ## need to copy previous array elements + nestlings.f[ncr,is,1] <- round(nestlings[T.count]) ##JUV[T.count] + N1.f[ncr,is,1] <- N1[T.count] + N2.f[ncr,is,1] <- N2[T.count] + N3.f[ncr,is,1] <- N3[T.count] + N4.f[ncr,is,1] <- N4[T.count] + N5.f[ncr,is,1] <- N5[T.count] + N6.f[ncr,is,1] <- N6[T.count] + Nterr.f[ncr,is,1] <- Nterr[T.count] + N2wild.f[ncr,is,1] <- 0 + N2released.f[ncr,is,1] <- 0 + + for (fut in 2:PROJECTION){ + + fut.survival[ncr,is,fut] <-min(imp.surv[fut,is]*mean.phi.terrvis[2],1) ### invalid parent error if survival>1 + + ### probabilistic formulation + nestlings.f[ncr,is,fut] <- (fut.fec[fut,ncr] * 0.5 * Nterr.f[ncr,is,fut]) ### number of local recruits calculated as fecundity times number of territorial pairs + N1.f[ncr,is,fut] ~ dbin(min((imp.surv[fut,is]*ann.phi.juv.telemetry),1),round(nestlings.f[ncr,is,fut-1])) ### number of 1-year old survivors + N2wild.f[ncr,is,fut] ~ dbin(min((imp.surv[fut,is]*ann.phi.sec.telemetry),1),round(N1.f[ncr,is,fut-1])) ### number of 2-year old wild survivors + N2released.f[ncr,is,fut] ~ dbin(min((imp.surv[fut,is]*ann.phi.capt.rel.first.year),1),round(capt.release[fut,ncr])) ### number of 2-year old captive-released survivors + N2.f[ncr,is,fut] <- N2wild.f[ncr,is,fut] + N2released.f[ncr,is,fut] ### sum of the N2 cohort derived from wild and captive-bred juveniles + N3.f[ncr,is,fut] ~ dbin(min((imp.surv[fut,is]*ann.phi.third.telemetry),1),round(N2.f[ncr,is,fut-1])) ### number of 3-year old survivors + N4.f[ncr,is,fut] ~ dbin(fut.survival[ncr,is,fut],round(N3.f[ncr,is,fut-1])) ### number of 4-year old survivors + N5.f[ncr,is,fut] ~ dbin(fut.survival[ncr,is,fut],round(N4.f[ncr,is,fut-1])) ### number of 5-year old survivors + N6.f[ncr,is,fut] ~ dbin(fut.survival[ncr,is,fut],round((N5.f[ncr,is,fut-1]+N6.f[ncr,is,fut-1]))) ### number of 6-year or older (adult) birds + Nterr.f[ncr,is,fut] <- round((N4.f[ncr,is,fut] * 0.024) + (N5.f[ncr,is,fut] * 0.124) + (N6.f[ncr,is,fut])) + + + } # fut + + for (fut2 in 1:(PROJECTION-1)){ + lambda.t.f[ncr,is,fut2] <- Nterr.f[ncr,is,fut2+1] / max(Nterr.f[ncr,is,fut2],1) + loglambda.t.f[ncr,is,fut2]<-log(lambda.t.f[ncr,is,fut2])## for calculating geometric mean of overall population growth rate + } # fut2 + + + #### FUTURE POPULATION GROWTH RATE ######### + fut.lambda[ncr,is]<-exp((1/(PROJECTION-1))*sum(loglambda.t.f[ncr,is,1:(PROJECTION-1)])) # Geometric mean + + } # end scenario of imp.surv + + } # end scenario of n capt released + + } # END of the JAGS model + ",fill = TRUE) +sink() + + + + + +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# SET UP MCMC PARAMETERS AND RUN MODEL +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +## Parameters to be estimated ('monitored') by JAGS +paraIPM<-c("mu.fec","lambda.t","b.phi.age","b.phi.capt","b.phi.mig","ann.phi.capt.rel.first.year", + "ann.phi.juv.telemetry","ann.phi.sec.telemetry","ann.phi.third.telemetry", "ann.phi.terrvis", #"breed.prop4","breed.prop5", + "mean.phi.terrvis","mean.lambda","fut.lambda","Nterr", "Nterr.f") + + +# MCMC settings +nc <- 3 +nt <- 4 +ni <- 50000 +nb <- 10000 + +# specify initial values +initIPM <- function(){list(mean.p.terrvis=runif(1,0.5,1), + mean.phi.terrvis=runif(2,0.75, 1), + sigma.obs.count=runif(4,0.1,100), + mu.fec = runif(1,0,1), + z.telemetry = z.telemetry, ### matrix with initial values for the true states of tracked individuals + mean.phi.telemetry = runif(1, 0.9, 0.999), + base.obs.telemetry = rnorm(1,0, 0.001), + base.fail.telemetry = rnorm(1,0, 0.001), + base.recover.telemetry = rnorm(1,0, 0.001))} + + +### RUN THE MODEL ### + +EGVU.IPM <- jags(data=INPUT, + inits=initIPM, + parameters.to.save=paraIPM, + model.file="C:/STEFFEN/RSPB/Bulgaria/Analysis/PopulationModel/vultures/EGVU_IPM_Balkans.jags", + n.chains=nc, n.thin=nt, n.burnin=nb, n.iter=ni,parallel=T) + + + + + + + + \ No newline at end of file