search

lizzieinvancouver / ospree

Budbreak review paper database
3 stars 0 forks source link

pop Up! Try to fit a 'jointish' model #406

Open lizzieinvancouver opened 3 years ago

lizzieinvancouver commented 3 years ago

We want to fit something similar to the traits joint model, but in the traits joint model we need to use a model to estimate a trait value per species ... in this case we have observed range 'traits' (if you will) per species. So we can just send that data in....

lizzieinvancouver commented 3 years ago

BTW, this is a Pop-Up Power Rangers project or a PUPR (pupper).

lizzieinvancouver commented 3 years ago

I posted two versions of what I think that we think is the same model ... here's the one where climvar is a vector of just 20 values (because we have 20 species):

+ mod1 <- stan('stan/jointish_climvar_emw1.stan', data = faker1,
+            iter = 5000, warmup=4000, chains=4)
+ )
   user  system elapsed 
 38.410   1.604 719.527 
> summary(mod1)$summary
                            mean      se_mean          sd          2.5%           25%           50%
alphaForcingSp[1]     -0.8721852 0.0048071262 0.123691088 -1.115785e+00    -0.9519442    -0.8733018
alphaForcingSp[2]     -0.8727996 0.0032577173 0.091305587 -1.045541e+00    -0.9352620    -0.8742532
alphaForcingSp[3]     -0.9259078 0.0038824060 0.101677990 -1.127473e+00    -0.9969916    -0.9287131
alphaForcingSp[4]     -0.8627635 0.0053787094 0.137279104 -1.140315e+00    -0.9506067    -0.8640570
alphaForcingSp[5]     -1.0042306 0.0034588516 0.096029225 -1.184241e+00    -1.0708773    -1.0071405
alphaForcingSp[6]     -0.8770218 0.0034926379 0.096595825 -1.059097e+00    -0.9416645    -0.8780712
alphaForcingSp[7]     -0.7753613 0.0029251913 0.086442874 -9.519398e-01    -0.8320421    -0.7749118
alphaForcingSp[8]     -1.0210610 0.0053405094 0.131948077 -1.288403e+00    -1.1053015    -1.0195811
alphaForcingSp[9]     -0.9771594 0.0046934695 0.120246076 -1.216709e+00    -1.0537158    -0.9723179
alphaForcingSp[10]    -1.0642241 0.0043241040 0.111709921 -1.275959e+00    -1.1420967    -1.0673766
alphaForcingSp[11]    -0.8455084 0.0033082482 0.096412044 -1.042572e+00    -0.9063058    -0.8424194
alphaForcingSp[12]    -0.8403149 0.0038011240 0.103133107 -1.035058e+00    -0.9106108    -0.8425737
alphaForcingSp[13]    -0.7801135 0.0023517993 0.080883257 -9.292549e-01    -0.8343221    -0.7833504
alphaForcingSp[14]    -0.9355176 0.0030136590 0.086944509 -1.108178e+00    -0.9915678    -0.9339033
alphaForcingSp[15]    -0.9454989 0.0049940600 0.125184569 -1.201980e+00    -1.0238539    -0.9417964
alphaForcingSp[16]    -0.8457470 0.0046018019 0.118114190 -1.077593e+00    -0.9225361    -0.8447717
alphaForcingSp[17]    -0.7761088 0.0040038333 0.105815635 -9.835414e-01    -0.8466749    -0.7756105
alphaForcingSp[18]    -0.7200880 0.0051333181 0.132529477 -9.631874e-01    -0.8087734    -0.7246042
alphaForcingSp[19]    -0.8128534 0.0035416935 0.099842475 -1.001318e+00    -0.8793698    -0.8165481
alphaForcingSp[20]    -0.9208422 0.0043364925 0.114834127 -1.154357e+00    -0.9914855    -0.9202166
muForceSp             -0.8840126 0.0038128721 0.086680322 -1.058448e+00    -0.9410007    -0.8825551
sigmaForceSp           0.1222828 0.0011250530 0.036080089  5.361277e-02     0.0986272     0.1205141
alphaPhenoSp[1]      149.2230283 0.0255756114 1.424192666  1.463217e+02   148.3294446   149.2674845
alphaPhenoSp[2]      149.3755948 0.0312695363 1.417600426  1.465035e+02   148.4259221   149.4228659
alphaPhenoSp[3]      148.6386246 0.0361732715 1.459066596  1.458168e+02   147.6090624   148.6909364
alphaPhenoSp[4]      150.3038899 0.0264505340 1.441114629  1.475401e+02   149.3226870   150.2887776
alphaPhenoSp[5]      148.9555724 0.0427734984 1.610891150  1.457203e+02   147.8735738   148.9993455
alphaPhenoSp[6]      149.2422833 0.0313174860 1.440968884  1.463511e+02   148.2971733   149.2797735
alphaPhenoSp[7]      152.1274943 0.0466251746 1.572159481  1.491596e+02   151.0261902   152.1016964
alphaPhenoSp[8]      149.0433717 0.0355365371 1.524981827  1.459860e+02   148.0121051   149.0866352
alphaPhenoSp[9]      150.8144128 0.0286814016 1.492101077  1.480043e+02   149.7878382   150.7499936
alphaPhenoSp[10]     146.9904982 0.0675939976 1.952925271  1.430548e+02   145.6548902   147.0714787
alphaPhenoSp[11]     152.5525735 0.0421904519 1.561295418  1.496449e+02   151.4602838   152.5039276
alphaPhenoSp[12]     149.1601407 0.0312593665 1.473283550  1.461488e+02   148.1890234   149.2427947
alphaPhenoSp[13]     149.5237928 0.0400822780 1.541821039  1.463714e+02   148.5276466   149.5651693
alphaPhenoSp[14]     150.4855769 0.0358346348 1.449136246  1.477705e+02   149.4816762   150.4389702
alphaPhenoSp[15]     150.5036692 0.0262353889 1.406086481  1.477652e+02   149.5533713   150.4694236
alphaPhenoSp[16]     151.0595517 0.0276864654 1.437780682  1.483355e+02   150.0921178   151.0087554
alphaPhenoSp[17]     151.7768208 0.0401151720 1.496433286  1.489846e+02   150.7200896   151.7551249
alphaPhenoSp[18]     151.5751020 0.0444720228 1.692882430  1.482960e+02   150.4268677   151.5388497
alphaPhenoSp[19]     150.1347069 0.0284387806 1.462757764  1.471439e+02   149.2082990   150.1431283
alphaPhenoSp[20]     151.5992569 0.0299203180 1.480782045  1.488962e+02   150.5957380   151.5469973
muPhenoSp            150.1582621 0.0222623699 0.794861736  1.486061e+02   149.6420037   150.1497421
sigmaPhenoSp           2.0331897 0.0286312342 0.699892789  7.374330e-01     1.5678989     2.0032555
betaTraitxPheno       -0.2077395 0.0003062803 0.007209969 -2.214449e-01    -0.2124526    -0.2077646
sigmapheno_y           1.9884894 0.0003084689 0.022897896  1.943486e+00     1.9727075     1.9881545
betaForcingSp[1]      -3.9918760 0.0012575665 0.070900841 -4.128438e+00    -4.0385082    -3.9940223
betaForcingSp[2]      -2.2890288 0.0015696377 0.070786970 -2.423812e+00    -2.3361927    -2.2909741
betaForcingSp[3]      -2.7002408 0.0018066361 0.073099636 -2.838362e+00    -2.7502992    -2.7015891
betaForcingSp[4]      -4.7380730 0.0013161876 0.072216548 -4.883424e+00    -4.7844410    -4.7383949
betaForcingSp[5]      -2.1554048 0.0021189163 0.080211829 -2.306919e+00    -2.2110870    -2.1581596
betaForcingSp[6]      -2.4862788 0.0015571857 0.071823663 -2.622257e+00    -2.5348594    -2.4866962
betaForcingSp[7]      -1.5760769 0.0023226246 0.078320214 -1.730676e+00    -1.6284530    -1.5739194
betaForcingSp[8]      -4.3531238 0.0017756875 0.076135102 -4.501170e+00    -4.4053751    -4.3553995
betaForcingSp[9]      -3.8986942 0.0014412290 0.074584288 -4.052787e+00    -3.9461304    -3.8946375
betaForcingSp[10]     -2.1608840 0.0033425492 0.096850576 -2.338568e+00    -2.2307544    -2.1650087
betaForcingSp[11]     -2.5056514 0.0021120535 0.077934219 -2.665129e+00    -2.5569833    -2.5032401
betaForcingSp[12]     -2.8630261 0.0015704649 0.073863900 -3.003464e+00    -2.9140688    -2.8685262
betaForcingSp[13]     -1.0881068 0.0020166173 0.077454041 -1.231684e+00    -1.1401902    -1.0901436
betaForcingSp[14]     -1.9962023 0.0017964653 0.072158361 -2.140743e+00    -2.0440948    -1.9939391
betaForcingSp[15]     -4.1819461 0.0013091574 0.070067569 -4.321816e+00    -4.2268293    -4.1800115
betaForcingSp[16]     -3.8734796 0.0013711216 0.071641804 -4.016060e+00    -3.9212947    -3.8716264
betaForcingSp[17]     -3.1175392 0.0020217039 0.074597006 -3.262985e+00    -3.1690368    -3.1161369
betaForcingSp[18]     -4.2452320 0.0022366237 0.084838183 -4.421026e+00    -4.3007669    -4.2433240
betaForcingSp[19]     -2.6250536 0.0014249988 0.073418227 -2.763123e+00    -2.6744599    -2.6253486
betaForcingSp[20]     -3.6455362 0.0014797159 0.073538316 -3.800270e+00    -3.6910881    -3.6427957
lp__               -4738.8314475 0.3235320824 6.757529278 -4.752056e+03 -4743.1616053 -4738.8836706

With a Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.

... I think just 20 species could be too small for this model?

lizzieinvancouver commented 3 years ago

I think my other model has some issues ...

+ mod2 <- stan('stan/jointish_climvar_emw2.stan', data = faker2,
+            iter = 5000, warmup=4000, chains=4)
+ )
2021-03-25 21:22:56.194 R[1605:245785] _mthid_copyDeviceInfo(288230376989705376) failed
    user   system  elapsed 
  38.738    1.831 2858.989 
Warning messages:
1: There were 4000 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded 
2: Examine the pairs() plot to diagnose sampling problems

3: The largest R-hat is 2.74, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat
MoralesCastilla commented 3 years ago

It seems like the fake data code ´joint_climvar_3param_db.R´ using less observations per species and betaTraitxCue set at 0.8 estimates parameters correctly. It still asks for a larger sampling size as some Rhats = 1.02, but may work just fine adding more samples. Parameter estimation:

                        mean      se_mean          sd       2.5%
betaTraitxForcing  0.7987747 0.0002883781 0.005397887  0.7884098
betaTraitxPhoto   -0.7970646 0.0002295233 0.005045735 -0.8067002
betaTraitxChill   -0.8045046 0.0002570714 0.005491033 -0.8153085
                         25%        50%        75%      97.5%    n_eff
betaTraitxForcing  0.7951064  0.7987083  0.8024428  0.8092751 350.3669
betaTraitxPhoto   -0.8004965 -0.7971257 -0.7936497 -0.7869484 483.2760
betaTraitxChill   -0.8082450 -0.8044602 -0.8008208 -0.7937969 456.2476
                      Rhat
betaTraitxForcing 1.023319
betaTraitxPhoto   1.004419
betaTraitxChill   1.003448
lizzieinvancouver commented 3 years ago

Ospree data are running! It's quite quick....

> threeparam_jntsum <- summary(threeparam_jnt)$summary
Warning message:
package ‘rstan’ was built under R version 3.5.2 
Warning message:
package ‘rstan’ was built under R version 3.5.2 
Warning message:
package ‘rstan’ was built under R version 3.5.2 
Warning message:
package ‘rstan’ was built under R version 3.5.2 
> threeparam_jntsum[grep("betaTraitx", rownames(threeparam_jntsum)),]
                         mean      se_mean         sd        2.5%         25%         50%
betaTraitxForcing -0.03031993 0.0004577382 0.02539916 -0.08046285 -0.04729785 -0.03022923
betaTraitxPhoto    0.03510949 0.0004925764 0.02743391 -0.01897828  0.01708794  0.03536639
betaTraitxChill    0.05943293 0.0020155715 0.07475655 -0.09044616  0.01046504  0.06038009
                          75%      97.5%    n_eff     Rhat
betaTraitxForcing -0.01347184 0.01954228 3078.963 1.000091
betaTraitxPhoto    0.05302144 0.08861849 3101.903 1.001731
betaTraitxChill    0.10931032 0.20255568 1375.631 1.002820
> threeparam_jntsum[grep("mu", rownames(threeparam_jntsum)),]
                mean    se_mean       sd       2.5%        25%        50%        75%      97.5%
muForceSp  -6.632500 0.01193242 1.010539  -8.560249  -7.314854  -6.636042 -5.9785624 -4.5884189
muPhotoSp  -1.545179 0.01492896 1.074353  -3.716318  -2.240851  -1.517326 -0.8279111  0.5249435
muChillSp -10.912947 0.03459044 2.929700 -16.689187 -12.803654 -10.954031 -9.0253547 -5.1561508
muPhenoSp  28.646315 0.02476051 2.416847  23.933472  27.058989  28.631369 30.2210094 33.5111501
             n_eff      Rhat
muForceSp 7172.148 1.0000641
muPhotoSp 5178.867 1.0005522
muChillSp 7173.547 1.0003320
muPhenoSp 9527.500 0.9999936
dbuona commented 3 years ago

I got similar answers with zscoring the climate variable too. It's always fun when your models directly contradict you hypotheses

dbuona commented 3 years ago

and by similar I mean in terms of sign not magnitude. Currently running on STV and range area too

lizzieinvancouver commented 3 years ago

@dbuona Hmm, Nacho, Faith and I though these agree with our hypotheses?

lizzieinvancouver commented 3 years ago

Higher temporal GDD variability reduces forcing cue and increases photo and chill cues.

lizzieinvancouver commented 3 years ago

Anyway, you never pushed your code so I have added files with commit 0460aefed97bf7af7b733d47a0186b6ea1b0dfb2

lizzieinvancouver commented 3 years ago

popUpmodels We also worked through how to plot results, this figure also pushed.

lizzieinvancouver commented 3 years ago

@MoralesCastilla updated fulljoint_fakedata_db.R and got the fake data parameters back.

lizzieinvancouver commented 3 years ago

He just reduced sample size so it runs faster.

dbuona commented 3 years ago

but isn't a "stronger cue" more negative? so it should decrease it i thibk

lizzieinvancouver commented 3 years ago

Here's my answer on STV (did not save/push this code:

> bb.3param <- with(bb.stan, 
+                       list(yPhenoi = resp, 
+                            forcingi = force.z,
+                            photoi = photo.z,
+                            chillingi = chill.z,
+                            species = latbinum,
+                            N = nrow(bb.stan),
+                            n_spec = length(unique(bb.stan$complex.wname)),
+                            climvar=unique(bb.stan$STV.cent)
+                       ))
> 
> 
> threeparam_jnt = stan('popUP/stan/joint_climvar_3param_osp.stan', data = bb.3param, # this stan code is similar to joint_climvar_3param_emw.stan but with a more reasonable prior for the intercept mu
> threeparam_jntsum[grep("betaTraitx", rownames(threeparam_jntsum)),]
                         mean      se_mean         sd        2.5%         25%         50%
betaTraitxForcing -0.03031993 0.0004577382 0.02539916 -0.08046285 -0.04729785 -0.03022923
betaTraitxPhoto    0.03510949 0.0004925764 0.02743391 -0.01897828  0.01708794  0.03536639
betaTraitxChill    0.05943293 0.0020155715 0.07475655 -0.09044616  0.01046504  0.06038009
                          75%      97.5%    n_eff     Rhat
betaTraitxForcing -0.01347184 0.01954228 3078.963 1.000091
betaTraitxPhoto    0.05302144 0.08861849 3101.903 1.001731
betaTraitxChill    0.10931032 0.20255568 1375.631 1.002820
> threeparam_jntsum[grep("mu", rownames(threeparam_jntsum)),]
                mean    se_mean       sd       2.5%        25%        50%        75%      97.5%
muForceSp  -6.632500 0.01193242 1.010539  -8.560249  -7.314854  -6.636042 -5.9785624 -4.5884189
muPhotoSp  -1.545179 0.01492896 1.074353  -3.716318  -2.240851  -1.517326 -0.8279111  0.5249435
muChillSp -10.912947 0.03459044 2.929700 -16.689187 -12.803654 -10.954031 -9.0253547 -5.1561508
muPhenoSp  28.646315 0.02476051 2.416847  23.933472  27.058989  28.631369 30.2210094 33.5111501
             n_eff      Rhat
muForceSp 7172.148 1.0000641
muPhotoSp 5178.867 1.0005522
muChillSp 7173.547 1.0003320
muPhenoSp 9527.500 0.9999936
lizzieinvancouver commented 3 years ago

Which seems oddly identical ... so hoping something is wrong! And @dbuona maybe! I guess I would have to plot it before committing either way.

lizzieinvancouver commented 3 years ago

Just added notes to my stan code, note that jointish_climvar_db_fj.stan is very similar (we think) to jointish_climvar_emw2.stan and that these three models:

  1. jointish_climvar_emw1.stan
  2. jointish_climvar_emw2.stan
  3. jointish_climvar_db_fj.stan

Are all different ways to do the same thing.

dbuona commented 3 years ago

Here's what i got for z-scored climvars: GDDLF:

>                mean    se_mean       sd       2.5%        25%        50%        75%     97.5%    n_eff
>muForceSp  -6.610749 0.01724754 1.028159  -8.546165  -7.316989  -6.625479 -5.9452291 -4.549839 3553.579
>muPhotoSp  -1.553275 0.02169065 1.085228  -3.700667  -2.260165  -1.525772 -0.8216678  0.551102 2503.209
>muChillSp -10.787781 0.04305072 2.918009 -16.346341 -12.765340 -10.795988 -8.8129717 -5.027160 4594.225
>muPhenoSp  28.632628 0.03627260 2.423646  23.800835  27.069804  28.636783 30.2547055 33.443286 4464.583

    >                   mean    se_mean       sd       2.5%        25%       50%        75%     97.5%
> betaTraitxForcing -1.185603 0.02535942 1.037302 -3.1861268 -1.8698276 -1.198662 -0.5235083 0.9113758
> betaTraitxPhoto    1.436883 0.02812626 1.115323 -0.7876327  0.7009272  1.405563  2.1968015 3.6413928
> betaTraitxChill    2.120400 0.11281205 2.908758 -3.4800977  0.1919739  2.131598  4.0485119 7.8638045

STV

      >          mean    se_mean       sd       2.5%        25%        50%        75%     97.5%    n_eff      Rhat
>muForceSp  -6.905474 0.01824638 1.014055  -8.882092  -7.592387  -6.910892 -6.2165099 -4.950378 3088.656 0.9993248
>muPhotoSp  -1.433797 0.02367997 1.101669  -3.670309  -2.166297  -1.411062 -0.6762489  0.665697 2164.412 1.0011698
>muChillSp -10.073937 0.04514473 2.804129 -15.494010 -11.985613 -10.080327 -8.2229745 -4.600555 3858.174 1.0011802
>muPhenoSp  28.689617 0.03653590 2.531763  23.593150  27.105327  28.696208 30.3199310 33.660097 4801.826 1.0000508

    >                    mean    se_mean        sd      2.5%         25%        50%         75%    97.5%     n_eff     Rhat
>betaTraitxForcing -0.6852567 0.02430261 0.9274110 -2.550505 -1.30301282 -0.6833394 -0.04838725 1.100520 1456.2594 1.001903
>betaTraitxPhoto    0.6102142 0.02137575 0.9368385 -1.158233 -0.01522165  0.5817210  1.23983561 2.508698 1920.8201 1.002291
>betaTraitxChill    3.5228529 0.10781632 2.7384754 -1.675302  1.62367921  3.4613933  5.37473195 8.932030  645.1322 1.015001

Range area:

>muForceSp  -5.983058 0.02779696 1.124715  -8.181468  -6.740190  -5.996958 -5.208986 -3.7128855 1637.158 0.9999861
>muPhotoSp  -2.216582 0.03027692 1.180744  -4.499631  -3.006336  -2.234730 -1.432503  0.1376175 1520.856 1.0001183
>muChillSp -10.237172 0.06793814 3.214857 -16.478756 -12.422126 -10.246302 -8.146321 -3.7823688 2239.217 1.0005190
>muPhenoSp  28.689687 0.03986042 2.504622  23.666482  27.008795  28.698919 30.374677 33.6254614 3948.213 1.0004273

     >                  mean    se_mean       sd       2.5%       25%       50%       75%    97.5%     n_eff     Rhat
>betaTraitxForcing  2.276480 0.04521120 1.400154 -0.3904778  1.316071  2.257980  3.234274 5.087126  959.0906 1.004410
>betaTraitxPhoto   -2.213859 0.04867302 1.574264 -5.2366863 -3.247098 -2.220223 -1.127718 0.881044 1046.1127 1.000540
>betaTraitxChill    1.281155 0.15528139 3.687943 -5.8831203 -1.237428  1.309937  3.798297 8.429928  564.0654 1.003017
lizzieinvancouver commented 3 years ago

@dbuona So range area matters to forcing! Huh ...

Also, totally fascinating how big that chill effect is for STV ... are you just dying to break the data down by continent as a quick check?

lizzieinvancouver commented 3 years ago

@dbuona I think it would be handy to have a two dataframes:

1.

2.

Then we could plot the lame JPG above and also compare how much the species-level effects change....

I also would be interested to include estimates from models where you run the NA and European species separately.

With these dataframes in hand I'd be exited to pop UP again!

dbuona commented 3 years ago

I made a data frame that combines all the relevant thing above, and a subsequent R script that subset them appropriately (though its rather hacky).

the csv is called betasandmorefromPOPUP.csv and the script is called quickplots_jointish.R

Here is what I found using a model where gdd2last frost is z.scored:

the species level cue estimates are are similar to our main ospree model

The slopes in gdd2lf appear to be in the opposite direct of our predictions

lizzieinvancouver commented 3 years ago

@dbuona Nice work! What are the axes in `The slopes in gdd2lf appear to be in the opposite direct of our predictions'? I don't understand the mean versus slope (y versus x).

dbuona commented 3 years ago

Mean and 50% credible intervals for Europe and North America for gdd to last frost as range trait. Note N. America model still has 10 divergent transitions.

North America: muForceSp -2.393017 (-3.576335 ,-1.2064065) muPhotoSp -2.122897 (-3.559914, -0.6803292) muChillSp -8.058155 (-11.897566 ,-4.3217652) muPhenoSp 24.635504 (21.680406, 27.4517953)

betaTraitxForcing -3.0992045 (-4.086134 , -2.1282060) betaTraitxPhoto -0.6363018 (-1.727646 , 0.5509935) betaTraitxChill -5.7735011 (-8.611913, -2.8715750)

Europe: muForceSp -3.05978212 (-5.757523 , -0.3348938) muPhotoSp -0.06132316 (-2.686767 , 2.5561890) muChillSp -5.76995515 (-9.254344 , -2.2917936) muPhenoSp 32.85736910 (30.726224 , 34.8665522)

betaTraitxForcing 7.5546654 (3.020477 , 12.210167) betaTraitxPhoto 0.9698243 (-3.607834, 5.392932) betaTraitxChill 3.1113403 (-2.423276, 8.569316)

This makes me feel we should probably just run all models separated by continent.

lizzieinvancouver commented 3 years ago

@dbuona I agree, these seem pretty crazy different. I suspect folks may want to see these models and the original ones though (continents lumped). In case it helps I pushed edits to your quickplots_jointish.R that show how to plot the JPG drawing from the Friday popup a couple weeks ago. I can help make more of these (and or better automate making them) if needed. See commit d85b787fa8fa8ca222d6163ea8b32c221729d99c

lizzieinvancouver commented 3 years ago

@dbuona I tried to put continent into the current model ... I am not sure I am thinking through it correctly, but it still seemed worth a shot. Check out my commit 02507a3ccd68815d653e02d3c78e4557445ef397

Here's the output:

> goobsumCont[grep("mu", rownames(goobsumCont)),]
                mean    se_mean        sd       2.5%        25%        50%         75%      97.5%     n_eff     Rhat
muForceSp  -7.566131 0.08229148  1.968881 -11.270203  -8.897317  -7.602597 -6.26467569  -3.655785  572.4388 1.009428
muPhotoSp  -1.340405 0.08065795  2.002290  -5.399411  -2.669522  -1.298353  0.06686449   2.413889  616.2538 1.005991
muChillSp -12.487485 0.11177645  3.751143 -19.739605 -14.988102 -12.536954 -9.87217753  -5.136650 1126.2287 1.003587
muPhenoSp  55.503028 1.77527979 28.262064  31.916088  36.242200  40.316101 74.79982998 119.174506  253.4394 1.009326
> goobsumCont[grep("betaTraitx", rownames(goobsumCont)),]
                       mean   se_mean       sd       2.5%        25%       50%       75%     97.5%    n_eff     Rhat
betaTraitxForcing -3.597192 0.1521513 3.782222 -10.816759 -6.2451897 -3.587035 -1.089241  4.065519 617.9348 1.007985
betaTraitxPhoto    1.712991 0.1486610 3.804397  -5.863439 -0.8361332  1.745826  4.320981  9.129367 654.9037 1.004633
betaTraitxChill   -1.472822 0.2008851 5.815275 -12.620075 -5.4356028 -1.508736  2.414328 10.180568 838.0020 1.002234
> goobsumCont[grep("betaFS", rownames(goobsumCont)),]
         mean       se_mean            sd          2.5%           25%           50%           75%         97.5%         n_eff 
  3.127290855   0.205229187   4.832428947  -6.615523285  -0.007941371   3.159068336   6.397933990  12.451705099 554.437636677 
         Rhat 
  1.008177163 
> goobsumCont[grep("betaPS", rownames(goobsumCont)),]
       mean     se_mean          sd        2.5%         25%         50%         75%       97.5%       n_eff        Rhat 
 -0.4260576   0.2002695   4.8506858  -9.5280198  -3.8509977  -0.5584549   2.7844262   9.5583633 586.6469108   1.0051739 
> goobsumCont[grep("betaCS", rownames(goobsumCont)),]
       mean     se_mean          sd        2.5%         25%         50%         75%       97.5%       n_eff        Rhat 
  5.3205996   0.2711819   7.3728968  -9.3696001   0.3383355   5.4387297  10.4015572  19.6005753 739.1877839   1.0041221 
>

I coded Europe as 0 and NAM as 1 ... So a negative effect of forcexTrait for Europe, but that is basically erased for NAM? Does that make sense in comparing to the models where you have you Europe and NAM running separately?

dbuona commented 3 years ago

That isn't what happens when I run the two continents separately, though it was what happened when I ran our previous "sequential models" with continent as an interaction. Definitely something to dig in to. I will talk a look at the model and try running it.

lizzieinvancouver commented 3 years ago

@dbuona Okay, I must have something off then ... I will try to think more. If you have any ideas, let me know.