From 681050944a715131456dec52fc3d77f6b779f8a6 Mon Sep 17 00:00:00 2001
From: Georges Kunstler <Georges.Kunstler@gmail.com>
Date: Fri, 20 Dec 2013 17:14:11 +1100
Subject: [PATCH] survival and growth model with lme4
---
R/analysis.tar.gz | Bin 0 -> 20770 bytes
R/analysis/_region_.tex | 5 +
R/analysis/glmer.run.R | 192 ++++++++++++++++++
R/analysis/lmer.output.figs.R | 12 +-
R/analysis/lmer.run.R | 23 +--
.../model.glmer/model.glmer.LOGLIN.AD.Tf.R | 4 +
.../model.glmer/model.glmer.LOGLIN.E.Tf.R | 4 +
.../model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R | 4 +
.../model.glmer/model.glmer.LOGLIN.ER.Tf.R | 4 +
.../model.glmer/model.glmer.LOGLIN.R.Tf.R | 4 +
.../model.glmer.LOGLIN.nocomp.Tf.R | 4 +
.../model.glmer.LOGLIN.simplecomp.Tf.R | 7 +
.../model.lmer/model.lmer.LOGLIN.ER.AD.Tf.R | 5 +
.../model.lmer/model.lmer.LOGLIN.ER.Tf.R | 4 +-
R/analysis/run.local.R | 109 ++++------
R/process.data/test.tree.CWM.R | 26 ++-
launch_all_glmer.bash | 22 ++
launch_all_lmer.bash | 16 +-
18 files changed, 344 insertions(+), 101 deletions(-)
create mode 100644 R/analysis.tar.gz
create mode 100644 R/analysis/_region_.tex
create mode 100644 R/analysis/glmer.run.R
create mode 100644 R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R
create mode 100644 R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R
create mode 100644 R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R
create mode 100644 R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R
create mode 100644 R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R
create mode 100644 R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R
create mode 100644 R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R
create mode 100644 R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.R
create mode 100644 launch_all_glmer.bash
diff --git a/R/analysis.tar.gz b/R/analysis.tar.gz
new file mode 100644
index 0000000000000000000000000000000000000000..77bf3dabd28d0d9ede0c9c015b591c20137d630a
GIT binary patch
literal 20770
zcmV)8K*qlxiwFQi*t1dq1MEF%bK6F;{)}HSft(c(gW#dVn>cVSWl^@?%94swHv5W?
z<$_3%LP7!z0P0A}`R%Wd8C(P@QIhT5ieZ-s40`&So}QkW2H*zn<Tm!>`m;}WG~jso
zatD8#FLxUGUv)fdwpu&eFLzqoJFRDp=2ok@^-S!1LY~rN9w%-j#Is-LL7YsyXc^pH
z_y0=|SGE5y?l`u?c`}<PcK2y%U<2>&ZZB#7=GN{mwf{C)$yQ?<;F>$T*#3>DWuyA|
zi?@Fq&ZB{6=-qmiCDcaq!0zg%rfrCLKApOeAA7Ojn_?73!ktV6RuI#0=uKj6<WB(L
z$C{Y<abnYzp)-7aGW8<e5|KBUN3nn9Ip^J<I+hTl`D9`T?i4D@qM4f{UK9vN=+Q6+
z`hWqEMi+iCM645I;?8E1TZ1TA8FO~tjuZD{;+dv)9r=l858^8WZdd~AAZp1L)~H8b
zJf9>n0s)J^&BFvi4%$a&z)loiEB?3$6rwot+#z8uAuZyIlgJOoaXUV62T2&^@oWgt
z%GeXKg1GjROEH-EQx`ZGgi}D~0f-xja3nl;a4Eb&7<pqq478ZU0$3ZHA<qrTXoDIr
zQ|-^9aNxzUHw1okfVZ)B%+--ga}}EfO3(cv5}Ubk;u*y|v6ooI65w_KreSI*gZ;jJ
zj#U$GydtDJK_Kh!c)#;+r%cv5`03a|L{{M1(Uhzu#e8}n%%^K%9=eIkP*MyeL4~6T
z?7{Gd4%vmFgflJp<ESgRI8_&?WW)v4;Wa7@e^WIn!y>C7dtG3{y0~?8f5_dUhqbBm
zIgWmG^7iQP*lzFZ!s<VdC*jPqlgRZ6<Gq3FElo>8do~;y(EkSsn3_GCgb9NXC*0eO
zH%Eu>?SHkqjt<g-1!)2fedK2yss=X@6Bym4=Z`NFu&HLl0y7_iojc<1=0j+mE~9n&
z6tudXdkEI;w4O?~61{5SnBJZ8Y1-IJ?4ft%yTO3WRgpH(jw}JDF4!c%S)CoV0a<|l
zmpFkkq3q)Lr=ugdsmgQaV*Www`G}pcK9>j`iTQ8TN1lVnxrni}VehCiuz$z#hbCbd
z+DICZ<E@^rU?fQ=#~8RQC6nHlb*IZKE~(QLHLO=>y4pvnt^LC{4^h-{`_pR2je$*Y
z1OC9cLrf4r?^A&vNLW=0qu)myWFj~CH5Mr8l{=ZkNQQ$q4JXW$8w|xH9Jr_eON?Q3
zj^X;+ixaLO4Bz%_ae)Oy))bdug6VywNE-?LCu$S<GHOyeQa3$yNrXv-bk-MQL(F4O
z1l|}la^(p>AOs=~MRMsTN;RNPkULIQ5#!K|3P){_UTn}adsRAP3G?cykl7Nfvy#4G
zgo2BF_cs4P$w<6;fB~0^ispg;+uSpx7EG$kSRF>Ks7eb`l<1dD25e@v6dX%cvcHsZ
zWgsT;4F|fzGKmW7S;t|KI(W>0y8|Q7nL2-$JMVTdAWAP<%KI|cdS2fO`z^R`!d(Mv
zvqTdVd^>S3lmdk`TKW^DtW3OM4DlFWn;8V^z!(?}OM<>Ye0be}-{f|NVddsz(-M>9
z7C9Qa@g+rZgOD^oj=WpN86G!^yCO!(H2zrP#=Sw7bRR0-=FTTBo;EG;2Z1A_Uw};<
zJ<v@>fW6a=A48Zm^dvkDMn6b`i0Fc)9>?>Gv<*<xYD|eu*1<IKZW06a2LUMHb6QZ_
z4*>7BDQVFkun)t>j+j!BYF39oGO*3i2^3qc1V=hTBU^3W#Xb`E{66;};S_C>ese}{
zl!%)P#BY)KZ*w<DAV@cEV7kN33!BVtfI*Q}HSF(qDk{%HKS*K&9Q<vs1uV@5IP^Dn
zO|bwPd*i8>0yhDePPY1H8M09IqX2Qe1$Cm!tk5J01Jf_#Q)-iFo2IC%PLy+Sy~<Hx
zVNd|egOIsjsnN;p7+lT3@S#7qEdDF}7uk>Okhky;J7vD(nfF&8xxO!{L&*Rra^t^3
zlN#$2%PEMd8ZCGjDh&_-i$}#MAP0#D3<~)H0qxPshJl9F1Z`+*@b}5t<EAd~*lILP
z^FK<BnUe&zmUERKZky;BLEA#c4r-kRtKbZus>RDf2q^?hryqIl$UeriGtYy0fg8to
z@pBjsfwds^6C79d_wJ3&^DShb2ladaQZ-?EXy}o14ooNtuXC{j3&9;PfWAY6H3t!H
zWK2f@O@yhM#UQ}*)QzaPWworPg~32Lb*W_A!hcOunAHO<;>q!}#CBN2MgPO@DOHAw
z3#y~4sus=>1z+krPBUM`hDo7=o<l?d#Fh#mnrfiSxs~*Eg(MY%<eN1=EtQeTUEZM>
zAR@|4m<8$bvH`9dPP*<P+sejyyhuwtN-CMIuhC>Emc+MGtdH7983M|Mlog#GorPLC
zDXzlSkp29^iPLOl%O<-f&34rCSZOASMXWq5W6vmyV~@s3{nU-zK{&^V+`yp``ZP}{
z<s{tfuZMmG9_^{u|I<}}T%5NfT<(8r5_jzXw3=IumfZi_-g(*D#r>aVWBcp=&*wNc
zmLJ+3_s~DQIX^r(7Uw6T3%=s$eWxorXXl6S+gPkk{ENtqZVfQ4ZBwqbrJEPrd$9}U
z$7);-7Z+~q*+U=4EZwe`w=e1(LF2TNm(dltg_ZTX4!<)8ypaV{P1=QmdvO)4gfq=v
zt=Y#26V0E!pJFZ_Q!Fq_ptQ$jjZ{rF3_=&r{CO5~6l05gXkn${NG=*NEu8<;hL<0V
zEn(X>36DcGdkq+-b3p+K03OeuO*~t&Cj|)uxw>yuvn?kL)Bt&+&LRh?s_vk)xX#5J
zf=kF231x;N*Nzq~2O9y#pS)ItF-4W+i17%NfElK(TfwO@wM53%5~eCt!AoVUFjwo~
zQ>j2?QqU%nn+z^(C^FJwQz?oTH8YHvo>I^}zXkj7R+&hqmJR1A4j4PT^FZMEWSGJL
zrf2TpS9k0Q@+%|^s5=wIqD#pyo{|_(#ZAt%Y#KA0sOUjq2c2Q?Ztm2Uael|{>ej;4
zS=_sR_mTBGUE0xKkfXJx96^Y?(~o4Pz5hrAu|EY5BsqaWW<jj;(F0ltDm}K55QY8N
zh9W*H+6u}a*j$g63|j2bIC709CDb(ZItfg4Y=6W9i#>p0a=ti~2R|2i$vg@SJ&8Qe
zhCo?->j&Qo-iap=@;``ifa5ojc7vfEyH~|6QSco=%e96h0<th+0+UDBm!3gnNDln$
zHntOg>KWuD37t%V!9*8yOcIfkBBZSriM88#f0U(C00u7@E~1@9evlY=tHzXbrJUYE
zk(46M6J(UAyMLCyVWFKkS^kHT38D`*Z=#%<bQs6Fd98@$em`}wi;Yjm)DLW!-Z;1`
z9lO}b(AK|u&T)HLjr`DI3Le5D@3*-hc`|Oq5}Box=L<J0arK20HHf$x^(LV^lp`45
z*>J=uXn}_cQj%OVVa>?o2+SCQ(2(@Oq@JTm6;!Z9itNZ#v~+}Ip15ICgX2AIpgxE^
z7XsC+)ynBo@nv}Zvl|7NspDieB?d0ct0Os5lBja>nA4h0<seAD`%>7P@~vB$2s%lg
zhh<gh{}-f_eQD&W-m#y}B5&qK5NK1TOej^lmL;t+uYzr@fna%yhz-F$MYv-BfHAIc
zMwO1rkeZN+q2@UGl^5tKsp}Y6*DYL9OF?tY>nfTSj1IR^t2ihqpG+I0BMx!6$R8P1
zHFL{&)t)hRVr=_#PZrtI(=IwniSr}wc|3x^a>U@Jo+fg`;7u!a^9r}O%;POOs>KvL
znc#>ua}^0)Aw-se_@iTF2#pUE_em&Zw_}sTj*k31lEUgfT^LLs+`lWkwJ6So7I!xn
zkNmp9l2rrfFdG}>s7Oofbx|TR@FZN!>Esqp*9}UIs~0NYokyuUN6|ADkW?5<Y50(_
zP3%Dom1L+??6nXv5o*i5i0wB&zqgYh8Cf$g8rT7@JXm|}^OJK+wZZ+pDb~o)I1T&z
z24UH2+i#uP4qUi5e9_!#fbuuc$4<j5oSJ!_bct;(WsvtkFK|+uG0g>z@r@tbft%gv
z^>jq7_xq+P{&qj6tB~w72D6g8-@PcjAAc2my;i!A*+`|^KPv6qm-{l;ufMmyFYe=`
z3~Uo3o_Pb`lgY$E2yGR|U^;4Pswh;3xk5vEs{vIl;dZeeSrzKm!o_dNR*^4IkUyPI
z1+Ir9`)g}Do`D`H_OE1$=L3`jZj7iS8Vj{?@$W8(*RpT)&7362w0jX=p%fzP1`m-V
zs}0h;F7+$0o<-&6Yd{<?r$Z?~Hv92{M2)`*{B@%txYnV4<rY2nO*W@8ZfcmY3`J%=
z9^wyey@P3uJxukBeAHH)4KGLs*QSz7&rWZru1y6FOMDXh9=5NdgM$C5Ga>ge3CC|$
z__abV1K_IM#0!ND1Okx#hoE9-p-Cf_Qh`ttH;z}>kCxNH3TQQ%r7E<|wUQ98XXg=M
z9jI+Bu#hcuLogs@oOvri*Ri;`71RhCOfk=c^U+>lLpv1RLdgkyx7f|O_Fhn`$90z1
z7%kLTiJ!O|`lFFdmqCg(i3cGTu&LVl4#J81XiX0ODGqpO8CbkuEY2Fj$xe$Pvh^Gl
zC>L&)TGJ_AmxvWab&8kEK%s@4!c`IQeB_jrT!CmWaOCge&5^u0s<_Fw_JVBHWraTo
z!mlpjt4sJxxP-q_um9xJSQ*NC*Jz_WciZ4p{A-vNi3<<&^@q5n<qBS2ri2nVdWzj2
zZDGoWR#oM_c<z>41}#kEHf8H?;sSy{GQN>%IkbFTk(idxOXGVhHjrT&6DQ%!&;zI>
z5rZ}_XrB)%3&tu^oK%>{CqUOC!zgqNcxBiZb2zPW|IMCZlM8}StXrEX=cc7YZ;_ku
zHlWCRiiLXe=CZmKHAcAZn**T^jZiwB*@!^P#n0vXcuOxX$E3Ilk^8DC_PWz*^0NEw
zb7+LavvYB9(iQ*HetQNsAPkk-N3#(AyKD?FUpH~s!*P^4hvf3^@cW+QD|Y2|S^BZH
zCQDfaum)G?bfn8pY4HNulV?QzwT^Vd4nrM#^ab+28aq2p%KvILwzsyL^!|sP#?H&H
z`Cp&oP&rx=w!bMhE|X;TL%n`|eQi$%H8P4Y8rPwp(RsSZrjxV8VJt{~#HZmRr}0`G
zfAE8$fAy+dD#yKratw)<-Nfx-5AVJL;&7{%qUo#X>w`ahC2!?1^+o=QfFRy;Cq2nt
zAA`ELuNKSpE6Uo<s=^kaoVrQm->7mcW4KT8$p^lErEOe!7@qg~JSZ<M2Rj^NIk+{5
ze_*tQSxEc`5tnMIHLdv!#u%sh8U}T^J@95w4{N_F!Zj3N<8^PVPwx--uvdU-*w_sv
z*bYI({&?!H24k4JA*em>0p?C-moACSNe*C8tI`f1bSy<E#!+~U3z0;b0?f1X_OV!v
zP1AT2L9Jy*FhvZff|^X_Rt20spIR&~{PXt1Ub7eU#h+prgwfRCcELxpr0(IAHwpGy
z6>zN#uBG6%D&V#<xUD`um;wiBIuF~ha!p{W%zt70B6&~EneVC1J(5$}a4)ca>7zF&
z86!Wpqn;@EC^eSp{mB2-oA{RydRFOLDbZI-r2d5=GhsY@)ANbw5T;YJC>+iQiJ_&F
zHKX|=*ZX?Dl^3O-?zk4j$R=f&7+<I$!Od4JjSNKG3Cumb@WYsk`5bv-EoYtw@b%T!
z)eCI=;;rw{>7Y0nNUH**-7LavRl&6KFsvkX#U@I9vR87ASY3j>V(DmAMHvZ8#+IUV
zz)f3f(|dKqLKX4O_h9H5T>6;Ek$F6k16Op&p$l&`0*|Sc?6lzevyPUR`Ikkj$8u8{
z=lZggL^%$F+U32(*${&;z@aqZHS}5Jhmj_2idbu)n)!?-FX6S(v>Tn;wpDAH+?CWm
zPDTZaOMS2`al_1%i8P7&B}{{&3=O|$244jb50#oHtmdsuv(m{-`c)KIpC3Wa4xkA~
z*f(-72trT)QPX8K0UZ@Zq*wOPl%9+I0EDx1PgNC`3TYpg9(LDDBI`~x4`5o|#Az1}
zmQ(k}X=X_lI&-eTwipbbXe<Lt?^%SDBN)NjJ7wr1@f_^-0$y?wYrxYg-g0o6NO0{&
z)A>v|;sRz#H9ufo%Je2YM<+>*@Bv>U<Irt-oeY`<C9J7Np%QEbhxA<CD8!}@4dApe
zil;iv14JFFU~uXBL5$A<X&Z_k%3SnrpveXjH|9Bk5ysBI;7z+2XRdSHdyit%(kHVI
zA9U05gTV($;Q(gT`~Vyd1pYL|(=uVm%V=%MENBSSGLW8S)u?jc!XL+=iy9aJU?K)s
zZj>@D%29Zp^}C)#1dZ^eBl(9N>HLo}0dZ~R=_i19y#KAy+-;@r|7&c~{12|a&i_8g
z!Q<E2o9^N1IX=g8-Z?uLd!2)mZU-;VI{U}s<XF5p`swVv)72_o7neOaMvw6#y7QC$
z6AY!hg!1Q;?mrPnTm0PqUo3d@uJh)f;<Ve@Kji0!+Hc#3$3XXF@2u1PS7)Co0|VXm
z;koD>9CY5CQ-+4jR$U{No3iC;bZ|D$j@p~};E%I;5aIKpn^_`{v-$3DbJ=3;=O;b(
zg8739b58vik-3nK|7BRn9!vNl<ljAs@*iYvoR<IjNvBT)e16j92MN9a3x($@zDT=2
znf|pkbGh$9>H*o^YdkzqaV*<!D<*hn_$X<P_X}m_Cq7P*Pw1_)X!Uo8br83Xxa-w6
zST8(ZftcAi9*UW3FJkR<AA_~iUBuYle+<U<ehDM11g!G$It3&p<Y|X{oLUzVvsxd4
zn0#;vDQor-NU3)(aQbMp%5Rh~YnmH`SFXVgV6mz!){MnOt+aKmRFhV!J1gbP%5iDp
z0(=bp^2W|O+sIE6ZzxHbF95(#uSAVc{xO?9|F`(|yKe#d#v7^IDYhWYVPbX^CV3@l
zd{V|FZ>N|#-5e$sM`4mzqQ)m>O!Bs-rSD*I=v6EDgciNqO}f?RlrLdrtXl)`#=&=;
zXqc)H8;u{DI~M-=t{`l0lPYfZEu5@wM&59<0yAZB6qQ`YZ`7I%^9OLR>g2nE07M7T
z+(l{hG?Yg7kSR@Uc^Ex5KqZ+%R-~$x=z`^DHnE;yXz9Jp^HG@qC(yn<df(NmW>DkL
zHHAVF<0bLMI2T}yKUrkCHuWbBF4w01N#e`3xp?vr^C33yWDW79S3JHR`J=duBOi~$
z_f9|IGT`d<zg@Vsqu&3s)!KRatkK-vZG6rD_$<f9H}!cO)i3;@?gdw(8xJCXmT1w`
zz$I{J5*wJ_ykym{^pS`RY@Ya&xK0VLa!BK6?qsm3(QA3sR8xho^xXPpE{3%U3<@wf
z-DEVFj>w;42uJksiQFC-mD^On+0pc7%i^1nQ_5Dlp+rT%nJJ$bQrG-|s(keRIgax~
zdk}h~k=-#dahl1c88PdV%P8@y9fn@wg5EF(Aj2=;ko0uew-~#&f1uu@#+&K%z^<Xw
zCpL7eiUK{X>n66Zl3`KvR!-|q?i@&D6bhd`$Q9$ml(|dweS<V1Gc_FkAcjBS%-{c_
z?_OlG;Kplig3yEsf(4cZr#qJT-V)pRB5Co0f=z(E2yZAXUH*-XKNr4XaT!Mb?_rR*
z6Mk2~#EqR6hOjoitr_Cjzuyn_>BCTe{lW)E?o;}li~52Ky;mO{opOD}Zb!hm9<l0U
zAB3HgyF>03m~Kpa?3@nZ2ZC~sFTyCC%LklcRzpzLk2P2Hkz5y{Y-YwSdrG;ZsP|7u
zDuoWB%CKvw^U6*!92xwj3!*m&uP{j%<)BwK5=@p4q>|*b9Qjg4d_HBqkR4A?Ad3}!
z@a57zDB17zSTdip$ZguC$x_k#{Xa83+@r%-Hxu{5VU+rN%)C(TiPBR(nr``wUf}y$
zR&s<1(@)9JDe{6}X8Qm3zO=n<97%Y;M!y2hWJY8R(IO?;lb&YpI*HTid2`I;xMyY`
ztu_xW(Kgevq#@<R?%woof2#@r30}HvD;GIE4g~^*LR~0a1;VC)$VdD$3f|k((^FFz
zYLaD??|C5R67#~{=R>m@j+(qlLe;beN-+BTI~Bf1Bv|511mp+tfdW!+ut=g0vL<Mo
zo9SF+B3G~CQUe_^`P%<(bU7P`Lm#JnG-!BcS@{~cedP={s!1Iq@x<U`$^98H(!acz
z#0tK<pf?o(byZasS)`=r4gq0{{T_+yv;0`=W1&Pp%n^^P_KNo5#x=-*QlN%fyqU-Z
zyr&Y|w1|~WyJ=Yi>#vgV<TFr6sZ1>^Lj0@ToR3gGLPI>&#)R%}DW6=8QErY?t~4jO
z4Wp6GqbH{>hFj`mhziuNw2DPmSk;EY2IO#?9IfKDmt5L~s|FZMs-XffX@v`Z>~8we
zrm9BDq0ORE=?GT?Kt!ItX2;JD*};?JUk+ZdgI5PH{`<|(Z`kW+lzepf=9d@8Z*o;J
zO)J0Ce#pCJl=2m3Vu0<Ld4Xcm5hdXQOS3SGl*^jTYvj^1vo0X519m@GZNAQ?3?L@&
zOS6Y~;ae2?$<e@!mL;4@vWas(!3a7DdE`e<a*K_i(q^O*<$0N(xYx2gFhXAPG6j|5
zE-#DID4C*|uu@6lr{`7Bx2c*tnSPy36(ZKe!*+PTLK$-{3{(Ak10^Mun+5Fe_X}b(
z7$AbnGp(~4y3$kDo<TzfV_MRpR|j<3OQ%aB1*J5pkde=a9{E;OP$JU~qQW9Fus^f4
zBqu3suS};!WLQ?d*yaUnUmc~-MU;<RJHk83#Pwf`TDDus1`F9+`Spjk#q?T0APU)Y
zRh1Nxf|gIEEfvzE)7-fQ(0K)-T!EJJLZ$$0r7ll_OXv@3rk*#YG2~4GB;7zU&!UdQ
z3!t>2Snp5O4GT74v9L=8btvG0nVSkiv+iKOs-~(A_D)RbXJ<gaH1qjX>1|?_7BsHP
zf-2)y1h&A$Wixx(pH`W@$r+q6fI7%s%$c(^vdv<>QeKu7iV5M*Qp0jOlE>Tn5rf)f
zTCR&x6VjSWzVM>YvX2#@2Ns>{$I5j?C)1O-vdDK{K4&;U6sh$TX#zSF+^83@$>qgK
z5aG-fWDBDAbc-kqPnG+nH;H7XyjN7&yyUG^PynqX_af<smd>>(f#H=SM14$X0Bl1{
zTO7H=JFI0WYk_9oSuF>2Zywf?;lXCrBp`Gxl)#Jb!KQU&A!Vz|PnZl5bZ*{8i3Zmg
z$Tc!^ofy0Xc}wiqVSWti60pM~?b2HYS`>uya}f4JjFzH*YZSdARA)tVKeokrZ=AAA
z5a`2?aO8R(Xf-@mgW!?oi|$s-(61=yd2KQ+RJ3g7rZ`<<S?Vw!;ny##F4xf?9yRII
zZTn-QlR`tq{d}XmY&C~Xp1qAWJ69)XKMS&?vSbAoUj$pCU^b2BI#2KjW3u;z7F;kz
z%-4!urO=OzT)d$p+ncECu*4i_iI2lsm6Nw@c8cdI$#XMnxcN<SUKf<az-MWPw;GEU
z&}^-%P1k0p%X#zV!7=^)JN%^vHm!lx|8O1-&+Tj{@<2#cX*~m+5@)DFa=~?7B(Zz4
z@W{FtyhupaCZLj5WDt_CeZj;*3)`k@X|Q5vL)v2Td6Fr?Ei)0mDrpg#=t9-RJT;li
zYwwcz_6E}?%2l{L$T|G{<uh^N{Je`*$*pwMH&E2op;)Ve5riF-;Gq3htcRMw3t0zC
z%%5VYO0SQdUaxxBuu6BP;P&7S-MrZ2T_MLqF)OdNXUBagi#$wm+D%?t@?#6{*h5?<
z+BT8ZEeGjQQ9THtlwQqRAfnblXwA^#(39KZyeu6EuGAo1I|R(4BB1T3e0b;tS6lM^
z3HcE#^(jK73O!36`tU~JhdI{-Pf-BIA85HbZlGC_?+dC)y`chfrK^L%)yy=7+3tsY
zV727GyRNm1&i#V#7B8u#&;RnhgByaBBdVjUFFzLX>vLzowIXU>Dhnth2dQC{pPRh2
z3|_~knam2oNEYq#A9)sBvW8~~2ZDB1zj}R4o2q|)^;5E?f^Os>CQ`^Vx<|~9ejQCe
zOqkgDl-y|7oI!m-1R$j}*i_iY<mu3+@+J3wl_4ueTw)LgR=y<5gH(M8ny)&`1sSb|
z?-!PlZFvx$%}v}~V9MCkWmw55)mnJM44x|-$P!`WQjBn~Fl{^|YH4)J$r@G*Iw1tZ
z_bu-GxNo4#_qYS*5)3C8puFKpy$T&N37@!sUhmEfycm#G6tNoWVwN<qkc$uf88TGu
zm3W*uqcQoO6dDUR?%^f@d521dewN<%SIW*(dGTm3%yJ=I0g&8Ib1SeiJ92IXqP~&m
zHt4o=7g1*CkXU9v)k2i2ogqnxQlFiL?}LdRN=l8&F~6u*zqLWT%CFQ;s1dC|o9-bJ
zGc?KU{)$tG<S$^%5Ddz}<$nAGif3R-*=1g(4{iC70Z4c}ldqjsi_g!YS?AL+4d*lm
zm2k6u$JjKl4#n9OibL?vpO2rjr$0Y?c6fC7>X;oL9lm+~`o&ZB9DX;{wl;k)RwJhU
zJ%+Yv#b}f}Lr9k&6>+k2dbyc9Jip-O30GA@1m~wf{)$-DS>qXH6ojnAK~6=cl2Dmb
z+q{t#mTU>!Gg30p_PjP9R6x%1>xU{v@`1W7;BK`=n~txQrdFB~9By4{v=TC|xVXW4
zn~Jp!%>wr_9M}C(ty-T@wSN2rswD=8N3Ph^D=?>5Igy+k-@uc-I37v1j(0e7>Az&_
zIO-Q2ei8?JgV?^ZuQM8z(;ZE#G+le|NA?x{z2@SA=h10qSH%K7F4tMFQhmMbPjTIZ
zsnyr-6FLLhvz<-6hD3EQcxltcWwo0(9!QRJ#C$e9C#KW>b=<k0{+J@bx%jblT7{>k
zM9pHq;~O9=CO#bQ>{)#ng-<Sl9^o)c(PkpY%>gdsSUb)@teI*+9NMKPTMO{<O^0vy
zHmN{^A}bIS#0mtp#RbFSy{~Q>$wJ2&@xj1SAVppYa!vUWIt<PBWYc8|p1+s4sua{e
ziQ)pWlP4uM6TDBLZ`#R@b+k(sBZnBn^rX{7UwWhe$jX|WWtX(Lk8xCK$(XveU^2%u
z8Pii~xKDz4fNOx~{)9~~$74CEBae<{0cU#ysNsip;y;L!WEF_0OgK*q{g$=6Hu(Ir
zQg!(Us}~nG$Vi}?n9vG=9YoPIvQKiCzG>A+V@tVte1TIYkXwpF9GqPMePfd?*Jgq+
zzzr(7<e08BM)1a)`N%#*jP!#_dh1rNW@qhcNO%^vsFf06J8j#!?Fw(&!u82|!xpYx
zVb&NTY$56BYYaE_H3&ppkpI=qp`4Z2mg^<(?K=Y{*Hvm?rRrV+A#F4+CZ4JbV8nNE
z>x}Aqkp`~#m7dC!fLhK#yJdPd4N(M_HzXX0cscO@^kVC^8xv}rD?pQo!#_KYNTl%e
z9~a2c>8Fz;IR>u<>5CVnjlIwt3(<{*r-FR|cAi_Z6J(2mnT*~ovvV9id+NUX1h56f
z9$yAtZm6OQVc7?FK(`<&j}g%T?}$Kwcl3w+O_>ufQSLu2^)d@WU&TVaZ3;I=%kTs1
zvGFq)w$l!!JW<n7Z-Z9x14U*_M;#f5a-sr|R!(b3%;T!jd~h{`WJMcPArZ;lW9r^9
zDlT*m#Yg*;>>gp_4gfpd6a`~oij5+}HHeH7NE+8gri<gYGqszHOYnLYQB<Fsk_!%g
zZnGMG5Sg_R@FW}tPp0E(6q^%oG7jDY$b1=5C<#aatsI)vUM?tNaom7F$u4&}jP)U9
zdf`(_m-Y`icz3xpq3XgE0-@PI_#*Gh5Kw>5jU+0jY$OTsV7XvO|6TPUBtxmvxf>Fu
zYzhAagX)>*UYn2S*Ip<49n4RNhez_I7`vDUS93E(ipAk5@H(7sd}MgpKoL%kC|OVf
zUjfEZV7dh`@4P(h6olN{a!6N{FpOIZLX}Dr?wywvxJd;vxZtBoeA%BF)G4@DsCOKk
z1(T84oKInJOPgs#idQadIC<tddzBp2rd9X}aHBbCirW#~wCaz^!LTx8P}$d28H6v)
z6nut?s(QBDQPoeWqRVy5&~(e-EUR2eaFP|2Y$YjH6KZ!R&7VMIoejlxGoj2*NoHxc
zo_t=5ZGv8gHtS34DsoBE;+9HhrIhCcvrPVQl0GYcW#X5HvMZa!sCKJ}@A&}(hEcz%
zS-rW%yMLdFXtwOh>z98&z)MA*A2N73czJwy^oITOFj--J`tv`2ek#5nJbv@zQ@kjb
z4sv{$ipWj*g~IIJ<W=8x-CHxexp}5it6MP};AuP(<J6Q<*)iUlBApJ>@-^UatwhJ9
z+h8cS!Qh&G(TguoP;@2G{y4m<?nR~1sc^opDw#MRq&HuIODKk<JY&47-9cwPJ5e%8
zRojWr%TAoucA~ZucWx(^T8gzsXP>pDxVf1?Og<YMv1klem=L}r`>@r4+)X9D@*L#B
z)8a`+bv8sBIpJ~ClXyDD%>$zluXV)e^>XS4??ahQCvkwsg)lT}lR}doxyi5w1PB9@
zHu;BeG^Z&QaZtePS*s2zB{rVazXc;9Cbyq?fV6%@0q~^&{7rFy0nec(Cl;EPE;kE-
z)9A`;;+c676Nq1r!waaa)wbJCk7w>VEqTOVxajP}{L4OPP7nULbek$oeg1>y`Q0AC
zz%1~>owmz+PDzecjufd%Rmf7u-nSZB4ARA-v8p6k<!ql$0-f!DQC#nQ!CW6^4hd&j
z4}bS=h4-M=l^m<&SEc4I)6$m*e=lWIx0+4hq<vPbx$?xJ^3~b2GEGa)b+1BXUaBTV
z!skt6TTy@i3#NXFY3%V%mgZk_>s0iBQn{+=1IsjQZ>b^ttSGU8W|8knPPBIE#6pEz
zBVukztld?WQ>4Jl(PD+wd`E;u^UXt374nj$UoHT8$s%z@KY;NCSITipIm$9Dgqz&g
zWo3+wHmE_P(-q{=W%7}d$w%cfsZX@7Bb13;(x7A=$y^3DKOG>|bXhTtqh9%#4W!sj
z1lxAF?dD{9Yu&|NN#ANx;x3Sq3dg;clGUW7TqqtN9KSxUWc%ZZCvP>Z<Qf53*Pr`V
zdDgaF>an^v=3CLh-du63SgYI!=pfe^QWXJwTO4+P!A{2#dlF5$T#`3W#{Tfv=87YL
zwMN0^Y+u#pDi(t$Q?qzo7UyO)?()jS<U+TXPhLJb8(&=5ber3)*0K2hKX}{@-T&9w
z@AcX)q&x6$ufG3p8_#`qJY^@K_|r+ypgTY73p&50CwBwvJ0HyWZP$S+zQ5nA*8g@J
z^}o~W?e|<ze9(Wl)BO+DyKO|uo^M?L%lI!|)$GLKjeCP7_P^WPb8`F#X{qi1R-Q%!
zZxSpyHf;UgxYxKxUSH$hJqq&xZk3Tpg~1{D{(Tls(BV?WAPo~je{9py40OeGkmWv-
zX@rwG9^;+O#W4uOi!LYOsp)Wx=+UADBhRV!(~NjLoyq85c91LWQfYU^($U(b>H5dV
zxnej;ycFs1?848N1dkBZLzG>#`X)+OmhCRB*B{T$Z&Dt;^99ANk5WKVpPnMn`0oQq
z*Bt%%)TXP{FesbAi5rid>JTtFUh^FHctMtyk&x>R;)<O-Vb5xWw6M7+muILX_iy=W
zu&OWYESi3ppX>2{{&w)wo67tWoX9TKF$3>K_-il@&+*2?vXbdZJ?HyZb~qV@?=65c
zo}T@L7Zd!!Mw4lDVb0HkseML!8T?DD@$(b>2-YQ>%%bVYgzToXBd$87BeVUG13Y=Z
zPCOe}V!MO-DnsYkbaz88Cu5DIHrhyy+{l)|1ICY?T-rlG;whG4I|99(2m=XTkQnr7
zV44Ghn-3Nfb{Xr~^Dv0~Xm}nDef{R@B#*pVP}Iu=LL#4_{6tY+X0aN;LulG0CA|<|
zJc<dVrx&XsiK69d(^pDP_IWdt$-OHnedXSZ#BgWWGU4xt<Qav0s2_#yLx!D<e#+0#
z=bj-1B!5m{MO7a~AUUzz641nwAb~>B)>(vTN*r}AzRNclQ!%dyQnwZ#bvC}5PFm;7
zS;>ae2~MYA{|KnhqHr2D0F5?<NU__8tuE$z??)8$nbFZy%s{cAXVVWR85rJjEoxeW
zeqxiz=#sewd)u=bs%_bVl!T=l>Cdv#NM|{5qAhG?xSFM7+sb^63SV<9Eqre<;w=({
zNOOLTTeFhN8458)@{dP?>tn%n2L)bA_;`?qDg=H^GMN&AbRi`IWylgiPyu8j(TN-#
zDNTbR@{|E-d}R08As~*oGM@(w{b_kU1T(CAP|2ru_;Sxu|Nr3Wo%{dYPOsqqcRIEI
zzl~=@|Nr3Wt@!^pcm8$nZ!OP#Gq?Q)KKq6a`;ENy6<zhRvfFUqKL=;M_SU~BZ@qTb
z*~-p3F6IffN?di8_s&~6e|%L(|Dk^po4E?C+L!<Qefio~|6+ahyLeWZ|95XX09fMx
z_u6~=S^wXKL4i~I|J!&rp8t1mYXRVvmH=+L2yoM7fPDE)5ZH7j;ihW|H(5=%`Fg_U
z!<4nhD7b5M{Y8c-x5!XfbO`E!Og)gPmmR*KWrrP>Rp5N_dYR%29NMhCOp%p=dSq0O
zjJ~;N`S|bVn}L_E|G0(qpS|vWJ^tIq^A+y@s&{bJJGkl{T=h_@9!jk|lw$QBt#5Ad
zQ;%-y(aj3}KZfE0;r}fd+0QrJ3bxGuZ|`?<>pyOLuRi~=m1jf$|9QRttKR=r@Bb=y
zL%wid|Arg1avXg~VLsqXX3G5`g%Q9g0(iI1H}1UGUi4?<Mb}<(?InLx&r<*Y)h{nz
z+;9nCssG<`^83GBr{4dym1jf$AEjdT5kM9yKI=LDSVDSnOTBh2?`pn!T3{6yy#Ri>
zV_qn-)E&nXcjkY8a*w-qxa-OFQr4~HpWm6EeE)Nuw_iAOUokO{K6LH+FL(WG_rB=f
z*WUd%^epxNH=Xp?p5^Di-QxN0djHQ>o_hYjp8v1s|8MoFha>g*@6TeZ&z+wCt`C2I
zDPx>^d81z5sC*LtpUko}B}9kY3j>TqUm^al7XN#%*J;=Bf4A{0;6J5^vEW^<|G1s~
zeK*U0x4rLn>+#<<o)vq_vTrQF8ycd3?B_U~{K`&(`G+8wFd;@T41{q(kiAcYIt#pJ
zRrELNz6s$={baNMC&|oQXSe-m*KIqUMy4`(xdoW*OqH_RL+GM*tuF*m;w(sVwkea-
zXuvqg79DN70RtiV?lh=ViGzQel!4PHd4aL#w4wVtS(X9oK84o7fJwD`SI|LrF}@uk
z3^^hLCEyw;s9^DDAfL$Kv@i+}0Bb=BtAU=3_yEJ$hiUrw0`Kauk=AsR$Aw7{{(uY*
zRR$`n<q)tkSriO=G9|Ym4{L`q;^+c0g+rodGWNf+n6nGR2w<bxf?s%X$b}D{s3Ty`
zqF{t2AT0{7%vT2-+RcfYl+9f&TW7Uw-PN-7R?D`RX8WO4m1%!TCVi}i3BRYF{=0lC
zK6&E7piDUF!22EKzYpUeoX>c80`7CRm|Idz$`TgB0#>#}r?jFytf}3ls#uF&u`Ijq
z*zUgDg1_xOfrog-k2|j2-IHZXQ0o>e0>5p4?-79L?fr#HSPwHv2+0C(j+u__xciXW
z-dFjw6qzksU9Y5Frqgng@}gL_&|W#!6Jl%RuWIIgS-r?PN5X)*Tk5w!p&>?7H0EhE
zW|K4XG#m$>Veiu2JHqO{75z1QHRP2$d+p>F&PODSnj{ZR=X;m$^A7LA-yYxN`xa*%
z{O{yCntt$(kA69{8k6DJkLeO1$8C)L>*?j(yx@%B;}wlmpp6q>q9ROlIKwNvt;TP{
z`l30*7c(!nJTH+vFK|3B7(CB9ysYR1265Ns5nVp}I0`-lV2e(o!2dOtHzA?H7RcZ@
zqpb{f*W&o4&E+?DpuUJ^?>IY7t84XJI5r<n@z!pdC<Ta7@ZLP3sxZ;)wj6Fe#z{}Z
zvbgp^0BUu2;QLW9!q1KH!Utkpc|O1PPWW$Lfh@n_r9c8Fel_+_JYzmF__cUDHGl@;
z#4|E;6sj{d#1&BjS1yR)Nm^hVCiCfKWcvJs7it{`@A0;!Hupy(;y!^g;FF<nPRhA!
zV_UdpYE4Yt4KdwY8&h%WK}^xf_i!=<#K40A*W+Xg>F7pGreX!wr(SN&=Im|QoX(A!
za}q!&FTE3V#1Ph+h<7$bytg*u7nkFCI2(g_lQX%9A8WRuw_zK)x6+2;cnX-$r?<h4
z?y4=B`4O1k2&{9P)BiL(<e(iARdCrm$e(EWxWk1JvlgPJ-&(e?gL->(sP|Fbw2$gy
zA9Z8n*ic<9(*G#yg5Q;O)Ed9I?RX{LvbKK-Ms^i$pZddaJbsH~1g*_(+U7lMGYRG?
zIvnUsIMza?+eFdU*D|_$(C003G<T$<$>-NIXyVu?AW<Cg8aPfvQ)Qk5{l@ca&y`S5
zc*dbE06i@Os>YF}U^PzAw04)_sfr*9CNSX45T(c{`kvgOobt$V^N!=4jbh53;hn?@
z$Jf*imoBi78F3dT_^(Bo*<VR!0QiXoe|&2tc_377O$@M!x4`~wvql>qm9l0tMmVic
z#t3TAm(2)#y-h|qtrZyo8NYQlXr}5G*+7QrQ?TK-%u=RT%JfO#O@!sWk^J)c^+z<L
zD%gCdy%A3Ti20*c?I2~oJ>C-rU!D?52DPeLv@t_k2~$21PZSDY9#@Ddu4GD?Y23Jr
zB7(Hm;E5z+%@t+KcB;9uv8wER4s5yY&P(DN+wY<7X=LU{$rjIj>+0VG2I-(v;E@}u
zs`05|tzB8?CY9k}@x!@4XJJgMxVU2{oSfm5+9Ws=*^BgNAacoi?td~G<A-w}XF%D1
zHmJ8hKKgHV{F?pp>gnOpo8yC5PucU={}feZkr8_v0r|_JSV)d%KD>ou5!sKz7=Y49
z{`7J@mdmHO4ob_9Zd<I@rYohOEG_y>i=~lVtX=SZxa9k2$@hkDyR?(>;jnQcUWX4)
z8bg&jYK+9|3A~b`)9xwCcM?uUCY0w;kV6R$S#dTbjAAp%nJ#X`K(l!ep@eisR?&G9
zd-jNc+1VE5H@UgF7MHu=sbs*KLx~Qc!q1`97)#n6HT@G@+41$o$rNXWY)3TTrWKYG
z^@7VOI*l?*<Jc@BDVttGr6v{URF*&c<&~E%81rjMV)pyh!|Pp>emoROAH^q6@_>N4
z8$oR?06Pn!c$R@h6>*f)xz9Z2v^krdhxpdxxI+P#1U#^>j~;Q5-!X2OC;FX?FVW0Y
z0CZD`Zf`w^Io=WX;QtmMnkRyh_zJwVedMTpLKC=CoCA)6nMpHwl0R6$vN$CZ^h~-{
zWW8o0uN=qBPvdpuzH+&8D7OcO7o^)}tAm<^K79*2xK<A*q<hrUaZH8J7WQRi4xvST
zz~V<L--7Og{TapW2mHAMKL^;hB{_s?1W3v%j_{Sbu9+Xt1Jj7XC<WA?q8{&{m1{AB
z$1yji(18FXY!P-8v1zt&vi5tNyiH#pLYAzY>eHpmOSY<t%3ZGP>A^pKe#)M_etL-8
zO*30+q-7pX;Uye~Fv^M>Q5)sxuZ%?wyx6_JvOuK-4oOw7abN6_Y0NQ!Jc!2JO<wPR
zm`1;%PLXLY-^+qBDZLfkBaU!6n}8IMg5bN5+fpBfu`#$M1Ndy3QF{k4eK>c}9eHR1
zt~-F|!#Ub0uyf;JWUI;Qxx?q6HDk_@t*8{#49DQ?T2BBD(;@&S==*%J>r0yEs|j=2
z9_MBaVq-@(2{%+wYg*8TswOzAv>|DP!qpH{z3gwVNfiDz4Vp3Z>7&sYyKnyl+eN+R
zQ~!N9f^YFK3TN}%mYT7frIjZs<xy8DDCPB_2Xp(}L6V?#Awq%43on9LbPD*{&3z91
z;d%PCgI@)dha&d^2EzBQqQ5B95@oE_1RC=8NF=313K*J5CeT=br=0vbf*O+p6rYE)
zbW<ACkwp&75+<T=Go2v0)&AphRFY#^;|o(X>knr~ptB>8*b&r$p5XH00-XdZBHc()
zgj3Q3@c8B+f8tmch>zI+m+b#kFE<bOU3iw9|Lb%dC%6B%=fLNB|K}E-``_(e#?dZ#
zqPxN5Jv&l%Vx5P*JUsg8kUcv(eD&lxbW+C$;0nyoZLp0O=$??jXMobpF6X;Rq9HQP
zqUkV*<7}oR(MX)%+4vg0?HFw55IXCl-Gbz{z{|6b&?!ptT+YKW_{Dg@r0m0mKZ{F0
z452>;t)_1{oWV1a%!lzu9L|HdLEQpsB_zLv<hMzq-r^@S+j<;Ao(1!_FoyW|W$+fU
zc{`j=qfwkWa2W@4&agP1ybSTvccV)&Pzpw-@KNw{zYjSFg9q~{2p+tcPDG(>sz0|Z
zcK;NH7c7AhH}26Hu2Vd_fwh5!qdEZwyq1B+^IHsnm45c;^f?$#hhSgIvdB-AErvc~
zauy%N&wL!a#@@4o7jF*9wYj1Z7LJ8%>?6&wKLci)c;qyu7F(Q^(@!+rqdZpQLj=sR
zhw*#x@)-_{Y~+KPl`;-}=y+@sA8i2zix^M`^jkJns&EsC1^$SDZU8agXdsdVO;YvK
z;W3kvCOYevK6Fn1NDtHI;ocoSh=z><4mJD~%%2FXfEDxD%R@uf^D%T*C10i!JbHs{
z&;ZJDZZ0!dSsuJRIM#psU1wzm*s5&B>{U!f!Vk}BPoffC=z{;|!$=Ev2h5}Cj5gE3
z`#GL#Jqgd)3^zQ3p!TmGfOvaeld?B?Q$lFcv|3Y@ZGvO)(x1^#A~hQGsSP706C9B7
z$Og+}CkUgw`lFFO4$cA`?$NGG+`|HccuphywAluZ^uZ%E2+mu8g9eC(`_Ni3&Szc^
zJt!Le3qQ3@5(l2zuK|;e0-O6M<A6+=t&-5B3FX6NF0+PlV3UgBxvLY~DOXn=5SRbj
zRELclIy9CNkcF(2#ouDiJRhM!sbx?VLIt0?BmF|3q!&FKU*L{^>_Tt5jwR=R;76DI
zzwUm|?e<*s|G*in{l9HIH#z_F?3Y(hj(>jr3TJwcetE@ydI63gJA8BeGgx#;Pkc8J
zt*-QfY%rNoaE2o#^pI@yLkeEdiH{~HjBvzj@8*5NU4s@m;>#216{+kT3mdxz{5$gu
zuo?!=+Qf7G7MJCf`7GMqjjjlazg&8j$dQ$n=ZFf3ip}<D$oQ!$$(cHK<fWkqhNki6
z#Q~b_e@>?(=m#cX_^%Cay!5YZVYK1Li@-m%U%|_pAiypnj-d^>ltG3TdT!2a+or|?
zcBVK%6{@pOG}}fE5Hsw7JiDhF0^?255VIU@`o@Jzs8{M0F-S<Yq^1nP(puw+YO@uk
zf&yYtE;+~q=o@sZFq0-2waUaS<Q&8>;cf}0c%UK0ssgDJQV=Z~iy)WN;I2r((pJcW
z9zz(4Tu#g+*;0g}RhATDN~a{Nj&DXjf<+>dlI6a48xJ@xN0$>u9h7M?6R@87!(V;y
z;K;C0Yk=NqTR8oc2;-}mO}Qm>G}*F6ujz;rCi<XMWTMadyNarG2t&(~P$g?XIJz-}
z!y_8s7MtU)-JHo3XS8mN(ZSOjgNVZmFh~*`Tv;!)?)Y>genGLDvkZaQH)k8+QpP#>
zurc#)p@G1^TY#dmP&TyDK$8R`rMVrRieASaLLYjWok5Y4Cs8oJj3%ama}ahoGT5KO
z$)B`AgD_4AafqF}sN09(48<|h0pb+S<=E$=QyVjJ8mNFfA}sA2a2(0JzqyWWoXs}L
zgbH}VX~D0;Zb?vH0+J>!B-~T~A{<|PFs7Pbgg#klh9kiO>=}&T=O%t?iR!gfH36*G
zFP7l2T}MyfWWMCbm8xmrqR5dz5CffKI}+#BWRFkZbQ-bhJ|eNk6(O1MO!M>vmjIxI
zk30v~{}y<!4ysD%-HavU#6-crFT*I1J|?Ef$vhpe@f3&j9uwyosq$oIJqdp%NX+hg
z?C(DAH^5$nXhq42Kv+g3qZ=s*BI<g8B9@^+RV=FHNoiEsoeHD`ixc^yG^`$vK!I1A
zaxV&omr)$P55#&si_fPY{^>`!t0S2L!NNm-ZW@v5R;jVn2t#*j4N8ZW=o(E<c7nLU
z(>54|5MF{7#{WzVb%Yi~K9M0JJu2QT3TA#JW*Nw_5ma4Enih#zM!wbnP|XyE3u$O(
zn9rU*6A+1{uYaXENim=I3hFjEE{K_`-ByMgByH&2$jCGDsjy9n5fS1hK591VENi4d
zXK%;jL8qtD$P0am6{s=*MZMJ~wB+!PU>c&RR+8Xl+)L@i-q^ncxAK`ne?DbW1_k#}
z>=B<E*04mtGkwj*!d@!lEIW{3E5OkQ-4#P_iasY6L$D>n3SpC|tDTFmRV5SD$)G#O
zT_aO?`7E{4jhvmxKq(QKoIRDDe;k05F;8JB_oU!b5NM7^*q;3Jr9Gd_Px&l}hV}%<
z?fmh<@#|x*n6SES!4fGLuWdgPQ9d3Y*gtu=9R|07&=04FyE^Zq#krWi$8k*{*9%@9
zG!g`9=hfv24hd7FMS^F37smF)Pe1hg#t7JC3<j3P{=6E}Wk8+^2v#X}zXj2BHT?1*
ze|%&=vNDxp1%hEoz&)yCG6h1$F>ZR0vo6^hQ^+Og>nw$KR1l(Aj#!Qj=YvLqRu(Vd
zWO#9T!9XSm{k6#u&wvi(xt1a&yjhxFH&01T%tC3L*IJh2M=7|pnL)J#k@iof@39pu
zmC;Swk(P$m{FoDZ_U97f>LtMMl+JvTA@$t@Ja;ep#yR};qY+(F_}<Uj)vI(wU|zLw
zgDXa#*^Ni|kK9>>b1NcxA5be!tv}X-(=PM#z)n7=hR2p8j0%e(EQ1iRtoeZd5>+SV
zMrF*3YoeeSPv@zuVNs(gOlMIy8);lj1ac{kaW9aFhSsPBxz*r`isS-tmK05wuN43a
z;5r!a0t)G@EcCu=yYMvw;>GObnh{kSbbZ({-XDgWowJn7s`bO7;)>FbC%JN>Om&IV
zLWza&kJasFLRF~7mW0!wOVM)BRB2u0#bP~*f`@;y{R>KlhkvrkDI>h}HJe3p&y)XV
zKjygU<$mWX%QV%?{;bLh9P%>1vKfz0y&REC06m^~^55deQ~B|<<fHiXc#_TnEYk4F
zv{v$3$-hh`|CaRl&PRgM-3&gUjEd3EeKyCVGaq1F!A=5piF?>WoDdKbm?}=qEv!8>
zii4hPbiLFbiP4WXwU!wOH6{C-I0_F>&F^FwMw%hjRw+wXBDwgsU;z^b0`Yu0GmQxp
zl7K0~4N?UPGXa$)Fey=ohg<Gh3X$P9ZRMNA!sl%W1*XWsWNzdX6C@kFiS2B11Jp$>
z0^5cNm#9odp`83!RkozWDb6KoSxBR4^CG?32%w=kx%u+=$;eJ=B(O5>^`b%Qx&u=(
zv$!*pz4-afF?;s<h$7E7L@%fIGp^nnuGj-wa*CF0U0Rozy{WhCd<IYHqRcJVv1<Hh
z-2c-N`+vI5ejWdL8_$M|a%<JGWdDD=<Ldl({-f=7y4`yJ*H)hL+kcp;me=Tp9_)EN
zqgo#HC7;Zygj()2X`u6mvvYj9gTCyJ>CnIbm%5`Cr2SvcMYVRw#iRAdq=Lqxx4KVR
z^I2~L#AbCY(|>o*P4&OM5Bl$P_j<MdZ{yii{}1n8?RVt3FsJ^Xnqptqv*y;Gjp|sc
z|ND-n|9Jm{)9dZk`oE25Q~f7-*opyAb4#@uxQl1aU56Xhv0VQ<ss6h?)c-xlsr7#w
z&&K`#(U$c8v~~pU^4VaM0qGDH@9^4nEY<&BFV+7J<b(J9y8qwG^LdZ|YmL9_XN~@U
z;~~J(`9CK;|JiZ&Q2*Qeb^MR5JU8h7w>AXOT!4B=aLdnnXO}jsW3BnWZimMIZMW9{
zZ9JRm{}$c<TFdY7N%UVk8+O|fAeOBEI&OO}AOCC5-LLh38&5_2ud+B>nuH{ww&?IF
zo#nzf3-l{F(IXG^<WHngAYT>BMI_BAh$#g)q~nGSdCpFTmJu-nv}iQCtYXDhL+ppg
zQ%_pCOmT5cb`;G0a2#7)znD~P&q?DdU-esi0}8+i5Mr==Z_xy)Q<PJ6c(SRrlh8nv
z@VsPG;*&=D6``xl%L$oNi5CG<Xi)hmh%LFbU;TbS(B;|A2dCi`g^&zRPwhhscTYZj
zrcP?oM#Cwc8#4`x&tvdxO)?~z9jeV_xbKi_sCwFa2KEE7-%}7uOUC0o+D{n{Dc#Ss
zm}Xd#s8!qaPxe%6_wCiO#Q%e?A#eZpx}Dnp+sZ@n8&9>}uR<UKvFTNwe{K{a`OqQ1
z+(GyB*gPKmLsLSB2ItZB6O2M)Vi`0B>1WP@D4IrTc<$kJGP2}$!B@Y$cwylwu>0)E
z`4n$)^VOMKieXA|q~-mO1V(HsEbS!<CxRNnv1JoH>*TQ*qorQ)<|w=Xd(>}Engf0x
zzd>f+xF<Z%Blm)NNAbTe{RoUOIJ_5B1!Q^f(vQwUAYA_Z-;q14D)I!H8X`W$QU^aj
z0i573)t>z}3f_l!sUOB$l&>YZGW~>yzfPCrdNUu{hs2Ed-_hJ0h5niOemH$kS9M#N
zg5uCrd3RhOG-S+@T2b|W3C-!4-#XMr=o=BdP%umL!sjZo7cy2+eqqKX#D}sUD+HVT
zj;!)Ac~n<rkekKx8f(hUxws1)58;7RW@<$+`1-vl_3g}z|Kw$U0Vq%s<eToSc%=xN
zEK`ZNbtD+aL3(eAtQGY?xs)W?M2s$yQ1lE61-ZMI9vBhFE9PcAkK)-lM4@9CB}}kD
zL|g3loxz~w>~V7R_?~dfz(14mRYn-7+ZbVVZbL>+U@8CtGbm;PGo(X=uX3-q!SO9H
z7=VoXExdr+8i)pwDL!H#ym;TdYE9hGk0(UBXsg&4u{dK8C%f%Y5&anKB}s=l^e1>*
z?)Y*Ph&Hz{ay7-7j@fiFzLw?r36z^mp)5O_PDd<+(L2Jt2(MHPU~9!?;08MP@~kU;
z0+8uPV6I`YXGgDJvh(?T_WkbeyLcR&4$lMo!!#Pr0A~yXda>)aoz8B%vkQMCw*~%b
zOPrc$T@aoU$(E>M7l>5U?F{R~<o=uc7{^o`>O{Q-+^FeEyHBo@V@G&;SPy}Y-$p?k
zo=s$Ib0wPp@S|r~9iGHWhYLc<LyP_C;h&(pGs&=py!o%DZaRzjME+znz0eT~uI6-;
zG~knDOhbihSRDBxIt&a^-dtY9rt<sA?N?qu<a;uWBET#iX1uf5I}r5(UciuWbm@<i
zdiBci8(n5HpGP`WBOV&o|GFupC;oVNImU3d7t%ztqn2bCGc)wxoloalbqP90G+shN
z?_yQYc4TXJ6e3nxb;H!QSJt?brFSI)bRkTpZfiLzk1^Ep(f_(JUIz|NloT8avM-$$
z(iQAxL<IJvJ|)mlB4@ImpfY*lh#iponTQf({N!K8F(7Sl)e(9}7fT}#Xx*Q8G~|GY
zhrL~hk|1rXjtZoL-KCCx1(kzbx2kVRZ@GfTfv)}6Z1&(^vZuT?ZPRhM4SmbuR@Asm
zh0AnKgLV_GOTet4ccoR{hUTS6tfPAxIvZ$TTIX%)U%|s$Lj#ktrJ5-GKj9)|#%k0u
zZJ@?ZBOo(%`tbGXvqsgVPey^&RFggvgvXVkP@W=lx;5z40e2tD%Yia8MwXh0qp~y=
z2d<Kqg^L*^J;g|z-}FqKB`WW83SViiJ|mWC0tV$y)-aZoR27!AYG+$WB!p;44_G67
zQaa=H@=~CE=ldQVxds+AvTk#6{9kd4!mZ9RESdl5b@mGPf9!Yb{r_8e8d-}%v7O@i
zNDd<CoHH#ECbwDWhed_`zz|R5b~i*csAR|}BhA8)l&<`ua01CF2^a5Wz}ktjNSSDa
zTue8&%~L2MVnZb(R>~P%K_xKwq^YDNz8qabaXJ`eLs5y)e2XJXWr`+0`dkYgNhucl
z&SS>Ei~(MjsT$qDQEH>&LXdKDQS5)JD~t+PhG?JZr7wPB&BsUoO?8*j8nYp}DZAb-
zOb|^g;is7=J?18wQYB18TPsRRdU;DKDM{K5BqcdrOk-^EgJeeie0*TP`TJKbK?xZ@
ztppK)<upc8RW593QYlZwhN&3S@`?=+X_|mYp6OEA@5_&xfM&(0R6w~brE2_0B-A7}
zdw)eSf(xOG=nw(u+`XsGYhps0QN}~TU?^&^r2@TY%3xBi6p1bcQOV5XX$SptVH^l4
zwFs)11OEkAGw`6rT}g>(8l&kj-o=Z}PRG*^yD&zHTXakxC4RBUHnQ`7Y5nOOA~ay4
z_GTE8Kf+0pi<;TkE;O+(+NuKp2R(eR%8Q_673DI;fWehz)ZALO%;@CYj_`IDS&rTd
zA#+e8?dU=Vgk@TJ|0JAZ>;r|Xe%T07rvIIpPZHrX6mK3tX>j`x=ts##mghT(O$p~I
zj{om#^BRk-!0r2JALlSyz&X9Zp8ypu!EX$Q@J)tI(;LxGe#mT|=L`gHdC@EpB5I0?
z8p)1dKYjf@d-D9?)lY|bqf|I~OEw!Yk-TT|6bF;wfzDz0Eyim|-Q{wqRS20<;E80Q
zReO&n2stBFmmt0$421ILf&&Ur@1+!%a51^rrZDde)JfS3J(-AkVk5cra4Ir^VNjLi
zT@QtlB&Oj|M3N;dfI-9iAi^N{D5W?#q74_W9!Li4s<#IT@hsMtaav>GrVJ#Qfk*?l
zGrV8JNYF5Ndd);bADLw4onA>rmsGtH_kbz&&n-Z`QW!yQmSmM_`K(FMP`YsP%5=D5
ze^L^%Sh*^4?Lr2AYKQwXd(!dW<$QYEfj{K4V*HP8r{_XG9RGElUayY-v5hC;|68G6
zZ|!+-b<1_E82<zLzt`<LIR5XoyY={g8_)lM`H#Ui{?5Le&_t8Awv1QJ*+(c4*2l*3
zj4-O1o}OYv*iO4=LAN;?;k5YEXS^AXnjCMT@iGpNMg7QwQUAkvI6Sw}tW(071s(Sa
z6%uAir7f!SQ`7THdvg2#?;`xu8u6bxJ@`cTzjfQS{olq@#{cxS>{gxi>sV(0yZgld
z&K{lw+3(Q#&%OHm*H)f2&VQ5(gnkI5{0zu=diK*F&PH*cDqyQ9glp9y%hYkEGS6!I
zf0NUn%k{t4QSl!--Tgh6&VTLI`+v9cY^eWb$3GWp|DBxvT-Q7-L1LMOP^f%O$5B_*
z{~MhCTtWX`TL10R60nH><aFxzk6U>*)&K42{$vv;VBP9;&CgYxRQ%OC2CM1+O`^Z7
z;Q#L>`rmE$DE^DvtNs72JR9nNA`o|F1h(3_e~XU4d{*}VZ{Pqd_y3*j{EzFn_4t1~
z&*uLBww?daz}^1}?Zj&Oe<Sr@uK!Lu_5b%M{-d*3$A8?)v!VWPPyIiY|6dF|^EIe@
z;aN@pZxa1^mGOTs)qj`9|8BR||E)Y5>i@PT|32Mh-;&9_uf-UwI{$Z*6Co@1|JwRr
zySG=zf7;5kssDe6%YK_X{Pnco+Rs|^e|HrBvD?e*zq8-2`~R&xH=6&u6C<#_RsXNF
zOJD0*+5f+R1F*dR@8<jez5ROqZ!6En{{Njj0GoUMwcEe?v-14^UFyGE@c;LndjHQ>
zo{i`Kw=e+Md<yV3EWo!h54f8CKi_u$zhB3H-_En4{y)F_{vSKwQFwa#X|zAB{p;0z
zd{$lmyUY2%?*4wk|99*4|E)ZmuK(S&`rmrDUoZQWK5L)<*<c-TY5(8W_W#lT-;UF+
z=l{3z-1z*@ov#Caf|b8-$O0t#|IR;)KS}(5r|WjRh4}yNI{w>M9u?SsHolH5{o0qy
zc{q-DNk8zS7ZLGa6Ozv*;(r=erwX|%j&aH;`25gp+wC4_ZQH>=^!AAU^eEuDuy!HF
z$nn?}1}=eDjBS4J=-i*n&?8q59{=*w8&S~%6?2grn8O%QB=mbvDVlsHBP7u$Et34d
zFZ~HnImkf_FJfB+FzsWwGnwBQq|KmzMNRQ|I+GEF;6r$7Qs9;0bSwh9R6({p&(Mk*
z){kvv@vEMdu+?D>coyGWnD??rMi-TlT|_bKW4@twbc3ea(Oonapx#|`gPSxrxRd4v
zcgfuwwKuw(21mEu;&1|{Wf;unn<;Ou_^q`$5BJs@ZQZN1)*_fzwnO3U-YP~LFT1#>
zMWf0f%NY>5?52uru&h#?%zS@`MRNw2U`5jp>M}FBvcg)@OnW*WZKJ7J&$gQT<z2-g
zFW@L3U*JyKfK{+g^a6DG4D?`ZVt_Y;$9>)#fZ94%R&{=(;BA!WV@X*NJZ?#Gm8%Ru
z5^ru9j5K2z7(FoEGDw-cWx&#acFRC%Ie||_NGN~RFfIjcDN`99xP-m1EaoUH_<U)>
z=j8>D4wn^#9OVVw<t1I6NO_fuIq8G<8z%+vE>lrmq)ipNlg4mPG%ism8Zl7pbDKJx
zoGX5feGM1iDWtB7;rhEk9vUq>4$9dVfBOcRAY(A_JT!}G`2lu(f*^5lc7ZoHs9H+P
zwljtKqO6ayS)p6@7U4t#0Gn<LFr+1zO~YmQWb0tF;W_i1JyCofT%oOLVIMH|PoU@Z
zk{cGIx%a5e*)X^QL)^Capy(8Aa)&d2^iK37ylq*H85q&|LRzYU%f~YA*&qS7tW<hJ
zg~=j75hab|WsnL#$|__{YoHeTH@#Y*Pcgrsu;;-a_Gb|kOak&$i_1L@nOVrTCmbo6
zCgj``C=3b>lStcTvZ(45eVc5hjSz%KI2?y~TSn4;yxsE2i=SWG|2R0(!HK#BDPZIn
zQtV1fV^_dW_?47m*OZVU6@Yz^$`n+XUVy@&YPy!<l=?=*AQ^%F;MlH}aNf*9f%L||
zpwcunE(JnWcwl5)5gM1t*tkSE6N_*mQk>1NXXv3f$4$=0^J~1db>zqA!N?#TF~Dd+
zFL4b-@`h&1;jV@44G|>4p~idohhwWf{KJtv9{#cDmW%OpI%`lgY7hw^oQyS$je!&u
zS5`C;6NGqR0WIRRR!FAyB>)n&fu&ng0ZZIZp-HlS_uZ9tm1~Y{eXWv2W<eCQ^+`A@
zYnc>M4+RXNXm}--K*+f!YC;V(In<>pH-ibC%}0!DE0Zdv-ImvX=~}NE?-?%-jtz{p
zo>Z8b1`xBiNvQ0hbcqtbh3no&E(&fo0ImUoM8G6A1WenI#*B1zTJe_tOu%c!?6PM`
z0@>mx456HU5|E|!_9OG+6pB(s1cQtG7OHbZ$u4)f!yQPOUig&K<@i4b{Ph$86wDAs
z()pg?q<p|H@8lN(o$cH!_>nO=UWx7T@@7a#NZfxn(4{2qhl8azN6{nA7sdgWwhLdc
z)ob>KNa(~vbW{qEFQ`jPQDJciq70oD-ulM!B?Oop!Lp!I14@t_1*uzr^lkw{DJ=vT
zHw7c}twf_tM$##YLcdU_jt*B;su+xEZpMk{CwlHm)iRm}KI(?br8VZO(l8R`3Og*M
l>Nd3}SMN8VLZ5n#q<-qBe(I-w>gQ{F{y)v?-;@9_0sx@V=0E@d
literal 0
HcmV?d00001
diff --git a/R/analysis/_region_.tex b/R/analysis/_region_.tex
new file mode 100644
index 0000000..39ec4cc
--- /dev/null
+++ b/R/analysis/_region_.tex
@@ -0,0 +1,5 @@
+\message{ !name(lmer.output.figs.R.tex)}
+\message{ !name(lmer.output.figs.R) !offset(305) }
+ddply(DF,'id',.fun=function(data) data$d[which.max(data$g)])
+
+\message{ !name(lmer.output.figs.R.tex) !offset(-5) }
diff --git a/R/analysis/glmer.run.R b/R/analysis/glmer.run.R
new file mode 100644
index 0000000..b9739f1
--- /dev/null
+++ b/R/analysis/glmer.run.R
@@ -0,0 +1,192 @@
+###########################
+###########################
+### FUNCTION TO RUN GLMER ESTIMATION
+library(lme4)
+
+
+get.ecoregions.for.set <- function(set, base.dir = "./output/processed/"){
+ sub(paste(base.dir,set,"/",sep=""),"",list.dirs(paste(base.dir,set,sep="")))[-1]
+}
+
+run.models.for.set.all.traits <- function(set,model.file,fun.model, traits =
+ c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass"),type.filling, ...){
+ for(trait in traits)
+ run.multiple.model.for.set.one.trait(model.file,fun.model, trait, set, type.filling=type.filling, ...)
+}
+
+run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), ...){
+ for (m in model.files)
+ try(run.model.for.set.one.trait (m, fun.model,trait, set,type.filling=type.filling, ...))
+}
+
+
+run.model.for.set.one.trait <- function(model.file,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), ...){
+ fun.model <- match.fun(fun.model)
+ for (e in ecoregions)
+ try(fun.model(model.file, trait, set, e, type.filling=type.filling,...))
+}
+
+
+#=====================================================================
+# Run glmer() (in package lme4) for one ecoregion, one trait, one model
+#=====================================================================
+model.files.glmer.1 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.E.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.R.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.ER.R")
+model.files.glmer.2 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.AD.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.R")
+
+model.files.glmer.Tf.1 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R.",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R")
+model.files.glmer.Tf.2 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R",
+ "R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R")
+
+
+fun.test.if.multi.census <- function(data){
+return("tree.id" %in% names(data))
+}
+
+fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
+ Start <- Sys.time()
+ glmer.output <- glmer(formula=formula,data=df.lmer,family=binomial)
+ end <- Sys.time()
+ print(end - Start)
+ print(summary(glmer.output))
+ saveRDS(glmer.output,file=file.path(path.out, "glmer.results.rds"))
+}
+
+run.glmer <- function (model.file, trait, set, ecoregion,
+ min.obs=10, sample.size=NA,
+ type.filling) {
+ require(lme4)
+ source(model.file, local = TRUE)
+ model <- load.model()
+ #= Path for output
+ path.out <- output.dir.glmer(model$name, trait, set,
+ ecoregion,type.filling=type.filling)
+ print(path.out)
+ dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
+ cat("run glmer for model",model.file," for set",
+ set,"ecoregion",ecoregion,"trait",
+ trait,"\n")
+ df.glmer <- load.and.prepare.data.for.glmer(trait, set, ecoregion,
+ min.obs, sample.size,
+ type.filling=type.filling) # return a DF
+ cat("Ok data with Nobs",nrow(df.glmer),"\n")
+ #= Run model
+ fun.call.glmer.and.save(formula=model$glmer.formula,df.glmer,path.out)
+}
+#========================================================================
+
+
+output.dir.glmer <- function (model, trait, set, ecoregion,type.filling) {
+ file.path("output/glmer", set,ecoregion,trait,type.filling,model)
+}
+
+
+#============================================================
+# Function to prepare data for lmer
+#============================================================
+load.and.prepare.data.for.glmer <- function(trait, set, ecoregion,
+ min.obs, sample.size, type.filling,
+ base.dir = "output/processed/"){
+ ### load data
+ data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.csv"), stringsAsFactors = FALSE)
+ fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling)
+}
+
+fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs,
+ min.perc.neigh=0.90,min.BA.G=-50,max.BA.G=150){
+## remove tree with NA
+data.tree <- subset(data.tree,subset=(!is.na(data.tree[["dead"]])) &
+ (!is.na(data.tree[["D"]])) )
+## remove tree with zero
+data.tree <- subset(data.tree,subset=data.tree[["D"]]>9.9)
+## select species with no missing traits
+data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
+ !is.na(data.tree[[BATOT]])),]
+# select species with minimum obs
+data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in%
+ names(table(factor(data.tree[["sp"]])))[table(factor(data.tree[["sp"]]))>min.obs])
+# select obs abov min perc.neigh
+data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh & !is.na(data.tree[[perc.neigh]]))
+return(data.tree)
+}
+
+fun.center.and.standardized.var <- function(x){
+return((x-mean(x))/sd(x))
+}
+
+### get variables for lmer
+fun.get.the.variables.for.glmer.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
+dead <- data.tree[["dead"]]
+logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
+species.id <- unclass(factor(data.tree[["sp"]]))
+tree.id <- unclass(factor(data.tree[["tree.id"]]))
+plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
+#= multiply CWMs by BATOT
+sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
+sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
+sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
+sumTnTfBn.diff <- sumTnBn-sumTfBn
+sumBn <- data.tree[[BATOT]]
+return(data.frame(dead=dead,
+ logD=logD,
+ species.id=species.id,
+ tree.id=tree.id,
+ plot.species.id=plot.species.id,
+ sumTnTfBn.diff=sumTnTfBn.diff,
+ sumTnTfBn.abs=sumTnTfBn.abs,
+ Tf=data.tree[[tf]],
+ sumTnBn=sumTnBn,
+ sumTfBn=sumTfBn,
+ sumBn=sumBn))
+}
+
+fun.get.the.variables.for.glmer.no.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
+dead <- data.tree[["dead"]]
+logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
+species.id <- unclass(factor(data.tree[["sp"]]))
+tree.id <- unclass(factor(data.tree[["tree.id"]]))
+plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
+#= multiply CWMs by BATOT
+sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
+sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
+sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
+sumTnTfBn.diff <- sumTnBn-sumTfBn
+sumBn <- data.tree[[BATOT]]
+return(data.frame(dead=dead,
+ logD=logD,
+ species.id=species.id,
+ plot.species.id=plot.species.id,
+ sumTnTfBn.diff=sumTnTfBn.diff,
+ sumTnTfBn.abs=sumTnTfBn.abs,
+ Tf=data.tree[[tf]],
+ sumTnBn=sumTnBn,
+ sumTfBn=sumTfBn,
+ sumBn=sumBn))
+}
+
+#============================================================
+# Function to prepare data for lmer with new CWM data
+# that will be used in a simple linear model
+#============================================================
+fun.data.for.glmer <- function(data.tree,trait,min.obs=10,type.filling='species') {
+if(! trait %in% c("SLA", "Leaf.N","Seed.mass","SLA","Wood.density","Max.height")) stop("need trait to be in SLA Leaf.N Seed.mass SLA Wood.density or Max.height ")
+# get vars names
+CWM.tn <- paste(trait,"CWM",'fill',"log",sep=".")
+abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',"log",sep=".")
+tf <- paste(trait,"focal",sep=".")
+BATOT <- "BATOT.log"
+perc.neigh <- paste(trait,"perc",type.filling,sep=".")
+data.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs)
+#= DATA LIST FOR JAGS
+glmer.data <- fun.get.the.variables.for.glmer.no.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
+ return(glmer.data)
+}
+
+
+
diff --git a/R/analysis/lmer.output.figs.R b/R/analysis/lmer.output.figs.R
index f044784..dc5591a 100644
--- a/R/analysis/lmer.output.figs.R
+++ b/R/analysis/lmer.output.figs.R
@@ -506,7 +506,7 @@ names(models) <- c('Absolute distance','Effect/response')
pdf('figs/R2.boxplot.two.pdf',width=16,height=5)
fun.plot.panel.lmer.res.boxplot.compare(models=models,
traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
- var.y='R2c.simplecomp',ylim=c(-0.015,0.06),
+ var.y='R2m.ratio',ylim=c(-0.015,1.5),
col=c(rgb(87, 157, 25,maxColorValue=255),rgb(255,102,51,maxColorValue=255)),
cex.axis=1.7)
dev.off()
@@ -537,9 +537,15 @@ names(models) <- c('Effect/response','Absolute distance')
pdf('figs/R2.MAP.two.pdf',width=10,height=7)
fun.plot.panel.lmer.res.x.y(models=models,
traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
- var.x='MAP',var.y='R2c.simplecomp',ylim=c(-0.015,0.065))
+ var.x='MAP',var.y='R2m.ratio',ylim=c(-0.5,3))
dev.off()
+
+fun.plot.panel.lmer.res.x.y(models=models,
+ traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
+ var.x='MAP',var.y='R2m.nocomp',ylim=c(-0.015,0.06))
+
+
pdf('figs/R2.boxplot.pdf',width=8,height=6)
fun.plot.panel.lmer.res.boxplot(models=models,
traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
@@ -550,7 +556,7 @@ dev.off()
### plot parameters
-models <- c('lmer.LOGLIN.ER','lmer.LOGLIN.E','lmer.LOGLIN.R','lmer.LOGLIN.AD')
+models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.E.Tf','lmer.LOGLIN.R.Tf','lmer.LOGLIN.AD.Tf')
names(models) <- c('Effect/response','Effect','Response','Absolute distance')
list.params <- list(c(Response='sumTfBn',Effect='sumTnBn'),
c('sumTnBn'),
diff --git a/R/analysis/lmer.run.R b/R/analysis/lmer.run.R
index 216a991..15d94cc 100644
--- a/R/analysis/lmer.run.R
+++ b/R/analysis/lmer.run.R
@@ -51,7 +51,7 @@ return("tree.id" %in% names(data))
fun.call.lmer.and.save <- function(formula,df.lmer,path.out){
Start <- Sys.time()
- lmer.output <- lmer(formula=formula,data=df.lmer)
+ lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
end <- Sys.time()
print(end - Start)
print(summary(lmer.output))
@@ -67,6 +67,7 @@ run.lmer <- function (model.file, trait, set, ecoregion,
#= Path for output
path.out <- output.dir.lmer(model$name, trait, set,
ecoregion,type.filling=type.filling)
+ print(path.out)
dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
cat("run lmer for model",model.file," for set",
set,"ecoregion",ecoregion,"trait",
@@ -86,10 +87,6 @@ run.lmer <- function (model.file, trait, set, ecoregion,
}
#========================================================================
-output.dir <- function (model, trait, set, ecoregion) {
- file.path("output/jags",set,ecoregion,trait, model)
-}
-
output.dir.lmer <- function (model, trait, set, ecoregion,type.filling) {
file.path("output/lmer", set,ecoregion,trait,type.filling,model)
}
@@ -106,8 +103,8 @@ load.and.prepare.data.for.lmer <- function(trait, set, ecoregion,
fun.data.for.lmer(data.tree.tot,trait,type.filling=type.filling)
}
-fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.genus,BATOT,min.obs,
- min.perc.genus=0.90,min.BA.G=-50,max.BA.G=150){
+fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs,
+ min.perc.neigh=0.90,min.BA.G=-50,max.BA.G=150){
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["BA.G"]])) &
(!is.na(data.tree[["D"]])) )
@@ -120,8 +117,8 @@ data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
# select species with minimum obs
data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in%
names(table(factor(data.tree[["sp"]])))[table(factor(data.tree[["sp"]]))>min.obs])
-# select obs abov min perc.genus
-data.tree <- subset(data.tree,subset=data.tree[[perc.genus]] > min.perc.genus & !is.na(data.tree[[perc.genus]]))
+# select obs abov min perc.neigh
+data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh & !is.na(data.tree[[perc.neigh]]))
return(data.tree)
}
@@ -186,12 +183,12 @@ return(data.frame(logG=logG,
fun.data.for.lmer <- function(data.tree,trait,min.obs=10,type.filling='species') {
if(! trait %in% c("SLA", "Leaf.N","Seed.mass","SLA","Wood.density","Max.height")) stop("need trait to be in SLA Leaf.N Seed.mass SLA Wood.density or Max.height ")
# get vars names
-CWM.tn <- paste(trait,"CWM",type.filling,"log",sep=".")
-abs.CWM.tntf <- paste(trait,"abs.CWM",type.filling,"log",sep=".")
+CWM.tn <- paste(trait,"CWM",'fill',"log",sep=".")
+abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',"log",sep=".")
tf <- paste(trait,"focal",sep=".")
BATOT <- "BATOT.log"
-perc.genus <- paste(trait,"perc.genus",sep=".")
-data.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf,perc.genus,BATOT,min.obs)
+perc.neigh <- paste(trait,"perc",type.filling,sep=".")
+data.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs)
#= DATA LIST FOR JAGS
if (length(table(table(data.tree[["tree.id"]])))>1){
lmer.data <- fun.get.the.variables.for.lmer.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R
new file mode 100644
index 0000000..42cf802
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.AD.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTnTfBn.abs"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R
new file mode 100644
index 0000000..efa145f
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.E.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTnBn"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R
new file mode 100644
index 0000000..4173585
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.ER.AD.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R
new file mode 100644
index 0000000..04663a4
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.ER.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTfBn+sumTnBn"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R
new file mode 100644
index 0000000..f6429cb
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.R.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTfBn"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R
new file mode 100644
index 0000000..e1a7144
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R
@@ -0,0 +1,4 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.simplecomp.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)"))
+}
diff --git a/R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R b/R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R
new file mode 100644
index 0000000..7a12a81
--- /dev/null
+++ b/R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R
@@ -0,0 +1,7 @@
+load.model <- function () {
+ list(name="glmer.LOGLIN.simplecomp.Tf",
+ glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn"))
+}
+
+
+
diff --git a/R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.R b/R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.R
new file mode 100644
index 0000000..ca20de2
--- /dev/null
+++ b/R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.R
@@ -0,0 +1,5 @@
+load.model <- function () {
+ list(name="lmer.LOGLIN.ER.AD.Tf",
+ lmer.formula=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs"),
+ lmer.formula.tree.id=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|tree.id)+(1|plot.species.id)+sumTfBn+sumBn+sumTnBn+sumTnTfBn.abs"))
+}
diff --git a/R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R b/R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R
index 29e481a..040d4bc 100644
--- a/R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R
+++ b/R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R
@@ -1,7 +1,7 @@
load.model <- function () {
list(name="lmer.LOGLIN.ER.Tf",
- lmer.formula=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTnBn+sumTfBn"),
- lmer.formula.tree.id=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|tree.id)+(1|plot.species.id)+sumBn+sumTnBn+sumTfBn"))
+ lmer.formula=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTfBn+sumTnBn"),
+ lmer.formula.tree.id=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|tree.id)+(1|plot.species.id)+sumTfBn+sumBn+sumTnBn"))
}
diff --git a/R/analysis/run.local.R b/R/analysis/run.local.R
index 2d541f5..f17aefb 100644
--- a/R/analysis/run.local.R
+++ b/R/analysis/run.local.R
@@ -1,88 +1,65 @@
-#####
+##### SCRIPT TO TEST BEFORE TO SEDN ON CLUSTER
source("R/analysis/lmer.run.R")
-sets.inv <- c("NVS","US","Sweden","France","NSW","Canada","Spain","Swiss") #
-sets.trop <- c("BCI","Fushan","Paracou","Mbaiki","Luquillo")
-#==========
-# Test lmer
-model.files <- c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.R.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.E.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.AD.R",
- "R/analysis/model.lmer/model.lmer.LOGLIN.HD.R")
-library(multicore)
-options(cores = 5)
-
-run.models.for.set.all.traits('France',model.files,run.lmer,type.filling="species")
-
#### TODO START TO WORK ON A WAY TO CHECK PREDICTION AGAINST OBSERVED TO SEE TRAIT EFFECT
-run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.R",trait='SLA',set='France',ecoregion='AB',type.filling='fill')
-run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R",trait='SLA',set='France',ecoregion='AB',type.filling='fill')
-run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",trait='SLA',set='France',ecoregion='AB',type.filling='fill')
-
-df.lmer <- load.and.prepare.data.for.lmer(trait='SLA',set='France',ecoregion='AB',
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.R",trait='SLA',set='France',ecoregion='HI',type.filling='genus')
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.AD.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.HD.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.AD.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+
+df.lmer <- load.and.prepare.data.for.lmer(trait='SLA',set='France',ecoregion='HI',
min.obs=10, sample.size=NA,
- type.filling='fill') # return a DF
+ type.filling='species') # return a DF
-simple <- readRDS("output/lmer/lmer.LOGLIN.simplecomp/SLA/France/AB/fill/results.rds")
-nocomp <- readRDS("output/lmer/lmer.LOGLIN.nocomp/SLA/France/AB/fill/results.rds")
-ERcomp <- readRDS("output/lmer/lmer.LOGLIN.ER/SLA/France/AB/fill/results.rds")
+## simple <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.simplecomp/results.rds")
+## nocomp <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.nocomp/results.rds")
+## ERcomp <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.ER/results.rds")
+## ADcomp <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.AD/results.rds")
+simple.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.simplecomp.Tf/results.rds")
+nocomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.nocomp.Tf/results.rds")
+ERcomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.ER.Tf/results.rds")
+ADcomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.AD.Tf/results.rds")
+ERADcomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.ER.AD.Tf/results.rds")
+
+
+anova(simple,nocomp,ERcomp,ADcomp,simple.Tf,nocomp.Tf,ERcomp.Tf,ADcomp.Tf,ERADcomp.Tf)
res.fix.simple <- fitted(simple)+resid(simple) -simple@pp$X %*% fixef(simple)
res.fix.no <- fitted(nocomp)+resid(nocomp) -nocomp@pp$X %*% fixef(nocomp)
res.fix.ER <- fitted(ERcomp)+resid(ERcomp) -ERcomp@pp$X %*% fixef(ERcomp)
- plot(df.lmer$sumBn,res.fix.no)
+par(mfrow=c(2,2))
+ plot(df.lmer$sumBn,res.fix.no,cex=0.1)
lines(0:15,0:15*fixef(simple)['sumBn'],col='red')
- plot(df.lmer$sumTnBn,res.fix.simple)
+ plot(df.lmer$sumTnBn,res.fix.simple,cex=0.1)
lines((-10):10,(-10):10*fixef(ERcomp)['sumTnBn'],col='red')
- plot(df.lmer$sumTfBn,res.fix.simple)
+ plot(df.lmer$sumTfBn,res.fix.simple,cex=0.1)
lines((-10):10,(-10):10*fixef(ERcomp)['sumTfBn'],col='red')
-
-lapply(c(sets.inv),run.models.for.set.all.traits,model.files,run.lmer,type.filling="species")
+ plot(df.lmer$sumTnTfBn.abs,res.fix.simple,cex=0.1)
+lines((-10):10,(-10):10*fixef(ADcomp)['sumTnTfBn.abs'],col='red')
-mclapply(c(sets.inv,sets.trop),run.models.for.set.all.traits,model.files,run.lmer,type.filling="species")
-mclapply(c(sets.inv,sets.trop),run.models.for.set.all.traits,model.files,run.lmer,type.filling="genus")
-mclapply(c(sets.inv,sets.trop),run.models.for.set.all.traits,model.files,run.lmer,type.filling="fill")
+ plot(df.lmer$Tf,res.fix.ER,cex=0.1)
+lines((-10):10,(-10):10*fixef(ERcomp.Tf)['Tf'],col='red')
-
-## # load results
-## mod.AD <- readRDS("output/lmer/lmer.LOGLIN.AD/SLA/France/AB/results.rds")
-## mod.HD <- readRDS("output/lmer/lmer.LOGLIN.HD/SLA/France/AB/results.rds")
-## mod.nocomp <- readRDS("output/lmer/lmer.LOGLIN.nocomp/SLA/France/AB/results.rds")
-## mod.simplecomp <- readRDS("output/lmer/lmer.LOGLIN.simplecomp/SLA/France/AB/results.rds")
-## # summary
-## summary(mod.AD)
-## summary(mod.HD)
-## summary(mod.nocomp)
-## summary(mod.simplecomp)
-## # predictions (back-transforming the log)
-## pred.AD <- exp(fitted(mod.AD)+var(residuals(mod.AD))/2)
-## pred.HD <- exp(fitted(mod.HD)+var(residuals(mod.HD))/2)
-## pred.nocomp <- exp(fitted(mod.nocomp)+var(residuals(mod.nocomp))/2)
-## pred.simplecomp <- exp(fitted(mod.simplecomp)+var(residuals(mod.simplecomp))/2)
+### GLMR
-## # observations
-## str(mod.AD)
-## obs <- exp(mod.AD@frame$logG)
+source("R/analysis/glmer.run.R")
-## #= R2 function
-## r2.calc <- function (obs,pred) {
-## SCR <- sum((obs-pred)^2,na.rm=TRUE)
-## SCT <- sum((obs-mean(obs))^2,na.rm=TRUE)
-## R.square <- 1-SCR/SCT
-## return(R.square)
-## }
-
-## # r2
-## r2.nocomp <- r2.calc(obs,pred.nocomp) # -0.48, Very bad model due to log and absence of competition
-## r2.AD <- r2.calc(obs,pred.AD) # 0.60 !
-## r2.HD <- r2.calc(obs,pred.HD) # -0.52
-## r2.simplecomp <- r2.calc(obs,pred.simplecomp) # -0.52
+#### TODO START TO WORK ON A WAY TO CHECK PREDICTION AGAINST OBSERVED TO SEE TRAIT EFFECT
+run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
diff --git a/R/process.data/test.tree.CWM.R b/R/process.data/test.tree.CWM.R
index e6b3509..650ce4d 100644
--- a/R/process.data/test.tree.CWM.R
+++ b/R/process.data/test.tree.CWM.R
@@ -280,8 +280,10 @@ var.name.fill <- paste(traits.name,"abs.CWM.fill",sep=".")
var.name.perc.genus <- paste(traits.name,"perc.genus",sep=".")
var.name.perc.species <- paste(traits.name,"perc.species",sep=".")
perc.non.missing.sp <- lapply(1:length(var.name.sp),fun.perc.non.missing,var.name.sp,var.name.perc.species,data)
-num.non.missing.sp <- lapply(1:length(var.name.fill),fun.num.sup.0.9.sp,var.name.fill,var.name.perc.species,data)
-num.non.missing.genus <- lapply(1:length(var.name.fill),fun.num.sup.0.9.sp,var.name.fill,var.name.perc.genus,data)
+num.non.missing.sp <- lapply(1:length(var.name.fill),fun.num.sup.0.9.sp,var.name.fill,
+ var.name.perc.species,data)
+num.non.missing.genus <- lapply(1:length(var.name.fill),fun.num.sup.0.9.sp,var.name.fill,
+ var.name.perc.genus,data)
perc.non.missing.genus <- lapply(1:length(var.name.sp),fun.perc.non.missing,var.name.genus,var.name.perc.genus,data)
perc.genus <- lapply(var.name.perc.genus,function(i,data) mean(data[[i]],na.rm=TRUE),data)
perc.species <- lapply(var.name.perc.species,function(i,data) mean(data[[i]],na.rm=TRUE),data)
@@ -366,15 +368,21 @@ mat.perc <- data.frame(lapply(mat.perc, function(x) (unlist(x))))
write.csv(mat.perc,file=file.path(filedir, "all.sites.perc.traits.csv"), row.names=FALSE)
mat.perc <- read.csv(file=file.path(filedir, "all.sites.perc.traits.csv"))
-mat.num <- mat.perc[,c(1:3,9:13)]
-mat.num[,4:8] <- mat.perc[,4:8]/mat.perc[,3]
+mat.num.g <- mat.perc[,c(1:3,9:13)]
+mat.num.g[,4:8] <- mat.perc[,9:13]/mat.perc[,3]
-names(mat.num) <- c('set','ecoregion','P obs total','P Leaf N','P Seed mass','P SLA','P Wood density','P Max height')
+names(mat.num.g) <- c('set','ecoregion','P obs total','P Leaf N','P Seed mass','P SLA','P Wood density','P Max height')
library('pander')
-pandoc.table(mat.num,caption="Number of tree radila growth observation per data sets and ecoregion.",split.tables='Inf')
+pandoc.table(mat.num.g,caption="Number of tree radial growth observation per data sets and ecoregion.",split.tables='Inf')
mat.perc <- read.csv(file=file.path(filedir, "all.sites.perc.traits.csv"))
-mat.num <- mat.perc[,c(1:8)]
-mat.num[,4:8] <- mat.perc[,4:8]/mat.perc[,3]
-pandoc.table(mat.num,caption="Number of tree radila growth observation per data sets and ecoregion.",split.tables='Inf')
+mat.num.sp <- mat.perc[,c(1:8)]
+mat.num.sp[,4:8] <- mat.perc[,4:8]/mat.perc[,3]
+pandoc.table(mat.num.sp,caption="Number of tree radila growth observation per data sets and ecoregion.",split.tables='Inf')
+
+par(mfrow=c(2,2))
+plot(mat.num.sp[,4],mat.num.g[,4])
+plot(mat.num.sp[,5],mat.num.g[,5])
+plot(mat.num.sp[,6],mat.num.g[,6])
+plot(mat.num.sp[,7],mat.num.g[,7])
diff --git a/launch_all_glmer.bash b/launch_all_glmer.bash
new file mode 100644
index 0000000..10ffefc
--- /dev/null
+++ b/launch_all_glmer.bash
@@ -0,0 +1,22 @@
+#!/bin/bash
+
+export LD_LIBRARY_PATH=/usr/lib64/R/library
+
+mkdir -p trait.workshop
+
+for site in "'BCI'" "'Canada'" "'France'" "'Fushan'" "'NVS'" "'Paracou'" "'Spain'" "'US'" "'Swiss'" "'Sweden'" "'NSW'" "'Mbaiki'" "'Luquillo'" "'Japan'"; do
+
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.1, run.glmer,type.filling='species');\n\"" > trait.workshop/glmerspecies1$site.sh
+ qsub trait.workshop/glmerspecies1$site.sh -l nodes=1:ppn=1 -N "glmerspecies1$site" -q opt32G -j oe
+
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.2, run.glmer,type.filling='species');\n\"" > trait.workshop/glmerspecies2$site.sh
+ qsub trait.workshop/glmerspecies2$site.sh -l nodes=1:ppn=1 -N "glmerspecies2$site" -q opt32G -j oe
+
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.1, run.glmer,type.filling='genus');\n\"" > trait.workshop/glmergenus1$site.sh
+ qsub trait.workshop/glmergenus1$site.sh -l nodes=1:ppn=1 -N "glmergenus1$site" -q opt32G -j oe
+
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.2, run.glmer,type.filling='genus');\n\"" > trait.workshop/glmergenus2$site.sh
+ qsub trait.workshop/glmergenus2$site.sh -l nodes=1:ppn=1 -N "glmergenus2$site" -q opt32G -j oe
+
+done
+
diff --git a/launch_all_lmer.bash b/launch_all_lmer.bash
index c809594..adbe387 100644
--- a/launch_all_lmer.bash
+++ b/launch_all_lmer.bash
@@ -6,17 +6,17 @@ mkdir -p trait.workshop
for site in "'BCI'" "'Canada'" "'France'" "'Fushan'" "'NVS'" "'Paracou'" "'Spain'" "'US'" "'Swiss'" "'Sweden'" "'NSW'" "'Mbaiki'" "'Luquillo'" "'Japan'"; do
- echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');sessionInfo();run.models.for.set.all.traits($site,model.files.lmer.1, run.lmer,type.filling='fill');\n\"" > trait.workshop/fill1$site.sh
- qsub trait.workshop/fill1$site.sh -l nodes=1:ppn=1 -N "$site-1" -q opt32G -j oe
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='species');\n\"" > trait.workshop/species1$site.sh
+ qsub trait.workshop/species1$site.sh -l nodes=1:ppn=1 -N "lmerspecies1$site" -q opt32G -j oe
- echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.2, run.lmer,type.filling='fill');\n\"" > trait.workshop/fill2$site.sh
- qsub trait.workshop/fill2$site.sh -l nodes=1:ppn=1 -N "$site-2" -q opt32G -j oe
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='species');\n\"" > trait.workshop/species2$site.sh
+ qsub trait.workshop/species2$site.sh -l nodes=1:ppn=1 -N "lmerspecies2$site" -q opt32G -j oe
- echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='fill');\n\"" > trait.workshop/fill3$site.sh
- qsub trait.workshop/fill3$site.sh -l nodes=1:ppn=1 -N "$site-3" -q opt32G -j oe
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='genus');\n\"" > trait.workshop/genus1$site.sh
+ qsub trait.workshop/genus1$site.sh -l nodes=1:ppn=1 -N "lmergenus1$site" -q opt32G -j oe
- echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='fill');\n\"" > trait.workshop/fill4$site.sh
- qsub trait.workshop/fill4$site.sh -l nodes=1:ppn=1 -N "$site-4" -q opt32G -j oe
+ echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='genus');\n\"" > trait.workshop/genus2$site.sh
+ qsub trait.workshop/genus2$site.sh -l nodes=1:ppn=1 -N "lmergenus2$site" -q opt32G -j oe
done
--
GitLab