From ae02ca27f2503c478e8001ef54a09e9ea9f8837f Mon Sep 17 00:00:00 2001 From: nicholasjclark Date: Wed, 20 Jul 2022 16:04:31 +1000 Subject: [PATCH] fixing some typos in stan code --- .DS_Store | Bin 0 -> 20484 bytes NEON_manuscript/.DS_Store | Bin 0 -> 14340 bytes NEON_manuscript/Case studies/.DS_Store | Bin 0 -> 10244 bytes .../mvgam_case_study2_files/.DS_Store | Bin 0 -> 8196 bytes .../Case studies/rsconnect/.DS_Store | Bin 0 -> 8196 bytes .../rsconnect/documents/.DS_Store | Bin 0 -> 10244 bytes .../documents/mvgam_case_study1.Rmd/.DS_Store | Bin 0 -> 8196 bytes .../mvgam_case_study1.Rmd/rpubs.com/.DS_Store | Bin 0 -> 8196 bytes .../documents/mvgam_case_study2.Rmd/.DS_Store | Bin 0 -> 8196 bytes .../mvgam_case_study2.Rmd/rpubs.com/.DS_Store | Bin 0 -> 8196 bytes NEON_manuscript/Figures/.DS_Store | Bin 0 -> 12292 bytes NEON_manuscript/Manuscript/.DS_Store | Bin 0 -> 10244 bytes .../Appendix_S5_IxodesJAGS_Hyp3.txt | 314 +++++++++--------- NEON_manuscript/next_todo.R | 251 +++++++++++++- R/add_stan_data.R | 2 +- R/add_trend_lines.R | 2 +- R/plot.mvgam.R | 4 +- R/plot_mvgam_fc.R | 35 +- R/plot_mvgam_resids.R | 46 +-- R/plot_mvgam_trend.R | 34 +- R/ppc.mvgam.R | 137 ++++++-- R/sim_mvgam.R | 2 - man/plot_mvgam_fc.Rd | 10 +- man/plot_mvgam_trend.Rd | 8 +- man/ppc.mvgam.Rd | 21 +- 25 files changed, 625 insertions(+), 241 deletions(-) create mode 100644 .DS_Store create mode 100644 NEON_manuscript/.DS_Store create mode 100644 NEON_manuscript/Case studies/.DS_Store create mode 100644 NEON_manuscript/Case studies/mvgam_case_study2_files/.DS_Store create mode 100644 NEON_manuscript/Case studies/rsconnect/.DS_Store create mode 100644 NEON_manuscript/Case studies/rsconnect/documents/.DS_Store create mode 100644 NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study1.Rmd/.DS_Store create mode 100644 NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study1.Rmd/rpubs.com/.DS_Store create mode 100644 NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/.DS_Store create mode 100644 NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/rpubs.com/.DS_Store create mode 100644 NEON_manuscript/Figures/.DS_Store create mode 100644 NEON_manuscript/Manuscript/.DS_Store diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..62ad0084f98a2516c3529def56b3e6759a7bba2b GIT binary patch literal 20484 zcmeHP33Oal8NPqIz4!7;UnsPbX(pL&DQQERG)t3~wrQ5O=|XGM6)Nc@GqfXiVOPv_r7^~@6AjXthR?U@0`h< zx&ON}@7?+D{l4%2|A`QSI!2-mLaY)(Xn1gInL@Oe3K2$Ik4N@kJ3n%cG3X&YX5%py zCR?G9{B(>9%hLt5tG4*JI0MBQ$ejTLU6aJcqFb!T-##%YA|fe9 zL>gLtSj5GU$n5DIP9{>xT{x2$ABl}UiPV6t@v+?KFX+AT-}1mSc}D-<^pJhVclw7C z@pPuX{&6Z9Gj`ng37TeT#zfG8*Xho zR2x3nP~Xy28(vu7bm)+-O`YAaq-W>g{sRZEJoK9P$%_|c?=YAa?{Sm5wb1}a$(ORX z_GI?fbfeUq7+lq}YPGqx&pLMxu6q4YCZ0;Bx5S3iX#1>5Y20bgKq|RmIL7v^1F6Kw zU^3ljl?_A^10#t@Cbm40=#3wUp>I-pU^taXY)GZ!{5IRSWn%j>Y|AliyCj_&&UDH* z>$9wpbZlb`?-kz@>%u#xx5e?1JkTQj$T?Fhrq7-`zo~6$=dnpv*}0R?GtEmdX~|6d z@_1}`^WJzgvwM3aJrGMqS{4tHEWH&b&r}^)w$SQqLoptzHJ~iIEed6 ztM2w-hh{yfUB;iiHj>V)ABrVMuWE&+ZD-5!XfzhpW@f+kneDaespnSoneFZN8JY-* z1)@u=6+5w54vKe(8^sav1#y@7qPR~yEFKq6i(iXBh-bxL#oxpWR6?aRk6rggN zMio>;3#gHrsGWLfHLay}w2>~Q9dsG>QIrNLNkjB%+DnJ&DtbM=f!VnoVvb}tS7Ituh(wqMLZ#~5uR)>(vn;9dvevv zRV#bUHD|9UKXHz7-}9$T4bQH*aKXZhmOJ0P@5}aliEby9ey1EuI$YMFq_-b~lTWJ% znN^$Nim`o#R)#f!olcWyPoGg$WwE!Zm!;KJCYFOPZF8$(E3qwg^I;hS&sVm_CbP;K z2PaqKTH8$18jDq+tV@={Iv)P69h+7j449o}x6Vt^5~s;w9}QiTA#Zu6{+zg}dK3(v0 zL-sy@?7)v+THdpI4WK7y|35Y0abU&;Dh`6=-41m9h_7%7f(<*-`%GrU7Lf#8ID79-gji)=|&umX`AJ@N6I2 z$crkknTr7tg*NIk>*ni#8-+I9Q?55QX}T4)(MEN(wxMfHc5_f4p)Xmg;p%IhJ;(m! zj!q35fNulJUe!u3+GIgS&VYk077LSh=d>j5zHRCV#ukfX=U+dBJf3x=N|@}OQ3QUKy#2U$qFK!IIFTY>dm?}mdTT}w(5&Q)_C}> zvM$k=0Wf&IDRW1tTDJ%^*`YdKFD&ST9};(nFN^!dL-4&Hf#oo8{0)Amf!$1z?52W3 zR81Fv;mo62Y6gh3!v7WnMpn>D$#!}HB3o!1Z3l=%DMq_!H@HrkGPFlBp4Wl#ya{0P zPWa+`=)GV*AE4WOzSpB(XaCYCbA8VT8`)EfZE*&QGaxw6KUPGc`ZFTp)%9XAY!|q6 zy>p(6Yn=`1zLDuwLt{O2t+pm8{g&ojs&z8z-nl^+$=aWbS4L?7cxB|JL8p~w<&Zem z?&kTyDjTUJqbM^b8!k3g2w>#Kj-CHH{n*4w=UBj*fWxGuc4r)%76=(0^aASOd|1MP zP{0jamPtu);?7oMQNSTm z3Jtk?NnlyP$D|}QvaJJG+D&F?d~1H_)w8|FCSo3+m@vYWHzy)N1c;MEAQ_rc0{lAy7tXuzSVIGvtsS z1~o_0f9G5DeQzGua(mP9CAefzMnx#gmHaNQn3W{r)&xm=g2CxCm_S*I_87rxokfCe ztDXy6iF{!M!TsU8))-|o6!EvlDXY=i#>&7NYg?J_TWhTGwzZ>^waaxZ!Sl2D!hZK* zQ24LF=N=P37C(i*{SLnN0{pC$%A}u_BREtEPg4ZyLTa%IRHxM8D`^9C_(tgPEeI4T zpNk<_#58^w1d3JoYh}oYRrodZE=iznp%2p$pV#FhP|h~!qfVs=b#;qg(Fx_fMbC%92#4A3Pd2F^;%Xl;;~v7oIs+}zfdBhOC;Gwj=QAx0iqZm=vj*18n* zsYHSEiSfK*DNnrF5tO;|+%QgBo?rb^<+*XH^4!Rg=Z29Z&+S#}mgiWhx$@jN-SQk_ zNg~M%(IhUy3fK&knmMkv1>s*7LzhDgC4nIOXg@II^)i~q>zD(-ACb|tyAb&O z2Hi{FrH2vtee$Ga`KYrVM~)k%kmE+j8Ff|ga{ioePzt_9Xo32&QDvuHNeQ9U)t zbi)!vblc%~oiY%_FF0_BCSi!cFuo8jYz9Z2-3YWA=rLk zgLE!OmN+?W6$WrvA=pge6n6fMU=5N9Z4y5^WiWR=vKpS$-gg8gcRB;>d4ENz1!RpOD;smup z-qG^@Phqr66F!*TZ_ZTCcIv@{O-;~&EeiTMq}be;TO@QN<_AZ?j4eyX6U>=QDi$ap z3aA32m&m|#DuZ5n(NNwHr?YwF))ATKnNtewi>??Y3Kuk8NlMPW1>P+ z7gar&ANsssA_r)T`oF0E2T}BxlgT*2&QjF>Mg4z9Durb(v8eyE&{b?-#xuY_q2UMr zgvQZ3E_Ql%mXmv904Q%?y~bRZcLzVU^!I6Gst#BFuglvZp$nC_Gg%E)2J_HTw#=vu z2HiPI6VzQtn>X9t`?Y zt&*0qtqWJ>rKP7QIB=r+_Stswase%s!Ga!C9IB+{rAS!z+ez{?jo8&*DnP*X?H`~I zAyLWo?RV2X^fg)E{wO`>)VKeUo}=g0vUCXQ#p-qTFZ~_~S-b?rvnop$mpY5xcpozEuagC| zT(qci?$@J$_Ezc6A4MVUo%C7y94g@Nfm45zxS*CR;2)=7s>^oVIi~{tRMN$#aaTQo z)8+h*sHp#)UlE}e=>)zgbIE5+|KsAus|#GkzyI$p{QrNlw^3}1Gf~{1f7rW&Hm?eHv1r literal 0 HcmV?d00001 diff --git a/NEON_manuscript/.DS_Store b/NEON_manuscript/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..eaf6690acf5b0a949262f677ed328c16df7e9293 GIT binary patch literal 14340 zcmeHN33Oal8UDXCZQi`THZP@}q|+?ZrfIsQS*Karrky5fn{=VI>6W(8PBKF~GMNdP zNkR#MfTAJ=7X(>EkVOv`6_FMd&?B1#Jt&7No1mh|aYIxNh|9tM-a9jyd6|?Hd!%~i zo%7zEd;fdi%)9g7`+eX4KO#a?pSOyrl87`uMQIvDnL3e^%fv(2#X0Jjz%Af2k5B$= z678ox>LnkA1L0A<(+%kvNYB9OGy@s@W5oK+8+ZD zjs-T;qiqfIET~(Q8Ra%JY}GTEE-D-_`&0KHW(XSp)H7+w6 zwH%b$c6NojyZLvARkPdCq?yOGv&DCB@Ps3qdVIl=U$sWlI)r7F*X#3Y#qPe>0Ww}zC%yY)AR!Ul73Hb(p&U*`UeO$ zOv7}{fE`6}Vh)N?i3TjiGOWdBd;sm(iVk$53qG8Ka}hxw25>$u!bkBjd>mKcYJ3hi z;TC)WU&c3ZFOK3qJcb|P$9M`a;AeORui-8H6aP}bHA6VDIL=f>?2301=SKdE?!1(9 zBW;Cwmbr=DCq6gwU$e!1OU_I+wa^CI$FMg@SJHKKm}llF-A51859t|tkzS=g(HlHB zf1|gNfeDz1NytVXW?~ldu>jSmfeUW5VI4MLBeq~S_TX&ng%{ljq6g>k>e5mcy!)LwSrG%B31A zU`E;HaBg=->+WbdtV%VUpoSBo-;Ah2ff`P5%I*ZWD%02W45VjZ^bAbprqSHLC(8YA zNzMI>B=?_p+H?QsrV7G;#`GEXTxZF`@|ya_rZ`u&(BmWqpGey!=TlT%QuG)FSJHrG z!L0lOqi7qGci;I|IT4c?W|pAo+Hui)_MD<3Q}AhIVkIjm292F8bMi`Ti8lM!`R-PVL$x%pd|MFlGqFCdKo^#^!}6hG(ICK z?9EK>Z^P}l17Bqddyk~~4>FZ~2oK{CJc=jr41R%El*5imWPdZ-d;8magMmBT_ode;hhtsi`hd&cjx0%*RS&oSmC5RNPUgRwbzHZ~`$`0Z=fw(2Rh_Q@&#xXw>d^$8ovD0y$f zcG#HVC@P-6xVmons^e35P3zP2>9MtCU31A2J|`;S#VC8OeWu3JJuc(bz4Gmaqb0oT zpKqUIk?``%OYEg15?;0})XPRByc}AgS0zYz+2qn|lO()sY|xikB)n{1VRxsJ@Y0rV zZ?Q=DsB)J(Qe2$ccu~)!awcxAw1-2oJ#>(D`&IM_x{+?DyXYQzfWA-1=m)G*#S(vk zwdzZ(R$r$-(_iTg06HeaV5K^nHT+zZV4+m1wWz~#td>i^4V!q`cgSV0)arTY#USg{ z%cPe7tkm$gN~L-f_u~=egz@FPF2O1JBOc{T$_2~L8H1#=Q#ohB>NRWMC+B?55dVrq znZn}IMU_jIt~|vUW6u@1p%;t~V?^&FXF)*}V?@<9w?LN&BP#RfvP~e2sFW>ci$WMt zt*X`wq6j0ZwRO54MHo?Cww$dBVMNuvimi+=ITB$`(g$zlILWu!|DL92+3SAA@uS!2 zZF+~JBiYF1*hoGKP>52Fk1RqN7NZvRGKSPFJx|1uwr~__8@6+t#KVy!FUOKZtfUuV zjwcP_Vvd(w%CK_?SK(6(I#)C7+=N?k*y4e077x5@ED!7swX`TwB!tTF6_Xx%c-7c_ zz2olNb6|Yivjx_&Xp2a*FsubALg55Jdpx#2-Y-DQ8yu{0l{u?wT>{XmDqXb#-0DUl znm{Z*#zeFm*l2moabT3*mH;Lv6Om4ZOYb4tY=x?vX}X>ONmjIqp_tseb23aObLQ0S zDXF1qR=!JiS&lNlOuwdA=rwu^nH)vUf`%NqPY6P~ z7!@36uHrcJN~}T?nwgf4>=7fBbV$ZhFXbrn4Km97d3*_9j}uU)pVBMogcvp2El-F} z4hE;v_thgWo*B-b&sB=$I*jyZzl zmii5$lT%#Mbn6uO)<9V5|vk?;88 zxX{Vsdlp{u89937=OQBR*$3m&D`aQwQ zm=BynEREfFL@?bQ!`$Q)XpEk@DG*JkFkKqohlH@Q=0?w$ZW!Z6D?}KrRBkkH zR4(`L<$$^oT$3As_(l9U^jLz8o(vdNkom5rlwqp~@j|Ks&8^M5Du z>lm6{liII~*0!!~GuEHJew`HeTOi{Grx@=3e7ej#&)+Js#7-#xgB=%qC)U`Qy^n=hCqtb`T8AXmUC2-;lKPi3**b(PIdOd*ldcrPJYALiu#5quNh;^4poGV1p&6Mk}^)#308TgG%fFXec^`L#2#> zjdG^6{#S3MSz^>_{V(1%OY8r%{!e@fGOhp9`ahXQiP}WfR4Ly7TRrys|EgE|nx295 z4E$ft05V(JTblUMVlNk1pHjSQ=j3xApW+bxhE*wt@Fn^>p8OzQ#jEEtez4swDqIe$ hQV!#~8unj%36S5BKQ-a*wEj=)f8oh|{EzDYe*>0q!-W6< literal 0 HcmV?d00001 diff --git a/NEON_manuscript/Case studies/.DS_Store b/NEON_manuscript/Case studies/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..aa69cddcecbbeee60b7e97d55540edd822209361 GIT binary patch literal 10244 zcmeHMYit}>6+Yj_ab_k-echzVuF1BOIB}EM33eREc8JMl?KE|hl2AK|)1=w#W0Fj} z-gS3o?KrN35&Wq76rvB5PzjX+f)pWvKmvqH1qcM?)uw7GDiT6~Xn(-Jsue{&_ug4Q zHVzOf)PiQLx#!+{?%A0+^PTgZJ0>C|%egipn}`fPr4&oJndgvtd0uuzf*Gj<+Y=2_ zin5d@kDSX^hZCU(LJx!<2t5#bAoRfd!2`_M@-P-_8MdJZLJx!q4tCIUI7 zrAi&#s4W1?8lF}M_jwvy%7iqjZY*?X&mqDc;5VRHht_(-z*>pCp&^ zys!5%W8O$6U(B%OPctKV*YjIjUqWO-&BEG6hG7}jQtOU^)5Wm?Kkb*jWV&=zjt-># z)0v<@SS;qJ`{~Rf*BMKR<-Ns%pLPpQNg2C2hnE~X>=tsxv1GAa$a(5_2DHL3Ov6lx z*u+FfXJ?Du*48%JVo$WSb+)(Io!hoePMXHbb!|QUCx*w*PTVtj?!10<1nGn7s$6m5 z>iv-SMtij%BI6Z5yj<~vX)TW4wC83kC9WOi$J;;RyTyWc$SHX&Q=-mem-J_gg##r= z>JMj&`SNhVONshyI-f1))4sDSpC53~I^K})obhGvcCT3Sd)2}rmu=--N{F)O+~%-> zJL>eYnI~XgTykB*4UOwJZRy;x>(bI|>O{S@%(|UDSMc4Fu2VXa_OebP=N3*K9CLI2 z>BFP0=VtQGU?Dy1nDvo!=B1@|*N7D>qtT&q#=X1j`s1R}Fw}2-AR4!h<@ids7IX5~ za)$UOuji}S(RlQz)by4dZgYXb+GvwC=<&@OBCQv0u$ryGyZF_b!qYRGt=8zld^+Rg zO_9^rX}3C}hex%?OtDu_c3In_L)qf+Fh7cC)}uvF_e6I_kIO6XOMCwQ5vL&U>L*&4 zP@E@Ru=wZ4*Ptn)t+w=_m zn4Y7b)34~)^a}lf{z!kKKhxjnZTbfYfQ1!Ui6~a12{&Q`HliI}=*EYz1G}*o`|uI$ z#}G~+jVy8)#wf;c78Cd+K8;847#_#x@D*Id)A%aBh412r`qm;9*xsVYl@>qZ-c^P# z!!`)t$np7L+pznlJ?{~0f3rkk_1b02ue-kS#^%J<&M&=9jGX0KG z`8p%=O-AEl)GIW`unN~R9@jG>o3RPo(8;LWj$ZUJD)-<3K8o8gh(kyP$aGLZ5hEzU z$0W|-WB52efd}v)9>S*-QlG$+xLAeF{t7C!#M4)g&C!8UAzv(E5cyr?q=ai`~)npOMcz&pL}>9-l8#WBs>S`(`o0yCs!*xMW`X zn#q{VGx0S|Dp30r16bCw;}By1HeIRx(I1DFc_&!*pWO z)3Ic)%o0^uD$v<`3}Y4p&E+3c$h}XD*3eekK?mpvS6pLEJde_2Bt!cf^ga3k{ephU zv?D|O_w*O~E4@i?(ckH9)S?c81AQ56ML5kI+zISLl7qXS!+Jmk_E9Dq2Pbiw1A7c- zFwT@S6VT`Jc|3(L;7d$7&*1xb7C*ytcoDDQw>p%k>E~8GuC#b{{viEV6ccBLNxrAR zox?P`PMUt9Rqt>%oc<}f%*A7cX*8qE;CXJrB^K!#My|dgzB(=>h&0ozjW@+bgmJEi z8{*BDSfDB>xtQyvxLBxaB{}3;$r3fHnvz4VlH#IP6-{#3W9_sA$(S=0-6Z2pmH=1K zRY<%d4!=xar)RnFd7kn2DwjR4(HmUw)H3o`FtylR@T|o;Y(^{Exa8T29_(ZkCYfS% z$s-BoDDJ>UDs6v5n!`*1(g%ft8#&f@|;izo2Ke~V(~Nt{(QQh{qI zLl1-=nDM{@-L>WS|J_&q|No3Z*b6-ndf@-g1E}d6=u7frRlYJ*KVgI?Y(7u$DG#}B zN=s$Jja(w<)#rE~$nT2lgSl2I<-z(yAg8oc<~*;=S?aI(|Mh2~&-|A; z^Z)bBO#WRWLbhb2i5iJW<0`3M&&_R$l*@WciG`dxreJ%bBAK+BM#!hpTeKra6oDuL zQ3Rq0L=lK0a3@57-)vd58t=Z6Mr{;pZ?U-Gi~;sH6@oA*ZwVZX<9?Q8vfVV`AI2KD)!uH)w|$1DcMmSM8W!9L3| z+)CChIffU^lT!js6Pg&*+-DGb95o*&F?E)o>j2T{dsT5 zbPUTG?yXpcKQh39bA0QlWfp~Kh@I+j@}uVEg^L>X_~Q7Is+H2cZma1{T7Riv9V=P> z3H6P5a;ebt^9K!H{chKgGa7*;9#1aUv;&ggSv0v-Y!P31KO5wGJoY49=9;xa=w8p` zyt~R&tb0%sy3uRr3#J{iHf&_p&c3p8H)L+!qKO8i`iI=nQTE4E?0B|Q)AWmn<+pd` zJ%7)b>0Dpcc1_zaB^`!g8d_5(r8Yl^KvO28M!I}z;*_E-w2cO6h#acW3-ltrM5pO2 zoul*gF?~*7&_()^zNRbmBmGRj(r@%T{Xu_WIy4B>qY+E73{6;z4cLfwY(fWiU?+Cr zQRHwChj182a1<`aa10*C@id;nvv?k_;0)fxTX+W_;6r?b3-|<|;wyZEAJhtCRq*as z<7$hm+^ca|8NA_PjUB!ncoW;UKa%}lzM4&`F*Vt8 z2wm8Q?dWDS=Fp3M=*NB(V8FyMMglDQD5HW&Jc*|mk?ZM~jIf5JlkMjR5Mpa$Q-T zv(~LGbL}M8BV1*XdXtIEGtpuqQuq;4pT*TAMd-zILrFvUUgx3gb5G$&-F zjO!SI7=arR5Zgy3lRPR=`P%&5AIdvkp0%IkhaB4vnwq{tVp`4g88a0{Q?xnSq1>QX z$puDG_OnK*U-Whx!C*dW_jsOjrETQ**;ZvhpWo@ZfnmE=DRj0?i%kyp+OFwUvR>IW z{cxU~R49t7r~`U(WTd5~tufuaxplNLJ<{B?c|&8mZPQ(&qpGr~p}D>LNb$tUkyE3m z&&lu*9T8MJA>8wm$LjdQoif(M^VL`{Rby4PS;_AEv>knR!Eqm?W*BR@zFtX+M*`LTBk&dX8SEae9Z| zr;q3h`jWn)Z|D;JM8D8)^auS(f6?EVfm-OO!#t$18f&l?8<9Z=vgkrL_8^BI>_b0H zSU8G7xTxSbPGAI&;VhoSQ+NjF@CshVYj^|i;RAe#3-}bD;amKGALTj|RpRcH{c4TN z+?r&QGI66#nm9ZeaZ~qh*RuaBalfA}*jhJt-uwm2>hD;$VdK`0D<5>iH@uckA%i4E zLh4kh5q`+)`N+IZF42~z_A5xXFHdun+~7bKG1m8n%JJ;6Me z-PNhJnm&zLtcJ$xwWgFloe3+oEn1tV*D#}{c9XUxrO#mQN^QHgP18y2IjPR3R%tr; z5SzfnDS7xNy-gSBbNZV3cbR^sU+H%cb5F-YEMw*^$4WF{J(|#r7PO)r+n9w}bYce{ zz%D$4f zrP#z5qp0!0OF&<|#^94i-;5GpjQ1!;)bQYg5BlPx(P-j-=FAo>5JHTxnsbsl|D6AF z=A3W7GiU#^0D$R|*$kiqfI=6SY(7<2DcsJlD@rWVa3fMYKmoFl2LpVtuJD!$iGV~v zA|Mfv2uK8G0s^#W^P_IvYNSyKbmvA97TPY8+B|yNG89= z*sPk_bLJ|Frf3VaLzxk;k_n8UfXrp*;XaX7I%7XVA!rz z41;acBA0{xwrhHov{!OXKkO$X42q&EYL>;v$6H$38ujLFtrLyunPg zs44g4Fr2wD3>?4QDKLmVQ^nvy6$4e97w^7H+dW|C zJvUhAusn0BOWzgzc;AwaP&cd7U4!jSa!e{U~du?e3xE?p)M%;oa z?7%d3VK?r>4EACl4xxz_9>EcGv4Y2N49D>yJcWLQGxA9$k56|I8 z_%VKk-{N;-oUtmQ_lS12#*1`!1B;pn9WCnE@YM*NxO115{#Ob8&3w+{x`m4tFIm~J zZqw#1J36Lb^t5+-In6?53F?HDs!}F=!GGl`@jbp=TanmL%(ISN;X{?5rF05O%?FlxA zI4!L0+D?rDUyH(=POQ}!(noI^L$Art*We8}2cN(f#LJ8D1N;QPAP^@RF2$9^$a-9j zYp@BE*o-aMitV_Q*qO#o+>Q6(Uc8qWI)sOC7!RX?_o0IY^l&uf=LtNC58}i42tGtrnvnCS@25(Cuv`Trf){{MfbPeYDRA|Mg?8wg-c zSEegX9aPT~ruW6q+I6}P)5Qz-n-nq(p^Ce3eRv#C4Sue;F485rE3(;aCQuq;4pT*TAMd-zILrFvqoj12F0g)@B6a&(o^=Vx)Y)|Bg#+^zMcKjp*bNV zWn9My#0cDofY?4Nm1vmCl)pBA_lNV2muKxK`C-TQgXZS%keF68ea1{h(G+ctb}%>O zRdRt5l>Mwx8W6o*Mlh6*+P$9VTxlEmy|z^u)aQ43ZeZB1RSKPL(_)hYeYR_Qm8@5G zO+TC`Cl!jKD(aw~935?KZEsAsY-$^8OpmrSZ(84&Zr^y<*qEv;TG`UkbGUf?#OTSf zQ)gv(h>i%Voe=K1$zyf=kuDi);`wT8_FZOe_+c$6RzBuQ$T3yrT&F>qreLL@1`wYKe zxu)$N-BYp6U}!&s=LYr>+bXGQU1D9YYZR@E3l}ZclS`6It6tLfs;{CqDE;NU{Ycpk zj>~r>)Ag#Nzk9&s74Gm%F{2r3lF9T+MHvv|T_uZaIjkjDuVDwH*XMYnWj5WZs=D9% zFn?j`sus;x)X1_wWHe#07kc&+sjNz>jjBi7Ih-$bPlP zWo}KfNtw9OCQTfgjJT4D) zWp{OIt)@?77OSCglh&Nlr!!%twpDA_^crTg)Na%^r}P=jU8(KRwrVwE0jVG9mPvZsV<2k&Dm+(5yPvBoqm3dO) zt;zW}k}J86=N@HkmAThD#=V;wfoXDei$(r_&(y#F-_$Y0i;EG65%_lqpr$+5o#i>E zI>kk~;n^I-RLIvt6JR zoA_cBH9q(g&=-GW@JXX@Mu{)R|0qV(@Zf_F`r@O}XyUnZXA3RVN0l%)x%a#Go^xjI zIrGiEduE9U*^=2vw3vuAt`fy8Zmv_LT-Iw!BATgU1ltq&g|a_XaJ<4by75FDfj9zj z1mXz95r`x3|A+ve*|KO;ocl@|w{ZmG2;7trko!ZLtHex%V+zBqgPSk~V3a99Sh&x5 zKxmVWFcaaJ!qAkjsj>$|rifAuD07ylr8;9K!ZC%)oS@7J(VY=xC6T4*jE?HsqLq!Uoks>woE$kdditz#3(?j< z)sx&jH_=DOA8uFQC=$wsl;bb4#@_jCu3Bhi%_3IM)8W-*365?H=1( zw#{JR00+$t?4!0-)b*Ogx-K_AXkA>mXt78xNiMBeDcvi!nyAxyN(K9&k{z5-_eiE| zbxpi;(Bu{F^h|k0Gn6Eg>6Mz+E9cvb7Po5JNUmPP2Kg?Z{f(Bn{&rm#e%FKiWR~F+RRh;b`LAJqsi7gHN6J)J^kL`Ap7I1>14KD)5J%I<-K?0{b1LS z_~>1Xzbg`*gN zi!zSm1V-=(&fp0=iKlTEFXI)wir4Wj-p2>HfKTu#zQOnSL9H`UA?!{yuC%z!y{leI zCT#RkEphmIgiYPO-N^o@g#B)|WNFRZdGi-6tGjL8hK*a>#vkOQXLmJ?LgpxHgp{dJ zB7A@odB}WDE-{v-_A;}q<5zgBcwJ$+m?=4=N%kjiP1UA^B##=dNUchV1oK-BSEtq* zVj8no50%#&4Jk35iK?VcMvEb;n8`}I$=H$-GnlhV+G=bw1j%r&MbTi&od9t;RU>i*KmFk?>Z}tQ;0Vw@@_a+ zbREw<#rA!F%UnCn z^$1s4q~4gq(1e>X9cK*FamMjK47s0XRcVuna7-c3aQW9i1jOh6`0vj6{0}L3Q#OAC D42j)( literal 0 HcmV?d00001 diff --git a/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/.DS_Store b/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..b4aee0ce0ed1e09f28641b139c9a7116759ca70b GIT binary patch literal 8196 zcmeHMU2GIp6h7ZlU}gY2<$no_TM9&=)RwmL6Dr$nl~h5HEwoUe%E3(;aCQuq;4pT*TAMd-zILrFvqoj12F0g)@B6a&(o^=Vx)Y)|Bg#+^zMcKjp*bNV zWn9My#0cDofY?4Nm1vmCl)pBA_lNV2muKxK`C-TQgXZS%keF68ea1{h(G+ctb}%>O zRdRt5l>Mwx8W6o*Mlh6*+P$9VTxlEmy|z^u)aQ43ZeZB1RSKPL(_)hYeYR_Qm8@5G zO+TC`Cl!jKD(aw~935?KZEsAsY-$^8OpmrSZ(84&Zr^y<*qEv;TG`UkbGUf?#OTSf zQ)gv(h>i%Voe=K1$zyf=kuDi);`wT8_FZOe_+c$6RzBuQ$T3yrT&F>qreLL@1`wYKe zxu)$N-BYp6U}!&s=LYr>+bXGQU1D9YYZR@E3l}ZclS`6It6tLfs;{CqDE;NU{Ycpk zj>~r>)Ag#Nzk9&s74Gm%F{2r3lF9T+MHvv|T_uZaIjkjDuVDwH*XMYnWj5WZs=D9% zFn?j`sus;x)X1_wWHe#07kc&+sjNz>jjBi7Ih-$bPlP zWo}KfNtw9OCQTfgjJT4D) zWp{OIt)@?77OSCglh&Nlr!!%twpDA_^crTg)Na%^r}P=jU8(KRwrVwE0jVG9mPvZsV<2k&Dm+(5yPvBoqm3dO) zt;zW}k}J86=N@HkmAThD#=V;wfoXDei$(r_&(y#F-_$Y0i;EG65%_lqpr$+5o#i>E zI>kk~;n^I-RLIvt6JR zoA_cBH9q(g&=-GW@JXX@Mu{)R|0qV(@Zf_F`r@O}XyUnZXA3RVN0l%)x%a#Go^xjI zIrGiEduE9U*^=2vw3vuAt`fy8Zmv_LT-Iw!BATgU1ltq&g|a_XaJ<4by75FDfj9zj z1mXz95r`x3|A+ve*|KO;ocl@|w{ZmG2;7trko!ZLtHex%V+zBqgPSk~V3a99Sh&x5 zKxmVWFcaaJ!qAkjsj>$|rifAuD07ylr8;9K!ZC%)oS@7J(VY=xC6T4*jE?HsqLq!Uoks>woE$kdditz#3(?j< z)sx&jH_=DOA8uFQC=$wsl;bb4#@_jCu3Bhi%_3IM)8W-*365?H=1( zw#{JR00+$t?4!0-)b*Ogx-K_AXkA>mXt78xNiMBeDcvi!nyAxyN(K9&k{z5-_eiE| zbxpi;(Bu{F^h|k0Gn6Eg>6Mz+E9cvb7Po5JNUmPP2Kg?Z{f(Bn{&rm#e%FKiWR~F+RRh;b`LAJqsi7gHN6J)J^kL`Ap7I1>14KD)5J%I<-K?0{b1LS z_~>1Xzbg`*gN zi!zSm1V-=(&fp0=iKlTEFXI)wir4Wj-p2>HfKTu#zQOnSL9H`UA?!{yuC%z!y{leI zCT#RkEphmIgiYPO-N^o@g#B)|WNFRZdGi-6tGjL8hK*a>#vkOQXLmJ?LgpxHgp{dJ zB7A@odB}WDE-{v-_A;}q<5zgBcwJ$+m?=4=N%kjiP1UA^B##=dNUchV1oK-BSEtq* zVj8no50%#&4Jk35iK?VcMvEb;n8`}I$=H$-GnlhV+G=bw1j%r&MbTi&od9t;RU>i*KmFk?>Z}tQ;0Vw@@_a+ zbREw<#rA!F%UnCn z^$1s4q~4gq(1e>X9cK*FamMjK47s0XRcVuna7-c3aQW9i1jOh6`0vj6{0}L3Q#OAC D42j)( literal 0 HcmV?d00001 diff --git a/NEON_manuscript/Figures/.DS_Store b/NEON_manuscript/Figures/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..f3a3eff10beb84395b4303a250794201ba0dd758 GIT binary patch literal 12292 zcmeHNYitzP6+Yi$U}ikVV*;4nHTHr512JG@V?%hY9}o;Ckl28E1hea1V4L->v$M7d z4vEt=ZK48E`XrB5X_~5PNt*sh(l)6SwdsSTRUg7lO4LSG`cz+u+9ndJ6ngHxv&-%p z6Od>i?T$2a=iYPg%-pl*-t&FuED@o(KU7OJkBBrrMQNEt1qP9?kSiI7kc)HFF^*fn zCyyh(w}|#rKP4zkKI$9QJKvDkKwbk^lLlP;F=G7|@L5RfD3Y~)V_dVN_CQN88V?e= zP6zv<5j$B^^BP>^3dT>EsOd(bS>)Z=zUhX|X16tYAQB$j+?R~RWA?Ui!sfEuDz*>A z2itn$u`P+PkZ_* zU{c6NiMn>$@kFvkKCCOk&sV*yemi_qm>(J$2)FX1Z43H|Nu?FD=Pg{)xUy;8NU>Ei z*_>id+ddGnBf)5R>tH05Jg_5R_k?4iNNoSkNGKdjM)pO*3B3g0(K};--tfqDw`a1q z%v(OzD9x(Wy845WTl*u)+eRi^CB?p(UUPSdzhF~5B%T(sp3!tEnKRSt-6bR~i7;P7 znkJj)FYub3Hc!uP9Sz3JT=4FtncE2OJo%a)to?E^~N)<&ms z)oRn*)f4aS<*C>usCJs0TfE*q;?3;=JGr?p9OJJ#>Ae1W{;C~9))WefZ_+CLJBQgU zh${WnvZP8?5J{a>`u(b|vu~I|OQ@CFX)ka7Lv)mm(--JzdX~OHFVc@`gkGcH(VO%Z z{geJhXW)VX3m$ke1Ldef4OU?d{Aj`kY{PbRA%J~2fCxSW8wYU+L%1FH;Nv)oW4Is3 z@hBd{7x6fr#IyJop2PF_9!}vGcm=QGH~2IDg1_Q#cn4=)nybh)O?{|Kn7TM_R7L9Q zD!rUjD?||{p%~NQ<_VTzK9-;k^*p^co?bgP;U?_DZtOufLg+;deYh2a7{*<=8y~?D zd>Ws@gLsH1_ykYzDSRC-;3Qtc%NXG)zK%cQExe6))#s-tbwg@8s>E3?i>QfuDM;P4 zkK)oV_fUXrN{V&Oaz`Orj=#@-QofegKwbm?ml_cJgO5*9`kiuEm2&Igi+~V9E+B+F z_?~%yZ1Xv#-zkSxDTneHIh0E^6b{9J8qV!L(Va`bQx2k1p zrKxIF*?f)Jf#s4_v$|}N&K$rRr=+!JoyKV28n2`c<}#i6fK{NRE6i0I!+vYRh}Z3z zzQ%8w775y-tKL;*h5BWK@-NdVdV~JXApQ)!iwP)Xj60L@?P4sK@Vy1CSclEn!U(qm zTP4B`Aee!2ZHDoOaVKNl4>OoQBEkFvjC3Ey=kR$3^3@aa zE5l8r-A&nsj6aNR#RDCQSTr8n&pm|eu&9j=G8ztXl*&FIG7b*fj@8r2l-0kQV!WE6 z8AZm;9ZYIFl7VEu-5f~l616RXG+ooXt+ZrfnAUj~++=JXj1zS|l_ z+|#69FnDw7j9GK9tyy~An%0pbxugvvgQ@3|_|#gI%eBN+$jf0`>GUGE$6b~ZcMIsY?bgB$X9zcjj)v6k0S_&2t;)7`$kyUT7HS7`$kyU1GR17`$j|Xf)iZ z{Y5f((YRuz;ZAQe0Gh!=ecVx^M%gVEz zk?igB%yOKOtj_UNInGE_jRhWdThbXx>0+bSaYj;FZ!|d0NJ^I*D;#GerK=6UE*76f zV)eO=IUZd)qf4~P$LIt-M?at+u{FLzuhVZij&ugVM3hKN6hRgdMVf~NxE2f1fJU~- z>#!1Q*&VItu**htvO5xSrJZb-A!(NpBsjLzj{&yNyVxDwgOB1~$38#B?ug?{Io8;w z>ZvQ9xEO29aXzWR`7*D8yaw_b$ZH_4fh(+m@%;V^V@x{q5DPCbRII;pUz|`iyTRbr`$7J(=xHE8c^k#FHocevZ|%p=V!r)e3~h*O(>VIcpN)aC~!|2L4LNn^zJ zDn3xLz4{)^DC?gGB_Zr%`Gy_ z$4gZgY5L1y%PO|_CjA>EBYW4#*j^RWsYZ_NHF0E5kg3hMkt2J8Ozpv5PSOW)i_Boh zCCL}Dy*uzR+$UpuM>&gqT*meu!I$u5j_rMwY5&vdg#D$(_Aanp=TNLMV%b-tA6rsa zC&GHmm(^GMmex0p3hTY6ALGBKWLdW4*%!pO1MW|G#$Z|NqHO`Ac2{c@6w; zY5)bT9j(ngw^RZlcbX#C?&EVWpW+bxhE*wt@Wn~TOPNf?b6X8}kq6)V#SNFks^mL9 e-uFKNh=1bc{r~8{>E-?Z`6+Yj_NoH;na-B!h*=CbX9OvmI-r8~MG>`S#Y3n*AA$DTt>Fj!@$-sKn z?Ch?SIB5-45rWbO&(Ol3Qa~s`B~TGULZy^PrItczfl_Ftf+`jMfIxtXKydEd*~DEt zQE3E?#Edle%$&!}+_T@g-@SK<2$`~-B1#gG!KadP7B^EI%C62yM=aDRBRHO@LJk$l zCzlG7ydy=FfhYq}2BHi^8Hh4)C1ij(TOGz#T1IV@fhYq}2F_@jDAS=F*D?s@XX@1h!cc@u4CrvSCzU(PbSTHQ)ZqjjP6*$O zaD{^4?Hrd3+zDwdqc+Mwl!3_^P-QcVc2S;euDFTocaQwss6Ud+7jo?LWNswy`laUP zUm!NEZu*Rwh9L|wTO95^UKs5yWlLp0lPw-pce}HtxZ)Wp>jT3ays&Pc9_Ie|T5V-I4MRS8M7?`w)-UZEal8v(b(3N2SC$Ujc5I^|-K6suWWbb) zY%E$vt+a)9(|(5jDBVq;pik06bdJ7G-=-(%8Tv8(gr1`p=oR`M{efPmztcbHpO}F; zkf_HzB(WN6uojz;Mh7zJLO1rH7k$`=gRtS?HXMhCQJla@jNwk)jgR0S+>ZzFIXsNd zaI=IG)1O_!(Zruk?4us_?y4-&b3_#8($BSB7s`uCXIi!8dW^O(OHQ!uOZ6 z6lUw^&YQnral@*0*KOL;asFqm`6VZUD!`EjN&uiLK!Tt1ji_XFTXC@{v6pe?oVb8$ zjq)nysGN{NK%!j&6&gsA;H604IVFw#xK#==yCcH{ghGo61_^lr@w(1 zff83^G2?G3mZK5t(To&Y(2A|t&REQ#6Fc!v?83Vkj|XuWM{pEbyc>B8qkxeBji+%t z?!X7}A$*uIc`wc~Dj&pW@L7BT=W0;cUB#c4`0}MtSm`Z#`GR+heH4lX?cq(~3U3MP zZDEgTCUu1oOKWR;Q?ji^F-v>P#-`-PmbXSMQ}*>ukV=o3a@fCof|*e4lbA+gGMm)T z6=DKw)KWvWe1f{G$9-mi>)HPqn=D>pi8s_Nmzo)j=?zOPYhu|OMxb7~$`Zi{AhV%k z(>fu7fq)B7ho-I*mL4FaZR)VLO~R@U6tX8&-mpbj=SLhK0z!4V!?J3|A07|%b9V4o zP2duvrg|@Mb5O?XdA)Ko9b(9j&`Bnhd+1(zkRGN-=uvu%en?LTwDK&yNWY=i=uh-# z`V0M){sjXIun^b4!ZI%C4P4Ug*cK8?4|a#dau`DatsLVbK7mtADyMN*Kq()^$GCXU z;Vbwm9>F*9T|9yB;rn<7&*9g)Y%fVF7lZC42_?RR8RhMP?*E8Mt`&(Y=xd8XGmN1? zBxeTX5KqRJDISq}_>VWnS1KNnvA`o$>Q#=6gxr0*N!rui0-sHQjHq94#x^ep|HUZ&sDAL(@fGkHRz(045sF|n+~ zs$gQX4sATSX~$-4!4BMrn|XrM!_%9+=*IvKFcxi|;v9pEB7BrlVNAXc?`K4QfQjXP zC`^8w3FVahZGl06T z-mVOPUUk=@_Df@|Nb-4)oZj_7C2ix&{qxKb7g}P=X)xq|3D95!7!dV}L bv(&r!|MfouqV+#79HRAqbz}MeasB^qa=HcT literal 0 HcmV?d00001 diff --git a/NEON_manuscript/Manuscript/Appendix_S5_IxodesJAGS_Hyp3.txt b/NEON_manuscript/Manuscript/Appendix_S5_IxodesJAGS_Hyp3.txt index a55a3db1..fbbff53d 100644 --- a/NEON_manuscript/Manuscript/Appendix_S5_IxodesJAGS_Hyp3.txt +++ b/NEON_manuscript/Manuscript/Appendix_S5_IxodesJAGS_Hyp3.txt @@ -1,157 +1,157 @@ -JAGS model code generated by package mvgam - -GAM formula: -y ~ s(siteID, bs = re) + s(cum_gdd, k = 5) + s(cum_gdd, siteID, k = 5, m = 1, bs = fs) + s(season, m = 2, bs = cc, k = 12) + s(season, series, m = 1, k = 4, bs = fs) - -Trend model: -RW - -Required data: -integer n; number of timepoints per series -integer n_series; number of series -matrix y; time-ordered observations of dimension n x n_series (missing values allowed) -matrix ytimes; time-ordered n x n_series matrix (which row in X belongs to each [time, series] observation?) -matrix X; mgcv GAM design matrix of dimension (n x n_series) x basis dimension -matrix S2; mgcv smooth penalty matrix S2 -matrix S4; mgcv smooth penalty matrix S4 -vector zero; prior basis coefficient locations vector of length ncol(X) -vector p_coefs; vector (length = 1) of prior Gaussian means for parametric effects -vector p_taus; vector (length = 1) of prior Gaussian precisions for parametric effects - - -#### Begin model #### -model { - -## GAM linear predictor -eta <- X %*% b - -## mean expectations -for (i in 1:n) { -for (s in 1:n_series) { -mu[i, s] <- exp(eta[ytimes[i, s]] + trend[i, s]) -} -} - -## latent factors evolve as time series with penalised precisions; -## the penalty terms force any un-needed factors to evolve as flat lines -for (j in 1:n_lv) { -LV[1, j] ~ dnorm(0, penalty[j]) -} - -for (i in 2:n) { -for (j in 1:n_lv){ -LV[i, j] ~ dnorm(LV[i - 1, j], penalty[j]) -} -} - -## shrinkage penalties for each factor's precision act to squeeze -## the entire factor toward a flat white noise process if supported by -## the data. The prior for individual factor penalties allows each factor to possibly -## have a relatively large penalty, which shrinks the prior for that factor's variance -## substantially. Penalties increase exponentially with the number of factors following -## Welty, Leah J., et al. Bayesian distributed lag models: estimating effects of particulate -## matter air pollution on daily mortality Biometrics 65.1 (2009): 282-291. -pi ~ dunif(0, n_lv) -X2 ~ dnorm(0, 1)T(0, ) - -# eta1 controls the baseline penalty -eta1 ~ dunif(-1, 1) - -# eta2 controls how quickly the penalties exponentially increase -eta2 ~ dunif(-1, 1) - -for (t in 1:n_lv) { -X1[t] ~ dnorm(0, 1)T(0, ) -l.dist[t] <- max(t, pi[]) -l.weight[t] <- exp(eta2[] * l.dist[t]) -l.var[t] <- exp(eta1[] * l.dist[t] / 2) * 1 -theta.prime[t] <- l.weight[t] * X1[t] + (1 - l.weight[t]) * X2[] -penalty[t] <- max(0.0001, theta.prime[t] * l.var[t]) -} - -## latent factor loadings: standard normal with identifiability constraints -## upper triangle of loading matrix set to zero -for (j in 1:(n_lv - 1)) { -for (j2 in (j + 1):n_lv) { -lv_coefs[j, j2] <- 0 -} -} - -## positive constraints on loading diagonals -for (j in 1:n_lv) { -lv_coefs[j, j] ~ dnorm(0, 1)T(0, 1); -} - -## lower diagonal free -for (j in 2:n_lv) { -for (j2 in 1:(j - 1)) { -lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1); -} -} - -## other elements also free -for (j in (n_lv + 1):n_series) { -for (j2 in 1:n_lv) { -lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1); -} -} - -## trend evolution depends on latent factors -for (i in 1:n) { -for (s in 1:n_series) { -trend[i, s] <- inprod(lv_coefs[s,], LV[i,]) -} -} - -## likelihood functions -for (i in 1:n) { -for (s in 1:n_series) { -y[i, s] ~ dnegbin(rate[i, s], r[s]) -rate[i, s] <- ifelse((r[s] / (r[s] + mu[i, s])) < min_eps, min_eps, -(r[s] / (r[s] + mu[i, s]))) -} -} - -## complexity penalising prior for the overdispersion parameter; -## where the likelihood reduces to a 'base' model (Poisson) unless -## the data support overdispersion -for (s in 1:n_series) { -r[s] <- 1 / r_raw[s] -r_raw[s] ~ dnorm(0, 0.1)T(0, ) -} - -## posterior predictions -for (i in 1:n) { -for (s in 1:n_series) { -ypred[i, s] ~ dnegbin(rate[i, s], r[s]) -} -} - -## GAM-specific priors -## parametric effect priors (regularised for identifiability) -for (i in 1:1) { b[i] ~ dnorm(p_coefs[i], p_taus[i]) } -## prior (non-centred) for s(siteID)... -for (i in 2:4) { -b_raw[i] ~ dnorm(0, 1) -b[i] <- mu_raw1 + b_raw[i] * sigma_raw1 -} -sigma_raw1 ~ dexp(1) -mu_raw1 ~ dnorm(0, 1) -## prior for s(cum_gdd)... -K2 <- S2[1:4,1:4] * lambda[2] + S2[1:4,5:8] * lambda[3] -b[5:8] ~ dmnorm(zero[5:8],K2) -## prior for s(cum_gdd,siteID)... -for (i in c(9:12,14:17,19:22)) { b[i] ~ dnorm(0, lambda[4]) } -for (i in c(13,18,23)) { b[i] ~ dnorm(0, lambda[5]) } -## prior for s(season)... -K4 <- S4[1:10,1:10] * lambda[6] -b[24:33] ~ dmnorm(zero[24:33],K4) -## prior for s(season,series)... -for (i in c(34:36,38:40,42:44,46:48,50:52,54:56,58:60,62:64)) { b[i] ~ dnorm(0, lambda[7]) } -for (i in c(37,41,45,49,53,57,61,65)) { b[i] ~ dnorm(0, lambda[8]) } -## smoothing parameter priors... -for (i in 1:8) { -lambda[i] ~ dexp(0.05)T(0.0001, ) -rho[i] <- log(lambda[i]) -} -} +JAGS model code generated by package mvgam + +GAM formula: +y ~ s(siteID, bs = re) + s(cum_gdd, k = 5) + s(cum_gdd, siteID, k = 5, m = 1, bs = fs) + s(season, m = 2, bs = cc, k = 12) + s(season, series, m = 1, k = 4, bs = fs) + +Trend model: +RW + +Required data: +integer n; number of timepoints per series +integer n_series; number of series +matrix y; time-ordered observations of dimension n x n_series (missing values allowed) +matrix ytimes; time-ordered n x n_series matrix (which row in X belongs to each [time, series] observation?) +matrix X; mgcv GAM design matrix of dimension (n x n_series) x basis dimension +matrix S2; mgcv smooth penalty matrix S2 +matrix S4; mgcv smooth penalty matrix S4 +vector zero; prior basis coefficient locations vector of length ncol(X) +vector p_coefs; vector (length = 1) of prior Gaussian means for parametric effects +vector p_taus; vector (length = 1) of prior Gaussian precisions for parametric effects + + +#### Begin model #### +model { + +## GAM linear predictor +eta <- X %*% b + +## mean expectations +for (i in 1:n) { +for (s in 1:n_series) { +mu[i, s] <- exp(eta[ytimes[i, s]] + trend[i, s]) +} +} + +## latent factors evolve as time series with penalised precisions; +## the penalty terms force any un-needed factors to evolve as flat lines +for (j in 1:n_lv) { +LV[1, j] ~ dnorm(0, penalty[j]) +} + +for (i in 2:n) { +for (j in 1:n_lv){ +LV[i, j] ~ dnorm(LV[i - 1, j], penalty[j]) +} +} + +## shrinkage penalties for each factor's precision act to squeeze +## the entire factor toward a flat white noise process if supported by +## the data. The prior for individual factor penalties allows each factor to possibly +## have a relatively large penalty, which shrinks the prior for that factor's variance +## substantially. Penalties increase exponentially with the number of factors following +## Welty, Leah J., et al. Bayesian distributed lag models: estimating effects of particulate +## matter air pollution on daily mortality Biometrics 65.1 (2009): 282-291. +pi ~ dunif(0, n_lv) +X2 ~ dnorm(0, 1)T(0, ) + +# eta1 controls the baseline penalty +eta1 ~ dunif(-1, 1) + +# eta2 controls how quickly the penalties exponentially increase +eta2 ~ dunif(-1, 1) + +for (t in 1:n_lv) { +X1[t] ~ dnorm(0, 1)T(0, ) +l.dist[t] <- max(t, pi[]) +l.weight[t] <- exp(eta2[] * l.dist[t]) +l.var[t] <- exp(eta1[] * l.dist[t] / 2) * 1 +theta.prime[t] <- l.weight[t] * X1[t] + (1 - l.weight[t]) * X2[] +penalty[t] <- max(0.0001, theta.prime[t] * l.var[t]) +} + +## latent factor loadings: standard normal with identifiability constraints +## upper triangle of loading matrix set to zero +for (j in 1:(n_lv - 1)) { +for (j2 in (j + 1):n_lv) { +lv_coefs[j, j2] <- 0 +} +} + +## positive constraints on loading diagonals +for (j in 1:n_lv) { +lv_coefs[j, j] ~ dnorm(0, 1)T(0, 1); +} + +## lower diagonal free +for (j in 2:n_lv) { +for (j2 in 1:(j - 1)) { +lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1); +} +} + +## other elements also free +for (j in (n_lv + 1):n_series) { +for (j2 in 1:n_lv) { +lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1); +} +} + +## trend evolution depends on latent factors +for (i in 1:n) { +for (s in 1:n_series) { +trend[i, s] <- inprod(lv_coefs[s,], LV[i,]) +} +} + +## likelihood functions +for (i in 1:n) { +for (s in 1:n_series) { +y[i, s] ~ dnegbin(rate[i, s], r[s]) +rate[i, s] <- ifelse((r[s] / (r[s] + mu[i, s])) < min_eps, min_eps, +(r[s] / (r[s] + mu[i, s]))) +} +} + +## complexity penalising prior for the overdispersion parameter; +## where the likelihood reduces to a 'base' model (Poisson) unless +## the data support overdispersion +for (s in 1:n_series) { +r[s] <- 1 / r_raw[s] +r_raw[s] ~ dnorm(0, 0.1)T(0, ) +} + +## posterior predictions +for (i in 1:n) { +for (s in 1:n_series) { +ypred[i, s] ~ dnegbin(rate[i, s], r[s]) +} +} + +## GAM-specific priors +## parametric effect priors (regularised for identifiability) +for (i in 1:1) { b[i] ~ dnorm(p_coefs[i], p_taus[i]) } +## prior (non-centred) for s(siteID)... +for (i in 2:4) { +b_raw[i] ~ dnorm(0, 1) +b[i] <- mu_raw1 + b_raw[i] * sigma_raw1 +} +sigma_raw1 ~ dexp(1) +mu_raw1 ~ dnorm(0, 1) +## prior for s(cum_gdd)... +K2 <- S2[1:4,1:4] * lambda[2] + S2[1:4,5:8] * lambda[3] +b[5:8] ~ dmnorm(zero[5:8],K2) +## prior for s(cum_gdd,siteID)... +for (i in c(9:12,14:17,19:22)) { b[i] ~ dnorm(0, lambda[4]) } +for (i in c(13,18,23)) { b[i] ~ dnorm(0, lambda[5]) } +## prior for s(season)... +K4 <- S4[1:10,1:10] * lambda[6] +b[24:33] ~ dmnorm(zero[24:33],K4) +## prior for s(season,series)... +for (i in c(34:36,38:40,42:44,46:48,50:52,54:56,58:60,62:64)) { b[i] ~ dnorm(0, lambda[7]) } +for (i in c(37,41,45,49,53,57,61,65)) { b[i] ~ dnorm(0, lambda[8]) } +## smoothing parameter priors... +for (i in 1:8) { +lambda[i] ~ dexp(0.05)T(0.0001, ) +rho[i] <- log(lambda[i]) +} +} diff --git a/NEON_manuscript/next_todo.R b/NEON_manuscript/next_todo.R index 1bdc01fc..f9e89f1d 100644 --- a/NEON_manuscript/next_todo.R +++ b/NEON_manuscript/next_todo.R @@ -6,14 +6,263 @@ dat$true_corrs mod1 <- mvgam(formula = y ~ s(season, bs = 'cc') + s(series, bs = 're'), data_train = dat$data_train, - trend_model = 'AR3', + trend_model = 'AR1', family = 'poisson', use_lv = TRUE, n_lv = 2, use_stan = TRUE, run_model = T, burnin = 10) +mod1$model_file summary(mod1) +plot_mvgam_factors(mod1) +plot(mod1, type = 'residuals') +fake <- dat$data_test +fake$y <- NULL +plot(mod1, 'trend', data_test = fake) + + + +plot_mvgam_series = function(data_train, data_test, series, n_bins, + log_scale = FALSE){ + + if(series == 'all'){ + n_plots <- length(levels(data_train$series)) + pages <- 1 + .pardefault <- par(no.readonly=T) + par(.pardefault) + + if (n_plots > 4) pages <- 2 + if (n_plots > 8) pages <- 3 + if (n_plots > 12) pages <- 4 + if (pages != 0) { + ppp <- n_plots %/% pages + + if (n_plots %% pages != 0) { + ppp <- ppp + 1 + while (ppp * (pages - 1) >= n_plots) pages <- pages - 1 + } + + # Configure layout matrix + c <- r <- trunc(sqrt(ppp)) + if (c<1) r <- c <- 1 + if (c*r < ppp) c <- c + 1 + if (c*r < ppp) r <- r + 1 + oldpar<-par(mfrow=c(r,c)) + + } else { ppp <- 1; oldpar <- par()} + + all_ys <- lapply(seq_len(n_plots), function(x){ + if(log_scale){ + log( data.frame(y = data_train$y, + series = data_train$series, + time = data_train$time) %>% + dplyr::filter(series == levels(data_train$series)[x]) %>% + dplyr::pull(y) + 1) + } else { + data.frame(y = data_train$y, + series = data_train$series, + time = data_train$time) %>% + dplyr::filter(series == levels(data_train$series)[x]) %>% + dplyr::pull(y) + } + }) + + for(i in 1:n_plots){ + s_name <- levels(data_train$series)[i] + + truth <- data.frame(y = data_train$y, + series = data_train$series, + time = data_train$time) %>% + dplyr::filter(series == s_name) %>% + dplyr::select(time, y) %>% + dplyr::distinct() %>% + dplyr::arrange(time) %>% + dplyr::pull(y) + + if(log_scale){ + truth <- log(truth + 1) + ylab <- 'log(Y + 1)' + } else { + ylab <- 'Y' + } + + plot(1, type = "n", bty = 'L', + xlab = 'Time', + ylab = ylab, + ylim = range(unlist(all_ys)), + xlim = c(0, length(c(truth)))) + title(s_name, line = 0) + + if(n_plots > 1){ + for(x in 1:n_plots){ + lines(all_ys[[x]], lwd = 1.85, col = 'grey85') + } + } + + lines(x = 1:length(truth), y = truth, lwd = 3, col = "white") + lines(x = 1:length(truth), y = truth, lwd = 2.5, col = "#8F2727") + box(bty = 'L', lwd = 2) + + } + + } else { + s_name <- levels(data_train$series)[series] + truth <- data.frame(y = data_train$y, + series = data_train$series, + time = data_train$time) %>% + dplyr::filter(series == s_name) %>% + dplyr::select(time, y) %>% + dplyr::distinct() %>% + dplyr::arrange(time) %>% + dplyr::pull(y) + + layout(matrix(1:4, nrow = 2, byrow = TRUE)) + if(!missing(data_test)){ + test <- data.frame(y = data_test$y, + series = data_test$series, + time = data_test$time) %>% + dplyr::filter(series == s_name) %>% + dplyr::select(time, y) %>% + dplyr::distinct() %>% + dplyr::arrange(time) %>% + dplyr::pull(y) + + plot(1, type = "n", bty = 'L', + xlab = 'Time', + ylab = 'Y', + ylim = range(c(truth, test)), + xlim = c(0, length(c(truth, test)))) + title('Time series', line = 0) + + lines(x = 1:length(truth), y = truth, lwd = 2, col = "#8F2727") + lines(x = (length(truth)+1):length(c(truth, test)), y = test, lwd = 2, col = "black") + abline(v = length(truth)+1, col = '#FFFFFF60', lwd = 2.85) + abline(v = length(truth)+1, col = 'black', lwd = 2.5, lty = 'dashed') + box(bty = 'L', lwd = 2) + + if(missing(n_bins)){ + n_bins <- max(c(length(hist(c(truth, test), plot = F)$breaks), + 20)) + } + + hist(c(truth, test), border = "#8F2727", + lwd = 2, + freq = FALSE, + breaks = n_bins, + col = "#C79999", + ylab = 'Density', + xlab = 'Count', main = '') + title('Histogram', line = 0) + + acf(c(truth, test), + na.action = na.pass, bty = 'L', + lwd = 2.5, ci.col = 'black', col = "#8F2727", + main = '', ylab = 'Autocorrelation') + acf1 <- acf(c(truth, test), plot = F, + na.action = na.pass) + clim <- qnorm((1 + .95)/2)/sqrt(acf1$n.used) + abline(h = clim, col = '#FFFFFF', lwd = 2.85) + abline(h = clim, col = 'black', lwd = 2.5, lty = 'dashed') + abline(h = -clim, col = '#FFFFFF', lwd = 2.85) + abline(h = -clim, col = 'black', lwd = 2.5, lty = 'dashed') + box(bty = 'L', lwd = 2) + title('ACF', line = 0) + + ecdf_plotdat = function(vals, x){ + func <- ecdf(vals) + func(x) + } + + plot_x <- seq(min(c(truth, test), na.rm = T), + max(c(truth, test), na.rm = T)) + plot(1, type = "n", bty = 'L', + xlab = 'Count', + ylab = 'Empirical CDF', + xlim = c(min(plot_x), max(plot_x)), + ylim = c(0, 1)) + title('CDF', line = 0) + lines(x = plot_x, + y = ecdf_plotdat(c(truth, test), + plot_x), + col = "#8F2727", + lwd = 2.5) + box(bty = 'L', lwd = 2) + + } else { + plot(1, type = "n", bty = 'L', + xlab = 'Time', + ylab = 'Observations', + ylim = range(c(truth)), + xlim = c(0, length(c(truth)))) + title('Time series', line = 0) + + lines(x = 1:length(truth), y = truth, lwd = 2, col = "#8F2727") + box(bty = 'L', lwd = 2) + + if(missing(n_bins)){ + n_bins <- max(c(length(hist(c(truth), plot = F)$breaks), + 20)) + } + + hist(c(truth), border = "#8F2727", + lwd = 2, + freq = FALSE, + breaks = n_bins, + col = "#C79999", + ylab = 'Density', + xlab = 'Count', main = '') + title('Histogram', line = 0) + + + acf(c(truth), + na.action = na.pass, bty = 'L', + lwd = 2.5, ci.col = 'black', col = "#8F2727", + main = '', ylab = 'Autocorrelation') + acf1 <- acf(c(truth), plot = F, + na.action = na.pass) + clim <- qnorm((1 + .95)/2)/sqrt(acf1$n.used) + abline(h = clim, col = '#FFFFFF', lwd = 2.85) + abline(h = clim, col = 'black', lwd = 2.5, lty = 'dashed') + abline(h = -clim, col = '#FFFFFF', lwd = 2.85) + abline(h = -clim, col = 'black', lwd = 2.5, lty = 'dashed') + box(bty = 'L', lwd = 2) + title('ACF', line = 0) + + + ecdf_plotdat = function(vals, x){ + func <- ecdf(vals) + func(x) + } + + plot_x <- seq(min(truth, na.rm = T), + max(truth, na.rm = T)) + plot(1, type = "n", bty = 'L', + xlab = 'Count', + ylab = 'Empirical CDF', + xlim = c(min(plot_x), max(plot_x)), + ylim = c(0, 1)) + title('CDF', line = 0) + lines(x = plot_x, + y = ecdf_plotdat(truth, + plot_x), + col = "#8F2727", + lwd = 2.5) + box(bty = 'L', lwd = 2) + } + + layout(1) + } + +} + +dat <- sim_mvgam(T = 100, n_series = 4, n_lv = 2, + mu_obs = c(4, 6, 10, 14), trend_rel = 0.3, + seasonality = 'hierarchical') +plot_mvgam_series(data_train = dat$data_train, + n_bins = 20, + series = 'all', + log_scale = TRUE) # Good for testing model files without compiling stanc(model_code = mod1$model_file)$model_name diff --git a/R/add_stan_data.R b/R/add_stan_data.R index 1159ee57..90725caf 100644 --- a/R/add_stan_data.R +++ b/R/add_stan_data.R @@ -299,7 +299,7 @@ add_stan_data = function(jags_file, stan_file, use_lv = FALSE, sub("\\)", "", sub(".*\\:", "", stan_file[b_raw_indices[i]-1]))))), - ') {\nb[i] <- mu_raw[', i, '] + b_raw[i] * sigma_raw[',i, + ') {\nb[i] = mu_raw[', i, '] + b_raw[i] * sigma_raw[',i, '];\n}') min_beta[i] <- as.numeric(sub("for \\(i in ", "", sub("\\:.*", "", diff --git a/R/add_trend_lines.R b/R/add_trend_lines.R index 685d2c21..fe430e93 100644 --- a/R/add_trend_lines.R +++ b/R/add_trend_lines.R @@ -280,7 +280,7 @@ add_trend_lines = function(model_file, stan = FALSE, paste0('row_vector[num_basis] b_raw;\n\n// latent factor AR1 terms\nvector[n_lv] ar1;') model_file[grep('// dynamic factor estimates', model_file) + 6] <- - paste0('LV[2:n, j] ~ normal(ar1[j] * LV[1:(n - 1), j], sigma[s]);') + paste0('LV[2:n, j] ~ normal(ar1[j] * LV[1:(n - 1), j], 1);') model_file[grep('model \\{', model_file) + 2] <- paste0('\n// priors for AR parameters\nar1 ~ normal(0, 0.5);\n') diff --git a/R/plot.mvgam.R b/R/plot.mvgam.R index 08ef3bb2..fa85cf39 100644 --- a/R/plot.mvgam.R +++ b/R/plot.mvgam.R @@ -44,11 +44,11 @@ plot.mvgam = function(object, type = 'smooths', # Other errors and warnings will propagate from individual functions below if(type == 're'){ - plot_mvgam_randomeffects(object) + plot_mvgam_randomeffects(object, ...) } if(type == 'pterms'){ - plot_mvgam_pterms(object) + plot_mvgam_pterms(object, ...) } if(type == 'residuals'){ diff --git a/R/plot_mvgam_fc.R b/R/plot_mvgam_fc.R index 007f914c..3bfb9338 100644 --- a/R/plot_mvgam_fc.R +++ b/R/plot_mvgam_fc.R @@ -18,8 +18,10 @@ #'\code{realisations = TRUE}. Ignored otherwise #'@param hide_xlabels \code{logical}. If \code{TRUE}, no xlabels are printed to allow the user to add custom labels using #'\code{axis} from base \code{R} -#'@param ylab Optional \code{character} string specifying the y-axis label +#'@param xlab label for x axis. +#'@param ylab label for y axis. #'@param ylim Optional \code{vector} of y-axis limits (min, max) +#'@param ... further \code{\link[graphics]{par}} graphical parameters. #'@param return_forecasts \code{logical}. If \code{TRUE}, the function will plot the forecast #'as well as returning the forecast object (as a \code{matrix} of dimension \code{n_samples} x \code{horizon}) #'@details Posterior predictions are drawn from the fitted \code{mvgam} and used to calculate posterior @@ -30,8 +32,8 @@ #'@export plot_mvgam_fc = function(object, series = 1, data_test, realisations = FALSE, n_realisations = 15, - hide_xlabels = FALSE, ylab, ylim, - return_forecasts = FALSE){ + hide_xlabels = FALSE, xlab, ylab, ylim, + return_forecasts = FALSE, ...){ # Check arguments if(class(object) != 'mvgam'){ @@ -168,6 +170,10 @@ plot_mvgam_fc = function(object, series = 1, data_test, ylab <- paste0('Predicitons for ', levels(data_train$series)[series]) } + if(missing(xlab)){ + xlab <- 'Time' + } + pred_vals <- seq(1:length(preds_last)) if(hide_xlabels){ plot(1, type = "n", bty = 'L', @@ -175,13 +181,13 @@ plot_mvgam_fc = function(object, series = 1, data_test, xaxt = 'n', ylab = ylab, xlim = c(0, length(preds_last)), - ylim = ylim) + ylim = ylim, ...) } else { plot(1, type = "n", bty = 'L', - xlab = 'Time', + xlab = xlab, ylab = ylab, xlim = c(0, length(preds_last)), - ylim = ylim) + ylim = ylim, ...) } if(realisations){ @@ -223,20 +229,29 @@ plot_mvgam_fc = function(object, series = 1, data_test, time = data_test$time) } + # Show historical distribution in grey + last_train <- (NROW(data_train) / NCOL(object$ytimes)) + polygon(c(pred_vals[1:(NROW(data_train) / NCOL(object$ytimes))], + rev(pred_vals[1:(NROW(data_train) / NCOL(object$ytimes))])), + c(cred[1,1:(NROW(data_train) / NCOL(object$ytimes))], + rev(cred[9,1:(NROW(data_train) / NCOL(object$ytimes))])), + col = 'grey70', border = NA) + # Plot training and testing points points(dplyr::bind_rows(data_train, data_test) %>% dplyr::filter(series == s_name) %>% dplyr::select(time, y) %>% dplyr::distinct() %>% dplyr::arrange(time) %>% - dplyr::pull(y), pch = 16, col = "white", cex = 0.75) + dplyr::pull(y), pch = 16, col = "white", cex = 0.8) points(dplyr::bind_rows(data_train, data_test) %>% dplyr::filter(series == s_name) %>% dplyr::select(time, y) %>% dplyr::distinct() %>% dplyr::arrange(time) %>% dplyr::pull(y), pch = 16, col = "black", cex = 0.65) - abline(v = NROW(data_train) / NCOL(object$ytimes), lty = 'dashed') + abline(v = last_train, col = '#FFFFFF60', lwd = 2.85) + abline(v = last_train, col = 'black', lwd = 2.5, lty = 'dashed') # Calculate out of sample DRPS and print the score drps_score <- function(truth, fc, interval_width = 0.9){ @@ -306,13 +321,13 @@ plot_mvgam_fc = function(object, series = 1, data_test, dplyr::select(time, y) %>% dplyr::distinct() %>% dplyr::arrange(time) %>% - dplyr::pull(y),pch = 16, col = "white", cex = 0.75) + dplyr::pull(y),pch = 16, col = "white", cex = 0.8) points(data_train %>% dplyr::filter(series == s_name) %>% dplyr::select(time, y) %>% dplyr::distinct() %>% dplyr::arrange(time) %>% - dplyr::pull(y),pch = 16, col = "black", cex = 0.65 ) + dplyr::pull(y),pch = 16, col = "black", cex = 0.65) } if(return_forecasts){ diff --git a/R/plot_mvgam_resids.R b/R/plot_mvgam_resids.R index e7791fc7..3a0379f1 100644 --- a/R/plot_mvgam_resids.R +++ b/R/plot_mvgam_resids.R @@ -320,13 +320,14 @@ x <- sapply(1:length(idx), # Plot plot(median_preds[1:length(object$resids[[series]])], object$resids[[series]], - main = 'Resids vs Fitted Values', + bty = 'L', xlab = 'Fitted values', - ylab = 'Residuals', + ylab = 'DS residuals', pch = 16, col = 'white', cex = 1, ylim = range(resid_probs, na.rm = T)) +title('Resids vs Fitted Values', line = 0) rect(xleft = x[seq(1, N*2, by = 2)], xright = x[seq(2, N*2, by = 2)], @@ -358,8 +359,8 @@ for (k in 1:N) { y = c(resid_probs[k,5], resid_probs[k,5]), col = c_dark, lwd = 2) } -abline(h = 0, lty = 'dashed', col = 'white', lwd = 2.85) -abline(h = 0, lty = 'dashed', col = 'black', lwd = 2.5) +abline(h = 0, col = '#FFFFFF60', lwd = 2.85) +abline(h = 0, col = 'black', lwd = 2.5, lty = 'dashed') # Q-Q plot coords <- qqnorm(series_residuals[1,], plot.it = F) @@ -379,13 +380,15 @@ pred_vals <- coords$x[order(coords$x)] pred_vals <- pred_vals[complete.cases(cred[1,])] plot(x = pred_vals, y = cred[5,][complete.cases(cred[1,])], - main = 'Normal Q-Q Plot', + bty = 'L', xlab = 'Theoretical Quantiles', ylab = 'Sample Quantiles', pch = 16, col = 'white', cex = 1, - ylim = range(cred, na.rm = T)) + ylim = range(cred, na.rm = T), + tck = -0.04) +title('Normal Q-Q Plot', line = 0) polygon(c(pred_vals, rev(pred_vals)), c(cred[1,][complete.cases(cred[1,])], rev(cred[9,][complete.cases(cred[1,])])), col = c_light, border = NA) @@ -399,7 +402,7 @@ polygon(c(pred_vals, rev(pred_vals)), c(cred[4,][complete.cases(cred[1,])], rev(cred[6,][complete.cases(cred[1,])])), col = c_mid_highlight, border = NA) lines(pred_vals, cred[5,][complete.cases(cred[1,])], col = c_dark, lwd = 2.5) -qqline(cred[5,][complete.cases(cred[1,])], col = 'white', lwd = 3) +qqline(cred[5,][complete.cases(cred[1,])], col = '#FFFFFF60', lwd = 3) qqline(cred[5,][complete.cases(cred[1,])], col = 'black', lwd = 2.5) # ACF plot @@ -426,16 +429,16 @@ cred <- sapply(1:NCOL(resid_acf), probs = probs, na.rm = T)) cred <- cred[, -1] clim <- qnorm((1 + .95)/2)/sqrt(acf1$n.used) -plot(1, type = "n", +plot(1, type = "n", bty = 'L', xlab = 'Lag', - ylab = 'ACF', - main = 'ACF', + ylab = 'Autocorrelation', xlim = c(1, N-1), xaxt = 'n', ylim = range(c(cred, -clim - 0.05, clim + 0.05))) axis(1, at = seq(1, NCOL(cred), by = 2)) +title('ACF', line = 0) N <- N - 1 rect(xleft = x[seq(1, N*2, by = 2)], @@ -468,10 +471,10 @@ for (k in 1:N) { y = c(cred[5,k], cred[5,k]), col = c_dark, lwd = 2) } -abline(h = clim, lty = 'dashed', col = 'white', lwd = 2.85) -abline(h = clim, lty = 'dashed', col = 'black', lwd = 2.5) -abline(h = -clim, lty = 'dashed', col = 'white', lwd = 2.85) -abline(h = -clim, lty = 'dashed', col = 'black', lwd = 2.5) +abline(h = clim, col = '#FFFFFF60', lwd = 2.85) +abline(h = clim, col = 'black', lwd = 2.5, lty = 'dashed') +abline(h = -clim, col = '#FFFFFF60', lwd = 2.85) +abline(h = -clim, col = 'black', lwd = 2.5, lty = 'dashed') # PACF plot pacf1 <- pacf(series_residuals[1,], plot = F, @@ -497,16 +500,16 @@ cred <- sapply(1:NCOL(resid_pacf), probs = probs, na.rm = T)) clim <- qnorm((1 + .95)/2)/sqrt(pacf1$n.used) -plot(1, type = "n", +plot(1, type = "n", bty = 'L', xlab = 'Lag', - ylab = 'Partial ACF', - main = 'pACF', + ylab = 'Autocorrelation', xlim = c(1, length(sorted_x)), xaxt = 'n', ylim = range(c(cred, -clim - 0.05, clim + 0.05))) axis(1, at = seq(1, NCOL(cred), by = 2)) +title('pACF', line = 0) rect(xleft = x[seq(1, N*2, by = 2)], xright = x[seq(2, N*2, by = 2)], @@ -538,11 +541,10 @@ for (k in 1:N) { y = c(cred[5,k], cred[5,k]), col = c_dark, lwd = 2) } -abline(h = clim, lty = 'dashed', col = 'white', lwd = 2.85) -abline(h = clim, lty = 'dashed', col = 'black', lwd = 2.5) -abline(h = -clim, lty = 'dashed', col = 'white', lwd = 2.85) -abline(h = -clim, lty = 'dashed', col = 'black', lwd = 2.5) - +abline(h = clim, col = '#FFFFFF60', lwd = 2.85) +abline(h = clim, col = 'black', lwd = 2.5, lty = 'dashed') +abline(h = -clim, col = '#FFFFFF60', lwd = 2.85) +abline(h = -clim, col = 'black', lwd = 2.5, lty = 'dashed') layout(1) } diff --git a/R/plot_mvgam_trend.R b/R/plot_mvgam_trend.R index 8feb2f1f..894a3be9 100644 --- a/R/plot_mvgam_trend.R +++ b/R/plot_mvgam_trend.R @@ -12,10 +12,14 @@ #'\code{realisations = TRUE}. Ignored otherwise #'@param hide_xlabels \code{logical}. If \code{TRUE}, no xlabels are printed to allow the user to add custom labels using #'\code{axis} from base \code{R}. Ignored if \code{derivatives = TRUE} +#'@param xlab label for x axis. +#'@param ... further \code{\link[graphics]{par}} graphical parameters. #'@export plot_mvgam_trend = function(object, series = 1, data_test, realisations = FALSE, n_realisations = 15, - derivatives = FALSE, hide_xlabels = FALSE){ + derivatives = FALSE, hide_xlabels = FALSE, + xlab, + ...){ # Check arguments if(class(object) != 'mvgam'){ @@ -114,6 +118,10 @@ plot_mvgam_trend = function(object, series = 1, data_test, c_dark <- c("#8F2727") c_dark_highlight <- c("#7C0000") + if(missing(xlab)){ + xlab <- 'Time' + } + if(derivatives){ .pardefault <- par(no.readonly=T) par(.pardefault) @@ -122,10 +130,10 @@ plot_mvgam_trend = function(object, series = 1, data_test, mai = c(0.8, 0.8, 0.4, 0.4)) plot(1, type = "n", bty = 'L', - xlab = 'Time', + xlab = xlab, ylab = paste0('Estimated trend for ', levels(data_train$series)[series]), xlim = c(0, length(preds_last)), - ylim = range(cred)) + ylim = range(cred), ...) if(realisations){ for(i in 1:n_realisations){ @@ -157,9 +165,11 @@ plot_mvgam_trend = function(object, series = 1, data_test, if(!missing(data_test)){ if(class(data_train)[1] == 'list'){ - abline(v = length(data_train$y) / NCOL(object$ytimes), lty = 'dashed') + abline(v = length(data_train$y) / NCOL(object$ytimes), col = '#FFFFFF60', lwd = 2.85) + abline(v = length(data_train$y) / NCOL(object$ytimes), col = 'black', lwd = 2.5, lty = 'dashed') } else { - abline(v = NROW(data_train) / NCOL(object$ytimes), lty = 'dashed') + abline(v = NROW(data_train) / NCOL(object$ytimes), col = '#FFFFFF60', lwd = 2.85) + abline(v = NROW(data_train) / NCOL(object$ytimes), col = 'black', lwd = 2.5, lty = 'dashed') } } @@ -171,11 +181,11 @@ plot_mvgam_trend = function(object, series = 1, data_test, probs = probs, na.rm = T)) plot(1, type = "n", bty = "L", - xlab = 'Time', + xlab = xlab, ylab = '1st derivative', xlim = c(min(pred_vals), max(pred_vals)), ylim = c(min(cred, na.rm = T) - sd(first_derivs, na.rm = T), - max(cred, na.rm = T) + sd(first_derivs, na.rm = T))) + max(cred, na.rm = T) + sd(first_derivs, na.rm = T)), ...) if(realisations){ @@ -222,10 +232,10 @@ plot_mvgam_trend = function(object, series = 1, data_test, } else { plot(1, type = "n", bty = 'L', - xlab = 'Time', + xlab = xlab, ylab = paste0('Estimated trend for ', levels(data_train$series)[series]), xlim = c(0, length(preds_last)), - ylim = range(cred)) + ylim = range(cred), ...) } if(realisations){ @@ -258,9 +268,11 @@ plot_mvgam_trend = function(object, series = 1, data_test, if(!missing(data_test)){ if(class(data_train)[1] == 'list'){ - abline(v = length(data_train$y) / NCOL(object$ytimes), lty = 'dashed') + abline(v = length(data_train$y) / NCOL(object$ytimes), col = '#FFFFFF60', lwd = 2.85) + abline(v = length(data_train$y) / NCOL(object$ytimes), col = 'black', lwd = 2.5, lty = 'dashed') } else { - abline(v = NROW(data_train) / NCOL(object$ytimes), lty = 'dashed') + abline(v = NROW(data_train) / NCOL(object$ytimes), col = '#FFFFFF60', lwd = 2.85) + abline(v = NROW(data_train) / NCOL(object$ytimes), col = 'black', lwd = 2.5, lty = 'dashed') } } diff --git a/R/ppc.mvgam.R b/R/ppc.mvgam.R index 8fa1dadb..62448306 100644 --- a/R/ppc.mvgam.R +++ b/R/ppc.mvgam.R @@ -16,7 +16,11 @@ #'number of bins returned by a call to `hist` in base `R` #'@param legend_position The location may also be specified by setting x to a single keyword from the #'list "bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" and "center". -#'This places the legend on the inside of the plot frame at the given location. +#'This places the legend on the inside of the plot frame at the given location. Or alternatively, +#'use "none" to hide the legend. +#'@param xlab label for x axis. +#'@param ylab label for y axis. +#'@param ... further \code{\link[graphics]{par}} graphical parameters. #'@details Posterior predictions are drawn from the fitted \code{mvgam} and compared against #'the empirical distribution of the observed data for a specified series to help evaluate the model's #'ability to generate unbiased predictions. For all plots apart from the 'rootogram', posterior predictions @@ -41,13 +45,15 @@ ppc <- function(x, what, ...){ #'@method ppc mvgam #'@export ppc.mvgam = function(object, data_test, series = 1, type = 'density', - n_bins, legend_position){ + n_bins, legend_position, xlab, ylab, ...){ # Check arguments type <- match.arg(arg = type, choices = c("rootogram", "mean", "hist", "density", "pit", "cdf", "prop_zero")) + optional_args <- list(...) + if(class(object) != 'mvgam'){ stop('argument "object" must be of class "mvgam"') } @@ -194,6 +200,15 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', pred_props <- pred_props[-which(pred_props < lower)] } obs_prop <- length(which(truths == 0)) / length(truths) + + if(missing(ylab)){ + ylab <- 'Density' + } + + if(missing(xlab)){ + xlab <- paste0('Predicted proportion of zeroes for ', levels(data_train$series)[series]) + } + hist(pred_props, lwd = 2, xlim = c(min(min(pred_props), min(obs_prop)), max(max(pred_props), max(obs_prop))), @@ -202,8 +217,9 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', length.out = 15), border = "#B97C7C", col = "#C79999", - ylab = 'Density', - xlab = paste0('Predicted proportion of zeroes for ', levels(data_train$series)[series])) + ylab = ylab, + xlab = xlab, + ...) abline(v = obs_prop, lwd = 3, col = 'white') abline(v = obs_prop, lwd = 2.5, col = 'black') box(bty = 'L', lwd = 2) @@ -212,14 +228,17 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', legend_position = 'topright' } - legend(legend_position, - legend = c(expression(hat(y)[propzero]), - expression(y[propzero])), - bg = 'white', - col = c(c_mid, - 'black'), - lty = 1, lwd = 2, - bty = 'n') + if(legend_position != 'none'){ + legend(legend_position, + legend = c(expression(hat(y)[propzero]), + expression(y[propzero])), + bg = 'white', + col = c(c_mid, + 'black'), + lty = 1, lwd = 2, + bty = 'n') + } + } if(type == 'rootogram'){ @@ -276,14 +295,23 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', repped_x[k] + min(diff(xpos))/2 else repped_x[k] - min(diff(xpos))/2) + if(missing(xlab)){ + xlab <- expression(y) + } + + if(missing(ylab)){ + ylab <- expression(sqrt(frequency)) + } + # Plot the rootogram plot(1, type = "n", bty = 'L', - xlab = expression(y), - ylab = expression(sqrt(frequency)), + xlab = xlab, + ylab = ylab, xlim = range(xpos), ylim = range(c(data$tyexp, data[,13], data[,5], - data$tyexp - data$ty))) + data$tyexp - data$ty)), + ...) rect(xleft = x[seq(1, N*2, by = 2)], xright = x[seq(2, N*2, by = 2)], ytop = data$tyexp, @@ -335,6 +363,15 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', pred_means <- pred_means[-which(pred_means < lower)] } obs_mean <- mean(truths) + + if(missing(ylab)){ + ylab <- 'Density' + } + + if(missing(xlab)){ + xlab <- paste0('Predicted mean for ', levels(data_train$series)[series]) + } + hist(pred_means, xlim = c(min(min(pred_means), min(obs_mean)), max(max(pred_means), max(obs_mean))), @@ -343,8 +380,9 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', breaks = seq(min(pred_means), max(pred_means), length.out = 20), border = "#B97C7C", col = "#C79999", - ylab = 'Density', - xlab = paste0('Predicted mean for ', levels(data_train$series)[series])) + ylab = ylab, + xlab = xlab, + ...) abline(v = obs_mean, lwd = 3, col = 'white') abline(v = obs_mean, lwd = 2.5, col = 'black') box(bty = 'L', lwd = 2) @@ -353,6 +391,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', legend_position = 'topright' } + if(legend_position != 'none'){ legend(legend_position, legend = c(expression(hat(mu)), expression(mu)), @@ -361,6 +400,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', 'black'), lty = 1, lwd = 2, bty = 'n') + } } @@ -384,11 +424,20 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', ymax <- max(c(max(cred), max(true_dens$y))) + if(missing(ylab)){ + ylab <- paste0('Predictive density for ', levels(data_train$series)[series]) + } + + if(missing(xlab)){ + xlab <- '' + } + plot(1, type = "n", bty = 'L', - xlab = '', - ylab = paste0('Predictive density for ', levels(data_train$series)[series]), + xlab = xlab, + ylab = ylab, xlim = c(min_x, max_x), - ylim = c(0, ymax)) + ylim = c(0, ymax), + ...) polygon(c(true_dens$x, rev(true_dens$x)), c(cred[1,], rev(cred[9,])), col = c_light, border = NA) @@ -407,6 +456,8 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', if(missing(legend_position)){ legend_position = 'topright' } + + if(legend_position != 'none'){ legend(legend_position, legend = c(expression(hat(y)), 'y'), @@ -416,6 +467,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', lty = 1, lwd = 2, bty = 'n') + } box(bty = 'L', lwd = 2) } @@ -436,21 +488,32 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', } bin_lims <- range(c(truths, as.vector(preds))) - delta <- diff(range(preds)) / n_bins - breaks <- seq(bin_lims[1], bin_lims[2] + delta, delta) + #delta <- diff(range(preds)) / n_bins + breaks <- seq(bin_lims[1], bin_lims[2], length.out = n_bins) xlim <- c(0, max(max(density(preds[1,])$x), max(density(truths)$x))) ylim <- c(0, max(c(max(hist(truths, breaks = breaks, plot = F)$density), max(hist(preds, breaks = breaks, plot = F)$density)))) + + if(missing(xlab)){ + xlab <- paste0('Count') + } + + if(missing(ylab)){ + ylab = '' + } + hist(preds, breaks=breaks, lwd = 2, main='', - xlab = paste0('Predictive histogram for ', levels(data_train$series)[series]), + xlab = xlab, + ylab = ylab, ylim=ylim, xlim=xlim, border = "#B97C7C", col = "#C79999", - freq = F) + freq = F, + ...) par(lwd=2) hist(truths, breaks=breaks, @@ -467,6 +530,8 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', if(missing(legend_position)){ legend_position = 'topright' } + + if(legend_position != 'none'){ legend(legend_position, legend = c(expression(hat(y)), 'y'), @@ -476,6 +541,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', lty = 1, lwd = 2, bty = 'n') + } } @@ -497,11 +563,20 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', function(n) quantile(pred_cdfs[,n], probs = probs)) + if(missing(ylab)){ + ylab = paste0('Predictive CDF for ', levels(data_train$series)[series]) + } + + if(missing(xlab)){ + xlab = '' + } + plot(1, type = "n", bty = 'L', - xlab = '', - ylab = paste0('Predictive CDF for ', levels(data_train$series)[series]), + xlab = xlab, + ylab = ylab, xlim = c(min(plot_x), max(plot_x)), - ylim = c(0, 1)) + ylim = c(0, 1), + ...) polygon(c(plot_x, rev(plot_x)), c(cred[1,], rev(cred[9,])), col = c_light, border = NA) @@ -528,6 +603,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', legend_position = 'bottomright' } + if(legend_position != 'none'){ legend(legend_position, legend = c(expression(hat(y)), 'y'), @@ -536,6 +612,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', 'black'), lty = 1, lwd = 2, bty = 'n') + } box(bty = 'L', lwd = 2) } @@ -563,8 +640,10 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', barplot(pit_hist, lwd = 2, col = "#B97C7C", xlab = paste0('Predictive PIT for ', levels(data_train$series)[series]), - border = NA) - abline(h = 1, lty = 'dashed', lwd = 2) + border = NA, + ...) + abline(h = 1, col = '#FFFFFF60', lwd = 2.85) + abline(h = 1, col = 'black', lwd = 2.5, lty = 'dashed') box(bty = 'L', lwd = 2) } } diff --git a/R/sim_mvgam.R b/R/sim_mvgam.R index 3adcf36d..9d2db088 100644 --- a/R/sim_mvgam.R +++ b/R/sim_mvgam.R @@ -181,8 +181,6 @@ sim_mvgam = function(T = 100, trend_alphas <- runif(n_lv, 0.75, 1.2) trend_rhos <- runif(n_lv, 6, 12) - cat('trend_alphas = ', trend_alphas, '\n\n') - cat('trend_rhos = ', trend_rhos, '\n\n') # Generate latent GP trends trends <- do.call(cbind, lapply(seq_len(n_lv), function(lv){ set.seed(runif(1, 1, 100)) diff --git a/man/plot_mvgam_fc.Rd b/man/plot_mvgam_fc.Rd index e73e011f..a1af8d2b 100644 --- a/man/plot_mvgam_fc.Rd +++ b/man/plot_mvgam_fc.Rd @@ -11,9 +11,11 @@ plot_mvgam_fc( realisations = FALSE, n_realisations = 15, hide_xlabels = FALSE, + xlab, ylab, ylim, - return_forecasts = FALSE + return_forecasts = FALSE, + ... ) } \arguments{ @@ -42,12 +44,16 @@ empirical quantiles of the forecast distribution are shown} \item{hide_xlabels}{\code{logical}. If \code{TRUE}, no xlabels are printed to allow the user to add custom labels using \code{axis} from base \code{R}} -\item{ylab}{Optional \code{character} string specifying the y-axis label} +\item{xlab}{label for x axis.} + +\item{ylab}{label for y axis.} \item{ylim}{Optional \code{vector} of y-axis limits (min, max)} \item{return_forecasts}{\code{logical}. If \code{TRUE}, the function will plot the forecast as well as returning the forecast object (as a \code{matrix} of dimension \code{n_samples} x \code{horizon})} + +\item{...}{further \code{\link[graphics]{par}} graphical parameters.} } \value{ A base \code{R} graphics plot and an optional \code{matrix} of the forecast distribution diff --git a/man/plot_mvgam_trend.Rd b/man/plot_mvgam_trend.Rd index 87215f03..ee0c10bc 100644 --- a/man/plot_mvgam_trend.Rd +++ b/man/plot_mvgam_trend.Rd @@ -11,7 +11,9 @@ plot_mvgam_trend( realisations = FALSE, n_realisations = 15, derivatives = FALSE, - hide_xlabels = FALSE + hide_xlabels = FALSE, + xlab, + ... ) } \arguments{ @@ -34,6 +36,10 @@ estimated 1st derivative for the estimated trend} \item{hide_xlabels}{\code{logical}. If \code{TRUE}, no xlabels are printed to allow the user to add custom labels using \code{axis} from base \code{R}. Ignored if \code{derivatives = TRUE}} + +\item{xlab}{label for x axis.} + +\item{...}{further \code{\link[graphics]{par}} graphical parameters.} } \description{ Plot mvgam latent trend for a specified series diff --git a/man/ppc.mvgam.Rd b/man/ppc.mvgam.Rd index 27c8f3a7..72e264f7 100644 --- a/man/ppc.mvgam.Rd +++ b/man/ppc.mvgam.Rd @@ -4,7 +4,17 @@ \alias{ppc.mvgam} \title{Plot mvgam posterior predictive checks for a specified series} \usage{ -\method{ppc}{mvgam}(object, data_test, series = 1, type = "density", n_bins, legend_position) +\method{ppc}{mvgam}( + object, + data_test, + series = 1, + type = "density", + n_bins, + legend_position, + xlab, + ylab, + ... +) } \arguments{ \item{object}{\code{list} object returned from \code{mvgam}} @@ -28,7 +38,14 @@ number of bins returned by a call to \code{hist} in base \code{R}} \item{legend_position}{The location may also be specified by setting x to a single keyword from the list "bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" and "center". -This places the legend on the inside of the plot frame at the given location.} +This places the legend on the inside of the plot frame at the given location. Or alternatively, +use "none" to hide the legend.} + +\item{xlab}{label for x axis.} + +\item{ylab}{label for y axis.} + +\item{...}{further \code{\link[graphics]{par}} graphical parameters.} } \value{ A base \code{R} graphics plot showing either a posterior rootogram (for \code{type == 'rootogram'}),