From ff27d5d3fa6d14ecdd22d781acc2c525892d7c8a Mon Sep 17 00:00:00 2001 From: Gerda Wyssen Date: Wed, 31 Jan 2024 15:07:48 +0100 Subject: [PATCH] updated downloadable files, fixed typos --- datavisualization.qmd | 37 +- datawrangling.qmd | 2 +- docs/datavisualization.html | 42 +- .../figure-html/unnamed-chunk-4-1.png | Bin 25504 -> 26043 bytes docs/datawrangling.html | 2 +- docs/downloadable_files/faux-stub.Rmd | 939 +++++++----------- index.log | 2 +- index.tex | 86 +- 8 files changed, 468 insertions(+), 642 deletions(-) diff --git a/datavisualization.qmd b/datavisualization.qmd index 1bb64b9..91f27ef 100644 --- a/datavisualization.qmd +++ b/datavisualization.qmd @@ -20,19 +20,17 @@ d <- read.csv("data/DatasaurusDozen.csv") |> d_summary <- d |> group_by(condition) |> summarise(mean_x = mean(x), - sd_x = sd(x), - mean_y = mean(y), - sd_y = sd(x)) + mean_y = mean(y)) ``` A graph contains always: -- **data** +- ***data*** -- **geoms**, visible forms (*aesthetics*) such as points, lines or boxes. +- ***geoms***, visible forms (*aesthetics*) such as points, lines or boxes. -- **a coordinate systems / mapping** describes how data and geoms are linked, also colors or grouping variables are specified here +- ***a coordinate systems / mapping*** describes how data and geoms are linked, also colors or grouping variables are specified here Further components could be: @@ -42,25 +40,25 @@ Further components could be: - coordinate functions -- **facets** +- ***facets*** - scales -- **themes** +- ***themes*** -(we will only cover the contents in **bold**) +(we will only cover the contents in italics) ::: {.callout-tip appearance="simple"} ## Good to know -- It is easiest when your data is in *long* format. +- For plotting with `ggplot()` it is easiest when your data is in *long* format. -- What variables do you want to plot (categorical? continuous? ...) +- What variables do you want to plot (categorical? continuous? ...) affects which `geoms`can be used. You can try out what is suited with the `esquisse`-package below or find ideas [here](https://www.data-to-viz.com). ::: ## Data, geoms and mapping -We start with entering the current data frame and add geoms and mappings (`aes()`) with arguments such as +We start with entering the current data frame and add `geoms` and mappings (specified with `aes()`) with arguments such as ```{r} ggplot(d_summary, # data @@ -80,7 +78,7 @@ Examples of available geoms: - boxplots: `geom_boxplot()` - violin plots: `geom_violin()` - + ## Facets @@ -101,13 +99,14 @@ ggplot(data = d, y = y)) + geom_point() + ggtitle ("Title") + - xlab("Value 1 [a.u.]") + - ylab("Value 2 [a.u.]") + + labs(title = "Title", + x = "Variable A [a.u.]", + y = "Variable B [a.u.]") + theme_minimal() # also theme_classic and theme_minimal are nice ``` ::: {.callout-tip appearance="simple"} -## \`esquisse\`-package +## `esquisse`-package With this package you can use the data frames in your current environment or load a new one to try out which geoms might be useful @@ -120,9 +119,11 @@ esquisse::esquisser() ``` ::: + + Further helpful ressources: -- [`ggplot`documentation](https://ggplot2.tidyverse.org/) +- [`ggplot2`documentation](https://ggplot2.tidyverse.org/) - [Website PsyTeachR: Data Skills for reproducible research](https://psyteachr.github.io/reprores-v3) @@ -132,4 +133,4 @@ Further helpful ressources: - [Here](https://www.data-to-viz.com) you find ideas for plots. -- https://rstudio-connect.psy.gla.ac.uk/plotdemo/ +- [Demo of different geoms](https://rstudio-connect.psy.gla.ac.uk/plotdemo) diff --git a/datawrangling.qmd b/datawrangling.qmd index 3ce5c37..84d5654 100644 --- a/datawrangling.qmd +++ b/datawrangling.qmd @@ -75,7 +75,7 @@ Data visualization and analysis is often easier in this format. If your data is - each measurement has one column -- each entity (e.g. person) has one column +- each entity (e.g. person) has one row This data format makes it easy to spot missing values or outliers and count how many observations you have. diff --git a/docs/datavisualization.html b/docs/datavisualization.html index fe4c9d3..c0c5e04 100644 --- a/docs/datavisualization.html +++ b/docs/datavisualization.html @@ -238,26 +238,24 @@

d_summary <- d |> group_by(condition) |> summarise(mean_x = mean(x), - sd_x = sd(x), - mean_y = mean(y), - sd_y = sd(x)) + mean_y = mean(y))

A graph contains always:

Further components could be:

-

(we will only cover the contents in bold)

+

(we will only cover the contents in italics)

@@ -269,15 +267,15 @@

    -

  • It is easiest when your data is in long format.

    • -

  • What variables do you want to plot (categorical? continuous? …)

    • +

  • For plotting with ggplot() it is easiest when your data is in long format.

    • +

  • What variables do you want to plot (categorical? continuous? …) affects which geomscan be used. You can try out what is suited with the esquisse-package below or find ideas here.

4.2 Data, geoms and mapping

-

We start with entering the current data frame and add geoms and mappings (aes()) with arguments such as

+

We start with entering the current data frame and add geoms and mappings (specified with aes()) with arguments such as

ggplot(d_summary, # data
        aes(x = mean_x, y = mean_y)) + # mapping
@@ -298,7 +296,7 @@ 
<
 
  • violin plots: geom_violin()
    •
    @@ -314,7 +312,7 @@

    +

    4.4 Themes and labels

    ggplot(data = d,
@@ -322,9 +320,10 @@ 
                         y = y)) +
     geom_point() +
     ggtitle ("Title") +
-    xlab("Value 1 [a.u.]") +
-    ylab("Value 2 [a.u.]") +
-    theme_minimal() # also theme_classic and theme_minimal are nice
    + labs(title = "Title", + x = "Variable A [a.u.]", + y = "Variable B [a.u.]") + + theme_minimal() # also theme_classic and theme_minimal are nice

    @@ -335,7 +334,7 @@

    -`esquisse`-package +esquisse-package
    @@ -347,14 +346,17 @@

    +

    Further helpful ressources:

    diff --git a/docs/datavisualization_files/figure-html/unnamed-chunk-4-1.png b/docs/datavisualization_files/figure-html/unnamed-chunk-4-1.png index 2382eb24088ef1b39c00ee22542db7795d53912e..52ae9144d9aabf05cfa9528990030ee8b4eadf97 100644 GIT binary patch delta 20517 zcmb@uc|26{`!`&Z!Hh@@B3st6XJ4`$`;uksQmBL}QCSkjX+oAP*|Q(}AT5@XwGu^P z>>-3?5Q>p~JcqvD-*dm7`}Mr;fA0UxIoIcNKI^r-ulEJ|c)0%a;m1OU6%&N<36Ft< zz62U8vn#?$_{02QBSSNfp?^tZYwP=WrVWg3&sjg+X9KQW%nB5}cnl$83Pf4*Pv%F( zUpXHwxwyP{{O7N`F7uUT^DefY{(HF--Gi!k1MHtWH_i`_In``mza7S`Q~8y`be~!0 zsSbff$4CT&Jr#E4|MySu|0wZafnFoXLdkc0wDFFI?`%m+tGf)dPNZ`(J!KvA-$+9U zhn>&~2VAHdKg`g{%$FnE)&Z?I^_%Z+xIcGpFos8RRQvxs94BND1)N#*r1J*3=_!5d zH(Z)l+n$fM8PcEfNk&s9+9T z{AS76Qq^u0InzrQett?E8MNK~vvFF$Lg|^Ms1Cj5eG`Hv*M(ZG+S+7h9lty)NcXDd zq@wyMJuxn%8n+HOb(-aro)kOHmRr{qfiQd+Jdwwy)#&6ShGEQWcnmp&IqKj~B8FV> z&TU(;BbhlLlY$rn->i0iM&Hw9J*6iqPE*%BxB2sT_^9VbXBu_xZ>-DsjHHN8q+LI~ z-PG+p&D$Ln5_>)yz>L-}K4Qb> zqulD)cxIg+3heOaNz%a!sw0;NpiT*~{@?yn;(zZrNgB~O_p;}x_fwi{!Ud*SE_aOS zH61b1p+uy4k{BIJlDC6?QXEYXXevPn3E9gm$#o9CbEAHD+%V3i;Y7P|IV&S;^g~W$ zp9ysAwJ0nC*3T|9XYWpZG?G4-jemku_=0_+ga}@{d6%8svN+1Om;=5D-kzS5A@Elz zjquuuo)Ez>%A<>;fW7I6(JCN!>`az7{H2Oqvk(+Ay{d=lh1$E*c_99Yev5%EO05lbi6mQQ>xLGp8-Az5*MrTw2c|dP?Uvco?RBl;ZMWr&*Mo7 zXURYHwj99{zF47T;b?<7;hhtLZ?5>~Urix(ihP7x6l390};M2T6r<&W=PbWxF zGR2Pzi7*M+wirR3$hy1K$qKpV^_NwsXR>A-&D*g;+$Z_enw0cm6~VU6a$jX=+aI;Z+;# z_TCGu+bzCwRq)M)GQh2$99hmv-nn((?NWG&Y^f44p{d^xdM%o-^>uIi_lCt$SPXms zFl3SwYj(RD8D__|I5G%tcMMeac)DEUO^{A5hD;pP$^NDMr0(V1XWw$YUr+4a zS^BW>l81rfR-cK%8G6ca<6Oh$ri5$E zkKP90u9nl~SAl4Z**SyrgAWs8ZZnhb{oSd$tA51G5MLm^Gk13*!~H3ZV{%XNz$Kvq z9-5#x(}8bNsU z5fv8N1}DLhlkNDs`xnrsWp%i5SD|=TD<L=G$ts!!Uyn7d?f%{@NxvZdtXa8vFxT4~ zyo&PDB|o-=E^g=r@iQze^E|8-IobCZH%yhFrtf^dx6w5QK zKn&bQ$(U=KK(@Y)wvd!83l{Dd%=!|0&Q=v$Psol#>k!^qA30a z&jDV`XtU9Vs3NC^i-A!c0^(d2V6TX-y9AzS8Rn*_kqtnfxm0w9LtXDg9Kv0Kc6pt$@z8+NNoM73<= z{7k5k3gICHyu(4vasKf+Mq17n@FaZsT|aj}_s6a}3!cL9moCioJ;9P>4G=}Z@hmjo zh21_@y98}FGl@>JpUkm=;r9}%F;fWB!W)&j80-l&=B@iDvaMdZqPcmsyI*De_WX5U z_tY*GX+53uzltERS;DwxkR=-+q*igt2eI|cN7JZXazc3JJ-Y)AC3@-#XCRPj9s;CJ zvOZ9bA*PyZ6g~h}j~thgtTTHqF#OO&0g3jDI}6Q><-jAFxn-1=YKG}~?Del{(&ABH^Ou+4F@_xQK9ncm16(1Y?pu0(dy_3>=UOMrn(&Nkl%i{p?7 zbQ-~Du+JI- zfwuOi3l>j-m-gPL=jp54o=LkgU!{9`9^HQ!$px*f|1h9d^OIIidkR%|G)cR0v9kPo za?QkQ@zwa?DvbX%-p#zOjNOs8-Gxtx(hKTs*r53%U3G*U z6~!lXW%f{8v8;9hPsnS}!Rrh~(XzdB@JZUBhG4HYXq~^uYWxx~Gb{XL@X_@?F?fJT4mQc@ce59udk!}YtCcv>=^9h z4g^xkOQ>KB$z?>{ka9)%XE_{GF4Qh)Jl4p3Fz@#gBSr>dB(+eZGZBKK1A}AGe_m-2 z@P!;`*zeTcz=@C}K4NR@_-m?>VCmB&P`?0s&KyxP8Pr7}H4Nr~t} zS#8PiQ3~7TIm+s28V9H*I5Yjv&s0DWHlk<@eK{-q)_jFOktK!(-qTC|4=(8(E;M(w zdT?;xpPHXG#UlRb1=${G3)pVAgyoBU!GkS5!uyjXyl4)+q|=k^Rz0&L^4|w96LE}v zAy+VN)pJMh3P$80%!F~;c3bt8F(HK)7 zf2D=-7OA<7`GTWp%u!^9!HX~~?sOE&Yp6Nm{I!J?9rzc9cm#|M13mG=2rre#k1-zg zOqwITztTp|@1F-5&;@S=&y*p+=dqEgYVl+#};$Ce>mX zaoA;@IHhcq2@Qj29=PwvUi9bMOC%9BTB6fXU(%n(b4t4&MLr>CT1-ob*CwUPg z?T@gQW;n>VlHoRLfOGjJGnm6gl))8(`H0pDT@V6NbO~yoDu#q@N}!KAJF2oo(JxbB zI!Fu$_TG~Wge|<+S5N~{OD|vk4IQf{pWjmT#P`)xwnXRt)&5S&V^-7(7`>K6v&^3jUAiMK`#opJjb!GuO<>Oif|M z6!A!{j}yb!C2&{$1iqxertEEpwz7S8@Z)vbFyN8qF&5-VZSHe3h)WSZjk5OHBLw>D zuqS{_aL{kF!j4mc1+13kdVkJHb% z#GVE|?Fi(N6MEuTD)asflfRE;N(U_TDm#?udH&k)c?wh+i_?+^TXmyW+B8Sv+wL)}J}DKx^1Z!hN`a+9l_ z=tGF3OQ;4tx9`8xE^N1ZJbDNqn2V_x#a;1M5@#&CVCqrL;HMtje*R1}Q&Zzr)1(UE z^?3I=E{8ADA%n2`Gq-dw5HG5Wfz;6wh5R1BKs3z?bXTYlyYl~VmAZOWIm zeH_KafO@W?!1WC}?URn{*5O_P02QmVQocwbu9zQ zvro@KEji5M<8lNM1HqWrkHhUe1i1|lGh>&?w73jQIs=3gxGKB@iPAh$HiT5iwC61V zzDHMVJ`YC?E@^gK1VV2MP6D0=vnNn~aD-ET)*=$~ytdHqzuYT7R$0OCU>#5beB4*pW2{on=bR?`gQZo8k6JVE zeM$OHOW~liNW*_q(XATAPYAJ_GlmSe7;M@)&75`*?Q+6#ekat&Z?(cz-hpR0jmTE- z@9Mzru*Jw-W56g5o&E-E<--D`3t;ins%J4u_YoQh*|rUa-;ZlWFt6$pJ{S_*51|x> za#5bPU$C{8W4ME2VC-XbQx@9WtIs^MchMc>>O1sdl5exBImG?C`ybWm21!z|jGYL{ z#s-CPq8Kn~gfPz~)w|Gc_nAodH6PM);oTNYy@f9L z9j|obXX{j}k@T$S9nNz5NY8%iw_M)PYfn1$2(*?$q(oLaWdo*QBw@-zJN5_5WOp>! zw6fK=N1r!NTLiQe^5t#HH@{54LC0w}wKX46BNHy@?4-L8Uiy24lmNH<&oVU5g&?b-N&xi%mii8`k0O zi^%g(!I~He1BW$#-@Z8X`!eh_kU)jv=Dok^H^#=E`ZLd|$BKE5z4u|S2ywc&6T#zt zqkeqxFfzJ015k}Met1P-_gaZ8k7%#Qvl!LGu@kf!JEyl=n3@jNZ6_hZtu3G|BN+i{ zl1vEQqJ_j8$zRL*tyy1Bb`jKIbAWeOKY_(5up!ZHMh2S-BK+aIjfJW5NHAkUhZTNN zL!0pXL_4g23q3%K1xDECmLLpfv|&U%9@~-PsR&;R*GPR0HNP<_6|P=`y>u!_G`s}<=S~@{@LT(+oJNB|bV_RUy86x14 zvnRQ;yccvU)Tna0_$-w7#1=UpfzWpiaPXc)fV-641t#em80S0RZBPm)s)~UlI0*lX z>9f>%C~0ApRggh9)udg< z{pVQzwmG-7!!gPCOA!rgmEua^{|FTe8?-ixh_+0;E0=D$vCE>C|f`s=?413SF8P>tHV}JeM zR0l$r?F4cDeNPheOB43OQ}AYCQOkTp`tL1z)ZM{M8E}1TVM2Rf2Y4|kNaZN?l#H<^ zyfPvD)yvpiB6NJ8{&nAh9lu;(Y0knbemw=T7^$+>c!Ui$dTDzvO2l zuC5RNQLJs`^o;D?#%dNMD0}wVjss^qz7%Lo(oCp=)q>TqLNT)lufHH8R&}FP>+s#i zT9)wB7Sc$y&1?h@bWI?bnE}COn9RG=pF30E*uaAPZnXb$xAjZ_b<$p%dR!%Y+(>$R z;v$PH0IfWhdyudFIELuVE=&R+^ckV61-gDjufov>><-i; zY*iwYr(XT~Dn^Q~O_DHsoe;xr;kxxJUXvOY-D4sx#ZbLhJD4o=dZx{sVD};E!lhTS zSCUr&0=xk0w<x7PHBZ3MK=^K+07(U& z7terdn^OsVy$uNL-REr+EPI@P(2QIH#d(ovJ$?VUKmB5K32$)_s-FdmH*sWvo3IV2 zzPBeT70Zg55|5!32qp-r#>NDs73}h?tPQ#BF)f|U;!xnGoN6j4%vX!JBT-`o{aG=S zxap5t(;eNvESb_D1%lWMdy%kH^&@lw%W%!f1aHP30Q4? zeq9NKk#&#diKVx9iDYxK3Rx8K(b(AE`KmYxjeH9CovN15lo(G_W5rKOrIlhQ7ff@` zQjfiF4!kpYn2NChZIvs9NX0OGNj%IV@A%nJk!Iv0)C#7oetd1eT;CE!o}6g0aG{+* z+S+Cp(~S6Vu_8y}eo28-ri7COztTa}&y2TdT`F9b4L_nh(yZ+^{sz@5cs z8hPF(3VHg_AICi47Jk8?Nw9j08YewAAllz9K*Lp$@-#+iT@)F85J}GGrvvrWiou2$ zrT8>>^NG!AE_`eAD!Le9h4*qu$8N(v?(*^@b=Sp63J*RST7WXmx)aO&uJ?1o+8loH z_da7<(cAhi0GoLL@=cW>C0&-32c_Pa(Bv1?pCz46aQ1m?bH&ESi3`Do+c7JIO=VW| zIkIxXWjboc7n@BXDL$ts2(P#YfWD!rJ{+C5;oFNGpwhp_j+#xUkt11DLI+cNbd2}I z4fbAU;vuxc;ZneQEvdTU;w5yo>h6&maoqCmi{evmx&)Dn`cPF;|8OtJt%vB-!s#N| zIBoQyTUZk2GGB%xjoKSCHJ7ZuCvPdo`V#?RT6kHLeXmCmN_LH@E*9*3a@?pP;JL#wA=ROKJ{di`&Z`ec|RVc5{MB&aI1sE!-0klAPm78~N*6)tC)b_M2#xu#pNiXog zEuPIcOl=enagt-+hZT=YNxwDHxh0Hl_aq({hr%A#{>WD1&eOzZ@5Du z6OQO9K`oO*#d_zj*xfRSfTSn+t}D50Jr8L9otd^5EVxzm+v4>NkXz+QHcdLatMXG+ z_1KBCA{e~t^ifG?s6v8Eh$NcX84w^dKCA`OQvw2%ms>}NbN#HKF%_DTn=!(yOag53 z?eBV#BMtf5h~O82z5y)crJkQPch1mx_+Xj@-s#j0mj+H*(VM%9y9GVwkS6_Qi#<;M z%r~()rw`lBHLV2d0G+N!&MfG|i_N)tszML)%rOGLt*IfrPL*W=XOa~RAVdEg8Zc<< zV`s~B`%lY42C9hF4T11#c7l?w6I3F2-WX!+S`@>Hb}ucK$ZG6OPA)9J=GtwOYi5|p zbP$H|=W1kK`mK8!Sv_HwiO$O9Ebsu9>Eg_P-eci@2;K`i(KlIc)ZM8OD;AkUl|D{h zp-bS(m&E02wQ$3n(<`?bl<39r{Tg7^G7Hg{UigCR2asHM-;8p(6Z z0)+Yb)6x_ng@=^@-AC}=V_nXHH9}pPcHx+4fU?W>sy39{R6l$St8HnQ?M$Hrqyrq9 z2~EOfO_wz_sRbfOI2Uhjv0rxTd*W@Ot(-ZuJia+^@4XjXrBa5ns{<4Asl0OxK?uft z*vmql_RHfkQ-T^RG-nxt0ujww>s&+e+!z5i|JF#ST#@-gnzh&$b;qhkYb zKJ+eEpGWh(U8QxGd(%VV$5mm=kK7F=9nPXQ9Bq5WBpFXdrQ~~u)RpkmAhNxNKIP52 z%+T~UHb(@#B{z-}K=fP=PIjJV{!QwDTBO26Wr8k@K+cU#=I8}v!5b)SIJ#Cs(13 z*_=)#7gtQpAszDHT-h6w-{ywqo~_&2D^6 zv^gKH6xWv@$L|~&HZmBvz0QfjX85E0mb=fm?B+IaO-$5{MGZw4WRO;&5^gp{MC>?= zTyN|Y%JzGaE zprVX29BKe=SOCwvF9sdbLM{f7`r^(uYZvy@#`+aF=&(I8W`6JhE91a|*BRd`IHBrDH%dL_Q&O$gzIBw=cKv}pa3;28cCqGD$tDtRUarNA7L0mQ;4 zmFl#?n){7$*YXW_(tUYu@ROE!z$WLeX%Je$NWEi{C1S69-XH{v8y=kT#3W-wL_;T^ zpM_xku`@GQy?T7-oWN%;z4u8%q-m|s%*2P2oQzb{uR$SEQvveaA-()oTDa~byMwoQ z$yy=XOYW-!NR49(bIe@K#Hjp+7m-&@Ym~3 z_Sy+&%|}m+TmnAF(l{2}%ts8|^xOBz>c;ZoQ~zE@mLbuYvEl^X?vKNV=Y}Q8 z(A6x#H^x5=;e<#=NwC3%-7a4+cOwZAn}972B7_-frXB=kk?w@0mp{96P-cw^&rX;` z_PdtinWI2N_V1N=(U85N2D-8Yf|fsB(4}Nj@wmlWctRMQThb#iZBll%*Wm<|!Nb|y zvLlUhruaVj*&o*af~Wz7gd7CVY2KJUIb_-6OYh2{KHUdtLdXe|$*Qy=Lr*c%+% z1vAe~NFDQxNw{#W%$&}Opwwu?2t9J``yVF+p+U2dVa0i4XN&StmG5oDB#cei zWlP%@rc1rLevjEW{OMP4Y~LHlP3Q&P+@SU(b^h!c>-3=(8VW%3ep?d8*L%VX`kwfm zNFL*XDFL*9oi@NRm)hogFP86Q(K&{FbWs(B&iU$LIAMX<-Iy}$7m2qKQ9Iie$}ZHU zwV~$Fzq55^eGAWjoPjX-Gd6T;Dycu^%L?UWe&4zc%%c@1!KyTbSL{z9g03&fU?BM? zC!7=MXz_FIQAG^U>H%48< zL09zxJT{602|*gerE!+=RiRpXfsrYpah!VBUj9@iaY>8)N&>G=1Apa?Dhu{)2%7{) z0?z*C#(jIoe71AL{4B^iUHg(kA)g&pcKn}qhtBR0K*#BwE&ZKrVaKF?V0$BpbuIxY zw-yo6?$M(e5$RIzV4aVCX7Hh0vZqXlZlo1X#wIc$K(lM)5pJXjLHwf;1P5xY2#p2< zI7!?AXQOeD%-)@;b+Gd;H1oY~vnbosRS{0D1wKURSHbbBc((?Y5VA%6E))@JM1QmK z-ne-`Jiyr)s)eF*D`=tBYDA{Mzr+UU?;zmM+6USCTU(bZGJEWa~B zqfFSa{4W|KC3GCVdES|su?aCg0>Nk&Gd~WZ%IXWi^V`ce54NOtsPXS6sC$hquhljK zCI_t6`cA|V!g6rTL^A1b7qMZ>P8deVxp_ra@T)%;_xRrOd6htF`jlr~D*lg9an!J! zIX>bq((Nw1&RT^j1V-E{5X~Ht-ZSiNR&ygSFTBsTNfoUrtK+@&G3>EW_V zwj5@A@gC8T*H#j8cz|Mqx@>O z0~^?a`mu@B0kIYURGzV5KxyRh<9hNIGkVNK30zN%_Eg~r^knn+v6kjHX<{;yYd%jW zl7TZ@gE~6#czT=s;v7!1m6GTfG!V=>_il3A&_XHmOn`Dkm8DRUmYM*#ExpscB=gNh z^J_`ODli~;toof|Gk9dIzW&ZF3Bx9s3vf;`(ei(XUed56+5&R|Bjqf9-ub-h3LMPu z`#!9a7%c-IW5Q475uMjzBFjqTv)ePUcMH5r-uHQ2cBFS=e4`UCSV5{dsKtyORA$Xv zcXHd*gJ|shbdqha@DH!xn?k@B0u0YfdWk-GaMkpEcr2x8>g+}?U!v8zLcftT z{P``OlVUgrn5g5p7dmgW7e(oHi=_y^X;g zHFWhRGb%O#$O8vAs-}z7zbjDO7cY2`&LnEo<^Q=TV5C(XD(Y+n($t-C>R7#M`sPBH zws@|&Ta0&0zA@fC+yhD}nB^oVzfacf2G$#PX6&(~3(qw0ZxCizOQ!j)1Yl3efROr= zTY}UVjH8~LR;P8szeJc_BTgt3+rYk&7edNh%?fJXE7z6WNVo$lxEE@M!9?$XYl@Ky z=S4eN$3%KXbwNV@RTg~U;3Wi%eYPCxO$$jpX?#e6jz^57Z|a)y+}0l21%f$p`kr(> zHGxVd*Uz|=uU;2{1;lZ1bQ%9++qaq*IbRg6HI&k$B2Kfl65ujyb7oWh;tZ5#J6pVi zkDjy3#t?Ja4MUtbWAg#1}dIA=FvuZ87y%P;l~# z`0~Z20r68B1r;Vf9VQ~Z%^Q6F7MjPtk`3I@7a(91Mk0=!#5YZtXvM&B0wKPZB z961TX2hF3pNpqfVNBd#QwNoDMUVNcksuT=1=6&T2$1;XS&H0Zy$dSLU!D0KZ|9h;| ztJLsof-pK)IdJ5_M>4!xq4?BVkkjARqvX8Lad_z;o>icM!{<@P^(!$~`Yp3fd>csG zLTCkA7*WRy1Y2L`w5=3eo&TjwJPx{-_B~OF#K8IK$|vVhjU;7p@vR1l=u4ezo1FB>ubV(}O+3AME=atr26a) z{zw;vFJyQ!kr%dr45eLBMemyj8GS2Y#)^Yc)cpR}9?$Ect4+b%p4n|S@&A!LTiW$q z@opl(;9ue`&=@dn7q??&QmsfmVnLX4F!hgX;Aj#S$Cr#vH+fD^oOB}u^`u&Al_40p z@1QXeJ@bD3ZRDa-+N9ldFPs~S@mK6KU>Extf5?pRkVzdG-J**BCY*Sp1(=S`luqzD z;-`&Z1Q+AVfFm|t*cac~+!XG77sGM)U11+T{PYt)E^1mXAQWM1I2I$o;A{y4ItHeR z=Yr^72n+#_ggNfRZw6v?B#aN07MOR)R_;6XDlte4+$!*s; z*6jAfj|-2zyhwm839EXjQ5`blZo0(+S$LbDctnbN(*beqE8~k;F zM#-4fCQHb!?o5en9_-O3^McELL-)^r?8gThNTB?_eVRg!?97M3>KF%s_R&(`?{0D3 zT;R99gp97govGk@7MkNE`i$^?urpu;2+mlo_!aeF##;aWH@prImen+%%O=30`|vMc z!?@BI$NPvzoIdouV&OP(Iz8Zph>Y2}`K38JTL075FOn@;aS#WVr0_dJ%~Q*5v#h6Ogf`(5yN@A%Q-=8MF7Yv`gm;gvZ~979f-@{ukxA0E_p>XMyqexXkz zNLV@|*AM5rPk#hPJ+Q~oe!$#I@p4#1c1U>^#J{@~!x+hbdB8>#j( zPb;s*PuZfrH0GLwIS1@cxDEvBnZ=sRT(?9>eO25n>x-gLWBzu&e#v5Mo#hOl?zXbk z$S;{$(QxMdAn~5c)%f8MJN{jh@YDG~^n=$h)J98}aMMARJP1q1aq;Ng+CBc5gm4mm zs(MZm6f#EJwBfm|m**)k-Koba`_BVw<=Cr*DW8uiD;EWzlnWjTGSKR1*1fC9iqe-!rbZ4vS`Cd#5e5|Ka_~WEl|D zPN?U*R^9>?;)bR#iEa$n1U}ZXEM=@NEYK6aaAe?@1sW5TIFG?xDAz80-_dO9n#j<< zX`p^IUQ4$HXj%WfL&uIMf?~qL;ulxk;21n?-=nvgtyMW)#B!skN1@nC;(|yo2ffXS zcRMOHp}9@LzJD5S4n2_NIDY+eA)p9vS39xmR>ctu-SKd+^X245vllkg{y8-iPzkBmv#fzS47dj)HM)ZOP z;K)8zGggdecFBgF=;_+p`QfIS>%0i*(7%8^HHvbp&b_)6I9Wt#tC0$(76LtR6NhZM z50_u&;QXlyxt9J14|EDI%jO=#c~Pr#Exit&ScP3q(x6cWtzN4o2-X?Mb?n>Z0~VlT zHc2)7e|>N+`|FJwI7eK62Bx~hY*w1{*^15faMtSA(mur;zPNaQyWmQg&?X{Wz3@5& zW^?cJW)a$Uvmesn+<4Q(8RbOyW-!9YL#nBzRxTRG9!Q(*h6XGquulTdKQd8lVOIXq zWn|vcDGu)+11Rc!V|Pp8GQ?@ zCD$jQ--ofm_TUrn?jtbxj32!}jz{Kw0`9;T6Zdyc_BfeWknoCw4DO8G-X2`nRNL9v z++8FbdA_HVbXZ>u1U1E1D{>dXz}wx$fJPBL8(wE`Rf(&&FyPPJ=iZPZdnmeCHhYTn zv%HH>wv@O)@pfdO%>Av{do`-D_nrKuWb}rKU2p1H7cm4~Lgm!0a*LgS;j!nUsOAD- zz_d_n2FQFa&Ye?+U!d@EF;E8RU|#=TW8>|Kz3rRZN7_Nh4xc$G{qG7Qdce4f^#}3g zfJWpznu9--1u?JABSzaSn&7?7uyFQy0D7 z%0l=HX0`R7bVb;*rJodxXi5UKnwz!O)#Fd`$Si!kB=*xWI{UR&U}t@ng89&taV+I? z6xUD0B2Q1cw)%+9c{zWYt&P(`6hFAq;A8#X&doFRiOY>cmy$Q+xd5N6XRX7-iMWOC zcRd3qzFjUAAe9)pX=*aQ|Jat;pWVqEz3kF&7HW`a zt?Rftck&x~xtuFMNC9ucz@X)d9{ue;xGEwdbtvs@z#sqP_DvFT^Go)1U8}SHapEk6oUO%(>E2l;t z#3{#aop}ox$X&CBI7k_Y>VPlA!=QAH-XmY?{N&%0C4HT%Pq)Y%i7`Ja*Azt1HpiNl z-wW!8pm$TCJ(#Ef1LrUJQjN5z@)YdnOj zZhBC2wtEKl3kvHc9QCZ#BRCzMjo5P!q1q z2i{>>$nHmxj%={`i=M{}M zRj+uN@3)ugQ)OI4(DhLuwLbxO90IIYgjQ!Ch8o%aYnlLj)R;5r1hQE^t(3%tMZcxV zxLv&DI5=rx?pg*r6ncjKK|8AR*VvAjB>^XRP$`E{{kxYqoX1B$yvSRykob%bhQDM| z<{LCWNlslun56*B!YDS9KJ-rH+u~En0vL%T;K9N#Sb-DKo}W~(e`9Xpn7_L@@BQxwnqwEP*Ux$4|Jy@9+V7{rJ|o96bJ1<< zUNzqqE36nB5!si83g(d`D<8aZyci-EJS#L^g8GSeVhBN}MAleIVIv)33)6!2CKN>Y z;6o1r5GU(gQar_7{*OyOlX4X;G&%6@ZQcZ|$1Rl;#vdibI5LwH!hj?#31#2IKey$% zh^ebK`>bB%<>uxE@)?d_qJDATf=i0psN7B;j&VHl{!Q+7fUZ5G#BvA0=VFs8y>Mo6 zD{V8Z>(=n;>cZuP+N31QN!kL=l<+|W_^+qN=FfGvTc~atc}(BbVr0nY1xHAC%=SsZ z=nr?mH2=j^o-cQFob>LIousEMjCCVa#Bh@Egr|3%`fD=(o}NbeAD7P$P+jlE;630; zw{6!GJ*^@0URWp7a5+fUy>rdt)_V`za@p|-G19%n64`bnSHqknsLI{4KFw)6Wg(Hh zKcz+(pwanR4^9V;y}qp;8@e9;$oM_m6`KgH3cvl?3{o{kRbs`yiJsiD_SE4))&6PV zy?W&I=&40Kwvz}KF~P3gFZ+_6exlVyS3i;d6zF94!b`i$XX#HSdij(n7DylPtj`(w zo0!Lv4!;2{({ZzgU(B{TW8Wo8bQ4(*dPud(My191FKD~zFd*k0F0Rc|{?pFdzc}ZH zTb^9r8-zyFv_rnOfDG}=e_GkuWG7}8&4jT!PV?qn+>|7Ei+mUFYj}~KLwTbftpC+WZdK{Yu#kYWd{L;*!CF< zw8Fy!F8D%e7Zt(gt*&S=yt2AJ%Q`0$tFm895G$~9sZhpppI*`~5PsfLf}}bGSS&a@ z-rDYSBUMY87qGXbZ}ca9i|ml?O@lXQ0nQtV;BW1XKKI6fP)W-3&JRe(tn#q6R$)i) zHJuoWFg%!O8m^NM`I=LDSmVnxS&0^D+64BTiV{6Zsh)4@-{CSVxFb&%a@X{Oxy+%7 zN=KP<0&XuGW}Y*Q#|)Wc@W+kKV6+vMDQCgIh>sorOpZ>ZN?cd)AC01_OU-lO{g_KJ z=xR|AK1AsLwH%K(7uPXsE|=ELYwm-a#x^Z32*lW%kh3TwW*OuN13aKQ8ZNb}7zTxJ zTHEU`bgswaTq$VE&U|gU+9eeG&wZXb=2HhOlfDw=zSY!HMt!|4cYN9*u}{0uspYq| zRqIvL{qcCR;>lO8-J>cl#`gJE_kiUs$kVq(zK{KG>-N@+z54qSN!-0B|3A^c8sTB> zU#*CwECf&E6dShErF+yTFqqc|fSeixCq&Nf-;c0?Id^ZNG~lC5xQ5bZK<@|qBpJd< zN`E{r4_4|;{)e|uIlXFW18K%1xs_DIvtOwkbEoE%TTY~cAIo82Px8Qf<60CauV&#|(q01tZ zj$zG$c);677drXZ-L2(Da%L|-7!8R_LDQ&hXq@an>48o7VIBHAYAI>!YRPp(EB*Pl_JkW60^ zJK~&9QAEYB?{Li$xNqYqXurDa>oGqW{084bjA>FO+s=mO`C%sU=cylE{iF`mklrGH z@AWwrZCYK=TMMjpo;jqpV-uN9}{t>anu&pBujX4+tKV+k@TI)d;z|I53^ymVs_`1h0 zz2|Y-NYvUgu;FidU(_N0AB~JbRRe7563!E&Qf{ah#p;>R%-=0Ac=4apHp_|Y)&fN_ zS2WKfDp;fr4)ceL{)2BI{k5GUv>!JDIe%bmnePsqGZcOCtQkvi*0K|k6{4?JxqVl- z3}&N6zL|gegS7&&cqYZ`Z7BrMks1VE2;@{Xs;2%DF`TCA90M^>l%0!YCm6=V*ss?!85?SH*K4VkdK5oXQA{S%AJ_{8n-Y>~ zOY5Uo*7i`2ezcNvc(ov3jQR`TOV6z!9@xdh63s7pAZQ%(PhU-2mvVK2PKDF1OF_11 zu$lqZfXIPib4dU7DeIkc#F&2m)!Wye;ItX2-+eVtvJB#I`8qK)Wj2tqP`u}!XH`X_ zzpLl=y6`0#>!isGnx&|B>_>V$c;4C>W)=3q-e30LuXtLt=J?4RLh30^GlrZ8(Mm}^ zYO@zK#q$10Ij=L^4$x_6c)-F2Tb&}%1A+@gt@Qmn>T(|d1@!7kB!hV+9+`)Be(Ia5 zO|eX`1K?R6J(%S_WU_vhUcrDQpYQwZMcI4LFqh0z9OTnP($OR5-j?pe9Q*Z!1TxXw zc)6Sz+tA06`>54M_?w6wQ~*^rwc3Wov7{UnK7Gn4x9(p|`U2v_NwsBHMxSDFdjgd= z`_B1B+3h2zQ}dWW@ybt>&WF?n$ted-Gbo4Xaguno{g8xoj~aY*2G3^;i0QGXH=ip^ zw*lj5IO6ZFyzlh8H9`BGYU+8lQP>f+Z&J&xFNc4Ud!~F=A?a5)K|6by&IN5z$me~i zyu1Hc`dTW87lzN)LPeJkL2W|Hn7Csq^Qv$tH@tAXQ;_oBkd+H)bw4pYgpM;eS!T$_&Xz)e1P!O z9%4BkUcXIgY_V`drKs^=KQ(lhg}{_be;t@=#)Dtqf8d*HxYfwvTypFm4&91K`6~W( zLi7d?_{E6);#R_80d}7zF=q13<9oDvh?tK&7r^Q-LIrG7I{)ncpd>D;LXVbbH<8=J z=Xty_p8Blsyf$I0dsdC1zg`U~4zolnuMte>ZwE$bS^sQ~xBEgbDyRB=<<;F3@C8s{ z-_eR-ld3HMCgF9(;^+#TPbiZ8UY=uUGW+1I6lDMCsa)_xm5CnI6bztjUVVpI->?PA z17>6D^YqTu8}DiCG6!6^+-wok`m*<8DO`xnPp{oT*q9k*M#EUJ6}-S0a?Dl3ctQ|hAfQ{uBB)ySRos z0>J_QVb_}9+A;_kNs}anNpwqsL2O<#pQ`ls4Ua=MX#Fjk#q+EXY4d1r|JTImk;aE9 zg550pCi0Sh*G&nNGO2XHmrfmRmuZ;s&6MxIs#VqdEUHtHg9r+yzjS;06On?uXk`AM zm*#ZBCSOGrR-N5k0aQg(+JE$|=o4KX3oFY-3!KgZX#V2g0Kn{vaj1ccF2ToPN-BNs zQ81aF0RI<(YGdqs2e_L1;NCf8U(67~(^qrH(OKTbm#p6l)FEI(U`^e`QWDnm+Mhh&|$dWK~NN@ z1>LPNh}*e*HW|VAcAloHQoE2D+jCNN<&hMO>mDj+{nAk^8x8^fu=ftvC!QPD+hWEL z&EZ|$uTCT*(Kc`%gn9qJBT(R2$3}fy@qqnII&}6p`QEV+x(xr_QpG%9xOPLPAoAo<-Tkd( ze8EgmyPX#Z1q$+;u+G0Vfd=HzoKswI%!R?MEvtOx~NOyUQ7b-h^I(FCDkyCbw6kuncykvylQ${ zoW6R>?6A4({==S;iP@sgD9Y*IhzcL;?=+&lq3}LCke3b5$Ix68$Qv*Ryiot|?kcf3 z&N5+<*AVbkC})jH2(vr!pah8h@7`Y3k%B*Z;7^3X+ggFDQI18$|MzH=o8hIvXUYXN zP||$|dH5e<2!ie@s0t;z|EG_04NEFp<9Lx|nwL^XM~`XKDQ{`nRWT99FiTUh%<|UV zv8GH@O1$nN%}ZLDkXR~GrdD>+8J$c)ImwhcX+U05liQ>NIw}~5S7a7c&pBVt;oD|y z)?RD#?)B{TzW@IZ$8y*fYOLw9!@@XJ61F$Fw}|st^B5ctE@EHfr}OJt^LoK*(F($# zC8IaRsEaZcS-%zRUJeCC^@QPoWt_d*L~O1J%U*2hr5KLik=>p^3la5sTlM{=;1(~J z6L1~i45w&DczU)63ZX#(N76i_Vpol1Xol(c*!Ud{ymOVj%Fu*xzG+X+$T4tJ{U3M( zM+1<$q>pON=E9nP>NnTOeZa16^sueZaZ8NbvVk>UzXF!y3cpjHscUO7piZoI(=R7% zI=2ovR;cAGUsCI;dHjFTg$U0wU|t7Z%3r@R)|hANO5f0 z0TdGlx0UzVTWAb8A&Z!~SD?NOvZo;7@(!O}M-2kH|E?moIzX^%qv4?kszJIoN3PC_ zmcr*onaMz0>v7F7ay=O|20a0dA-QKG5` z@=B|+25o(Tzj$VPAswXGw806!>V&F;SLh0`4XLcwa7L0foDn$rYa=FH2hd8M!cytV zKS%i1T3;GyG`&Q7pj&Fvl?xT#Sx^}TGs@h%36TwDWO5@gty%VfUs466$+1m1=yFVw zu{;xlxO|5fn7*iG+!vfI@w$vyt6rR?BQ6m4%{DiKbiY;=+oTYFH4> zHw{ow8k8*fi!p9VI`ImTG9_#tSCAqKn4rTbg?(g>>8Dv`3?+AX!#6a3^N9+q#iRTk`2`zY8@S zn7VmIn-Uhff;9oH@fZrZ??q2~jOH=%64im2Ef(K7ebnr2dUy4&cZBBQE~NhEQF2D50t?DJ2C#6FJesGZB{E8=|--lXA%2 zy~D^Uaud5o6&(K1^G zx3|ljKE}0yA$N4Ot6?NYEb3pPmiGH!gq-ks%^K%=#4d7UWaeuXO=5%MOlu?_50hZY`M~8 ztJ%D`nUa80FlV~3-NVNXC{&PaXWB0!Eynigr{ZZO!mS!Bj@CKqVl4!lfS6=R+ zLkp8OwQL|=EWUVp9cAo&ldP+|=k>{Z=zWJS+?EZ@ymv{>op_;xtJ%vKwZp;)unmVA zk*h5!aUR)XZCkc?@#f}Tansph1e*`C*8tiS4t^P7rO!gd_b* zi5&q=C}&Ym7 z$;`W>{tYmtVv!tnu(ZQ9MDQ%fg@;&7J3s2k_7=Tt+1W3i990GrDO*&zXxlIgEk8m5 zFh)J<9CnumlguQkYHi=*kgTqZrI$oM z9NRC*Q+EQn?Isw3S}y3MRhi<6r?>DT{Qak4$Ox=;r)$LV&FpfP|5lGpmz)3T+#HN@ zk-m1x!F^d>W2x)c z=#qj&LyXD#4%VWynxXprPe=B9{~7~~>yhKm*GI&S%#%Pbt@rJZ6Ej|&#&1tNxP{gZ zv&h*0`M$8qtEX-(Ol;gOL1?)h6$D2~I*UK8lyju|BIkF9s2`n5@1dcb6}qHu5w*IS zs7T*a-1!%Ou?$KQfEg~8n2XQJqQ@^A*SV!8Xr^rQdCF6WK_dOGX#1rzE&z=-aJF40 zX5$IB&{@-HE|HXZGv%(k{%A)Cwaclr@xsNqVS8PY9u>WBnW8r_c4z;glq$~DBXrct zaM8D=zd@IK#r=_|$Oa7e(ZnUdJYDCcADlcsxE1b02VuFi!DDVSdA{iIIbPb7x7l&vyGgYZoH3ZnJ9>T2W@>_LO9C9|8ftVW35Z77rcng{}abu$Vv!EuuCON z$#4Ug9jZ?(S^ysvQF4rtBd;hMe5%UjU?I-Tz0QCV=3b)j%$veZ2PH8o<)@m2%~0aB z$#kTT2i#{l7zTVk$IGm`89%)pfa5kBgx=M`HT`>1>$0?Bx<5gveS;ny{qT8v&3imA z{%YLn;fuJDXToDx>TiZ{pHS!CzBk9xB2{`7kQ?_$uqNGgj7_!ZLR~lgfVzG*5Vdk- zA-{52xpw}9O#sfyQBK)Ns u{#9F`?EhttdW9fz!d;20U?yfms-*`CdlLo2>wK4h&pt1I%u~v5YJ+b|U+doiK&L7_yE%Yf2auB9irVP$5fW%M!+tT@oRZ$P&U> zi$aPFLPGZO9QuBK&-2IghZoMhoO{2o_w~Nsr?KZavGe$2;p68CB6!&2FN9bvLU%GC zoBeeC*s)`DxAbjqA3MhP@bKq&!tAQ+v10=41qHpBW^DnG}X+3@yTi0^; z>0>TxWn;HofxBo6b zcWMe6X>4q~zB%@b_Dtlps(HR+XCg1%Pg-#lO#1D32Ku1`3J7S;$&cHAm#bA&0>V@s zZvKUxi6rYrqk(;nhkyi5JL-%j?fPth>0|+#?#u&x2-A8}W2WZz*sqaCDo4rUIBnXl zNSH|EW22)#3(C+5Z{LJ2y<18+0a8|j>tv_4|Im2;5Q=10Gl1{>>71+Dsocx;jErO! z)CG#rm%}zuU*_pOg@hxSB@k2#^)x$ofraheiRGdFXz~qD9q9QPNqEq1|IaU-x*C)9 z*eu#JnB$W0ox8WDDZjs)Yts(gHvk!*ZOd0q)15IG5GPGDvpgiwp1EEnK+>Ew64j2v z9<>ZphE}_Up*2gug3B}n)x$r51Orgxv_q|Cp%GNnOoFukzf`-g-|YgO>K2K_FWZ>#OubPivZzySaa7!TC*4gm`(Zl^+rVJSm?^H{{FRS zVkGl3GBbJE=JFf=S%6Oiw_UeOBk4`~^}k_=gB#GE=>efJOZwx;|G&97nh&mmM`-V$ zn@SNpQlMRzIBU!UE_hM~4xw}`$;{C;Ne<4DbjuO@-)~^ZfGE2R?5_`%ZLWCMc?0CF z_D2MVm#6``Tzs0DTm{t*KA464X!ncy<^9=T9GY0)t_95d$h<4qmYP%KGQBU=XnxPB zty`OWEq~(@3Y)3IKaS&rFJ(EW;gmj#r-{P$6Ey6 zw?fwqj|siL=AU;(dq6JV*^WUli+6jWW&Zm9Sd)qd_r3zdqsOT?z2^B&*a{%UZI5R% zFc@6N855!Alb?{fv3)-SLn8}qwN=ChKcV&kym0PWGbIDr=}yirkup}}vPiFUYSets zuN*x`iZfe_agMjOIz2Rxi;rlWd;C2Qa;M@A<{l3*Hm9W*xJG-Cwv8gJF zhP$wl(e%-CaUR@ldu)3F)35l}6tO$z9H55U#Lk6;GkpFi#v9VVgr#hHu$7@DftUtd zM!4jaRpO^-;pod>H634TPIygi<2Q3wypwD=wE&tS!@rNr)H%CcgMgK*XAiKA^nMLK z2xOdMKKH(dNIoCc&aX+o$-TW4d91I86R{~rM*hnMcue8heCZqS^|9w^{WO9pW9r#p ztx$~_8)aQR#8}Z9Cf$u6T98W|V1zY&U47%w)A{!CwWFSXDpNz$Xozrm+LvJ!#+?9S z#MoLw{P(9*$Sh~!CqHQjC>?OLNFYB>b1CjxO%v4S=C8j)?s#j)&2I-2C;pl!UCN3Jn zI~)<6*Yzg#@|+fL1uTn*Cbm+oGp!qKrTK8Paw3rWXVX8y$}M94oRh$v89XTrw^MvZ zi0Fu3cNWJ95Ux2%HO8m0j1NoVJfHmjjx33I-Ykiph#G06_~UjuzYW%Tdj!&wb=s}L z-7=PdmcR-0u;ML{h?-V`#1+++7{Sl_2eAh7EaGrP!J+^Q+zx)lfQ_zZ6tmO8NtX?W5aa0tuO5B)sGFsB&C{#US+9tqeHFXX$kqoOROV^opi(1fAcz~!I4j+bTPE(HTv z)3dj^Sv9Bw!?@{dh-WNBBp<^3<=4atm>({h^2nV0OfZW}9F!wyyc3W?3{TljHY}Xq zKZQWH=9*wZ}a^LG4}S+ctj*E2L|$tpgcBTW^9PlV~-?i0Uo+{B>vs&GC470amhB2sS3BPy+l zKdA|is~FuzYZdT6OOF~@uu#rJj<$#hBnH8+-{fXh8lMhmbnC6fm2wJYWy9d%D=QBZ zJ9ewvRR~D>;bBbBudb{5L{AfgXPoov(EV^1$jJ_FMZD3)9zTBLZdi~r2er2{r$(d_ z2ebI-a62Hby3PXa6Eej~FFQEnTFnJFIE>*}+i&SzSOv~icU^N`5K8_=aq7>BQ?TQ` z4S>|B=sG8fVzUH124F5|SW)%zj-81CwBK77sNP2kGMg6I+&?-k%1z{q#4}oB&jGxS3Q#AJw7cgy zMUlp~So=Qj|9c=sHzK#Bhc&}y_PtQ?B}m(5ggX@5( z-=FmBfd|%tqz5PHF8b6ie#@8c4J&EdfYyT?lKd1+z!8mCw$fRPcNEPS1;Jdjdt`Sc z?EYpR&VH0?d|1x**5Uc|>RHe4$E>+}WGPOaI=yjMg0IcZCk(?9S4Fm@>8a#|l04D`&UE*hI4)eSq zBD$z8NQ?ONiCNb4Y5K=)&aN-eSu`|c{shf{qR3HZ@%nbx*a5Fk8$wI)bD0T8%A)c!ox7}(-c;ZF#~$10fDcR*h64H= zqS*nM^#iD;=cozG=cAF!-wg2vmEl{r_p*9xVJTR#R)=VF;<0jGzuO~YG=bh|+8Xkj zfEh7vcYB9r30` z^8WgjRL}~BBML*fE##9Jdc)M3ru2vgh4NNoYx9^Ogc85dL>5ffLWA)3>AMFz^p>5N zoo?8H#0Ae|!2b-wih6=CwC8f|?|2G92rj={lYpdEVmrDE+tQ7)9ySFWQ!t{loFimA z+Yk+g3_+NWh{YUD5g}HRgM63y;%;3n6DKbw<0kwy!!aZ|IXq%U|EF92AhEmM58H_D z97=ohDxduAOJ;}C7{ZU|m^g_}G!?RP za$A{3uhc%EoIOeXH&-JH`PurFKuR_!%3os8)s^?1OPq!veAbC5qN5gZjd4F+7RM+! z{pD;&V8LF5{&_3bcT$Xk5D1B!rB?#fa?pTbc>0@2Nt}sGcc(NC!@?a^`%_8Oq=4ha zC6GEkG#Aa9@AKw@E%viBtx19GJO2V0>pbfh>Fa_9Bhl+P%hNk*H`vK9#3EzrjEK<# zv_GV;<&KQecA6cdiaFnKJlSm+vNUX@NjFSs*}db@w;LjcbM9%dcNM2rvsO%mX6FEZ z(4d~mz0ue{c@}a5JZPDY@IA%)uX9v&2i#T>tgTN)RuwLrfbmdrH#|tt`b}%=TP}4r zmPrnrqft zn>gXy(QE9ue{7!%9|`Lv4haC}f4IxK@$}-P6P&KY?Bds-U#aIpD4oC%UsH^1-5t3+ z(N$>8{j7wk{Hr$!hEehzm+w7nQ_g)c5httx<XmlMih z08$SKxN*)0XjVQ3s^tXa_CG)4^{nwi6)FX14(J+QQeWHjMqgO2N zjxQK=VLnvogwSiuR&`$G8_K$-c5$7S1W9gQ4^R2%eNz~Tr$P)fwdHKHaoPv`3_u|D ziqJ_+B=P)rzGi3sIv3#$?#4mxs- z2+MccbKq6+lH5TAvRwUAn>1-)2fBSBt0UO?(y_Z$oLvy^ohc+fK=jzw?D&JL{CCot zmV)YNkm4EFVY5o2T&pLW_IgbZwwclg?#d_&wJ~OkZ^-q6aJ?Sa}4pkbgH5?fQ`T zFnVL?@h~^0)kQe}O(e}=b#u`MH3Au%g1%wwf=H_765^T_V0rv}U%f;HI{C_}0J|H4 zp^wX79^JkAQw)=$LNoJ$1KeXG0QW23K@YOP}KrpnBDxv z!*&gZ0H1!EuF;Lqz!B6!q%NH$RYWD95cv0WN|LOna4u4u``;1iqRra)h7dW2M%}G~d zubtGN5i49{ij^NqL-w(ZS#jM5-ULd^vdfm6C>?^4f%PtY{AMBOHc3}(3g5!Hc_DfjhbAfzN!Y_6cZ0f1c)G6Y4&VNeUW>kiP+OFz}d;@cy2;2XvGC=v|n1w8^+<0mIUEBaqB?V@%7h|`SM2`>ZRF*V9bdjM~daXEsbmOfN-CtBn>7F_~=GitzOnDKupV)W7wR9aTHq@fwEJ@5Gc7g zF#O@5)u2YSwQZBhk_wJ!t3eNa_h1bjS3l)ZhQyl|q_3PVKIznYqUtPZGQc%USk&^f zBz*jo&q6SbWMqUs`XrO!22&R(N3=W>|D^NR#-IVR?d7m9y*F^XPT0v$U4HW?Mbqd?7;Ael$E_(0hC9- znKrraLmW#jeEdN&Wq785v*whP-rG$jQ8lbW3=8?nSMrZu4w(ZzC=U{!TyuH^%2Iul z1rr2iL9ml43;We79Di-p+Hu5F=9TkD%pK>jI_bU0-l}YC&P_y{zKm_wsl+a;C#$Uz&03tuv`PqlV$Roz zknM3>+9&O=$O!PGu(62Yv?shC!YaUoNmiD^o6=Wg}u+X1&#d32%18kF9S zZ8VX<$;H8>v+GI)^YfZM(AB)Ryzn_-`T+0;rYV1>!#Gr#ESa)V&{4yQPy+2>aN~=H zuq`!3D(*ufGgr9#eHWc4JwfM zHz7j^2S@ED9FcDL6yC?{r8pb8XX-8`?BQnINhRL(^;1w=g)|j!iDUY8WwYyxEN@ELm~WfVGO)tA|rv%O3+fs7Y&g? zwhNdY0}Ym80}yS`UE@Fucza;MaX|Idq_UL9mIXhgC5#w87R^;Lo6ikfav}MGULiOl z;cDjkpj40~L;VbGQ6)lpz&bBMa+@hqlgQr)vD71~dquFDy6mh6 zXoMxim_0#DGYAGw3$UayCX9VFs9^yScdkgoJi(V<@YiV(&bE@vp}vAINS|-e(v<1QE^p7@=gPd9jC{y%3!_ zq=_qoZA>!V;5I+U48Gl)SfOs-`&m@4C&53EXeb2(va$;LQh!1tL;<;3ia!t~VR0vz z$ztA1)xT$%Oq#>~ivE@t)Ek~r_?ugzOHBGwu-2{?M`v_^mZr13isJ*Z5**~tqO)JO zb2_JlV4*%h@VBz&uLpk@_JhAHec|s2H+obhIE@$)A`Pe&xdfUN(Buo?%oseNOFZaj zQX`PrJJ)tX)Vp! zV-yU)(|7Ogoc|5%ExWCaQ|?EhbG(|v_GbtF#xVu-^WU7>66d{U7dt^xi0gW#3b$ip z0}^W5x~zB8_R@eU^s4c+0zZ1hSS+5-GSbR$-I zE${OsNXO7%b{Hf5di(zWDBg!Y`FdYgc%671M;i}2$5UFm76bf9a&jjsB3Ltm>gtaP z(Y2?D3s>f1W3M4^(KTx>UHG#~3Bz9^&dC=vp=90j`m+Uzoi4mYJPAwDdPC;m*M}r?2~WDzBWkMz2$1b&fY7J1SNWvcl) z-C~1NsN2pT`MuBOY2x?%c02SyMbgE`9`bkIXQ02tRbA5tQx*H(eTR=XE#MuFH(E+w%#ouqJ;|y=f_yvjv%1 z8DcW`zbO@(7t#by!#pn!#PPdtF0_!+=eBd6(!_{4=u4dSq(7J)e^UfJr&?Ep&brW8m(#Jn*Yy6=l(+KEE7vfN zw<HQRcYz3+2)3xK8yu@)TrY1}-s1UaK*RPA@;R8ccA}~jafVoaGd|}O9 zHNUhQgl+u($2hBH`Dy?MqdDFW`!z?bia6Bw9ffg@d>!|*BRg>Hn3JYtH-3;2@9nR|z3;|oY2U|Z^3B}vScQ2qKwK~Y@cSJr(v*#; zsyNY;4`@Ft&dgc%(5pbZe^K@&=R(+OHxai7IRg$Trx%XrZwApw0tLa^?ZvFcb9UB? z0)rtvSVnadVl(B@{*!N-yWUe&2WDGY=?-^RR&p6t!}|5H!D5#28eoGNDz2+)3d#&Y zOw)ubJHL8+=UM=^Nh6uY%}p(ub%%Er#FR^s=*HgeTTAr>odg#0``({x2S1i`_I}MC zp!QQtgt5M3P0fI}X#Lce12^IPCRRzkF9L~7SrBEX2`%GSyT5w(UDCq*H`m=!K2W1K z7oH;geuQ`*WD;WkniTW8LPj_Fqd0u$)G^n3APvA-3Iv0D*@jrT$Y&AT5IM!foAq&9 zJ_O)&TGYzP(k3I*?yPSX_c&lunlHPxrjEulORM0C*AT} zqF?`u^sC+;oH!3q+P?sorQ)w^H%DYk+;P}p5PWtHU%=DOQu@v)oEeD;Y)M43{uz~R zf3ecKa&UU@RK4b||7Mc6BEQ`76~P?JjA32l@7J;bj*}4KN{0ldWSW?G=CSmr@UxXB z1!pUbh~V!bs_P+wp{_}96?p+9KI$2*Uwzy3L1_blA7NSRzs(I7iJ9th-dbcK7d7*Q zhduw3M@{%GOt-$(uRv9LgQNv8bDui=7?CZt_1v+>%Q8eaH zt6-u0$fO?+Xf6IU4Xw;OA&$gHChL4&q_(K3pMzOBrHZl36W%VT@&I%-*l7&Bbei08 zEHq_JMN@rs)7x`7Ex; z(n-UHKVott_0>~>7;1tzjw$4Hu{kwmn2&tlOgl43bA$uSQtH9Y3SR($K^?L8Yo;Th z6PIYseERJ_B@Kkg>lyl{U_;>cV;< zc)+g9b!p-3NGH*GAFOS%j3*frc`?_u{$b)STGH?x>22yZBIRY~s-GeU}MITZmGpbgLcj%@E4_K>W3g|K)C} zd&LQ2JJF>>i=!~GsC zG?0@2w3%W!dBWJ`5(_iXo`I7gt#y?72W(AB$P5e&D9Pa9uhfo>&P-VcNxNrGw4rw%Hxkc1Fv0Z}a`ify-YIo?f!?%P$!Y3pj!wpuXey3gps{-|((a)1 zT_Ew^mOSqp_shi7K_l`Dr^(YSl&wwGZTC0}i4FH1p=C?@Pd##BCP4vPs*+ygfarQ+~cmaEoj*(XUd=dEnwc%dSBfE{m1T9kl0{C zF1qF2Tl#~_UHZ)@Oq-4G0*Rd8;e^?K0{Zx0pzPA?sK|!!_cU(?St>r1QWF^WzEr4O zrYsm<{+^@UAs`XW3&g|iLKLjms`($|Y;Li^oZO>Fg_TJfl(Pq!`ppmTfA8$pdT3o_ zWE7qkOLO#!BhBN$;E&r=;W<7gRlcP2v`QtpH*O=v;ibGn zAfQbV9ekpBul)yPJ9^Tt|K^4pfRdg1#P%n@H8m`+P7L9MURz3%7NwhSxG2Nm=K>$V zbgNB3Ysr3fhzP&j>@VICQeW)sIFDYK@_*6Iwe zd`X$6E(RVR>xvpUe1yAdoCIBqNd;g=WRoS_8cGY>ZTeTwV)XiJPp@-Wh3@Ee34Nsj zN8ktS&z4K^2IEjETY@mpntUqD0h+)!GLHV9bf0YrAH4D&E4 zHuCbU>05oU^ggTtbAXm0eYU>Bz}SLV_v23nMlVSr@gh4k1|as_qZ{6|=7-3`IK`fb zD7AP1QoOoz*ja+KwTO8-UzGDRO<-0U1j+~g<^WFg&K{G2kS~#l(!E4p>D)1WbK&xE36UnwD$oqpZ##~LF=zS{=UmPErmF z*aB}8ATdW3`XFE;ss9Al_agWc4aV1A8aJVW52Z!Y^UNFm_7)s@2MJN2Q=X4U0(Xhn z3%5InSkOO+ULSOgTNgnfKS1vOgFyx7__j!xj`0W=#Hk&*1ov3y2}=wYbM!CEb2qCF zvLLK*S1ll3@JNUtxo0fw90BMgo~YoCtD7%~4ChN8W&W?|HUH5i(-N?DenJs4OBtYL zCeKoA4o;S5l`#bAz0|W&I&DNX-|6?e&5$27@f&x^=W?bXR%Rc@6bW+*QGzU3R2N!t#a@!HmjtLt| zGlx=Y=xS8wp@h~2-i_RbiINe4!JAo1rhb|FIVT}4%0=R=Z4nw2fHy>W7_GfpK-BOT zue9Fg9^yhEAErbwvK41Uuj~EUSF=liZ49>NEOz605teMOBd;YiGp;sRiNP0&+6k9! z0Bk;vSXwWk#v;8xR)tFXOE|n_P=)n9W<<0>ysicCckiz0Sm(x_O*{V=ld(bKOTI|) zc){c6Hrzig=0}D)cm`Bj)6H{vS|wVj5>7vi^Rg^_2S?QGZe7WMXHJ;`W+eWem&(^9$ak&hFWI}QO+yeT zaYAF47rkyA?P}Q-|FNukiik~nX5#W)z#Q_NeYom1ri(Ieta=W3Jmu?;YCU^LFtqL_ zwAi3h5H1rf@!`Mzsc2M+WOqwY9THG{(56CF4_~rnV8eoWw))WvAPoNU@xoKP?tS|t z_iU(b>__fX@N5A}EEREJr}FU?qo9?`8Ie+&!_uX!5+sGO3V&`3_5U(HLq&Z0brFC` z4`($*uTRE7d(BW{V04QKk{Y?lkEs1|+Da`ZX5C5HeI5sRLspiqJ$Ansp>*q;1dbP5 zHe|+=|D=%H=&yUNR*Hzy3zY`gM*i9S<4SOD)NiCX;~dfFiT>4-C*q{M8Zg}qBsiT? z&G}Lo15PdfUARs1z&Ef6?PHV#*xsdiAsB%ciSPN!fUxBASByD(R~Irjc;qj-&GWsB zaX2BJ>TnE7u~I}-2p^^RC^o(-9wFWW8+gP$GaxUcV*Y2SN;5>wDGhKzQQVIBmK|X~ zDgybx?%S55)uM)JE72odNI}w~0Y$X?j#D~XEl%QtFC6mo1D-F$NqarN7O#)|wa7%P z5j5Ju;n0C>{QIQfJUf2b^{?}%Q6ZeHHSCDe^V6Rm2rzYFQIWJ~HqWh!UVkaqw~Gbs z=9iy#L_UD=;ez&>OfR4C;gqHv`uWrf6|Nsg ztZTD@;dRgTi1F@|Eg5f=TuT~@_&JALj9SdG#%R!qh#q>k8QKHuRaJl>{%vBM>3~Hp zTDIJVn7-6v7_4~hAUqdxs)eh;$5C1(I3i#Z)APC*bf6u`-UD`o8B+w#m{S1kA)?0T=^bGO=Uw&5H)D% z-$8FOp;Sn)meKtBHF;75_!`R{f0mQ>Vh!6sDY7_r8nZ*$n7C6Z~~}O3-epH;Apu^ov&BS`FVKRrGK6 z^W$jU=oli}2Rf`JORH1SQXDUi1eoX%&l`okEJ2ew4BP zPhrg{MdPUJ&&^hBp%&rUH_6!3JA%zXN8E_BzJ)?fGlzr`~zRW3%~I&-KFplKh-&|iVT>Xk~FZSq9NRZ8gbMyZuo-p@2L2$j{N~ur-X-D z-`qd_)S_UXog*d2v!vW&u&maUzzG|aTHtO~P!kYR2MTc)WpHymt8y18CmBvKOv$eC zUn1rfCMOr)>v}S~180GZh(wF5AK%{r$p&vyutRAE6}Qs$FALzVP_6x7J{h+2x3_(* z)seo#0rNsS2{r=hq#8^xf$W&BSjWL6W48v{bnXQ*Xx9>O{UO{ekEZ<6Hm@3fs)p zf6f&Xyb?rjS)UXUiJF+<*Usqm>5p}6*o(Vu|NY#(hgmRB1RN0wOw`X_;)iuY7f#L*&YhPjzEH~ zEB**VOnZp-tL>pM?5m>(2B%`~AIh1NHmDQ;{fPxkJej*oyzxA@+o2Y*;JUx|9Pt%DI{7%@)SAA^>oDx*jRzY_^MA8VMq&+9b&vj88fX6^%AZSonp3@&0%@*14 z(4n4SBer!yG+j5P7!gYa?HPiw#pnJxZQ9Dz!g#wf z%zBi6TOYjwx-n@X5mcgKlAgc=!4uQat*Jd$)QKSm$5M?u5f4^XzAPsDrk3yY^YdnH`|bL78I1eROX+?SFt@(owGn+F)401KJRN{MghNHV zP|5r6E!aG>CUsnnwAX!Ku9<$TFJhvfCp$c|Q46=B z{c_w1b7UHg(8IT%#P_ryRekg7qs`hc%@nzpFJ3yZBarOxez9#{eg6D=XcdO2UIfT* zbflrJbboey$XyX8acpkRWkV$Q57&Lwv+WLsRcn96Dg&P+b5lW?i+^_GB<0Hd-_PxNzsZ-P0^-x>polUdOl==~PKr zXprU4q#=Qjn`bXhjO|tJxXISN{ti6#BW>lytAv0J>y=BfA+rDNV~ys-+f_4sVCZW{Y%mSG`C9b);^d-Kqh-^&F`zW;?7Y_|sg)DW zGGR^kMv=rO7c9fp?Gh2Yr)xF|wU^Ild1bqg=rte8sw;}*&4@w2Lpe@MNMeGEViz{& z%nc#HoIt&Ko}aHh@Kgzt(CeF@rX+BloUD*GE*~#V5_`Xz=f3%OC@s60gG=XSV$Xy9 zESS{;WpDwqLVKILUjypZn%g11L7WydW%>{k`x4c&2kF6U8_&b-!UgVqWlP1>o z-)sF7RlBY&d|WZ^Rm=#?yE*%ThX?cW9s#|L2@>cnQ*2b(pj8sGT%6!jw^sBp{kr9^ z1+goTRv~g^{>&tQ|JF~u*3nUc(epRBrVECf36pv`5jB#~{$?AIc*zDw zU=z?{lmhhPp%T+9D=B_TgQS9Ct*Sf>s^^c*_ba&CS2{ebi$JVO!`Fkq1+K0-$r z*MG6K#X{WNVb_5t&wps}k3oh7^J&r4FQ9PfbkI$h7Bp@Tg}=+~_C6k1Rzcgdd6W z?fxk1yrS9E#0~!X&v%ACk^ZZpn2r3=0}FZI^3HwfS>MKMu8ToTU=SSDu8$p|ubGXW zAYKs1DMKTGWE2dPOElYu-8XSbyovv<2roSbjow=iBaPT=!}rjDqUfL)|9r38_ht;A zmm356eOUyk(+o)(^SlI`AmaX&07-}r6a}oke!@O(x+~$N=Sq$YrM58O(?EdCXHzoo z>s=GWxP7&TD1I6CiEjfrRp=3ic^5<4YGgyy7OByp zEhq5-qmMc964c?2-Nxb@Doa;7SQiA{Fe|>x&m<1HK6LD6fZ3-)bXEbApevoyi%WhH zovM)A2SRvP-!xkgw}nkz66uxh4%zyu0uNR?6-(Q4tPf@4P(zB)_df<4bgVUaFWu#) zS)u1hij?|L%m~;W$*e2(LrrdRThyDLTQ3@oM6ulZ5vV>$?OdxfK7h+GqdA)0&y^w9NHPo(K_} zY%{T=8qPjLvC8l*jz&OYWyEmRyM_Mw|CK9TFjAlD)^^pt4G#+GtinWKoN(#?D^*^M z5>YYip%bg=P}gqxV+)%wih5=9h%E)r0}a(;iAg2H3o2I0uXnCTtd8UkPaZAEQH!j7xa8Nn*JEGK&aL!| z4PNH<*%`Wr=Fy#buYWo3P+FR#i>S-9Cf%x7`BYT*QPU+l9%`mBavC}px zVz?Zc2$=K*d(-Gkc6qM@-(C4vkuWK%6F#I;+3fnt=_RK7YGBW0jyGJDS>hzY#EgO? zqyTwvocVd0o(!;w?Av_CFVwQ`2igtGlE?z)5|E@v90v#>ix<`VeOs0#7Mfd^v#uXM z_>)xppqYxQ;r4zi%VPnEr{;i0$DdBg&y9-%Ul06UX=dyYr$f(?ta_LurQ*E~BC(QrGa--i>jdqB(rN z%|n3R`Mc4(#_716i4g9gZ|{j#AMLwubZ(TN0|{#SnH z7o$T@)2rU}MV$BaIsI`@-6sF7bb<{~=R0GYcNjxcZX(ZUe35Elb#C2l(e3T;f@nOg zBG@;@D>I^=s}04HL3tR;h6Y+J$Tm!7tG|No7dkLtPbw`v%KVP5YM&4q_XCJ6D44s7$qY!x$I@ z@6Z2V9c2q2f6bg|a35}>f6wvtJS`IcVY;itcstzGt_RxshiDXG8SBjP*Sqh>`i{l9 zs2N;&<9^%z#hsYfH($5WUi*iZk;R;2GW9GdjbL8o?{%|#=jMpgR4*e9nC!<|7_^|s z8ORoZCem+gnX8wBTn=jS$k>}~HpM^1(50NN+e_#q+x;0Xp8Mz|=T2nEU_H*A%3p>{+H}%X`+toek)pOuPF>|rA@{zn514*Ry9u7Y? zd!i|JUfGy@^Him#J&z>-(IAUvhcuIbL77DfTC5&*Ae^gEIpcJ|4)S91(3?jmMqUld zn09z2N$3VHzG?vhePwvW_B&qQ)1VW|r8CYjD%mq+VjBCWNd2$MT}w^*p2d_XJA}4k$%r6BznFi)4Sg)U?5{;Fm-|0QJMl#zSQe(h--~I- zOkZFQG5Vz}4&NY3NS2|t#oqkXT?3xzWZL0PY@GA|R$HGUF&i)u2R?+Ua}pN5Bl$6jp9@Ri7E}^XsnxByxVq(*}=yH znCJiaC_ifDOLZyg$yc$-8;7vVrLsP4*+82}dZgQ1JBm@dzAW>PRh3-)(03*4yoG-m z6o7&4Psue>;|qS=7T9F9%wB=%4~;Js50)w@rF+N^Bum#_z=}QfUkMKk)zLbhz&TpFYv`Cx#`AsOr!;B20vvwnS0R_!15ZBSX4^%jKeedC$VkEKl&rY#cDlN}RMggS%TKtHJT8#Gu}| z9IdD28v3Y=mQO_TC~)VovQ1I(7qRqQOMvKO)m3Y2Mx^Lon=Y*W$x;D>PTuwri9eRt z2=kc?`+>G?RUsK@b+i6+xS_=3)((08L^N1lhtf5VI;#Lw>noN}cMYsPcncdqSz zEAhF;phvWDjSu0z?%3_7AJXFTwA#Dl?x=JWeFWhx24WH_(IT~ofX(^fhdJFWlK$w= z+RrzIVSHy#eZJQQ47|7xL~FJM{s@wNW5OSMnZs3i&&|)`>jX!ez`4?#XkBZYDYw^Z^nr~iL~cfR2L50$m#cI zFAnfWSVAPqwGmhf^RGXQK;cJ#zas6F$u-}1v!0uhS3m$;nfwSYX?IHxY=yQYj>w2@ zlDkftzgo>6k{LB|bNxBHt?SNP_mn7f=b z+rn}lcl)?yC7dt$Yx?hF!&jwpZz7>XSlf}%P`Z$*<45djudsCeNZrh!vq zkH0Aag?BkuZy8$iofYTBd{+{cIJ{W=y1%L@+{ognKXvzcE!zFWt~S^$dYtx&r2R$% z5jrT$9GH&=ixvXh(>&>~Oejz3loXhSHj6@-MQe3v%6S7pT?U=>_?F7?7CbKgCF`d=cWLGlYz{LiNDj(aS2 zp|bzo?s5+!vkUUZbBK|n2Tz^<*MJs)m$lSBOHY8#^HyOm3luK~5l7)3Nw~-=aoj&4 z0TlD5bWt2*`=1XBimJC0M*cVp1uy>>(tkm4)#tzm-eTzG5pK0VkEG2Y z#d#U)#|0|xMB$%K6$GK?&n(Z-C!l)EO^HAldh>c5?{<5tIGp?52%Vd2g#H;^o{H}b zhkYIrf6Q7eEd|o@Yk><;G#eCOSN$@2{fsy+>*$7oBKPQ##U@AKMb3!K#wb`431djI ze0dmnM?MDW>}lcL!WDR1=tBUcGf(v^nGs5u@2;$6tS67YNT35!1Nu!q{fA+Dhu;!_ zVnI^{N!zDun3NWGEN*DeU81XT!k$v$d~i6G+|FlU{=c->&}&$9vgk0Ff-ZAKDzquL z=|gn%u^#5)8~k=O8EQ0o?%Iw}L2ImhQu)7vSLq`3`%J+O@Mj;)CY_{4F8dT?*)5LG zIA8{NpkN|~>U$nWNWAAan@YALhU-e=b!}J37Lkf*4yu4cuAudh_Y$&P{7b&Y$Ed*n zA&#v539E5lOt_g2al^zzEsE+eS55PIgUq2cXG>;$dnW8@$I-)A&HaX%tn*jC!R=-X z-_ad8CSX;hsRT|s+hL4$S=POAmj8Juq1nADT+_#gLkv?vTLJp+fDm1U9F-}^eHhSu zI01u^cF;?Y@B}S1t#+)}y}4OWAK3^{Ve3qi9mT;LdqXQp9?lMZ zxz&gsRhgzNS~FSq_;0yC6ZRro@A$Rn>-1Pcm;sUNlQ!=Eys z`p!Z2i@nIcZvumKN6pNz6e|dTFqdi^-M=L)UCKJt3=Qg2{e&YY?qli^x<4H7pySnw z|0k!ZG)aD^8s$e_bN;_a+BrJ!2xnhofH{t0f<8;rrIL?3LsB^iQC(vL zV4if4^VADVcp`Ravcqg41tS!BssuT=nzXzQLrx%y!{z#+ht&;RWuX)ciHqSc_oCP1 zAi*>gPjbl3t}?r8l|bBVk$mny3UnM%Q(z4R+_KIOKu8xwzdo`H(R_89v1eU{pLoLU zEXyH{f$bx_Y7{_KSP8}VPc3gR$i9EX)#W-war*1yefVB(6BpQ!4KbnCJET6C;tUB% zTiT(oUK|q21_X8(6y7;6Npjz0oXduIL`wcQ-$P$H#KGSGe@eOfza-B!E)+x2WJle~ z&FT5Ld|M{fsXQvEGnbkjP-|xgGtJEDI=)plO3-IyO7`L_3kwB&#ZGgrgDdill#uU-L*e5`s)i8Ciax z-vXo1zqIw0vXFnnG+qU29HL{3dR0Vy;CX1W9L$tF8UJcMt#V7IA7MVB%sRJ zJm2;AulH$-Pid833`~@jIyyAPK1sTV!0g?ezWNi+y+e~}oRP59%Ba7zpPpk3MJ@tT z;alO4`v7b8ec$#?9}<&ZqQaiP8_tbk?yg54TtDzo6D@9t?>eGaFU>==A%TACcAb znj}IeTxY^bXRrH5I!>t{jbD<|+>yyeA1|z^nlx=M01P+c$<|y`CZ$b@Extm{A^P-w z%M^zJ{fe;^{e==Es^votAPmy(2|~}Ec#}#cK70R_&~(DQf+u5nZ}DW^g%pI%vBlRD zD!lf@Dbm8N=0_C0`M;r&i*NN$iVyFJZE~1Cy;e{1m|?Wi-u&=c2Q%81TXpu06*bq3 zRO}JvJ!r>CyE(VK)ZmW()obHk%p=Kr6#2Ijs7a)#3=6Dt3iRj5HLZe5nVVf?%NdV^ zB}U}Ncea}tnxbV@hKZ6RNn$11;`ou$*?IR=N0QH4j`!deVzJY3qcF-V3M(1vQdunRiFTVa@6C=NM@+Y<{2jVao#}~8{uj)3`@~ghz^B|CcUxhlt@jhDMcG+4Q{2YHN9y@>fN>wj9pp9y! zId!N=5kjvu7+&YuEm8d#PDdS_#dyejFz)Y|lBwaknGqs)ie>Gl)2H?L;Yn9t;0Lz& zQ1Im1yZGwe27LdjT#6Q|#bUC1>bm;t!;}Il(GKKj2dJs;13G+U3zNU>ZddLhX3tA! zXOw5~P>Juf24wpN7rG!Xz$q4T{Vt~ZrO^ME1C6+gC@;-DnlrRUq%4zqGNgvB?iS)t zL6|~(e1mrW>gBx`lie+w3Foybml*k-Vl1S$+v!Jx;q=4uZ$r+l>n9&9x_G36ivQ0) z$TcxQtEb#*MUUrF)Sp}WtFr_K-Lnzt#MpAL)XCbUaw*>jPtLAwgm7uwUo;mjo)A5* zc5%<-{$!0|WTf-ioGHof=j!r;u`j$3Yh$_Wc> z434nDFtyL-)m0iDW}@MKHqcdSUW3{d8$vI49y>g9SW+_1dU6HuqbtP2{-|A}U((Bc zY_5&Ev=)|@qgo*mO{9VL>wyqQ%*xS;n=-Pm3M;gLU1jo{WObg@Dgak^M9)XLa0yEi z(?JTLxe4#?qvuUwAAVU6)N=Z{SuK!N8PW}~gN6#aa-|ZF=#J%i#T_NqRRaiuC$C>S zA4@Zs1J&odT~Pk%@F$k2T}NwzBU~Pijs-w-O+))O17nghUjAwVvSWW5v;Zjc{0|7* z6!z0C273ojuOEArM%EgQx0_ z>T(VsK+EK?cCJ%sWvSI6vY)QDus|qRM7Kr|0u`95uw{!1LHOdpMbFj!@&Uh0Wc)q+ zd@7xIIS6bcSGMyy9e#02N@RUPf{yaJEawp6`xT#p3$Rq|EUH=qKfTQu;l zgds*$%eLt<<^Dmvi(A{lFE@yrmqmLtWB_y|b64{seB_0&vLwDcfep3aYTul0o=(Ov zLqW2Fbfj}4{Ks0sO=8ky0xYTPCL#&N@L24u3twJ98M9eKv;S^%naq^wyrNtgNH>%4 z_{i;O9JFxPwBsT?MccQ4?^3Ca%Fb;7lgkSk;+k%c*Xk1?*9rx$wyK3hw*%k. ## Setup -```{r setup, message=FALSE} -library(tidyverse) # for data wrangling, pipes, and good dataviz -library(afex) # for mixed effect models -library(broom.mixed) # for getting tidy data tables from mixed models -library(faux) # for simulating correlated variables - -options(digits = 4, scipen = 10) -``` +We'll be using 4 packages in this tutorial. -## Simulation +```{r libs, message=FALSE} +library(tidyverse) # for data wrangling +library(faux) # for simulation +library(broom) # for tidy analysis results +library(afex) # for ANOVA -### Random Factors +set.seed(8675309) # Jenny, I've got your number +``` -First, set up the overall structure of your data by specifying the number of observations for each random factor. Here, we have a crossed design, so each subject responds to each stimulus. We'll set the numbers to small numbers as a demo first. +A seed makes randomness reproducible. Run the following code several times. Change the seed to your favourite integer. If the seed is the same, the random numbers after it will be the same, as long as the code is always executed in the same order. ```{r} -sub_n <- 2 # number of subjects in this simulation -stim_n <- 2 # number of stimuli in this simulation - -dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) - -dat +set.seed(0) +rnorm(1) ``` -### Fixed Factors - -Next, add the fixed factors. Specify if they vary between one of the random factors and specify the names of the levels. +## Normal -Each subject is in only one condition, so the code below assigns half `easy` and half `hard`. You can change the proportion of subjects assigned each level with the `.prob` argument. - -Stimuli are seen in both `congruent` and `incongruent` versions, so this will double the number of rows in our resulting data set. +Let's start with a normal distribution using the base R function `rnorm()`, which returns `n` values from a normal distribution with a mean of 0 and a standard deviation of 1. ```{r} -sub_n <- 2 # number of subjects in this simulation -stim_n <- 2 # number of stimuli in this simulation - -dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) |> - add_between(.by = "sub", condition = c("easy","hard")) |> - add_within(version = c("congruent", "incongruent")) - -dat +rnorm(n = 10) ``` - -### Contrast Coding - -To be able to calculate the dependent variable, you need to recode categorical variables into numbers. Use the helper function `add_contrast()` for this. The code below creates anova-coded versions of `condition` and `version`. Luckily for us, the factor levels default to a sensible order, with "easy" predicted to have a faster (lower) reactive time than "hard", and "congruent" predicted to have a faster RT than "incongruent", but we can also customise the order of levels with `add_contrast()`; see the [contrasts vignette](https://debruine.github.io/faux/articles/contrasts.html) for more details. +You can change the mean and SD. Simulate a lot of values (1e5 == 100,000), save them in a variable, and visualise them with `hist()`. ```{r} -sub_n <- 2 # number of subjects in this simulation -stim_n <- 2 # number of stimuli in this simulation +x <- rnorm(1e5, mean = 30, sd = 5) -dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) |> - add_between(.by = "sub", condition = c("easy","hard")) |> - add_within(version = c("congruent", "incongruent")) |> - add_contrast("condition") |> - add_contrast("version") - -dat +hist(x) ``` -The function defaults to very descriptive names that help you interpret the fixed factors. Here, "condition.hard-easy" means the main effect of this factor is interpreted as the RT for hard trials minus the RT for easy trials, and "version.incongruent-congruent" means the main effect of this factor is interpreted as the RT for incongruent trials minus the RT for congruent trials. However, we can change these to simpler labels with the `colnames` argument. +## Multivariate normal +But how do you create correlated values? You can do this with `MASS::mvrnorm()`, but you need to construct the `Sigma` argument yourself from the correlation matrix and the standard deviations of the populations, and then you need to turn the resulting matrix into a data frame for many use cases. This isn't very difficult, but can be tedious with larger numbers of variables. -### Random Effects +```{r} +n = 1e5 # this is a large number to demonstrate that the result is as expected +mu = c(A = 1, B = 2, C = 3) +sd = c(0.5, 1, 1.5) +r = c(0, .25, .5) -Now we specify the random effect structure. We'll just add random intercepts to start, but will conver random slopes later. +cor_mat <- matrix(c(1, r[1], r[2], + r[1], 1, r[3], + r[2], r[3], 1), + nrow = 3) +Sigma <- (sd %*% t(sd)) * cor_mat +vars <- MASS::mvrnorm(n, mu, Sigma) |> as.data.frame() -Each subject will have slightly faster or slower reaction times on average; this is their random intercept (`sub_i`). We'll model it from a normal distribution with a mean of 0 and SD of 100ms. +cor(vars) |> round(2) +``` -Each stimulus will have slightly faster or slower reaction times on average; this is their random intercept (`stim_i`). We'll model it from a normal distribution with a mean of 0 and SD of 50ms (it seems reasonable to expect less variability between words than people for this task). +### rnorm_multi -Run this code a few times to see how the random effects change each time. this is because they are **sampled** from populations. +In faux, you can create sets of correlated normally distributed values using `rnorm_multi()`. ```{r} -sub_n <- 2 # number of subjects in this simulation -stim_n <- 2 # number of stimuli in this simulation -sub_sd <- 100 # SD for the subjects' random intercept -stim_sd <- 50 # SD for the stimuli's random intercept +dat3 <- rnorm_multi( + n = 50, + mu = c(A = 1, B = 2, C = 3), + sd = c(0.5, 1, 1.5), + r = c(0, .25, .5) +) +``` -dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) |> - add_between(.by = "sub", condition = c("easy","hard")) |> - add_within(version = c("congruent", "incongruent")) |> - add_contrast("condition", colnames = "cond") |> - add_contrast("version", colnames = "vers") |> - add_ranef(.by = "sub", sub_i = sub_sd) |> - add_ranef(.by = "stim", stim_i = stim_sd) +The function `get_params()` gives you a quick way to see the means, SDs and correlations in the simulated data set to make sure you set the parameters correctly. -dat +```{r} +get_params(dat3) ``` -### Error Term - -Finally, add an error term. This uses the same `add_ranef()` function, just without specifying which random factor it's for with `.by`. In essence, this samples an error value from a normal distribution with a mean of 0 and the specified SD for each trial. We'll also increase the number of subjects and stimuli to more realistic values now. +If you set `empirical` to `TRUE`, the values you set will be the **sample** parameters, not the **population** parameters. This isn't usually what you want for a simulation, but can be useful to check you set the parameters correctly. ```{r} -sub_n <- 200 # number of subjects in this simulation -stim_n <- 50 # number of stimuli in this simulation -sub_sd <- 100 # SD for the subjects' random intercept -stim_sd <- 50 # SD for the stimuli's random intercept -error_sd <- 25 # residual (error) SD +dat3 <- rnorm_multi( + n = 50, + mu = c(A = 1, B = 2, C = 3), + sd = c(0.5, 1, 1.5), + r = c(0, .25, .5), + empirical = TRUE +) -dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) |> - add_between(.by = "sub", condition = c("easy","hard")) |> - add_within(version = c("congruent", "incongruent")) |> - add_contrast("condition", colnames = "cond") |> - add_contrast("version", colnames = "vers") |> - add_ranef(.by = "sub", sub_i = sub_sd) |> - add_ranef(.by = "stim", stim_i = stim_sd) |> - add_ranef(err = error_sd) +get_params(dat3) ``` -### Calculate DV -Now we can calculate the DV by adding together an overall intercept (mean RT for all trials), the subject-specific intercept, the stimulus-specific intercept, and an error term, plus the effect of subject condition, the effect of stimulus version, and the interaction between condition and version. +### Setting r -We set these effects in raw units (ms). So when we set the effect of subject condition (`sub_cond_eff`) to 50, that means the average difference between the easy and hard condition is 50ms. `Easy` was coded as -0.5 and `hard` was coded as +0.5, which means that trials in the easy condition have -0.5 \* 50ms (i.e., -25ms) added to their reaction time, while trials in the hard condition have +0.5 \* 50ms (i.e., +25ms) added to their reaction time. +You can set the `r` argument for correlations in a few different ways. -```{r sim-dv} -sub_n <- 200 # number of subjects in this simulation -stim_n <- 50 # number of stimuli in this simulation -sub_sd <- 100 # SD for the subjects' random intercept -stim_sd <- 50 # SD for the stimuli's random intercept -error_sd <- 25 # residual (error) SD -grand_i <- 400 # overall mean DV -cond_eff <- 50 # mean difference between conditions: hard - easy -vers_eff <- 50 # mean difference between versions: incongruent - congruent -cond_vers_ixn <- 0 # interaction between version and condition +If all correlations have the same value, just set r equal to a single number. -dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) |> - add_between(.by = "sub", condition = c("easy","hard")) |> - add_within(version = c("congruent", "incongruent")) |> - add_contrast("condition", colnames = "cond") |> - add_contrast("version", colnames = "vers") |> - add_ranef(.by = "sub", sub_i = sub_sd) |> - add_ranef(.by = "stim", stim_i = stim_sd) |> - add_ranef(err = error_sd) |> - mutate(dv = grand_i + sub_i + stim_i + err + - (cond * cond_eff) + - (vers * vers_eff) + - (cond * vers * cond_vers_ixn) # in this example, this is always 0 and could be omitted - ) +```{r} +# all correlations the same value +rho_same <- rnorm_multi(50, 4, r = .5, empirical = TRUE) +get_params(rho_same) ``` -As always, graph to make sure you've simulated the general pattern you expected. +You can set `r` to a vector or matrix of the full correlation matrix. This is convenient when you're getting the values from an existing dataset, where you can just use the output of the `cor()` function. -```{r plot-dv, fig.cap="Double-check the simulated pattern"} -ggplot(dat, aes(condition, dv, color = version)) + - geom_hline(yintercept = grand_i) + - geom_violin(alpha = 0.5) + - stat_summary(fun = mean, - fun.min = \(x){mean(x) - sd(x)}, - fun.max = \(x){mean(x) + sd(x)}, - position = position_dodge(width = 0.9)) + - scale_color_brewer(palette = "Dark2") +```{r} +rho <- cor(iris[1:4]) +round(rho, 2) ``` +Notice how, since we didn't specify the names of the 4 variables anywhere else, `rnorm_multi()` will take them from the named correlation matrix. + +```{r} +rho_cormat <- rnorm_multi(50, 4, r = rho, empirical = TRUE) +get_params(rho_cormat) +``` -### Interactions +Alternatively, you can just specify the values from the upper right triangle of a correlation matrix. This might be easier if you're reading the values out of a paper. -If you want to simulate an interaction, it can be tricky to figure out what to set the main effects and interaction effect to. It can be easier to think about the simple main effects for each cell. Create four new variables and set them to the deviations from the overall mean you'd expect for each condition (so they should add up to 0). Here, we're simulating a small effect of version in the hard condition (50ms difference) and double that effect of version in the easy condition (100ms difference). +```{r} +# upper right triangle +# X2 X3 X4 +rho <- c(0.5, 0.4, 0.3, # X1 + 0.2, 0.1, # X2 + 0.0) # X3 -```{r sim-simple-main-effects} -# set variables to use in calculations below -hard_congr <- -25 -hard_incon <- +25 -easy_congr <- -50 -easy_incon <- +50 +rho_urt <- rnorm_multi(50, 4, r = rho, empirical = TRUE) +get_params(rho_urt) ``` -Use the code below to transform the simple main effects above into main effects and interactions for use in the equations below. -```{r sim-effect-calc} -# calculate main effects and interactions from simple effects above +## Factorial Designs -# mean difference between easy and hard conditions -cond_eff <- (hard_congr + hard_incon)/2 - - (easy_congr + easy_incon)/2 +You can use `rnorm_multi()` to simulate data for each between-subjects cell of a factorial design and manually combine the tables, but faux has a function that better maps onto how we usually think and teach about factorial designs. -# mean difference between incongruent and congruent versions -vers_eff <- (hard_incon + easy_incon)/2 - - (hard_congr + easy_congr)/2 +The default design is 100 observations of one variable (named `y`) with a mean of 0 and SD of 1. Unless you set `plot = FALSE` or run `faux_options(plot = FALSE)`, this function will show you a plot of your design so you can check that it looks like you expect. -# interaction between version and condition -cond_vers_ixn <- (hard_incon - hard_congr) - - (easy_incon - easy_congr) +```{r} +simdat1 <- sim_design() ``` -Then generate the DV the same way we did above, but also add the interaction effect multiplied by the effect-coded subject condition and stimulus version. -```{r sim-ixn} +### Factors -dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) |> - add_between(.by = "sub", condition = c("easy","hard")) |> - add_within(version = c("congruent", "incongruent")) |> - add_contrast("condition", colnames = "cond") |> - add_contrast("version", colnames = "vers") |> - add_ranef(.by = "sub", sub_i = sub_sd) |> - add_ranef(.by = "stim", stim_i = stim_sd) |> - add_ranef(err = error_sd) |> - mutate(dv = grand_i + sub_i + stim_i + err + - (cond * cond_eff) + - (vers * vers_eff) + - (cond * vers * cond_vers_ixn) - ) +Use named lists to set the names and levels of `within` and `between` subject factors. +```{r} +pettime <- sim_design( + between = list(pet = c("cat", "dog", "ferret")), + within = list(time = c("pre", "post")) +) ``` -```{r plot-ixn, fig.cap="Double-check the interaction between condition and version"} -ggplot(dat, aes(condition, dv, color = version)) + - geom_hline(yintercept = grand_i) + - geom_violin(alpha = 0.5) + - stat_summary(fun = mean, - fun.min = \(x){mean(x) - sd(x)}, - fun.max = \(x){mean(x) + sd(x)}, - position = position_dodge(width = 0.9)) + - scale_color_brewer(palette = "Dark2") +You can set `mu` and `sd` with unnamed vectors, but getting the order right can take some trial and error. + +```{r} +pettime <- sim_design( + between = list(pet = c("cat", "dog", "ferret")), + within = list(time = c("pre", "post")), + mu = 1:6 +) ``` +You can set values with a named vector for a single type of factor. The values do not have to be in the right order if they're named. -```{r table-ixn} -group_by(dat, condition, version) %>% - summarise(m = mean(dv) - grand_i %>% round(1), - .groups = "drop") %>% - pivot_wider(names_from = version, - values_from = m) +```{r} +pettime <- sim_design( + between = list(pet = c("cat", "dog", "ferret")), + within = list(time = c("pre", "post")), + mu = c(cat = 1, ferret = 5, dog = 3), + sd = c(pre = 1, post = 2) +) ``` -## Analysis +Or use a data frame for within- and between-subject factors. -New we will run a linear mixed effects model with `lmer` and look at the summary. +```{r} +pettime <- sim_design( + between = list(pet = c("cat", "dog", "ferret")), + within = list(time = c("pre", "post")), + mu = data.frame( + pre = c(1, 3, 5), + post = c(2, 4, 6), + row.names = c("cat", "dog", "ferret") + ) +) +``` -```{r lmer} -mod <- lmer(dv ~ cond * vers + - (1 | sub) + - (1 | stim), - data = dat) +If you have within-subject factors, set the correlations for each between-subject cell like this. -mod.sum <- summary(mod) +```{r} +pettime <- sim_design( + between = list(pet = c("cat", "dog", "ferret")), + within = list(time = c("pre", "post")), + r = list(cat = 0.5, + dog = 0.25, + ferret = 0), + empirical = TRUE, + plot = FALSE +) -mod.sum +get_params(pettime) ``` -### Sense checks - -First, check that your groups make sense. +You can also change the name of the `dv` and `id` columns and output the data in long format. If you do this, you also need to tell `get_params()` what columns contain the between- and within-subject factors, the dv, and the id. -* The number of obs should be the total number of trials analysed. -* `sub` should be what we set `sub_n` to above. -* `stim` should be what we set `stim_n` to above. +```{r} +dat_long <- sim_design( + between = list(pet = c("cat", "dog", "ferret")), + within = list(time = c("pre", "post")), + id = "subj_id", + dv = "score", + long = TRUE, + plot = FALSE +) -```{r mod-ngrps} -mod.sum$ngrps |> - as_tibble(rownames = "Random.Fator") |> - mutate(parameters = c(sub_n, stim_n)) +get_params(dat_long, digits = 3) ``` -Next, look at the random effects. +### Multiple Factors -* The SD for `sub` should be near `sub_sd`. -* The SD for `stim` should be near `stim_sd`. -* The residual SD should be near `error_sd`. +Set more than one within-or between-subject factor like this: -```{r mod-varcor} -mod.sum$varcor |> - as_tibble() |> - select(Groups = grp, Name = var1, "Std.Dev." = sdcor) |> - mutate(parameters = c(sub_sd, stim_sd, error_sd)) +```{r} +dat_multi <- sim_design( + between = list(pet = c("cat", "dog", "ferret"), + country = c("UK", "NL")), + within = list(time = c("pre", "post"), + condition = c("ctl", "exp")), + mu = data.frame( + cat_UK = 1:4, + cat_NL = 5:8, + dog_UK = 9:12, + dog_NL = 13:16, + ferret_UK = 17:20, + ferret_NL = 21:24, + row.names = c("pre_ctl", "pre_exp", "post_ctl", "post_exp") + ) +) ``` -Finally, look at the fixed effects. -* The estimate for the Intercept should be near the `grand_i`. -* The main effect of `cond` should be near what we calculated for `cond_eff`. -* The main effect of `vers` should be near what we calculated for `vers_eff`. -* The interaction between `cond`:`vers` should be near what we calculated for `cond_vers_ixn`. +Because faux uses an underscore for the separator, you have to set the `sep` argument to something different if you want to use underscores in your variable names (or set the separator globally with `faux_options`). -```{r mod-coef} -mod.sum$coefficients |> - as_tibble(rownames = "Effect") |> - select(Effect, Estimate) |> - mutate(parameters = c(grand_i, cond_eff, vers_eff, cond_vers_ixn)) +```{r} +# faux_options(sep = ".") + +dat_multi <- sim_design( + between = list(pet = c("cat", "dog", "ferret"), + country = c("Glasgow_UK", "Rotterdam_NL")), + within = list(time = c("pre", "post"), + condition = c("ctl", "exp")), + mu = data.frame( + cat.Glasgow_UK = 1:4, + cat.Rotterdam_NL = 5:8, + dog.Glasgow_UK = 9:12, + dog.Rotterdam_NL = 13:16, + ferret.Glasgow_UK = 17:20, + ferret.Rotterdam_NL = 21:24, + row.names = c("pre.ctl", "pre.exp", "post.ctl", "post.exp") + ), + sep = "." +) ``` -### Random effects - -Plot the subject intercepts from our code above (`dat$sub_i`) against the subject intercepts calculated by `lmer` (`ranef(mod)$sub_id`). +### Anonymous Factors -```{r plot-sub-ranef, fig.cap = "Compare simulated subject random intercepts to those from the model"} -# get simulated random intercept for each subject -sub_sim <- dat |> - group_by(sub, sub_i) |> - summarise(.groups = "drop") +If you need to make a quick demo, you can set factors anonymously with integer vectors. For example, the following code makes 3B\*2B\*2W mixed design. -# join to calculated random intercept from model -sub_sim_mod <- ranef(mod)$sub |> - as_tibble(rownames = "sub") |> - rename(mod_sub_i = `(Intercept)`) |> - left_join(sub_sim, by = "sub") +```{r} +dat_anon <- sim_design( + n = 50, + between = c(3, 2), + within = 2, + mu = 1:12 +) +``` -# plot to check correspondence -sub_sim_mod |> - ggplot(aes(sub_i,mod_sub_i)) + - geom_point() + - geom_smooth(method = "lm", formula = y~x) + - xlab("Simulated random intercepts (sub_i)") + - ylab("Modeled random intercepts") -``` - -Plot the stimulus intercepts from our code above (`dat$stim_i`) against the stimulus intercepts calculated by `lmer` (`ranef(mod)$stim_id`). - -```{r plot-stim-ranef, fig.cap = "Compare simulated stimulus random intercepts to those from the model"} -# get simulated random intercept for each stimulus -stim_sim <- dat |> - group_by(stim, stim_i) |> - summarise(.groups = "drop") +Faux has a quick plotting function for visualising data made with faux. The plot created by `sim_design()` shows the *design*, while this function shows the simulated *data*. -# join to calculated random intercept from model -stim_sim_mod <- ranef(mod)$stim |> - as_tibble(rownames = "stim") |> - rename(mod_stim_i = `(Intercept)`) |> - left_join(stim_sim, by = "stim") - -# plot to check correspondence -stim_sim_mod |> - ggplot(aes(stim_i,mod_stim_i)) + - geom_point() + - geom_smooth(method = "lm", formula = y~x) + - xlab("Simulated random intercepts (stim_i)") + - ylab("Modeled random intercepts") -``` - - -### Function - -You can put the code above in a function so you can run it more easily and change the parameters. I removed the plot and set the argument defaults to the same as the example above with all fixed effects set to 0, but you can set them to other patterns. - -```{r sim-function} -sim_lmer <- function( sub_n = 200, - stim_n = 50, - sub_sd = 100, - stim_sd = 50, - error_sd = 25, - grand_i = 400, - cond_eff = 0, - vers_eff = 0, - cond_vers_ixn = 0) { - dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) |> - add_between(.by = "sub", condition = c("easy","hard")) |> - add_within(version = c("congruent", "incongruent")) |> - add_contrast("condition", colnames = "cond") |> - add_contrast("version", colnames = "vers") |> - add_ranef(.by = "sub", sub_i = sub_sd) |> - add_ranef(.by = "stim", stim_i = stim_sd) |> - add_ranef(err = error_sd) |> - mutate(dv = grand_i + sub_i + stim_i + err + - (cond * cond_eff) + - (vers * vers_eff) + - (cond * vers * cond_vers_ixn) - ) - - mod <- lmer(dv ~ cond * vers + - (1 | sub) + - (1 | stim), - data = dat) - - return(mod) -} +```{r} +plot(dat_anon) ``` -Run the function with the default values (so all fixed effects set to 0). +You can change the order of plotting and the types of geoms plotted. This takes a little trial and error, so this function will probably be refined in later versions. -```{r sim-lmer-default} -sim_lmer() %>% summary() +```{r} +plot(dat_anon, "B1", "B2", "W1", geoms = c("violin", "pointrangeSD")) ``` -Try changing some variables to simulate different patterns of fixed effects. -```{r sim-lmer-changes} -sim_lmer(cond_eff = 0, - vers_eff = 75, - cond_vers_ixn = -50) %>% - summary() -``` -### Power analysis +## Replications + +You often want to simulate data repeatedly to do things like calculate power. The `sim_design()` function has a lot of overhead for checking if a design makes sense and if the correlation matrix is possible, so you can speed up the creation of multiple datasets with the same design using the `rep` argument. This will give you a nested data frame with each dataset in the `data` column. -First, wrap your simulation function inside of another function that takes the argument of a replication number, runs a simulated analysis, and returns a data table of the fixed and random effects (made with `broom.mixed::tidy()`). You can use purrr's `map_df()` function to create a data table of results from multiple replications of this function. We're only running 10 replications here in the interests of time, but you'll want to run 100 or more for a proper power calculation. +```{r} +dat_rep <- sim_design( + within = 2, + n = 20, + mu = c(0, 0.25), + rep = 5, + plot = FALSE +) +``` + +### Analyse each replicate -```{r power1} +You can run analyses on the nested data by wrapping your analysis code in a function then using `map()` to run the analysis on each data set and `unnest()` to expand the results into a data table. -sim_lmer_pwr <- function(rep) { - s <- sim_lmer(cond_eff = 0, - vers_eff = 75, - cond_vers_ixn = 50) - - # put just the fixed effects into a data table - broom.mixed::tidy(s, "fixed") %>% - mutate(rep = rep) # add a column for which rep +```{r} +# define function +analyse <- function(data) { + t.test(data$W1a, data$W1b, paired = TRUE) %>% broom::tidy() } -my_power <- map_df(1:10, sim_lmer_pwr) +# get one test data set +data <- dat_rep$data[[1]] +# check function returns what you want +analyse(data) ``` -You can then plot the distribution of estimates across your simulations. ```{r} -ggplot(my_power, aes(estimate, color = term)) + - geom_density() + - facet_wrap(~term, scales = "free") +# run the function on each data set +dat_rep |> + mutate(analysis = map(data, analyse)) |> + select(-data) |> + unnest(analysis) ``` -You can also just calculate power as the proportion of p-values less than your alpha. +### ANOVA + +Use the same pattern to run an ANOVA on a version of the `pettime` dataset. + +First, simulate 100 datasets in long format. These data will have small main effects of pet and time, but no interaction. ```{r} -my_power %>% - group_by(term) %>% - summarise(power = mean(p.value < 0.05), - .groups = "drop") +pettime100 <- sim_design( + between = list(pet = c("cat", "dog")), + within = list(time = c("pre", "post")), + n = c(cat = 50, dog = 40), + mu = data.frame( + pre = c(1, 1.2), + post = c(1.2, 1.4), + row.names = c("cat", "dog") + ), + sd = 1, + id = "pet_id", + dv = "score", + r = 0.5, + long = TRUE, + rep = 100 +) ``` +Then set up your analysis. We'll use the `aov_ez()` function from the {afex} package because its arguments match those of `sim_design()`. There's a little setup to run first to get rid of annoying messages and make this run faster by omitting calculations we won't need. -## Random slopes - -In the example so far we've ignored random variation among subjects or stimuli in the size of the fixed effects (i.e., **random slopes**). +```{r} +afex::set_sum_contrasts() # avoids annoying afex message +afex_options(include_aov = FALSE) # runs faster +afex_options(es_aov = "pes") # changes effect size measure to partial eta squared +``` -First, let's reset the parameters we set above. +This custom function takes the data frame as input and runs our ANOVA on it. The code at the end just cleans up the resulting table a bit. ```{r} -sub_n <- 200 # number of subjects in this simulation -stim_n <- 50 # number of stimuli in this simulation -sub_sd <- 100 # SD for the subjects' random intercept -stim_sd <- 50 # SD for the stimuli's random intercept -error_sd <- 25 # residual (error) SD -grand_i <- 400 # overall mean DV -cond_eff <- 50 # mean difference between conditions: hard - easy -vers_eff <- 50 # mean difference between versions: incongruent - congruent -cond_vers_ixn <- 0 # interaction between version and condition - +analyse <- function(data) { + a <- afex::aov_ez( + id = "pet_id", + dv = "score", + between = "pet", + within = "time", + data = data + ) + # return anova_table for GG-corrected DF + as_tibble(a$anova_table, rownames = "term") |> + mutate(term = factor(term, levels = term)) |> # keeps terms in order + rename(p.value = `Pr(>F)`) # fixes annoying p.value name +} ``` -### Slopes +Test the analysis code on the first simulated data frame. -In addition to generating a random intercept for each subject, now we will also generate a random slope for any within-subject factors. The only within-subject factor in this design is `version`. The main effect of `version` is set to 50 above, but different subjects will show variation in the size of this effect. That's what the random slope captures. We'll set `sub_vers_sd` below to the SD of this variation and use this to calculate the random slope (`sub_version_slope`) for each subject. +```{r} +analyse( pettime100$data[[1]] ) +``` -Also, it's likely that the variation between subjects in the size of the effect of version is related in some way to between-subject variation in the intercept. So we want the random intercept and slope to be correlated. Here, we'll simulate a case where subjects who have slower (larger) reaction times across the board show a smaller effect of condition, so we set `sub_i_vers_cor` below to a negative number (-0.2). -We just have to edit the first `add_ranef()` to add two variables (`sub_i`, `sub_vers_slope`) that are correlated with r = -0.2, means of 0, and SDs equal to what we set `sub_sd` above and `sub_vers_sd` below. +Use the same code we used in the first example to make a table of the results of each analysis: -```{r sim-subject-cor} -sub_vers_sd <- 20 -sub_i_vers_cor <- -0.2 +```{r} +pettime_sim <- pettime100 |> + mutate(analysis = map(data, analyse)) |> + select(-data) |> + unnest(analysis) +``` -dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) |> - add_between(.by = "sub", condition = c("easy","hard")) |> - add_within(version = c("congruent", "incongruent")) |> - add_contrast("condition", colnames = "cond") |> - add_contrast("version", colnames = "vers") |> - add_ranef(.by = "sub", sub_i = sub_sd, - sub_vers_slope = sub_vers_sd, - .cors = sub_i_vers_cor) +```{r, echo = FALSE} +# show the first 6 rows +head(pettime_sim) |> + mutate(across(5:8, \(x) round(x, 3))) ``` +Then you can summarise the data to calculate things like power for each effect or mean effect size. -### Correlated Slopes +```{r} +pettime_sim |> + group_by(term) |> + summarise(power = mean(p.value < 0.05), + mean_pes = mean(pes) |> round(3), + .groups = "drop") +``` -In addition to generating a random intercept for each stimulus, we will also generate a random slope for any within-stimulus factors. Both `version` and `condition` are within-stimulus factors (i.e., all stimuli are seen in both `congruent` and `incongruent` versions and both `easy` and `hard` conditions). So the main effects of version and condition (and their interaction) will vary depending on the stimulus. +The power for the between-subjects effect of pet is smaller than for the within-subjects effect of time. What happens if you reduce the correlation between pre and post? -They will also be correlated, but in a more complex way than above. You need to set the correlations for all pairs of slopes and intercept. Let's set the correlation between the random intercept and each of the slopes to -0.4 and the slopes all correlate with each other +0.2 (You could set each of the six correlations separately if you want, though). +## Non-normal Distributions +The newest version of faux has a new function for simulating non-normal distributions using the NORTA method (NORmal To Anything). The `dist` argument lists the variables with their distribution names (e.g., "norm", "pois", unif", "truncnorm", or anything that has an "rdist" function). The `params` argument lists the distribution function argument values for each variable (e.g., arguments to `rnorm`, `rpois`, `runif`, `rtruncnorm`). -```{r rslope-sim-stimuli} +This function simulates multivariate non-normal distributions by using simulation to work out the correlations for a multivariate normal distribution that will produce the desired correlations after the normal distributions are converted to the desired distributions. This simulation can take a while if you have several variables and should warn you if you're requesting an impossible combination (but is still an experimental function, so let Lisa know if you have any problems). -stim_vers_sd <- 10 # SD for the stimuli's random slope for stim_version -stim_cond_sd <- 30 # SD for the stimuli's random slope for sub_cond -stim_cond_vers_sd <- 15 # SD for the stimuli's random slope for sub_cond:stim_version -stim_i_cor <- -0.4 # correlations between intercept and slopes -stim_s_cor <- +0.2 # correlations among slopes +```{r} +dat_norta <- rmulti( + n = 1000, + dist = c(U = "unif", + T = "truncnorm", + L = "likert"), + params = list( + U = list(min = 0, max = 10), + T = list(a = 1, b = 7, mean = 3.5, sd = 2.1), + L = list(prob = c(`much less` = .10, + `less` = .20, + `equal` = .35, + `more` = .25, + `much more` = .10)) + ), + r = c(-0.5, 0, 0.5) +) +``` -# specify correlations for rnorm_multi (one of several methods) -stim_cors <- c(stim_i_cor, stim_i_cor, stim_i_cor, - stim_s_cor, stim_s_cor, - stim_s_cor) +The "likert" type is a set of distribution functions provided by faux to make creating Likert scale variables easier (see `?rlikert`). You may need to convert Likert-scale variables to numbers before analysis or calculating descriptives. -dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) |> - add_between(.by = "sub", condition = c("easy","hard")) |> - add_within(version = c("congruent", "incongruent")) |> - add_contrast("condition", colnames = "cond") |> - add_contrast("version", colnames = "vers") |> - add_ranef(.by = "sub", sub_i = sub_sd, - sub_vers_slope = sub_vers_sd, - .cors = sub_i_vers_cor) |> - add_ranef(.by = "stim", stim_i = stim_sd, - stim_vers_slope = stim_vers_sd, - stim_cond_slope = stim_cond_sd, - stim_cond_vers_slope = stim_cond_vers_sd, - .cors = stim_cors) +```{r} +# convert likert-scale variable to integer +dat_norta$L <- as.integer(dat_norta$L) +get_params(dat_norta) ``` -### Calculate DV -Now we can calculate the DV by adding together an overall intercept (mean RT for all trials), the subject-specific intercept, the stimulus-specific intercept, the effect of subject condition, the stimulus-specific slope for condition, the effect of stimulus version, the stimulus-specific slope for version, the subject-specific slope for condition, the interaction between condition and version (set to 0 for this example), the stimulus-specific slope for the interaction between condition and version, and an error term. +## Exercises -```{r rslope-sim-dv} -dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) |> - add_between(.by = "sub", condition = c("easy","hard")) |> - add_within(version = c("congruent", "incongruent")) |> - add_contrast("condition", colnames = "cond") |> - add_contrast("version", colnames = "vers") |> - add_ranef(.by = "sub", sub_i = sub_sd, - sub_vers_slope = sub_vers_sd, - .cors = sub_i_vers_cor) |> - add_ranef(.by = "stim", stim_i = stim_sd, - stim_vers_slope = stim_vers_sd, - stim_cond_slope = stim_cond_sd, - stim_cond_vers_slope = stim_cond_vers_sd, - .cors = stim_cors) |> - add_ranef(err = error_sd) |> - mutate( - trial_cond_eff = cond_eff + stim_cond_slope, - trial_vers_eff = vers_eff + sub_vers_slope + stim_vers_slope, - trial_cond_vers_ixn = cond_vers_ixn + stim_cond_vers_slope, - dv = grand_i + sub_i + stim_i + err + - (cond * trial_cond_eff) + - (vers * trial_vers_eff) + - (cond * vers * trial_cond_vers_ixn) - ) +### Multivariate normal -``` +Sample 40 values of three variables named `J`, `K` and `L` from a population with means of 10, 20 and 30, and SDs of 5. `J` and `K` are correlated 0.5, `J` and `L` are correlated 0.25, and `K` and `L` are not correlated. -As always, graph to make sure you've simulated the general pattern you expected. +```{r} -```{r rslope-plot-dv, fig.cap="Double-check the simulated pattern"} -ggplot(dat, aes(condition, dv, color = version)) + - geom_hline(yintercept = grand_i) + - geom_violin(alpha = 0.5) + - stat_summary(fun = mean, - fun.min = \(x){mean(x) - sd(x)}, - fun.max = \(x){mean(x) + sd(x)}, - position = position_dodge(width = 0.9)) + - scale_color_brewer(palette = "Dark2") ``` -## Analysis +### From existing data -New we'll run a linear mixed effects model with `lmer` and look at the summary. You specify random slopes by adding the within-level effects to the random intercept specifications. Since the only within-subject factor is version, the random effects specification for subjects is `(1 + vers | sub)`. Since both condition and version are within-stimuli factors, the random effects specification for stimuli is `(1 + vers*cond | stim)`. +Using the data from the built-in dataset `attitude`, simulate a new set of 20 observations drawn from a population with the same means, SDs and correlations for each column as the original data. -This model will take a lot longer to run than one without random slopes specified. This might be a good time for a coffee break. +```{r} -```{r rslope-lmer} -mod <- lmer(dv ~ cond * vers + - (1 + vers || sub) + - (1 + vers*cond || stim), - data = dat) - -mod.sum <- summary(mod) - -mod.sum -``` - -### Sense checks - -First, check that your groups make sense. - -* `sub` = `sub_n` (`r sub_n`) -* `stim` = `stim_n` (`r stim_n`) - -```{r rslope-mod-ngrps} -mod.sum$ngrps |> - as_tibble(rownames = "Random.Fator") |> - mutate(parameters = c(sub_n, stim_n)) -``` - -Next, look at the SDs for the random effects. - -* Group:`sub` - * `(Intercept)` ~= `sub_sd` - * `vers` ~= `sub_vers_sd` -* Group: `stim` - * `(Intercept)` ~= `stim_sd` - * `vers` ~= `stim_vers_sd` - * `cond` ~= `stim_cond_sd` - * `vers:cond` ~= `stim_cond_vers_sd` -* Residual ~= `error_sd` - -```{r rslope-mod-varcor} -mod.sum$varcor |> - as_tibble() |> - select(Groups = grp, Name = var1, "Std.Dev." = sdcor) |> - mutate(parameters = c(sub_sd, sub_vers_sd, stim_sd, stim_vers_sd, stim_cond_sd, stim_cond_vers_sd, error_sd)) -``` - -The correlations are a bit more difficult to parse. The first column under `Corr` shows the correlation between the random slope for that row and the random intercept. So for `vers` under `sub`, the correlation should be close to `sub_i_vers_cor`. For all three random slopes under `stim`, the correlation with the random intercept should be near `stim_i_cor` and their correlations with each other should be near `stim_s_cor`. - - -Finally, look at the fixed effects. - -* `(Intercept)` ~= `grand_i` -* `sub_cond.e` ~= `sub_cond_eff` -* `stim_version.e` ~= `stim_vers_eff` -* `sub_cond.e`:`stim_version.e` ~= `cond_vers_ixn` - -```{r rslope-mod-coef} -mod.sum$coefficients |> - as_tibble(rownames = "Effect") |> - select(Effect, Estimate) |> - mutate(parameters = c(grand_i, cond_eff, vers_eff, cond_vers_ixn)) -``` - - -### Function - -You can put the code above in a function so you can run it more easily and change the parameters. I removed the plot and set the argument defaults to the same as the example above, but you can set them to other patterns. - -```{r rslope-sim-function} -sim_lmer_slope <- function( sub_n = 200, - stim_n = 50, - sub_sd = 100, - sub_vers_sd = 20, - sub_i_vers_cor = -0.2, - stim_sd = 50, - stim_vers_sd = 10, - stim_cond_sd = 30, - stim_cond_vers_sd = 15, - stim_i_cor = -0.4, - stim_s_cor = +0.2, - error_sd = 25, - grand_i = 400, - sub_cond_eff = 0, - stim_vers_eff = 0, - cond_vers_ixn = 0) { - dat <- add_random(sub = sub_n) |> - add_random(stim = stim_n) |> - add_between(.by = "sub", condition = c("easy","hard")) |> - add_within(version = c("congruent", "incongruent")) |> - add_contrast("condition", colnames = "cond") |> - add_contrast("version", colnames = "vers") |> - add_ranef(.by = "sub", sub_i = sub_sd, - sub_vers_slope = sub_vers_sd, - .cors = sub_i_vers_cor) |> - add_ranef(.by = "stim", stim_i = stim_sd, - stim_vers_slope = stim_vers_sd, - stim_cond_slope = stim_cond_sd, - stim_cond_vers_slope = stim_cond_vers_sd, - .cors = stim_cors) |> - add_ranef(err = error_sd) |> - mutate( - trial_cond_eff = cond_eff + stim_cond_slope, - trial_vers_eff = vers_eff + sub_vers_slope + stim_vers_slope, - trial_cond_vers_ixn = cond_vers_ixn + stim_cond_vers_slope, - dv = grand_i + sub_i + stim_i + err + - (cond * trial_cond_eff) + - (vers * trial_vers_eff) + - (cond * vers * trial_cond_vers_ixn) - ) - - mod <- lmer(dv ~ cond * vers + - (1 + vers || sub) + - (1 + vers*cond || stim), - data = dat) - - return(mod) -} ``` -Run the function with the default values (null fixed effects). -```{r rslope-sim-lmer-default} -sim_lmer_slope() %>% summary() -``` +### 2b + +Create a dataset with a between-subject factor of "pet" having two levels, "cat", and "dog". The DV is "happiness" score. There are 20 cat-owners with a mean happiness score of 10 (SD = 3) and there are 30 dog-owners with a mean happiness score of 11 (SD = 3). -Try changing some variables to simulate fixed effects. +```{r} -```{r rslope-sim-lmer-null} -sim_lmer_slope(sub_cond_eff = 50, - stim_vers_eff = 50, - cond_vers_ixn = 0) ``` -## Exercises +### 3w -1. Calculate power for the parameters in the last example using the `sim_lmer_slope()` function. +Create a dataset of 20 observations with 1 within-subject variable ("condition") having 3 levels ("A", "B", "C") with means of 10, 20 and 30 and SD of 5. The correlations between each level have r = 0.4. The dataset should look like this: -```{r ex1} +| id | condition | score | +|:---|:----------|------:| +|S01 | A | 9.17 | +|... | ... | ... | +|S20 | A | 11.57 | +|S01 | B | 18.44 | +|... | ... | ... | +|S20 | B | 20.04 | +|S01 | C | 35.11 | +|... | ... | ... | +|S20 | C | 29.16 | + +```{r} ``` +### 2w*2w -2. Simulate data for the following design: +Create a dataset with 50 observations of 2 within-subject variables ("W1" and "W2") each having 2 levels. The mean for all cells is 10 and the SD is 2. The dataset should have 20 subjects. The correlations look like this: -* 100 raters rate 50 faces from group A and 50 faces from group B -* The DV has a mean value of 50 -* Group B values are 5 points higher than group A -* Rater intercepts have an SD of 5 -* Face intercepts have an SD of 10 -* The residual error has an SD of 8 +| | W1a_W2a | W1a_W2b | W1b_W2a | W1b_W2b | +|:--------|------:|------:|------:|------:| +| W1a_W2a | 1.0 | 0.5 | 0.5 | 0.2 | +| W1a_W2b | 0.5 | 1.0 | 0.2 | 0.5 | +| W1b_W2a | 0.5 | 0.2 | 1.0 | 0.5 | +| W1b_W2b | 0.2 | 0.5 | 0.5 | 1.0 | -```{r ex2} +```{r} ``` -3. For the design from exercise 2, write a function that simulates data and runs a mixed effects analysis on it. +### 2w*3b -```{r ex3} +Create a dataset with a between-subject factor of "pet" having 3 levels ("cat", "dog", and "ferret") and a within-subject factor of "time" having 2 levels ("pre" and "post"). The N in each group should be 10. Means are: + +* cats: pre = 10, post = 12 +* dogs: pre = 14, post = 16 +* ferrets: pre = 18, post = 20 + +SDs are all 5 and within-cell correlations are all 0.25. + +```{r} ``` -4. The package `faux` has a built-in dataset called `fr4`. Type `?faux::fr4` into the console to view the help for this dataset. Run a mixed effects model on this dataset looking at the effect of `face_sex` on ratings. Remember to include a random slope for the effect of face sex and explicitly add a contrast code. +### Replications + +Create 5 datasets with a 2b*2b design, 30 participants in each cell. Each cell's mean should be 0, except B1a:B2a, which should be 0.5. The SD should be 1. Make the resulting data in long format. ```{r} -# code female = -0.5, male = +0.5 ``` -5. Use the parameters from this analysis to simulate a new dataset with 50 male and 50 female faces, and 100 raters. +### Power + +Simulate 100 datasets like the one above and use `lm()` or `afex::aov_ez()` to look at the interaction between B1 and B2. What is the power of this design? ```{r} diff --git a/index.log b/index.log index 31df125..0c3f090 100644 --- a/index.log +++ b/index.log @@ -1,4 +1,4 @@ -This is XeTeX, Version 3.141592653-2.6-0.999995 (TeX Live 2023) (preloaded format=xelatex 2023.7.12) 31 JAN 2024 08:47 +This is XeTeX, Version 3.141592653-2.6-0.999995 (TeX Live 2023) (preloaded format=xelatex 2023.7.12) 31 JAN 2024 15:05 entering extended mode restricted \write18 enabled. %&-line parsing enabled. diff --git a/index.tex b/index.tex index fd9688b..1cbc613 100644 --- a/index.tex +++ b/index.tex @@ -190,7 +190,7 @@ \makeatletter \makeatother \makeatletter -\ifdefined\Shaded\renewenvironment{Shaded}{\begin{tcolorbox}[sharp corners, boxrule=0pt, frame hidden, enhanced, borderline west={3pt}{0pt}{shadecolor}, interior hidden, breakable]}{\end{tcolorbox}}\fi +\ifdefined\Shaded\renewenvironment{Shaded}{\begin{tcolorbox}[enhanced, frame hidden, breakable, sharp corners, boxrule=0pt, borderline west={3pt}{0pt}{shadecolor}, interior hidden]}{\end{tcolorbox}}\fi \makeatother \makeatletter \@ifpackageloaded{sidenotes}{}{\usepackage{sidenotes}} @@ -520,7 +520,7 @@ \chapter{Installing R and RStudio}\label{installing-r-and-rstudio}} \end{Highlighting} \end{Shaded} -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-warning-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-warning-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-warning-color}{\faExclamationTriangle} \end{minipage}% @@ -534,7 +534,7 @@ \chapter{Installing R and RStudio}\label{installing-r-and-rstudio}} \end{minipage}% \end{tcolorbox} -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-important-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-important-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-important-color}{\faExclamation} \end{minipage}% @@ -557,7 +557,7 @@ \section{Working with R projects and R accessible to your current R session. In short, R projects set the working directory. -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-tip-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-tip-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-tip-color}{\faLightbulb} \end{minipage}% @@ -573,7 +573,7 @@ \section{Working with R projects and R \end{minipage}% \end{tcolorbox} -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-caution-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-caution-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-caution-color}{\faFire} \end{minipage}% @@ -633,7 +633,7 @@ \section{CSV format}\label{csv-format}} \end{Highlighting} \end{Shaded} -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-tip-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-tip-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-tip-color}{\faLightbulb} \end{minipage}% @@ -657,7 +657,7 @@ \section{Data from Excel}\label{data-from-excel}} Importing data from Excel into R, you can use the \texttt{read\_excel()}-function. -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-important-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-important-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-important-color}{\faExclamation} \end{minipage}% @@ -692,7 +692,7 @@ \section{Useful functions}\label{useful-functions}} \hypertarget{common-import-mistakes}{% \subsection{Common import mistakes}\label{common-import-mistakes}} -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-tip-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-tip-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-tip-color}{\faLightbulb} \end{minipage}% @@ -710,7 +710,7 @@ \subsection{Common import mistakes}\label{common-import-mistakes}} \end{minipage}% \end{tcolorbox} -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-important-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-important-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-important-color}{\faExclamation} \end{minipage}% @@ -779,7 +779,7 @@ \chapter{\texorpdfstring{Data wrangling using the Functions of the \texttt{tidyverse} allow you to perform data wrangling easily. -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-tip-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-tip-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-tip-color}{\faLightbulb} \end{minipage}% @@ -879,7 +879,7 @@ \section{\texorpdfstring{Change the data format from \emph{long} to \item each measurement has one column \item - each entity (e.g.~person) has one column + each entity (e.g.~person) has one row \end{itemize} This data format makes it easy to spot missing values or outliers and @@ -897,7 +897,7 @@ \section{\texorpdfstring{Change the data format from \emph{long} to \begin{center}\rule{0.5\linewidth}{0.5pt}\end{center} -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-tip-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-tip-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-tip-color}{\faLightbulb} \end{minipage}% @@ -989,7 +989,7 @@ \subsection{\texorpdfstring{\texttt{pivot\_longer()}}{pivot\_longer()}}\label{pi $ value 32.33111, 51.20389, 55.99303, 55.38460, 61.41110, 83.33978, ~ \end{verbatim} -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-tip-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-tip-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-tip-color}{\faLightbulb} \end{minipage}% @@ -1248,9 +1248,7 @@ \section{\texorpdfstring{Creating graphs with \NormalTok{d\_summary }\OtherTok{\textless{}{-}}\NormalTok{ d }\SpecialCharTok{|\textgreater{}} \FunctionTok{group\_by}\NormalTok{(condition) }\SpecialCharTok{|\textgreater{}} \FunctionTok{summarise}\NormalTok{(}\AttributeTok{mean\_x =} \FunctionTok{mean}\NormalTok{(x),} - \AttributeTok{sd\_x =} \FunctionTok{sd}\NormalTok{(x),} - \AttributeTok{mean\_y =} \FunctionTok{mean}\NormalTok{(y),} - \AttributeTok{sd\_y =} \FunctionTok{sd}\NormalTok{(x))} + \AttributeTok{mean\_y =} \FunctionTok{mean}\NormalTok{(y))} \end{Highlighting} \end{Shaded} @@ -1258,13 +1256,13 @@ \section{\texorpdfstring{Creating graphs with \begin{itemize} \item - \textbf{data} + \textbf{\emph{data}} \item - \textbf{geoms}, visible forms (\emph{aesthetics}) such as points, - lines or boxes. + \textbf{\emph{geoms}}, visible forms (\emph{aesthetics}) such as + points, lines or boxes. \item - \textbf{a coordinate systems / mapping} describes how data and geoms - are linked, also colors or grouping variables are specified here + \textbf{\emph{a coordinate systems / mapping}} describes how data and + geoms are linked, also colors or grouping variables are specified here \end{itemize} Further components could be: @@ -1277,16 +1275,16 @@ \section{\texorpdfstring{Creating graphs with \item coordinate functions \item - \textbf{facets} + \textbf{\emph{facets}} \item scales \item - \textbf{themes} + \textbf{\emph{themes}} \end{itemize} -(we will only cover the contents in \textbf{bold}) +(we will only cover the contents in italics) -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-tip-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-tip-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-tip-color}{\faLightbulb} \end{minipage}% @@ -1296,9 +1294,13 @@ \section{\texorpdfstring{Creating graphs with \begin{itemize} \item - It is easiest when your data is in \emph{long} format. + For plotting with \texttt{ggplot()} it is easiest when your data is in + \emph{long} format. \item What variables do you want to plot (categorical? continuous? \ldots) + affects which \texttt{geoms}can be used. You can try out what is + suited with the \texttt{esquisse}-package below or find ideas + \href{https://www.data-to-viz.com}{here}. \end{itemize} \end{minipage}% @@ -1307,8 +1309,8 @@ \section{\texorpdfstring{Creating graphs with \hypertarget{data-geoms-and-mapping}{% \section{Data, geoms and mapping}\label{data-geoms-and-mapping}} -We start with entering the current data frame and add geoms and mappings -(\texttt{aes()}) with arguments such as +We start with entering the current data frame and add \texttt{geoms} and +mappings (specified with \texttt{aes()}) with arguments such as \begin{Shaded} \begin{Highlighting}[] @@ -1351,8 +1353,8 @@ \section{Data, geoms and mapping}\label{data-geoms-and-mapping}} \marginnote{\begin{footnotesize} -{[}Here{]}https://github.com/rstudio/cheatsheets/blob/main/data-visualization.pdf) -you can download \texttt{ggplot}-Cheatsheet( +\href{https://github.com/rstudio/cheatsheets/blob/main/data-visualization.pdf}{Here} +you can download \texttt{ggplot}-Cheatsheet \end{footnotesize}} @@ -1388,8 +1390,9 @@ \section{Themes and labels}\label{themes-and-labels}} \AttributeTok{y =}\NormalTok{ y)) }\SpecialCharTok{+} \FunctionTok{geom\_point}\NormalTok{() }\SpecialCharTok{+} \FunctionTok{ggtitle}\NormalTok{ (}\StringTok{"Title"}\NormalTok{) }\SpecialCharTok{+} - \FunctionTok{xlab}\NormalTok{(}\StringTok{"Value 1 [a.u.]"}\NormalTok{) }\SpecialCharTok{+} - \FunctionTok{ylab}\NormalTok{(}\StringTok{"Value 2 [a.u.]"}\NormalTok{) }\SpecialCharTok{+} + \FunctionTok{labs}\NormalTok{(}\AttributeTok{title =} \StringTok{"Title"}\NormalTok{, } + \AttributeTok{x =} \StringTok{"Variable A [a.u.]"}\NormalTok{,} + \AttributeTok{y =} \StringTok{"Variable B [a.u.]"}\NormalTok{) }\SpecialCharTok{+} \FunctionTok{theme\_minimal}\NormalTok{() }\CommentTok{\# also theme\_classic and theme\_minimal are nice} \end{Highlighting} \end{Shaded} @@ -1402,13 +1405,13 @@ \section{Themes and labels}\label{themes-and-labels}} \end{figure} -\begin{tcolorbox}[enhanced jigsaw, rightrule=.15mm, colback=white, colframe=quarto-callout-tip-color-frame, left=2mm, leftrule=.75mm, breakable, bottomrule=.15mm, arc=.35mm, opacityback=0, toprule=.15mm] +\begin{tcolorbox}[enhanced jigsaw, arc=.35mm, breakable, colback=white, opacityback=0, leftrule=.75mm, bottomrule=.15mm, toprule=.15mm, colframe=quarto-callout-tip-color-frame, left=2mm, rightrule=.15mm] \begin{minipage}[t]{5.5mm} \textcolor{quarto-callout-tip-color}{\faLightbulb} \end{minipage}% \begin{minipage}[t]{\textwidth - 5.5mm} -\textbf{`esquisse`-package}\vspace{2mm} +\textbf{\texttt{esquisse}-package}\vspace{2mm} With this package you can use the data frames in your current environment or load a new one to try out which geoms might be useful @@ -1424,11 +1427,21 @@ \section{Themes and labels}\label{themes-and-labels}} \end{minipage}% \end{tcolorbox} +\marginnote{\begin{footnotesize} + +If we do not want to load the whole library we can call the function by +writing the library name then \texttt{::} and the function afterwards +such as \texttt{esquisse::esquisser()}. This is helpful when libraries +have overlapping function names which can create conflicts or to easily +call a function without typing \texttt{library()} first. + +\end{footnotesize}} + Further helpful ressources: \begin{itemize} \item - \href{https://ggplot2.tidyverse.org/}{\texttt{ggplot}documentation} + \href{https://ggplot2.tidyverse.org/}{\texttt{ggplot2}documentation} \item \href{https://psyteachr.github.io/reprores-v3}{Website PsyTeachR: Data Skills for reproducible research} @@ -1443,7 +1456,8 @@ \section{Themes and labels}\label{themes-and-labels}} \item \href{https://www.data-to-viz.com}{Here} you find ideas for plots. \item - https://rstudio-connect.psy.gla.ac.uk/plotdemo/ + \href{https://rstudio-connect.psy.gla.ac.uk/plotdemo}{Demo of + different geoms} \end{itemize}