From 1274087e275e409ff28a7d1310da5d53831348e0 Mon Sep 17 00:00:00 2001 From: Nicholas Clark Date: Thu, 21 Jul 2022 09:16:54 +1000 Subject: [PATCH] move all variable declarations to top of blocks in stan; add plot_series function --- .DS_Store | Bin 20484 -> 0 bytes NAMESPACE | 1 + NEON_manuscript/.DS_Store | Bin 14340 -> 0 bytes NEON_manuscript/Case studies/.DS_Store | Bin 10244 -> 0 bytes .../mvgam_case_study2_files/.DS_Store | Bin 8196 -> 0 bytes .../Case studies/rsconnect/.DS_Store | Bin 8196 -> 0 bytes .../rsconnect/documents/.DS_Store | Bin 10244 -> 0 bytes .../documents/mvgam_case_study1.Rmd/.DS_Store | Bin 8196 -> 0 bytes .../mvgam_case_study1.Rmd/rpubs.com/.DS_Store | Bin 8196 -> 0 bytes .../documents/mvgam_case_study2.Rmd/.DS_Store | Bin 8196 -> 0 bytes .../mvgam_case_study2.Rmd/rpubs.com/.DS_Store | Bin 8196 -> 0 bytes NEON_manuscript/Figures/.DS_Store | Bin 12292 -> 0 bytes NEON_manuscript/Manuscript/.DS_Store | Bin 10244 -> 0 bytes NEON_manuscript/next_todo.R | 245 +-------------- R/add_base_dgam_lines.R | 26 +- R/add_stan_data.R | 46 ++- R/mvgam.R | 7 + R/plot_mvgam_series.R | 286 ++++++++++++++++++ R/summary.mvgam.R | 2 +- man/plot_mvgam_series.Rd | 51 ++++ 20 files changed, 402 insertions(+), 262 deletions(-) delete mode 100644 .DS_Store delete mode 100644 NEON_manuscript/.DS_Store delete mode 100644 NEON_manuscript/Case studies/.DS_Store delete mode 100644 NEON_manuscript/Case studies/mvgam_case_study2_files/.DS_Store delete mode 100644 NEON_manuscript/Case studies/rsconnect/.DS_Store delete mode 100644 NEON_manuscript/Case studies/rsconnect/documents/.DS_Store delete mode 100644 NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study1.Rmd/.DS_Store delete mode 100644 NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study1.Rmd/rpubs.com/.DS_Store delete mode 100644 NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/.DS_Store delete mode 100644 NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/rpubs.com/.DS_Store delete mode 100644 NEON_manuscript/Figures/.DS_Store delete mode 100644 NEON_manuscript/Manuscript/.DS_Store create mode 100644 R/plot_mvgam_series.R create mode 100644 man/plot_mvgam_series.Rd diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index 62ad0084f98a2516c3529def56b3e6759a7bba2b..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 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 diff --git a/NAMESPACE b/NAMESPACE index ba39516c..f62a0dc1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -29,6 +29,7 @@ export(plot_mvgam_fc) export(plot_mvgam_pterms) export(plot_mvgam_randomeffects) export(plot_mvgam_resids) +export(plot_mvgam_series) export(plot_mvgam_smooth) export(plot_mvgam_trace) export(plot_mvgam_trend) diff --git a/NEON_manuscript/.DS_Store b/NEON_manuscript/.DS_Store deleted file mode 100644 index eaf6690acf5b0a949262f677ed328c16df7e9293..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 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< diff --git a/NEON_manuscript/Case studies/.DS_Store b/NEON_manuscript/Case studies/.DS_Store deleted file mode 100644 index aa69cddcecbbeee60b7e97d55540edd822209361..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 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)( 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 deleted file mode 100644 index b4aee0ce0ed1e09f28641b139c9a7116759ca70b..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 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)( diff --git a/NEON_manuscript/Figures/.DS_Store b/NEON_manuscript/Figures/.DS_Store deleted file mode 100644 index f3a3eff10beb84395b4303a250794201ba0dd758..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 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 diff --git a/NEON_manuscript/next_todo.R b/NEON_manuscript/next_todo.R index f9e89f1d..6e62af76 100644 --- a/NEON_manuscript/next_todo.R +++ b/NEON_manuscript/next_todo.R @@ -1,19 +1,20 @@ library(mvgam) -dat <- sim_mvgam(T = 100, n_series=4, n_lv = 1) -dat$true_corrs +dat <- sim_mvgam(T = 100, n_series=3, prop_missing = .4) +plot_mvgam_series(data_train = dat$data_train, series = 1) -mod1 <- mvgam(formula = y ~ s(season, bs = 'cc') + - s(series, bs = 're'), +mod1 <- mvgam(formula = y ~ s(series, bs = 're'), data_train = dat$data_train, trend_model = 'AR1', family = 'poisson', use_lv = TRUE, - n_lv = 2, use_stan = TRUE, - run_model = T, + run_model = TRUE, burnin = 10) mod1$model_file +mod1$model_data + + summary(mod1) plot_mvgam_factors(mod1) plot(mod1, type = 'residuals') @@ -23,238 +24,6 @@ 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, diff --git a/R/add_base_dgam_lines.R b/R/add_base_dgam_lines.R index a5d797da..61d4d145 100644 --- a/R/add_base_dgam_lines.R +++ b/R/add_base_dgam_lines.R @@ -33,12 +33,16 @@ add_base_dgam_lines = function(use_lv, stan = FALSE){ } transformed parameters { - // basis coefficients - row_vector[num_basis] b; + // GAM contribution to expectations (log scale) + vector[total_obs] eta; - // dynamic factor loading matrix + // trends and dynamic factor loading matrix + matrix[n, n_series] trend; matrix[n_series, n_lv] lv_coefs; + // basis coefficients + row_vector[num_basis] b; + // constraints allow identifiability of loadings for (i in 1:(n_lv - 1)) { for (j in (i + 1):(n_lv)){ @@ -57,15 +61,12 @@ add_base_dgam_lines = function(use_lv, stan = FALSE){ } // derived latent trends - matrix[n, n_series] trend; for (i in 1:n){; for (s in 1:n_series){ trend[i, s] = dot_product(lv_coefs[s,], LV[i,]); } } - // GAM contribution to expectations (log scale) - vector[total_obs] eta; eta = to_vector(b * X); } @@ -98,12 +99,12 @@ add_base_dgam_lines = function(use_lv, stan = FALSE){ generated quantities { vector[n_sp] rho; - rho = log(lambda); vector[n_lv] penalty; + matrix[n, n_series] ypred; + rho = log(lambda); penalty = rep_vector(1.0, n_lv); // posterior predictions - matrix[n, n_series] ypred; for(i in 1:n){ for(s in 1:n_series){ ypred[i, s] = poisson_log_rng(eta[ytimes[i, s]] + trend[i, s]); @@ -130,11 +131,12 @@ add_base_dgam_lines = function(use_lv, stan = FALSE){ } transformed parameters { + // GAM contribution to expectations (log scale) + vector[total_obs] eta; + // basis coefficients row_vector[num_basis] b; - // GAM contribution to expectations (log scale) - vector[total_obs] eta; eta = to_vector(b * X); } @@ -167,12 +169,12 @@ add_base_dgam_lines = function(use_lv, stan = FALSE){ generated quantities { vector[n_sp] rho; - rho = log(lambda); vector[n_series] tau; + matrix[n, n_series] ypred; + rho = log(lambda); tau = sigma ^ -2; // posterior predictions - matrix[n, n_series] ypred; for(i in 1:n){ for(s in 1:n_series){ ypred[i, s] = poisson_log_rng(eta[ytimes[i, s]] + trend[i, s]); diff --git a/R/add_stan_data.R b/R/add_stan_data.R index 90725caf..1134a8ce 100644 --- a/R/add_stan_data.R +++ b/R/add_stan_data.R @@ -228,16 +228,39 @@ add_stan_data = function(jags_file, stan_file, use_lv = FALSE, jags_smooth_text <- gsub('##', '//', jags_smooth_text) jags_smooth_text <- gsub('dexp', 'exponential', jags_smooth_text) - K_starts <- grep('K.* <- ', jags_smooth_text) - for(i in 1:length(K_starts)){ - jags_smooth_text[K_starts[i]+1] <- gsub('\\bb\\b', 'b_raw', - gsub('dmnorm', 'multi_normal_prec', - paste0(gsub('K.*', - trimws(gsub('K.* <- ', '', - jags_smooth_text[K_starts[i]])), - jags_smooth_text[K_starts[i]+1]), ')'))) + if(any(grep('K.* <- ', jags_smooth_text))){ + K_starts <- grep('K.* <- ', jags_smooth_text) + for(i in 1:length(K_starts)){ + jags_smooth_text[K_starts[i]+1] <- gsub('\\bb\\b', 'b_raw', + gsub('dmnorm', 'multi_normal_prec', + paste0(gsub('K.*', + trimws(gsub('K.* <- ', '', + jags_smooth_text[K_starts[i]])), + jags_smooth_text[K_starts[i]+1]), ')'))) + } + jags_smooth_text <- jags_smooth_text[-K_starts] + } else { + # If no K terms then there are no smoothing parameters in the model + # (probably the only smooth terms included are random effect bases, which don't need + # smoothing parameters when we use the non-centred parameterisation) + stan_file <- stan_file[-grep('// priors for smoothing parameters', stan_file, + fixed = TRUE)] + stan_file <- stan_file[-grep('lambda ~ exponential', stan_file, + fixed = TRUE)] + stan_file <- stan_file[-grep('vector[n_sp] rho', stan_file, + fixed = TRUE)] + stan_file <- stan_file[-grep('rho = log', stan_file, + fixed = TRUE)] + stan_file <- stan_file[-grep('// smoothing parameters', stan_file, + fixed = TRUE)] + stan_file <- stan_file[-grep('[n_sp] lambda', stan_file, + fixed = TRUE)] + stan_file <- stan_file[-grep('vector[num_basis] zero; //', stan_file, + fixed = TRUE)] + stan_file <- stan_file[-grep('int n_sp; //', stan_file, + fixed = TRUE)] } - jags_smooth_text <- jags_smooth_text[-K_starts] + if(any(grep('b\\[i\\] = b_raw', jags_smooth_text))){ jags_smooth_text <- jags_smooth_text[-grep('b\\[i\\] = b_raw', jags_smooth_text)] } @@ -284,6 +307,7 @@ add_stan_data = function(jags_file, stan_file, use_lv = FALSE, '\n// random effect variances\n', paste0('vector[',n_sigma_raw,'] sigma_raw', ';\n', collapse = ''), '\n', + '\n// random effect means\n', paste0('vector[',n_sigma_raw,'] mu_raw', ';\n', collapse = '')) b_raw_text <- vector() @@ -300,7 +324,7 @@ add_stan_data = function(jags_file, stan_file, use_lv = FALSE, sub(".*\\:", "", stan_file[b_raw_indices[i]-1]))))), ') {\nb[i] = mu_raw[', i, '] + b_raw[i] * sigma_raw[',i, - '];\n}') + '];\n}\n') min_beta[i] <- as.numeric(sub("for \\(i in ", "", sub("\\:.*", "", stan_file[b_raw_indices[i] - 1]))) @@ -419,7 +443,7 @@ add_stan_data = function(jags_file, stan_file, use_lv = FALSE, stan_data$total_obs <- NCOL(stan_data$X) stan_data$num_basis <- NROW(stan_data$X) - if(any(grepl('## smoothing parameter priors...', jags_file))){ + if(any(grepl('// priors for smoothing parameters', stan_file, fixed = TRUE))){ stan_data$n_sp <- as.numeric(sub('\\) \\{', '', sub('for \\(i in 1\\:', '', jags_file[grep('lambda\\[i\\] ~ ', diff --git a/R/mvgam.R b/R/mvgam.R index c2a2b09b..2bc5efe3 100644 --- a/R/mvgam.R +++ b/R/mvgam.R @@ -858,6 +858,13 @@ mvgam = function(formula, parametric_tdata <- NULL } + # A second check for any smooth parameters + if(any(grep('K.* <- ', model_file))){ + smooths_included <- TRUE + } else { + smooths_included <- FALSE + } + if(smooths_included){ zeros <- paste0('vector zero; prior basis coefficient locations vector of length ncol(X)\n') } else { diff --git a/R/plot_mvgam_series.R b/R/plot_mvgam_series.R new file mode 100644 index 00000000..e05b34cd --- /dev/null +++ b/R/plot_mvgam_series.R @@ -0,0 +1,286 @@ +#'Plot observed time series used for mvgam modelling +#' +#'This function takes either a fitted \code{mvgam} object or a \code{data_train} object +#'and produces plots of observed time series, ACF, CDF and histograms for exploratory data analysis +#' +#'@param object Optional \code{list} object returned from \code{mvgam}. Either \code{object} or \code{data_train} +#'must be supplied. +#'@param data_train Optional \code{dataframe} or \code{list} of training data containing at least 'series' and 'time'. +#'Use this argument if training data have been gathered in the correct format for \code{mvgam} modelling +#'but no model has yet been fitted. +#'@param data_test Optional \code{dataframe} or \code{list} of test data containing at least 'series' and 'time' +#'for the forecast horizon, in addition to any other variables included in the linear predictor of \code{formula}. If +#'included, the observed values in the test data are compared to the model's forecast distribution for exploring +#'biases in model predictions. +#'#'@param series Either a \code{integer} specifying which series in the set is to be plotted or +#'the string 'all', which plots all series available in the supplied data +#'@param n_bins \code{integer} specifying the number of bins to use for binning observed values when plotting +#'a the histogram. Default is to use the number of bins returned by a call to `hist` in base `R` +#'@param log_scale \code{logical}. If \code{series == 'all'}, this flag is used to control whether +#'the time series plot is shown on the log scale (using `log(Y + 1)`). This can be useful when +#'visualising many series that may have different observed ranges. Default is \code{FALSE} +#'@author Nicholas J Clark +#'@return A set of base \code{R} graphics plots. If \code{series} is an integer, the plots will +#'show observed time series, autocorrelation and +#'cumulative distribution functions, and a histogram for the series. If \code{series == 'all'}, +#'a set of observed time series plots is returned in which all series are shown on each plot but +#'only a single focal series is highlighted, with all remaining series shown as faint gray lines. +#'@export +plot_mvgam_series = function(object, data_train, data_test, + series, n_bins, + log_scale = FALSE){ + + # Check arguments + if(is.character(series)){ + if(series != 'all'){ + stop('argument "series" must be either a positive integer or "all"', + call. = FALSE) + } + } else { + if(sign(series) != 1){ + stop('argument "series" must be either a positive integer or "all"', + call. = FALSE) + } else { + if(series%%1 != 0){ + stop('argument "series" must be either a positive integer or "all"', + call. = FALSE) + } + } + } + + if(!missing(object)){ + if(!missing(data_train)){ + warning('both "object" and "data_train" were supplied; only using "object"') + } + data_train <- object$obs_data + } + + 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), na.rm = TRUE), + 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), na.rm = TRUE), + 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 = 'Y', + ylim = range(c(truth), na.rm = TRUE), + 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) + } + +} diff --git a/R/summary.mvgam.R b/R/summary.mvgam.R index 28747c27..fccebbfa 100644 --- a/R/summary.mvgam.R +++ b/R/summary.mvgam.R @@ -106,7 +106,7 @@ rownames(mvgam_coefs) <- coef_names print(mvgam_coefs) message() -if(length(object$mgcv_model$smooth) > 0){ +if(any(grep('rho', rownames(MCMCvis::MCMCsummary(object$model_output))))){ message("GAM smoothing parameter (rho) estimates:") rho_coefs <- MCMCvis::MCMCsummary(object$model_output, 'rho')[,c(3:7)] diff --git a/man/plot_mvgam_series.Rd b/man/plot_mvgam_series.Rd new file mode 100644 index 00000000..82b0f3d2 --- /dev/null +++ b/man/plot_mvgam_series.Rd @@ -0,0 +1,51 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_mvgam_series.R +\name{plot_mvgam_series} +\alias{plot_mvgam_series} +\title{Plot observed time series used for mvgam modelling} +\usage{ +plot_mvgam_series( + object, + data_train, + data_test, + series, + n_bins, + log_scale = FALSE +) +} +\arguments{ +\item{object}{Optional \code{list} object returned from \code{mvgam}. Either \code{object} or \code{data_train} +must be supplied.} + +\item{data_train}{Optional \code{dataframe} or \code{list} of training data containing at least 'series' and 'time'. +Use this argument if training data have been gathered in the correct format for \code{mvgam} modelling +but no model has yet been fitted.} + +\item{data_test}{Optional \code{dataframe} or \code{list} of test data containing at least 'series' and 'time' +for the forecast horizon, in addition to any other variables included in the linear predictor of \code{formula}. If +included, the observed values in the test data are compared to the model's forecast distribution for exploring +biases in model predictions. +#'@param series Either a \code{integer} specifying which series in the set is to be plotted or +the string 'all', which plots all series available in the supplied data} + +\item{n_bins}{\code{integer} specifying the number of bins to use for binning observed values when plotting +a the histogram. Default is to use the number of bins returned by a call to \code{hist} in base \code{R}} + +\item{log_scale}{\code{logical}. If \code{series == 'all'}, this flag is used to control whether +the time series plot is shown on the log scale (using \code{log(Y + 1)}). This can be useful when +visualising many series that may have different observed ranges. Default is \code{FALSE}} +} +\value{ +A set of base \code{R} graphics plots. If \code{series} is an integer, the plots will +show observed time series, autocorrelation and +cumulative distribution functions, and a histogram for the series. If \code{series == 'all'}, +a set of observed time series plots is returned in which all series are shown on each plot but +only a single focal series is highlighted, with all remaining series shown as faint gray lines. +} +\description{ +This function takes either a fitted \code{mvgam} object or a \code{data_train} object +and produces plots of observed time series, ACF, CDF and histograms for exploratory data analysis +} +\author{ +Nicholas J Clark +}