From 34b85f256c80880f250b07e59d4117a1650c33c8 Mon Sep 17 00:00:00 2001 From: Boulangeat Isabelle Date: Wed, 29 Jan 2020 15:11:08 +0100 Subject: [PATCH] rebuild main analysis output document --- Main_analyse.Rmd | 17 +- Main_analyse.html | 324 ++-- .../figure-html/lambda, fig-TRUE-1.png | Bin 0 -> 101259 bytes .../figure-html/lambda, fig-TRUE-2.png | Bin 0 -> 97145 bytes .../figure-html/lambda, fig-TRUE-3.png | Bin 0 -> 111111 bytes .../figure-html/lambda_sign-1.png | Bin 0 -> 69990 bytes .../figure-html/plot_coeffs-1.png | Bin 0 -> 65205 bytes .../figure-html/plot_coeffs-2.png | Bin 0 -> 61082 bytes README.md | 1459 ++++++----------- 9 files changed, 679 insertions(+), 1121 deletions(-) create mode 100644 Main_analyse_files/figure-html/lambda, fig-TRUE-1.png create mode 100644 Main_analyse_files/figure-html/lambda, fig-TRUE-2.png create mode 100644 Main_analyse_files/figure-html/lambda, fig-TRUE-3.png create mode 100644 Main_analyse_files/figure-html/lambda_sign-1.png create mode 100644 Main_analyse_files/figure-html/plot_coeffs-1.png create mode 100644 Main_analyse_files/figure-html/plot_coeffs-2.png diff --git a/Main_analyse.Rmd b/Main_analyse.Rmd index 9fd34d8..af267fe 100644 --- a/Main_analyse.Rmd +++ b/Main_analyse.Rmd @@ -3,9 +3,9 @@ title: "Main_analyse.Rmd" author: "Isabelle Boulangeat" date: "23/01/2020" output: - md_document: - variant: markdown_github html_document: + keep_md: yes + variant: markdown_github editor_options: chunk_output_type: console always_allow_html: true @@ -97,7 +97,7 @@ res.tab.eval = as.data.frame(do.call(rbind, lapply(eval.models, function(x)unlis res.tab.vars = as.data.frame(do.call(rbind, lapply(eval.models, function(x)unlist(x[5:6])))) kable(res.tab.eval, caption = "model evaluation", digits = 2) %>% kable_styling(bootstrap_options = c("striped", "hover", "responsive"), full_width=FALSE) -kable(res.tab.vars, caption = "selected variables")%>% kable_styling(bootstrap_options = c("striped", "hover", "responsive"), font_size = 10) +# kable(res.tab.vars, caption = "selected variables")%>% kable_styling(bootstrap_options = c("striped", "hover", "responsive"), font_size = 10) ``` N.B. Some models do not fit well. The variables explain more or less depending on species but also on processes (extinction and colonisation). @@ -133,7 +133,7 @@ data.frame(unclass(table(coeffs.all$var, coeffs.all$process))) %>% kable_styling("striped") ``` -## Variables effect coefficients +## Variable effect coefficients This gives a summary of selected variables with significant impact on colonisation or extinction and their estimated coefficients. @@ -153,18 +153,15 @@ coeffs.all[which(coeffs.all$process == "extinction" & coeffs.all$Pr...z..<0.1),] ## Lambda (1-ext/colo) An overview of the performance of populations by site of by species -```{r lambda, fig-TRUE, fig.width=10} +```{r lambda, fig-TRUE, fig.width=10, warning=FALSE, fig.height=7} all.sites = unique(DATA$cover[,c("SemiUniquePLOT", "YEAR", "SUBSITE", "SoilMoist", "TotalCover")]) pDATA = proj.dat.shape(all.sites, DATA$climato) -str(pDATA) -summary(pDATA) selectVars = c("snowDays", "GSL", "GST", "SoilMoist", "TotalCover") pDATA.all = na.omit(pDATA[,c(selectVars,"SemiUniquePLOT", "YEAR", "SUBSITE") ]) pred.all = lapply(spDATA.fit, predictions, pred.data = pDATA.all, selectVars=selectVars) pred.tab = lapply(pred.all, function(x)data.frame(sp = x$sp, lambda = x$lambda, YEAR = x$YEAR, SUBSITE = x$SUBSITE)) pred.all.together = do.call(rbind, pred.tab) -head(pred.all.together) pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>% ggplot(aes(YEAR, lambda)) + @@ -180,7 +177,7 @@ pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>% ggplot(aes(SITE, lambda, sp)) + theme(axis.text.x = element_text(angle=90, hjust=1)) + geom_boxplot() + geom_hline(yintercept = 0) + - facet_wrap( ~ sp,scales = "free_x", ncol = 4) + facet_wrap( ~ sp,scales = "free_x", ncol = 7) pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>% ggplot(aes(SITE, lambda, sp)) + @@ -190,7 +187,7 @@ pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>% ``` Simplified version focusing only on the sign of lambda -```{r lambda_sign, fig=TRUE, fig.width=10} +```{r lambda_sign, fig=TRUE, fig.width=10, warning=FALSE} pred.all.together$lambdaSign = ifelse(pred.all.together$lambda>0, "+", "-") pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>% diff --git a/Main_analyse.html b/Main_analyse.html index 74c02c9..b616633 100644 --- a/Main_analyse.html +++ b/Main_analyse.html @@ -343,9 +343,9 @@ str(DATA$cover) ## $ GST_1989 : num NA NA NA NA NA ... ## [list output truncated] 
mostRecordedSP = DATA$cover %>% group_by(Name) %>% summarise(nobs = length(Name)) %>% arrange(desc(nobs))
-selection = as.vector(t(mostRecordedSP[1:22, "Name"]))
-selection = selection[-which(selection %in% c("Arctagrostis latifolia", "Persicaria vivipara", "Carex aquatilis", "Cassiope tetragona", "Luzula confusa", "Luzula nivalis", "Salix arctica", "Poa arctica", "Eriophorum angustifolium", "Petasites frigidus", "Vaccinium uliginosum", "Dryas integrifolia"))]
-

N.B. I have quickly selected 10 species among the 22 having the most observations, with no problem with zero-days growing season in the previous years before an observation, and no problem in the dataset… selection procedure to discuss

+selection = as.vector(t(mostRecordedSP[1:30, "Name"])) +selection = selection[-which(selection %in% c("Arctagrostis latifolia", "Persicaria vivipara", "Carex aquatilis", "Cassiope tetragona", "Luzula confusa", "Luzula nivalis", "Salix arctica", "Poa arctica", "Eriophorum angustifolium", "Petasites frigidus", "Vaccinium uliginosum", "Dryas integrifolia", "Stellaria longipes", "Saxifraga cernua", "Festuca brachyphylla"))] +

N.B. I have quickly selected 14 species among the 30 having the most observations, with no problem with zero-days growing season in the previous years before an observation, and no problem in the dataset… selection procedure to discuss

Prepare species specific datasets

@@ -409,7 +409,33 @@ spDATA.check = lapply(spDATA.list, check.nas) ## $`Rubus chamaemorus` ## ## AA COLO EXT PP -## 315 34 22 191 +## 315 34 22 191 +## +## $`Salix reticulata` +## +## AA COLO EXT PP +## 491 14 24 161 +## +## $`Geum rossii` +## +## AA COLO EXT PP +## 16 4 8 153 +## +## $`Empetrum nigrum` +## +## AA COLO EXT PP +## 448 18 23 133 +## +## $`Bistorta bistortoides` +## +## AA COLO EXT PP +## 43 20 22 96 +## +## $`Hierochloe alpina` +## +## AA COLO EXT PP +## 603 39 36 61 +

Here we can see how many transitions we have for each species.

Fit models

@@ -624,136 +650,95 @@ Rubus chamaemorus 0.50 - - -
kable(res.tab.vars, caption = "selected variables")%>% kable_styling(bootstrap_options = c("striped", "hover", "responsive"), font_size = 10)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-selected variables -
- -vars.colo - -vars.ext -
-Vaccinium vitis-idaea +Salix reticulata -GST:SoilMoist;GST;TotalCover:SoilMoist;GSL:SoilMoist;SoilMoist;TotalCover;GSL - -SoilMoist + +0.18
-Carex bigelowii + +0.75 -GSL;snowDays:SoilMoist;snowDays;TotalCover:SoilMoist;SoilMoist;TotalCover + +0.21 -SoilMoist;TotalCover:SoilMoist;snowDays;TotalCover + +0.73
-Ledum palustre - -SoilMoist;TotalCover:SoilMoist;GST;TotalCover +Geum rossii -snowDays + +0.25
-Betula nana + +0.70 -TotalCover:SoilMoist;snowDays;GSL;SoilMoist;TotalCover + +0.28 -TotalCover;snowDays + +0.86
-Eriophorum vaginatum - -SoilMoist +Empetrum nigrum -GSL;snowDays;TotalCover;GST + +0.16
-Salix pulchra + +0.77 -GSL;snowDays;SoilMoist;TotalCover:SoilMoist;TotalCover + +0.11 -TotalCover:SoilMoist;SoilMoist;TotalCover + +0.67
-Eriophorum chamissonis - -snowDays:SoilMoist;GSL:SoilMoist;TotalCover:SoilMoist;SoilMoist;TotalCover;snowDays;GSL +Bistorta bistortoides -snowDays:SoilMoist;GSL:SoilMoist;GSL;GST;snowDays;SoilMoist + +0.20
-Persicaria bistorta + +0.71 -SoilMoist;snowDays:SoilMoist;GST:SoilMoist;snowDays;GSL;GST + +0.04 -SoilMoist;snowDays;GSL + +0.59
-Dupontia fisheri +Hierochloe alpina -GSL:SoilMoist;snowDays;snowDays:SoilMoist;GST;GSL;SoilMoist;TotalCover - -snowDays;snowDays:SoilMoist;GSL;GSL:SoilMoist;TotalCover:SoilMoist;TotalCover;SoilMoist + +0.18
-Rubus chamaemorus + +0.77 + +0.28 + +0.77
+
# kable(res.tab.vars, caption = "selected variables")%>% kable_styling(bootstrap_options = c("striped", "hover", "responsive"), font_size = 10)
+

N.B. Some models do not fit well. The variables explain more or less depending on species but also on processes (extinction and colonisation).

Selected variables and coefficient, sumnmary table

@@ -770,31 +755,22 @@ Rubus chamaemorus return(coe) }) coeffs.all = do.call(rbind, coeffs.all.list) -str(coeffs.all) -
## 'data.frame':    79 obs. of  7 variables:
-##  $ Estimate  : num  -0.939 3.268 0.95 -1.496 -2.566 ...
-##  $ Std..Error: num  1.11 1.45 0.46 2.04 1.48 ...
-##  $ z.value   : num  -0.85 2.257 2.065 -0.733 -1.738 ...
-##  $ Pr...z..  : num  0.3955 0.024 0.039 0.4637 0.0822 ...
-##  $ sp        : chr  "Vaccinium vitis-idaea" "Vaccinium vitis-idaea" "Vaccinium vitis-idaea" "Vaccinium vitis-idaea" ...
-##  $ var       : chr  "GSL" "GST" "TotalCover" "SoilMoist" ...
-##  $ process   : chr  "colonisation" "colonisation" "colonisation" "colonisation" ...
+# str(coeffs.all)
-
-

Number of selection of variable per process

+
+

Number of selected variable per process, accross all species

data.frame(unclass(table(coeffs.all$var, coeffs.all$process)))  %>% 
-  mutate(
-    colonisation = cell_spec(colonisation, color = ifelse(colonisation > 5, "red", "black")),
-    extinction = cell_spec(extinction, color = ifelse(extinction > 5, "red", "black"))) %>%
   kable(.) %>%
   kable_styling("striped")
+ - @@ -802,90 +778,166 @@ extinction - + - + - + - + - + - + - + - + - +
+ colonisation + extinction
-<span style=" color: red !important;" >7</span> +GSL -<span style=" color: black !important;" >4</span> + +11 + +6
-<span style=" color: black !important;" >3</span> +GSL:SoilMoist -<span style=" color: black !important;" >2</span> + +5 + +3
-<span style=" color: black !important;" >4</span> +GST -<span style=" color: black !important;" >2</span> + +9 + +6
-<span style=" color: black !important;" >2</span> +GST:SoilMoist -<span style=" color: black !important;" >0</span> + +5 + +2
-<span style=" color: red !important;" >9</span> +SoilMoist -<span style=" color: red !important;" >6</span> + +13 + +10
-<span style=" color: red !important;" >7</span> +TotalCover -<span style=" color: black !important;" >5</span> + +11 + +7
-<span style=" color: red !important;" >6</span> +TotalCover:SoilMoist -<span style=" color: black !important;" >3</span> + +9 + +4
-<span style=" color: red !important;" >6</span> +snowDays -<span style=" color: red !important;" >7</span> + +10 + +10
-<span style=" color: black !important;" >4</span> +snowDays:SoilMoist -<span style=" color: black !important;" >2</span> + +7 + +4
+
+
+

Variable effect coefficients

+

This gives a summary of selected variables with significant impact on colonisation or extinction and their estimated coefficients.

coeffs.all[which(coeffs.all$process == "colonisation" & coeffs.all$Pr...z..<0.1),] %>%
   ggplot(aes(sp, Estimate, fill = var)) +
   geom_col(position = "dodge") +
   theme_bw() + facet_wrap(~sp,scales = "free_x", ncol=4) +
   ggtitle("Colonisation process")
-

+

coeffs.all[which(coeffs.all$process == "extinction" & coeffs.all$Pr...z..<0.1),] %>%
   ggplot(aes(sp, Estimate, fill = var)) +
   geom_col(position = "dodge") +
   theme_bw() + facet_wrap(~sp,scales = "free_x", ncol=4) +
   ggtitle("Extinction process")
-

+

+
+
+

Lambda (1-ext/colo)

+

An overview of the performance of populations by site of by species

+
all.sites = unique(DATA$cover[,c("SemiUniquePLOT", "YEAR", "SUBSITE", "SoilMoist", "TotalCover")])
+pDATA = proj.dat.shape(all.sites, DATA$climato)
+selectVars = c("snowDays", "GSL", "GST", "SoilMoist", "TotalCover")
+pDATA.all = na.omit(pDATA[,c(selectVars,"SemiUniquePLOT", "YEAR", "SUBSITE") ])
+
+pred.all = lapply(spDATA.fit, predictions, pred.data = pDATA.all, selectVars=selectVars)
+pred.tab = lapply(pred.all, function(x)data.frame(sp = x$sp, lambda = x$lambda, YEAR = x$YEAR, SUBSITE = x$SUBSITE))
+pred.all.together = do.call(rbind, pred.tab)
+
+pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>%
+  ggplot(aes(YEAR, lambda)) +
+  geom_point(shape = 1) +
+  geom_smooth() +
+  facet_wrap( ~ sp, ncol = 4)
+
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
+

+
pred.all.together$SITE =  unlist( lapply(as.character(pred.all.together$SUBSITE), function(x){
+  strsplit(x, ":")[[1]][[1]]
+}))
+
+pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>%
+  ggplot(aes(SITE, lambda, sp)) + 
+  theme(axis.text.x = element_text(angle=90, hjust=1)) +
+  geom_boxplot() + geom_hline(yintercept = 0) +
+  facet_wrap( ~ sp,scales = "free_x", ncol = 7)
+

+
pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>%
+  ggplot(aes(SITE, lambda, sp)) +
+  theme(axis.text.x = element_text(angle=90, hjust=1)) +
+  geom_boxplot() + ylim (-1, 1) + geom_hline(yintercept = 0) +
+  facet_wrap( ~ sp,scales = "free_x", ncol = 7)
+

+

Simplified version focusing only on the sign of lambda

+
pred.all.together$lambdaSign = ifelse(pred.all.together$lambda>0, "+", "-")
+
+pred.all.together[-which(pred.all.together$sp == "Rubus chamaemorus"), ] %>%
+  ggplot(aes(x=SITE, y=lambdaSign, fill= lambdaSign )) +
+  theme(axis.text.x = element_text(angle=90, hjust=1)) +
+  geom_bar(stat="identity")  + 
+  facet_wrap( ~ sp,scales = "free_x", ncol = 7)
+

diff --git a/Main_analyse_files/figure-html/lambda, fig-TRUE-1.png b/Main_analyse_files/figure-html/lambda, fig-TRUE-1.png new file mode 100644 index 0000000000000000000000000000000000000000..247bd28d33555ef6c16755801b473812f9e72d69 GIT binary patch literal 101259 zcmcfpWmp?;yZ(zqp~Z@o0xeps6xX6f3dOxx@#5}I0u(Rq?!`)ry9IZ*;%>o$6CC#N zJZt@bYkm9t_lv!cV;@IeBs0lmCNpzi=lwa)E5VBL5}4>j=m-c1m{O9TJ|iF?wIU!q z$wGYwEKyN?>5PDY{MuYpRMApYLe$pE)dVfS+HY^(ony}nR8H3In!lIK&yj52P`+e*06tAW%m5wKn;t&h zuK)EVA^IgEK}rb|A^VzRP~>+3hA97_EW&5m>Pn;xI)$=NxD+s4tF@Lz+d+D>e-`gE zi)+P6@!!~!8&q!n#P}YGcvPS|;DYFJ9ZfE@t4Rho*sf(0?qu<(TQ!{Sa1`_kH#gZ} zT$}j*R8`sHdm>KcVArMi^BcCZDZMN?i@jgln?DRvofosp?GE+tGp*aq0#23O$QjSV}3)$36U{T{MsKt2r~b zxH^Qw3+m5x%DBe#eS{m@uiQ5C`&VpNTANE}?R(9W6Mp2@;0CRf;~s|_Yvqg+=;rSG zlmU?Yr`zi$cS*RXLnbe52fTgipTLB2o0k~phu`V<-8)hXj96Y zIVJ~7i)CHMI)8ga<~bqvB)JU1-nGx?tcOpPp$^|lHJ<%e3JX`5W2)cXyHe^JCIn}K z1Gs+scsFBKqCI_=#2N%am-~W*YxGX|*)qY5S)m`+GxnFoFJFh~Wt$KZLt8(b@Kq%b zwKccL-S_tD`CZiwSJ4{<3HlOL5Zy);5lk2#{|?YL#1ryj{9S( zdv-)h_mB^A@K}-xJfp+w8EGpq-4EPTcTHg&pYv127I~{0v%WrX(93EoPZu#_)dW%v z(A_-Mf!ucMm22-Xax2xWy-sgXAHLGN+Ph_kGu&F_u6S>#AJ;ncXbv`3$ad^oeQ|VE zC9;0DPH;4q+!fqKdTM{Xp^h!?mC%;H^xEP|?X`UPr=l~zv!E(Q z2(9P4yKcykXTRr==dfp+&XZ{Vf%Ns}Rc1SKN4;*i#VL{IGKj$aRjuhxf@P631(g*( z8G%DP*Q$~}`w8uHR6`M^8J~pWy&Hl|mJ9(>EJO0Gf*b8-W;Mj%oAM`Ta`&pLDm!*@ zbfjvway-o+@xrJ@jl$x%R4(e&PqQiG6*fPc<0$hh;}w=qe#8!{%y-GMesXh5`lryQ z;`8pf;gaVQjR`)RJ(WyKD_0c z-OcN6a%Z8_u?Wiat(u!Zy=D2(X-ZJ|;WvV5xHn=}CRi0o8{-fGK^Q^m(?=CI#Qk)P zCIYdm9u$0J{3y^{Zqf2@B2lQh&vEFhe6*ZV%F3>$u#($ z@Xg31DJCyT@f*#`(jCA*{8acR8 zr{pJL8kcnrDeHnjd3mxIe#tvnr(o zH+(_vfBR;hfjFJ{rB4Xz%L)PBR@|mi2wZdHyI06JRR1@lxl*xB7DgcV+hXE4+akRI*uqL6!rVJ|NRiQK;A$o2viFJ_U<1)?Jgeh z>z&iRe@Xa{UkDe+XA-p+Q697SW#~UP_YHND_KCx{J!OkroW!A2?o<>qx7PjH%7Q?+ z(CsSrR$ij+FZ-qiq=5t`j)z;Q(5bGjY*G`t^oQ4wWjJiUIB@$|#Dd#q!E6PeRzdJ8 zSEJFzdV!7<&nE=UYQARZ!*M6x7JS~8+xu=$Z>U74k$y_!m+$h!{l!)|_tH;$4Z+5r zxI9~}93vkL)k?HqJLXr3Vhi1t%{DmyF60!lXKA^0n2nTh2R+3_q-PUqxCq-lfR7EtQy`suN$i)^?X2QQz2$;3al2wQ@rPs6fh@bh) z+cXyc2uxmVbWJ_%xvktxfvXi1&T$Oyg9wrMhF=Q^L9xmz%XvmK5$~FbDLNP|TNc24Mv@*&i`Za-3;i z+{{lBbv+8NHfV*+F)qsuds1F=m?T{^K?Q1*2BN4caTN zb?|Bgm2G;QG?5GTF}s=Ku$=L)@>!KBVsl-*IY2!=US``q*E-VCek+$^G*@G7xwf`` zP+AP!pg5)y@+?>B4*sn2Wk^CJF}3MBb+`dv@cV z(FC#h^ULZ!zeMUv^qzKvu;Qsu#+6lwomYB1=#s0;egCay3p13+Nva)ph9!;O{8uTm zAiO*>nZy3IGXrhY3S9-ydOzG*lLRY6htsVOIpf2C;WD|0-1`IZw5``^ z>J%Fd2|lhzH`+ResFfM-Z5y`N*C35{XqF#V1(tiI$xy||*e{>Uza@iqAbUE8?CRR~ z0iMmdmKD;~6YtPgP{-f->s&hY(Ih5$!e{4QlH*(fnaL|Uf*EhT|J{=wee71^qkw)p zGKcf7%0#K2&}q0WxWnmiX)(>d`}I8&?0ztr1LjG$Qhq$*dEn@>A)2N$8Z}j-!{5;e z#ml_(zgoriHrpD>w1we-GTvJz?2PBh3=5Ef_SyDk$iYW#CS4);TZ0@UoaC3SD1-IT z49kP&8yk=<5wgeGwmyTP>~ibz9dMP~3=7^^ahML1+%z!@^7fyLB0s_Dy&x-)*xmhI z%C*?w%;D8>z*M_nfW^BJWi{)x8LvJUxVG96H17~|;@A<0$u;7A9I98fxqyg5@|9sI z&N8mNnU%70(eIgo)$v5;eO0w}|4X7TbOY}y-#I?RCQF0d&YEVv?}T3OmDcvW`g7Jx z6**J!UEHBxuK#Vxrr;#ut3R9la_x=^`Fy0V3-()a%JNc-(-qu*B2aT>zg1Lt-oR&Z z!_UsfzsRQBU7c)hr8y3llM1%n!ECKxAMfcs%$>pJxyJ~!o#2R zyEab@qvhjt>*tqsUfp%EcLdh#-xlGetr$k|cB|pou?hk!=asAcY0GnC>@Qb? z<^!MV{<702b12n?yePGok7z8d2@geh$}LYBF5Em?=C?6mo2CcxsIJ!V5bv9|zbbV- zTzYnM4IGrA;CZ3KY=4`1L3;P*dt}o%SEONBo&CC%XdLyrEIfH_6MYHNrI+`7K{G zU;I(@n+RhX=LkRwTU^a$-$W~TUB4c|U&I09Z+-vL9mOoGeJ8JIV^Nff=n%POp%XFx zw%uptF7^7}?#S^Ck<+V+*5lT+`Lo+cBr*M^ zkq)u^r8y4KwmwnvG|ogVzY#Be0gc+>m_pWbRpppkajGPSD$lSu z@h$hq;#2MUe)O1q$_M{Q411>Q;kNF7d2+DI5*);Uw_Y;MhKu^~o02q7(sSp0M#{r+ zgaEm+@A#pHN>DK~PNtfV{gK2u6d~1g{Y3|_p@iA@eZKJl{o9XcYIF?~o=LLu_1#`A zWgsTUSupzf53$hwC`p4s61Bh}MYL|*&Y7M1HZ+rWL8(yDFM+~%WChqLlW36l{^97u zd3d_W+h+lW^4RePdTA4|2;L583ePHV0tDn#;$~UYm;w2=$Fq*isuH)#aF)PZz)eI* zziWdWI#1XesKth3&-~Y<3S{0xNW9TE@#yet*>7PM53145F6Wz;$_W>S~g0 zfsIoThSnf~XPlBsU2?uK%e1NB>1(a2QOTF_xpiT$*xH|y@QmK>D0o*yAmwmwkA;Mg z>kW4uf7i`fro>4>!0aq3GRS>-u1O?Oc%wvXcoxMV^A7aR?}y`FtW{fA^_Nsia8>g_ zxyX5=L&?Wd>gR||zW9NLkZF-@F|nA{#wzy43<}Mp!YcNli2gicf;IXKmG3$%F9Rsz z53xPBMQhZ9SXKO6PU1^vq&*68SiR`WQ9Gg-;v>NuR5j}HSL7wN{nNpCq_BYQ`mfCF zkfG-^hWbwl*s?!X{x}#e>Kd1SS2ucRdRe- zc-Jj~5EPdkJ*ec$^&K**EMr_9*AyG~kk!B7DK5i`oV*&#Y1W0tw7cP+W}If?;p3@) zhMXhwGirkJP%Q)M*%l@3QTOUZfyQ5#emQ%~voj7&UElQhMMk_UzBfQI5L1sj1H|7${TnYG08wu4GQ)v|oD+lEb8f%+b~T!s3^cOqjVqan5UwzTC<<>h z_4sNZrt$7YV^&kPjl>A+x#~&9UQ@3r2J=+cOdJMrEmK~Y%XFe+c%$0~#C(fXm!lYjh>ril-s#kbp%u?w5x+JYTF$u5Xr=b7@SGpQ?Y zESO)@W(TJbzPwlqPzlcJxf9bW92*#ABIvbe5snzE4CGkN$P7Wm_W`$$5_^OvGxQB4 z)=^!3so>*+xICZT;B{BRCgvKfQQTV3WFh05K0+I1?fQv)LW^%<6(^QwHYgHvA|EIF z`MyKSSe`jtONC^cw;k7nE%T?NR0ocG#xoB3UL{^&(q&DBjVaAj@Ej$NCL{VK<3<-W zwG3BtbG=03qPFGEx)BFkbm~je!23xyzJO|sFC^dT(0E?hyfWh?yJmVqfm$Z-w=C(0 zIo?ZjK}JimNgu8nR`sLE2~w}!v1UEl>-};C96OM;VM~T5=D2onDPA$& zkF#R5SR(f(@~n>bblgHWvmrwU4Q#k-h@LU^=a7# zwH_I&QEFwM5z<7>xpPeBs=XmgJ8ZQTpLddlN#h~}_~^>_Kf%}SQD*2)+=HD$>idEo zJMnp}C75kjd5V#02W&|6aGM!NHv82p{TF}YBWXHBUW<@!rX8jjtbNAR{~EfUNge(N z9^>{YGFE7vC|J6WU+Rj*qcnG!Cx7FLClv|`JM*vCHz(9}%hvfBq#uvUbOg?CCx_#^ zl5h2_hsm_#zFGbR+V_N22xN+{-ho}>j0&<(-cQ7RS&8P%`#9)29?Y8~#zQD@AcSvB z*pW|%kt1?T7kTHvRJh= z5*8?tYjF(-?PUcd#a1ZVG7fpoLzWbkgzNsti*|FlaZ-r-&t;{aiJiO(#gNkhXQ`3Q zNE-YEo<3JOZ8uU&FtS^6j})CdE!{)cPB?$Uhe=PPAd-UKJ6ruDT<6`O{$zMhVNF*Rh1RN=HCCFZLu-b|*tO9De> zW~+soIQBUmz2Wx^P)Jw(Z||;{%{81zyeJ{|o!P%S1R$yo?2)NUmsBoC0GW9+LUv8lPR-EY+;v zLYU{sf728(oG2u1Wa_85|0GAaKfRt1Jt3769nuIz7*oJS{q)ABOG&TVYNR0B6BkwD zjSa4-UgLki60RYP$m0warb9;m!@G$fVb)ztM2PnJh&B^V|4*Is-W2|GEolUT=9xNs zL%PJ5&*yY;SF5$|#`x&GEaY=hA}5tdg2k7mQcFWTW>cyR?PAV)Pl7a|%vym7s7 zixKz&w|VTid`;q~{L6Lo1d#a*nK7BM-n+tlhf>iXtgRvL<)U~8g_^|Fj>oxg%a9h6 z^RPv^K>tLdo)<5nOF+x`Psb_R`3anDbdK-Lgn$1;SWleikYurOGY3~rUsF)aF9Xvf zM-L_s_`bQBPs~5@)HhU4<|ocu2pqRY*JnF{`l1vhD?Fb3t_Pg4Z{^q}$UU97RC~gS zX;FA&MEU=Dsz%*VCxeli95$snjQd~mDD zuO%(XNZP{%d^KpLIm^Ieqep<1i@!fQ- zGdM6qE0E)2d#qzHjdz06e4LtA)~W3WwI_pGQ5n|J?bRuabvv3$Vy-K=qMw}~km?2Q zTgj&T6a<ef4BTL7}tlg_;N>RN<@hSEG&$h-@^jAA0h)b*+#<;!>b?>xQm@k#*#7lx{=(9ER( zxG|{O6x8Sf>rW-AfHb*<@FknnCkDcfm#}pb)MV-r(1e+foe4kRri2bWcy@2g8;Nf{O^yIMZw;mR(|Tb z$yK$U)y@9XD7x`+R`+~Xce~d88tq|07uh|aG^_5cjn2!NAjm@;qM3cK{lIL)F20+r`~tFuf%K&}&c zGK3E&L^pVEWm#{>JVF%FvbdAnTjlL4JkgG?ULN#}9-Qw@-~3`&QZSS{vAB8uB2$e% zu22RseXHQhFN!#Fo&Q}k@R`PVLuUZPI!;B04$1`*w-kvNThG#O6gHQAdgc3E>TWQY(yFQ+6r z4c@skZg(%d)^4i$cM15P{2v&* za;aSJDV(mNKhlzyV$%(!yyr0T8TGiBw-u^%gbsgI+SD3lG^jD`!KwsZ9DMbVVmJ%n zM9LtA%-7mt6DYjQi@s-F1VHIg)tqUjXOS*B!j(_JK@7D&q!%8hBdl0MJAVoby8r<+ z{!WzlrUF7jV*H?3KV3H<(4V>8sJ^LN!dp08Us`M%Sk$AYoF01yoy~EHNtKbl3Gh~o zfc+T|E0YGDm1P6O{!6(C!J@@tLp56VK%{Icm$mPn)pSX^vkI5>+{W^+AD(IC^jfLh zwk490-fYQ@hCLA`z}HTp!ADdpMtm)K{!kDj*q!tu{fW9~%{#%ER`kt@FuVg$6!fRn zQj`1e;dUfWlV+8fK9=C+%CE0_)-@6yfmV~zus<0v{Oc#aRxH6_5pzsTymoJp@&m^h!;u$P>;^(1JoSCphBT2t@P*R!{@~+y5YX zbdNB(%+J66eYOOFRf3qXzwHmsFMz%vwi@64FS+qQD+m}^I03~4CBrI_^0l-0NrBGyxTV&0z8r=Ct}B8UTD(s zdJk=h^lYIr?LZa4?SlEVwMY2&ie^|3z~hy6xYH420&}-F%MTHMlm^?B<$1fnmmTC> zBLI?QYcuSqtIKd(3oeIBgcGSo-QcTh>vuC%ZS(U}2Vsck-Ivdp~d7id;d{JD@L_``t( z*Kt4Cz$&uzqsr$4d&5Kp1FC4=@33R6UQ2iEYec@)x&Lh51HjX0GGEu5XqxIjRdOB{$1YemF!T<{dwvBd!FF2 z{wshmBAo^+fGLUd6+Z-GliRQk|Gc!98($UdIyo6&NdPAIY5tses2B>r4_$qP9bD{n@xoaj_T%0lHuYg%r^*7^y=q#a16}=AfX9 z_O;TxCcqJL@3_25Zn_41H*Vq#QbCBg-;YX2-OVMR1GD_4lB5Lc{ORcz5zY|7ePFKv z-|3}iI63Q%dpd0O!=haW<~JBPC~e%mdx?Efv_#`ln`D~3AhX9IX+FDjvveUSDiOq@Qgz2(YdKxQyG(dD z(<5l|Xcn5X1RD}DU5rA#$#x&PWBtVvsFKeJr16rMzlYQ{?Xk30xC>&=6SEokYt&e2 z414A{4=XDO{8{-&A9fsDi1O4IOXq0;>j zK=Jlv>0R}!%Pu#2TKfWaeDcqX z<*G>UG{HZ6X{A1kwJJ2V2%g4VwFF^vavRYz^tv52YM5iIe{d6l?^|yLz2iBgAJkLF{=*SpVD52K- zMA1k5KLKRRS2hR1{Axh&e!K{@Lo?B$C&b})^mXss zgL~hL39fZy&s7XP;ebD5xxsAE^wfc)G+9lhVEc@lqa*mo04D3oF?SrwHCF|6+}TYrv+%+(o#lxp3-+9~3^B z=Lg`9?1jzZP!kG!P3^9C4T@1>zWLnJ1^Z4tf2=HRQFEFaCcrr| z)3u{UqYocwDtMfuX@;N?6RX3E^V$VRtClD_kUgBV zfZ{du>}O6y#TT8Fd~;AJhw+{sJEkuARj;ID=c+`~V}f9-x!%!Q^mcTxH_M>2E4Ivh zmdpjp9YcUu3f32y$zZTuJmnsmL}F>!rkqbYa>^h4xNtKNVBca3_;khV9$T+XXHK9I z{0FRi+*YhTUfZHukVDVX`iJ|QAq?Bded!7TbyuRq*SunWYqke>;Tn}D)-_n`Cohld z&SD$~XYeo*LE&w1JDx<-CCg#>{ryEtYM1w=UvMQPIJw~=W3YvDH_Zu|(R~e#gK|CW z(Qrww(u$Sp0LhhJH1#|u<_a^aJUo|;pBUv`^{T{Ff?6fIAA$u>@9o|3AR7t~+|I+C zR;3yfCMK5oOM1WR>cjwX#=E=uTklKy@yG$&Q0k%WFwfod-x-$hMN@JuW`~s^p=fd)l|7R-rxw|8<)5$7ZLxMr}5 zFUP0Ooh5e*_oTlk_<&A|q6&```c#vRn*xbaVv2{gkl^!60O=3Q9H&V0Co@awRTx zmBcko;426IuvmAZ!My$g7LSM<*0c#ikIXB_7*6DrE+ezej#)CEZJtV#%?Hg1@7q|2 zqwGQ%kJA`sWNbtgPnKA?yCM8j@WYrddAj=R0$*OHfealJETJ{KM@-UKHZff5l7hj~ zGX^`A1G! zJT1(*UUh8r5f)_HH~8Ros4vC7W_G;VOq^$xmw@|z;2$Jhy5DkXJbHDipF|$903H}} zk@U`xsVo2w+g+b8w$wmZ&dS$?OkPcGMqEMBHrDv=PIg>yOV`CoOusPh-JTBDLdw;3 zJmd8#sE7jdwYj=>?=}hq>s=$d=Fbfa4B0|F&2vSvJXL(0)ej9z0V&+OB(MR29K^VVD}VoPdQLd{$1$Xf_UI|QrDLA`J55GX1e}%x3 z&S0=ax2JuYez|<$q}M6&a#c^W$BsIzN!JR^egxEY#;^|tTRER4jlJ&|Xy+@516_fwugp4Uzj)mrK*M$VJ@Fz177vdW#I3^*`5S=-A9GW1HH!0zy> z|7=y?gEP(Vqiau5fT#YXo0xJJBP^;JkL+X@yToCKwxV@cCNN1)HL6gEHI=9%jz&nwAmKGMejjP3l637fx;Zoq<;1R?4jPQ8bmqzUKEag=}}|J z+AwMxRHqn0T@wKA4l!4BQro0$<*E+XNyH!x@_H2L*s`_YdK-Lbx#Ly775||RCu6$v z=B9uoG1$*!dh8ZjDQ+RMr!AlbgtWXP=*c-|JYy_EAk`_3ZqF2>TrhYeYgZ5h^*?&9 zVb_lHl!zmZmHk`$uwGf7;)gmeqAehF;5v*RF?Oomak{30i=yd}6MLCF-M;Pj*+1at zt@|7!(I$~?$RE6L2FF#6?V&uO08yuSh}j(GYxd6Os09Jr;X%^33OsFWIUjoxe@5WR z^R(TbmwFF9yMouYSSQ&Y`aj()W78XLlp%G1}m@!_9>E zms`BCM5bm}k}3tJbiStJ`Ll<2Ar>~qQ7SdA&6x=X2kI9G(aBzBgea(i@Fd zIe1SnZ}`XfKFYv4(nCm@GM#!;-Bqw;_^NVnL8)KF6g0`b*e)$qPe-u_S?uPU5T`Zh(2>6F%Lr9O4s z?}2jHa!@VpQryq42XHItJdxrJGz(qQi$mO2t0<(XOS;XcR`WxAUQgax4isnBc^c%mc@0sYc7*L9=5H3gk{SnG z6~^&?i9OC1n!hg;a%s1h5ET57{sa6Axc=Ul^RKQ8G|7dHLF7GV~kx!g&0I5_sH_SBLNuSkQl1_~B zjgfz*z1G&A#iq+X|4~YkPT|H&YC{~N%sZky2oc@ete+U|tMs4RI}4CyaN>h>9kFo9 z67fwKXv%Ag<|1B*iELMdEZvu;L61un=CF%v=d4)HPIFkQYfE-eWl(#gKcoWv1~ zfH=M-vn8)2;rN%s$s|h5mt~n-Y<> z4gpg;%}6AlS%wvemJz!VO#+Uh;%_)x#E4?(uHJ;7xkl(%gRyZ<`7`}j#K=~YA9;Vu zvm1*_0U~@6uGT&{oXcz~5T=iJ^)GZg&Tuha2B`wJNXc^L{yyB5<(y;57;uH1hVD3L z8V?bI25HTwcz<`n-zkW)({SfT-=Fuv>B$eZCQO_@{YA5dk#|WSlitiyFDd>;)|`(? zZ@*toKnDFE=G-Jdz%%!BV%~}TedU?s{G1Rg)Ahb{K0TXINTm=BM0T&?@svqq9@41G zGabnYl)8NJ_m_2xWet;j6*T1lO88BE91fOUnxJ3@J7A?bo{dSEDmyRK**kozDr7Wl%5SRL+-WE7q#5@fNgcmI53v%^T7`K!WtDVE?|#d}0v56tbCuo)yrT zvtnJyoYTR=>|llfKe5_h_8}1Z2kt(eRr&Th zMPKwuSFTQ>^ZKJy5)qu>oN2=5r_eKd@&{e7ItRGuuET|bC~br2u5-h}R2L+sc%o|x zrv783g{w-QP^k9&P<6+Ph%Z?AUw0wF&0t<&RSb{m8B! zbnoQij&<`Ohb33OV8D7HTz$dWe0xKX;+#*%hWLWSzDF7}Q$Ye%N9hh&O!n`^o`DWt zHx0LU4o+P*C)PFt2&->z$@;{o#WG+My>8OHD-TukZzz4_c_Id*{(P>g)OCm5KqOd&P{ z{0&xo0J+eue4MygxdXn^Bes3akFx#$9K!@g4FfEpTRZ+R_oJ0IlbJI8Su^o-9mk>f zTVoRB)dJnXER8*skcF#u$>X;X|vIMg22?7PAfO&6+Z1yvR zOE_SHSq``|rGN$rDI+IsF0RiVN!d_-yzJ+N6Lpm!I`w(HZsR2dw13k%+M z{bAIu8_h^%GYnj`S*ROZnPgvZAI4F&I@OIgLeN%=b~W;miI4L zZ-!hDTz4n6*Oo7rx=>hA3s^Ok_;UX2`+Ioscq7P|>PXPK>{DFFDovab6x`bnbX(q@u#wk>xK z1ev0)*;Dm6aLV~uSRXL>1S}cTU%rPU|jcDZ@>YwM4-{j|w?M`()I-SupXRhEDHd&{YrC0!_jD{Ssh zkyb4;E$784jmHv@b!q;+9pr>i?`l@?BE2m8&GXmlX{zG_x$}EvJ?8P#P zpOK5O*kj(^b~1#Y^)&xP1}Z@YY49%f{S`Hlpxz7^+B=c$RDYA();wp0&EIw}5^Plp zw~7pG2EXXnH3IF#{VQ*v+gnkU(l4<4qgwT#5ALx1fl0HfES))fiPub_^07=ybl6Jf z+sWlDDQ8--oF5`oqT+WXUGeYa#(qtC$5a3Oqu@y|6|4+8Z`-o?;bOs&cNOLNE5$SV z;pgs;T@d@{XKBoUEn(X06k%}pKmEUdKQ|ji{_OfqCOZg87Nh&8QqJHbiVup8c*RC3 z(KS(-4QyUzRJK2YdCm1e4_Yn*McWDpjq#igUj^ko{a z;~>y60V8P_5tUeQwjuu06Q!|*H*o{7$Is`ykh#AHS`{H5KId<-N2G2mK2in!e8WxD z|LL3FDwYNcb06f6PW)|49!@ZGnM`s=ao2!8O5?YxAJ7%Hew8&$?ca{MswNVvzzo z`?nd{gx^K=(uPICM2{Sn-(Y1AC`?q%1Wl!ido(yXT6r?Z?{mA%`{|^mtF%`!TaV-w zSg|~CC(1}IS6Q<-Y0S+G2w1&IdC(zxn+M_qf35@f{O*IV;PPdv7}gWc!F$OmPGT}g zB!7?cO-M5ILmaFWu*x#$l0Ajh=nn!2f%;9DXn$)Me~`mbrF+))X5*AU-+UXhI2I}) zA}DkX>>^+E=*Wm;&yJHee5rctfx@jG9XSkk>u+JVP88h1)UXpU!5sd9o2QzC)l$GwpdUy|QLG1Nk;31`6PfaPCuMre`Pr zVBu6Q;5a)q6y5Ap#bTvfL(7y!r+B=mGxgB32X?Y(-m2WD$*Me9bQ+*WxKRM z5waNg9qB4rby;-#gzELhhx;oB`0pWjjqS3{B?e6eWUf-+#8wbZpQQqF)(_#dioEdp z_ubaa`B`f-!;0l%NL1h}FwT7(*sj`514~x}Zs}Mw#nQT$MV``OA)tGBTfHnS_5Qw& z+5p>o;cKiK_S=8fX37=3D5U*~(dK_Di0{IFG=L%A_O_Tz_OA`L(khyUezgXM;SWHY z$)oeN6(x#wysiSj-z#mg$?&>1=P>GxDvHRHN!(&w7zZ+^6TrB-?+*KiP61wW0;!`l zgk{hnWb2^ukPGNjDOzm4Sv=$dTnc-!P7P8wAYJ6r=Kw_j;NFC63TGdQ=h>LX;(B*@ zEK*utHc(amw|+MJfD91StyuoYu7IIO7Z!_upOwR9|I1^R=&_ij0KRvU00^E`Z8<$R zvh4&&2yy90i!bmHh#J8c`s0C?3_4|5AfEus0&;q1S*{HF?DJ&L7AcF|WU}ovrASP2(9GbsQ)s6Yxk%<}jVZCTTpQkUdK+2H#lx z58=Z6Q={Aj4ZS4c6QUm4Jz8@8HoCUo&p;|6>z82l4xeGSju0l_j|G;qU|FO632aCKZ zMK{HouqNfdMSLQ`&8Wvjl|auF^mCt=VXeSS*Q`{Y?79=6WeXQD^b#7JWIEuzrvIv3 z&guuMv9Ra@(do!v#S6glQ4OZ^j5_v^%S49l(21R;Oev34BMU7x~9`cUnwGal` zPy(jGjlNh+(VRpeR}{pwj^nrwsEYDaf7789DJ|Rgu=&~%UB@vo!ekV3@aX-sl?vm53h1Ne#EC$_ zR7D}N6XppRr+z}I_!4neS1IYGB(g7F6zr)4eE_iNe`yPyS$7T)36BQ&jG@$nToP;B5tM#$A$IBO z79zjCw!P*G|WW{}U$OX3GoW2*1PVN0eD6{B;=? z$Rx>QQ%Wr%w7aj>2T2YPb&CoU6~;3B+6IrQq?!LtCFyY~whyPgx3pSaMZX8SK2v~n zeb@!uumqz)5^~@mK7{wFz8&{(2Z!$O${3nB(IWvSX(*ZFR zUsG4-z?kVTi8S-C>25fP>;YVV1I%ZEY3PMjHq75GaEfj|v{r!o(|De-6Y;+e8^XCH zV7>zsgTlvSYqnk?OYtP)(U~Q z#b7?dCQj21`fd3%<09#?JvSH?I+$z-r5Z=l z-VJ#~g70s6^~sn|gI_z4)`}!y6cqnGf_Km=?qS)nNiNdHgrI53zT!ryODaSEYko#K z%8IssS7*@y<`RhiYb)PSjQ~(F{9`oYuT{dLZ-Ej)QBaa0_g9H{Ba@{&j5yp~g5`2~ zV3)Vnd_(Y9@y`OM`qsY6rEsGCBJ>15@EU241+!*0PIirf1IIkq(9|zTtGBNlytpl= z7mJ2w>FAAIjk?7_ z@oT{$2NtGz42^`(5l08uKz^{nJPy>!5$+IK*U%6GdN!LuRql#GwfltPRH}Lw2+^9& z`+bzm&s7tx1$w*dYs`USMZo*{dSq)TmG~Z*@*hTw`A-Go9?VrU{M1D!ycg-l_hgKy z@w{<41xBT-FZ4?i%Yc*2Ftz^jh;pq9_umRYxqNL_v^(CkzybG<)a|JhC%Cyx3QE)n zi@x~P2k!+@z%*(IZ9y|ef`IdsqrAX-LWb#0@uh*k!NC@2&#Ab5ud1#5HQo1KWl$Ot&E<-lOCWe{y)pZGmsNdSn^?A zd4(fCCp`;Tj5M=3IF&OPFw?i4pU+@ZNC6{Frql&$^s@?xr*A=lOf6^sPfn+NecEhm zWa2j&+2{T6m?h@=o{6st^@Mx^fYQk#Sra%YP$Pn;iC-sLf>Mz7|G&z1IoF_0F_bnL zWO(e-