lizzieinvancouver / localadaptclim

1 stars 0 forks source link

work on table formatting and inline estimates #24

Open lizzieinvancouver opened 1 year ago

lizzieinvancouver commented 1 year ago

@alinazeng I strongly recommend all tables in the supp and inline estimates of key numbers in the text. We should also work on what format we want the tables in the supp; the one table we have now reports stuff we don't want and possibly not all the stuff we do want. To start, can you:

alinazeng commented 1 year ago

@lizzieinvancouver

paste the 'summary' output of the spring event ~ MAT model below

summary(fitC_spring_mat)

Model Info:
 function:     stan_glmer
 family:       gaussian [identity]
 formula:      spring_event ~ (MAT_prov | species_garden)
 algorithm:    sampling
 sample:       4000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 671
 groups:       species_garden (20)

Estimates:
                                                        mean   sd    10%   50%   90%
(Intercept)                                           115.1    5.1 108.7 115.0 121.7
b[(Intercept) species_garden:Alnus_rubra_A]           -38.2    6.4 -46.3 -38.2 -29.9
b[MAT_prov species_garden:Alnus_rubra_A]                0.9    0.5   0.2   0.9   1.5
b[(Intercept) species_garden:Betula_papyrifera_K]     -18.0    6.0 -25.8 -18.0 -10.3
b[MAT_prov species_garden:Betula_papyrifera_K]         -3.4    0.7  -4.4  -3.4  -2.5
b[(Intercept) species_garden:Betula_papyrifera_L]       2.2    6.0  -5.5   2.1  10.0
b[MAT_prov species_garden:Betula_papyrifera_L]         -2.1    0.7  -3.0  -2.1  -1.2
b[(Intercept) species_garden:Betula_papyrifera_M]       5.6    5.9  -2.0   5.6  13.1
b[MAT_prov species_garden:Betula_papyrifera_M]         -1.4    0.7  -2.3  -1.4  -0.5
b[(Intercept) species_garden:Fagus_sylvatica_R*]       -0.3    7.5  -9.7  -0.2   9.4
b[MAT_prov species_garden:Fagus_sylvatica_R*]          -0.5    0.7  -1.4  -0.5   0.4
b[(Intercept) species_garden:Fagus_sylvatica_T*]        3.1    5.8  -4.2   3.2  10.5
b[MAT_prov species_garden:Fagus_sylvatica_T*]           1.0    0.3   0.6   1.0   1.5
b[(Intercept) species_garden:Fraxinus_excelsior_Q*]    41.9    6.9  33.0  41.8  50.8
b[MAT_prov species_garden:Fraxinus_excelsior_Q*]       -1.2    0.6  -1.9  -1.2  -0.4
b[(Intercept) species_garden:Picea_abies_S*]            2.8    6.1  -5.0   2.7  10.4
b[MAT_prov species_garden:Picea_abies_S*]               1.1    0.9   0.0   1.1   2.2
b[(Intercept) species_garden:Picea_engelmannii_B]     -12.3    5.2 -19.0 -12.2  -5.7
b[MAT_prov species_garden:Picea_engelmannii_B]          0.1    0.3  -0.4   0.1   0.5
b[(Intercept) species_garden:Picea_mariana_I]          25.1    5.2  18.4  25.2  31.6
b[MAT_prov species_garden:Picea_mariana_I]              0.6    0.6  -0.2   0.6   1.3
b[(Intercept) species_garden:Picea_sitchensis_D]        4.0    5.9  -3.6   4.0  11.4
b[MAT_prov species_garden:Picea_sitchensis_D]           0.3    0.4  -0.2   0.3   0.8
b[(Intercept) species_garden:Pinus_albicaulis_C]      -20.9    5.4 -27.8 -20.8 -14.1
b[MAT_prov species_garden:Pinus_albicaulis_C]           2.3    0.7   1.4   2.3   3.2
b[(Intercept) species_garden:Pinus_ponderosa_J]        51.2    6.8  42.5  51.1  59.7
b[MAT_prov species_garden:Pinus_ponderosa_J]           -0.1    0.5  -0.7  -0.1   0.5
b[(Intercept) species_garden:Populus_balsamifera_F]     2.7    5.5  -4.3   2.7   9.7
b[MAT_prov species_garden:Populus_balsamifera_F]       -0.3    0.6  -1.0  -0.3   0.4
b[(Intercept) species_garden:Populus_trichocarpa_D]   -44.2    5.3 -51.0 -44.2 -37.4
b[MAT_prov species_garden:Populus_trichocarpa_D]       -0.6    0.2  -0.9  -0.6  -0.3
b[(Intercept) species_garden:Pseudotsuga_menziesii_H]  13.4    5.9   6.0  13.4  20.8
b[MAT_prov species_garden:Pseudotsuga_menziesii_H]     -0.4    0.3  -0.8  -0.4   0.0
b[(Intercept) species_garden:Quercus_petraea_U*]       -3.4    5.9 -11.0  -3.2   4.0
b[MAT_prov species_garden:Quercus_petraea_U*]          -0.5    0.4  -1.0  -0.5  -0.1
b[(Intercept) species_garden:Quercus_petraea_V*]       -1.5    5.8  -9.1  -1.6   5.9
b[MAT_prov species_garden:Quercus_petraea_V*]          -0.8    0.4  -1.3  -0.8  -0.4
b[(Intercept) species_garden:Tsuga_heterophylla_E]    -30.2    6.9 -38.8 -30.2 -21.5
b[MAT_prov species_garden:Tsuga_heterophylla_E]         3.9    0.7   3.1   3.9   4.8
b[(Intercept) species_garden:Tsuga_heterophylla_G]     17.9   12.2   2.5  17.6  33.8
b[MAT_prov species_garden:Tsuga_heterophylla_G]        -0.8    1.3  -2.5  -0.8   0.8
sigma                                                   4.9    0.1   4.8   4.9   5.1
Sigma[species_garden:(Intercept),(Intercept)]         537.4  170.4 349.2 508.8 760.8
Sigma[species_garden:MAT_prov,(Intercept)]             -9.7   10.0 -22.2  -8.7   1.6
Sigma[species_garden:MAT_prov,MAT_prov]                 3.4    1.5   1.9   3.2   5.3

Fit Diagnostics:
           mean   sd    10%   50%   90%
mean_PPD 106.8    0.3 106.5 106.8 107.2

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
                                                      mcse Rhat n_eff
(Intercept)                                           0.2  1.0   579 
b[(Intercept) species_garden:Alnus_rubra_A]           0.2  1.0   916 
b[MAT_prov species_garden:Alnus_rubra_A]              0.0  1.0  7220 
b[(Intercept) species_garden:Betula_papyrifera_K]     0.2  1.0   804 
b[MAT_prov species_garden:Betula_papyrifera_K]        0.0  1.0  5819 
b[(Intercept) species_garden:Betula_papyrifera_L]     0.2  1.0   729 
b[MAT_prov species_garden:Betula_papyrifera_L]        0.0  1.0  5404 
b[(Intercept) species_garden:Betula_papyrifera_M]     0.2  1.0   739 
b[MAT_prov species_garden:Betula_papyrifera_M]        0.0  1.0  5664 
b[(Intercept) species_garden:Fagus_sylvatica_R*]      0.2  1.0  1266 
b[MAT_prov species_garden:Fagus_sylvatica_R*]         0.0  1.0  4530 
b[(Intercept) species_garden:Fagus_sylvatica_T*]      0.2  1.0   717 
b[MAT_prov species_garden:Fagus_sylvatica_T*]         0.0  1.0  6786 
b[(Intercept) species_garden:Fraxinus_excelsior_Q*]   0.2  1.0  1012 
b[MAT_prov species_garden:Fraxinus_excelsior_Q*]      0.0  1.0  6210 
b[(Intercept) species_garden:Picea_abies_S*]          0.2  1.0   777 
b[MAT_prov species_garden:Picea_abies_S*]             0.0  1.0  5548 
b[(Intercept) species_garden:Picea_engelmannii_B]     0.2  1.0   606 
b[MAT_prov species_garden:Picea_engelmannii_B]        0.0  1.0  6745 
b[(Intercept) species_garden:Picea_mariana_I]         0.2  1.0   605 
b[MAT_prov species_garden:Picea_mariana_I]            0.0  1.0  7027 
b[(Intercept) species_garden:Picea_sitchensis_D]      0.2  1.0   756 
b[MAT_prov species_garden:Picea_sitchensis_D]         0.0  1.0  6811 
b[(Intercept) species_garden:Pinus_albicaulis_C]      0.2  1.0   638 
b[MAT_prov species_garden:Pinus_albicaulis_C]         0.0  1.0  6452 
b[(Intercept) species_garden:Pinus_ponderosa_J]       0.2  1.0   903 
b[MAT_prov species_garden:Pinus_ponderosa_J]          0.0  1.0  6359 
b[(Intercept) species_garden:Populus_balsamifera_F]   0.2  1.0   654 
b[MAT_prov species_garden:Populus_balsamifera_F]      0.0  1.0  6217 
b[(Intercept) species_garden:Populus_trichocarpa_D]   0.2  1.0   632 
b[MAT_prov species_garden:Populus_trichocarpa_D]      0.0  1.0  4231 
b[(Intercept) species_garden:Pseudotsuga_menziesii_H] 0.2  1.0   749 
b[MAT_prov species_garden:Pseudotsuga_menziesii_H]    0.0  1.0  6144 
b[(Intercept) species_garden:Quercus_petraea_U*]      0.2  1.0   767 
b[MAT_prov species_garden:Quercus_petraea_U*]         0.0  1.0  6789 
b[(Intercept) species_garden:Quercus_petraea_V*]      0.2  1.0   712 
b[MAT_prov species_garden:Quercus_petraea_V*]         0.0  1.0  6040 
b[(Intercept) species_garden:Tsuga_heterophylla_E]    0.2  1.0   996 
b[MAT_prov species_garden:Tsuga_heterophylla_E]       0.0  1.0  5707 
b[(Intercept) species_garden:Tsuga_heterophylla_G]    0.2  1.0  2435 
b[MAT_prov species_garden:Tsuga_heterophylla_G]       0.0  1.0  3560 
sigma                                                 0.0  1.0  5493 
Sigma[species_garden:(Intercept),(Intercept)]         5.1  1.0  1124 
Sigma[species_garden:MAT_prov,(Intercept)]            0.3  1.0   938 
Sigma[species_garden:MAT_prov,MAT_prov]               0.0  1.0  1508 
mean_PPD                                              0.0  1.0  3469 
log-posterior                                         0.2  1.0   854 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
alinazeng commented 1 year ago

@lizzieinvancouver

Hi Lizzie,

give the mean slope of that model +/- 90% intervals here?

I am not sure if by "mean slope of the model +/- 90% intervals" you meant what's already included in the summary output above or something else. If it is something else, could you let me know how to look that up? what code? thank you!


> coef(fitC_spring_mat)$species
                           MAT_prov (Intercept)
Alnus rubra A            0.87315711    76.83859
Betula papyrifera K     -3.43686106    97.04248
Betula papyrifera L     -2.09609941   117.11437
Betula papyrifera M     -1.39615197   120.59769
Fagus sylvatica R*      -0.53129758   114.81497
Fagus sylvatica T*       1.04455939   118.24571
Fraxinus excelsior Q*   -1.19323612   156.82465
Picea abies S*           1.07532251   117.77697
Picea engelmannii B      0.05183523   102.84449
Picea mariana I          0.58266971   140.24542
Picea sitchensis D       0.26076261   119.04015
Pinus albicaulis C       2.28969754    94.27060
Pinus ponderosa J       -0.12447500   166.14255
Populus balsamifera F   -0.31024033   117.77528
Populus trichocarpa D   -0.64566448    70.88876
Pseudotsuga menziesii H -0.44995768   128.42657
Quercus petraea U*      -0.54415827   111.80105
Quercus petraea V*      -0.82197326   113.42651
Tsuga heterophylla E     3.94611951    84.86734
Tsuga heterophylla G    -0.81178943   132.62573
lizzieinvancouver commented 1 year ago

@alinazeng Okay, now that I see the output (thank you!) I think in general when we want the mean and 90% intervals we'll want them across all species. So for the above example, extract the posterior of each MAT_prov slope for each species, then divide by the number of species and use that posterior to calculate the mean and 10% and 90% quantiles. This is similar to what you did in issue #16 ... but do all species since we're summarizing the full model.

BTW, if you put the code in code chunks (not quote chunks) in markdown then they are easier to read. I updated them above.

lizzieinvancouver commented 1 year ago
> Estimates:
>                                                         mean   sd    10%   50%   90%
> (Intercept)                                           115.1    5.1 108.7 115.0 121.7
> b[(Intercept) species_garden:Alnus_rubra_A]           -38.2    6.4 -46.3 -38.2 -29.9
> b[MAT_prov species_garden:Alnus_rubra_A]                0.9    0.5   0.2   0.9   1.5
> b[(Intercept) species_garden:Betula_papyrifera_K]     -18.0    6.0 -25.8 -18.0 -10.3
> b[MAT_prov species_garden:Betula_papyrifera_K]         -3.4    0.7  -4.4  -3.4  -2.5
> b[(Intercept) species_garden:Betula_papyrifera_L]       2.2    6.0  -5.5   2.1  10.0
> b[MAT_prov species_garden:Betula_papyrifera_L]         -2.1    0.7  -3.0  -2.1  -1.2
> b[(Intercept) species_garden:Betula_papyrifera_M]       5.6    5.9  -2.0   5.6  13.1
> b[MAT_prov species_garden:Betula_papyrifera_M]         -1.4    0.7  -2.3  -1.4  -0.5
> b[(Intercept) species_garden:Fagus_sylvatica_R*]       -0.3    7.5  -9.7  -0.2   9.4
> b[MAT_prov species_garden:Fagus_sylvatica_R*]          -0.5    0.7  -1.4  -0.5   0.4
> b[(Intercept) species_garden:Fagus_sylvatica_T*]        3.1    5.8  -4.2   3.2  10.5
> b[MAT_prov species_garden:Fagus_sylvatica_T*]           1.0    0.3   0.6   1.0   1.5
> b[(Intercept) species_garden:Fraxinus_excelsior_Q*]    41.9    6.9  33.0  41.8  50.8
> b[MAT_prov species_garden:Fraxinus_excelsior_Q*]       -1.2    0.6  -1.9  -1.2  -0.4
> b[(Intercept) species_garden:Picea_abies_S*]            2.8    6.1  -5.0   2.7  10.4
> b[MAT_prov species_garden:Picea_abies_S*]               1.1    0.9   0.0   1.1   2.2
> b[(Intercept) species_garden:Picea_engelmannii_B]     -12.3    5.2 -19.0 -12.2  -5.7
> b[MAT_prov species_garden:Picea_engelmannii_B]          0.1    0.3  -0.4   0.1   0.5
> b[(Intercept) species_garden:Picea_mariana_I]          25.1    5.2  18.4  25.2  31.6
> b[MAT_prov species_garden:Picea_mariana_I]              0.6    0.6  -0.2   0.6   1.3
> b[(Intercept) species_garden:Picea_sitchensis_D]        4.0    5.9  -3.6   4.0  11.4
> b[MAT_prov species_garden:Picea_sitchensis_D]           0.3    0.4  -0.2   0.3   0.8
> b[(Intercept) species_garden:Pinus_albicaulis_C]      -20.9    5.4 -27.8 -20.8 -14.1
> b[MAT_prov species_garden:Pinus_albicaulis_C]           2.3    0.7   1.4   2.3   3.2
> b[(Intercept) species_garden:Pinus_ponderosa_J]        51.2    6.8  42.5  51.1  59.7
> b[MAT_prov species_garden:Pinus_ponderosa_J]           -0.1    0.5  -0.7  -0.1   0.5
> b[(Intercept) species_garden:Populus_balsamifera_F]     2.7    5.5  -4.3   2.7   9.7
> b[MAT_prov species_garden:Populus_balsamifera_F]       -0.3    0.6  -1.0  -0.3   0.4
> b[(Intercept) species_garden:Populus_trichocarpa_D]   -44.2    5.3 -51.0 -44.2 -37.4
> b[MAT_prov species_garden:Populus_trichocarpa_D]       -0.6    0.2  -0.9  -0.6  -0.3
> b[(Intercept) species_garden:Pseudotsuga_menziesii_H]  13.4    5.9   6.0  13.4  20.8
> b[MAT_prov species_garden:Pseudotsuga_menziesii_H]     -0.4    0.3  -0.8  -0.4   0.0
> b[(Intercept) species_garden:Quercus_petraea_U*]       -3.4    5.9 -11.0  -3.2   4.0
> b[MAT_prov species_garden:Quercus_petraea_U*]          -0.5    0.4  -1.0  -0.5  -0.1
> b[(Intercept) species_garden:Quercus_petraea_V*]       -1.5    5.8  -9.1  -1.6   5.9
> b[MAT_prov species_garden:Quercus_petraea_V*]          -0.8    0.4  -1.3  -0.8  -0.4
> b[(Intercept) species_garden:Tsuga_heterophylla_E]    -30.2    6.9 -38.8 -30.2 -21.5
> b[MAT_prov species_garden:Tsuga_heterophylla_E]         3.9    0.7   3.1   3.9   4.8
> b[(Intercept) species_garden:Tsuga_heterophylla_G]     17.9   12.2   2.5  17.6  33.8
> b[MAT_prov species_garden:Tsuga_heterophylla_G]        -0.8    1.3  -2.5  -0.8   0.8
> sigma                                                   4.9    0.1   4.8   4.9   5.1
> Sigma[species_garden:(Intercept),(Intercept)]         537.4  170.4 349.2 508.8 760.8
> Sigma[species_garden:MAT_prov,(Intercept)]             -9.7   10.0 -22.2  -8.7   1.6
> Sigma[species_garden:MAT_prov,MAT_prov]                 3.4    1.5   1.9   3.2   5.3

I think we should include the above for our model tables in the supplement. We don't need the other info (unless reviewers ask).

alinazeng commented 1 year ago

@lizzieinvancouver

I think we should include the above for our model tables in the supplement. We don't need the other info (unless reviewers ask).

Do I need to rename the first column or leaving it as is would be fine?

So for the above example, extract the posterior of each MAT_prov slope for each species, then divide by the number of species and use that posterior to calculate the mean and 10% and 90% quantiles. This is similar to what you did in issue https://github.com/lizzieinvancouver/localadaptclim/issues/16 ... but do all species since we're summarizing the full model.

does the below look correct to you? I used the posterior_interval() code

                                                               10%          90%
(Intercept)                                           108.65960039 121.68301469
b[(Intercept) species_garden:Alnus_rubra_A]           -46.33756308 -29.94556575
b[MAT_prov species_garden:Alnus_rubra_A]                0.22450184   1.51374102
b[(Intercept) species_garden:Betula_papyrifera_K]     -25.78773481 -10.26146234
b[MAT_prov species_garden:Betula_papyrifera_K]         -4.38011328  -2.47320549
b[(Intercept) species_garden:Betula_papyrifera_L]      -5.51734969   9.99250920
b[MAT_prov species_garden:Betula_papyrifera_L]         -3.04305379  -1.20312277
b[(Intercept) species_garden:Betula_papyrifera_M]      -1.99108091  13.12779721
b[MAT_prov species_garden:Betula_papyrifera_M]         -2.28904914  -0.49216861
b[(Intercept) species_garden:Fagus_sylvatica_R*]       -9.72365628   9.36060189
b[MAT_prov species_garden:Fagus_sylvatica_R*]          -1.41392466   0.38266225
b[(Intercept) species_garden:Fagus_sylvatica_T*]       -4.23224600  10.52696076
b[MAT_prov species_garden:Fagus_sylvatica_T*]           0.60962497   1.46435869
b[(Intercept) species_garden:Fraxinus_excelsior_Q*]    32.97972127  50.75911947
b[MAT_prov species_garden:Fraxinus_excelsior_Q*]       -1.90576091  -0.42452124
b[(Intercept) species_garden:Picea_abies_S*]           -4.98387445  10.42165002
b[MAT_prov species_garden:Picea_abies_S*]              -0.03404749   2.17788235
b[(Intercept) species_garden:Picea_engelmannii_B]     -19.00855980  -5.67483581
b[MAT_prov species_garden:Picea_engelmannii_B]         -0.38417879   0.49615796
b[(Intercept) species_garden:Picea_mariana_I]          18.37183368  31.62135308
b[MAT_prov species_garden:Picea_mariana_I]             -0.16715121   1.34155135
b[(Intercept) species_garden:Picea_sitchensis_D]       -3.56763043  11.40399726
b[MAT_prov species_garden:Picea_sitchensis_D]          -0.24728788   0.76744278
b[(Intercept) species_garden:Pinus_albicaulis_C]      -27.82496074 -14.10506772
b[MAT_prov species_garden:Pinus_albicaulis_C]           1.39587209   3.18532439
b[(Intercept) species_garden:Pinus_ponderosa_J]        42.53539477  59.68516359
b[MAT_prov species_garden:Pinus_ponderosa_J]           -0.72932300   0.45816103
b[(Intercept) species_garden:Populus_balsamifera_F]    -4.30241378   9.70927787
b[MAT_prov species_garden:Populus_balsamifera_F]       -1.00187063   0.41165546
b[(Intercept) species_garden:Populus_trichocarpa_D]   -51.04432212 -37.44416683
b[MAT_prov species_garden:Populus_trichocarpa_D]       -0.94897997  -0.33559242
b[(Intercept) species_garden:Pseudotsuga_menziesii_H]   6.03570741  20.78113223
b[MAT_prov species_garden:Pseudotsuga_menziesii_H]     -0.84951148  -0.04526081
b[(Intercept) species_garden:Quercus_petraea_U*]      -10.97329278   3.96477263
b[MAT_prov species_garden:Quercus_petraea_U*]          -1.01067725  -0.05979317
b[(Intercept) species_garden:Quercus_petraea_V*]       -9.05386154   5.94810094
b[MAT_prov species_garden:Quercus_petraea_V*]          -1.28889256  -0.35036696
b[(Intercept) species_garden:Tsuga_heterophylla_E]    -38.84578856 -21.52480433
b[MAT_prov species_garden:Tsuga_heterophylla_E]         3.11462470   4.79569273
b[(Intercept) species_garden:Tsuga_heterophylla_G]      2.53698491  33.78106314
b[MAT_prov species_garden:Tsuga_heterophylla_G]        -2.49960805   0.81169346
sigma                                                   4.75951935   5.12284244
Sigma[species_garden:(Intercept),(Intercept)]         349.21925741 760.79748493
Sigma[species_garden:MAT_prov,(Intercept)]            -22.22036486   1.57811210
Sigma[species_garden:MAT_prov,MAT_prov]                 1.88114370   5.29559177
lizzieinvancouver commented 1 year ago

does the below look correct to you? I used the posterior_interval() code

@alinazeng I don't think you posted the code you meant to .... if you did post the code you meant to, then, no I don't think it looks right as you need to extract the full posteriors for each species, sum them and then divide by the number of species, then calculate the interval on the one resulting posterior. Please check the code in issue #16.

lizzieinvancouver commented 1 year ago

Do I need to rename the first column or leaving it as is would be fine?

Good point! We should relabel that as 'parameter' for a header, and we actually only need the b[MAT_prov] estimates for each species so you could skip all the b[(Intercept) ones (but not critical if a pain).

Also, it occurred to me that we should add to the end of the first paragraph of the results something about the variability across species. Something like, 'Results were variable across species with some showing strong clines (e.g., Betula_papyrifera -- then give mean and 10% and 90% intervals for each garden).'

alinazeng commented 1 year ago

Hi Lizzie,

@lizzieinvancouver

ooo Thanks!! Yes I will follow issue #16

by "extract the full posteriors for each species, sum them and then divide by the number of species" you mean extract the full posteriors for each species x garden and divide by the number of species x garden, correct?

"then calculate the interval on the one resulting posterior" -> does this mean we would only want to report one line of interval per model rather than per species x garden? Since we sum and average things.

One more question: would increasing adapt_delta affect the estimates? if so I will wait till after we fixed the warnings of the models and then calculate the posteriors.

Cheers.

alinazeng commented 1 year ago

@lizzieinvancouver

Had more time than I thought this morning, so I tried out the following for the spring event ~ MAT mode. Hopefully this checks the second task of this Issue?


# following Issue # 16 to sum and average posteriors of one model

draws <- as.matrix(fitC_spring_mat)
all_spp_garden <- c("Alnus rubra A",
                "Betula papyrifera K",
                "Betula papyrifera L",
                "Betula papyrifera M","Fagus sylvatica R*","Fagus sylvatica T*",
                "Fraxinus excelsior Q*","Picea abies S*","Picea engelmannii B",
                "Picea mariana I","Picea sitchensis D","Pinus albicaulis C","Pinus ponderosa J",
                "Populus balsamifera F","Populus trichocarpa D","Pseudotsuga menziesii H",
                "Quercus petraea U*","Quercus petraea V*","Tsuga heterophylla E","Tsuga heterophylla G") 

#change spaces to underscores
library(stringr)
all_spp_garden <- str_replace_all(all_spp_garden," ","_")

# Get all the posteriors of each species x garden
all_spp_garden_draws <- matrix(data=NA,
                          nrow=nrow(draws),
                          ncol=length(all_spp_garden))
for (i in c(1:length(all_spp_garden))){
  all_spp_garden_draws[,i] <- draws[,paste("b[MAT_prov species_garden:", all_spp_garden[i], "]", sep="")]
}
# taking the average
all_spp_garden_post <- rowMeans(all_spp_garden_draws)

quantile(all_spp_garden_post, probs = c(.1, .5, .9))

10%         50%         90% 
-0.28666551 -0.10987028  0.07123643 
lizzieinvancouver commented 1 year ago

Had more time than I thought this morning, so I tried out the following for the spring event ~ MAT mode. Hopefully this checks the second task of this Issue?

Yes! This looks like what I was thinking. And yes to your previous queries (yes, I meant extract the full posteriors for each species x garden and divide by the number of species x garden, and yes, this means we only report one line of interval per model). Nice work!

lizzieinvancouver commented 1 year ago

would increasing adapt_delta affect the estimates? if so I will wait till after we fixed the warnings of the models and then calculate the posteriors.

It could (the reason we want to fix the divergences is that sometimes they can affect the model output, but it's not that common), but is unlikely to...

In general you want code that's easy to re-run, so aim to write code where it's easy to update the model, then update the estimates. Sweave helps a lot with this.

alinazeng commented 1 year ago

@lizzieinvancouver

Quick questions:

(1) is there a code to save posteriors so I don't need to run the codes everytime?

(2) to use \Sexpr{ }.... For example, I want to put in the median of the posterior of fitC_fall_lat by writing \Sexpr{format(median(evergreenspost_fitC_fall_lat),digits = 1)}

Where/how in the manuscript.rnw should I tell R sweave what evergreenspost_fitC_fall_lat is? or do I refer R to a separate script where I run to generate evergreenspost_fitC_fall_lat?

And what code should I use to tell the 90 UI of a posterior?

Thanks and I hope my questions make sense.

lizzieinvancouver commented 1 year ago

@alinazeng Good questions!

(1) You can save Stan objects (mystanmodel in the below) as RDS or Rdata files and read them in that way. That way you need run them very infrequently.

mystanmodel <- stan_lmer( ...)

(2) Usually I build -- as R objects -- all the things I was to use within Sexpr in an R file and then call the objects in Sweave so there's less R code that can go wrong. For example:

mypost <- rnorm(1000, 5, 3)
my90UI <- round(quantile(mypost, 0.9), 2)

Or something similar (be careful with the indexing) should work to give an 90% interval and then you'd use /Spexr{my90UI}

HTH!

alinazeng commented 1 year ago

@lizzieinvancouver

Usually I build -- as R objects -- all the things I was to use within Sexpr in an R file and then call the objects in Sweave so there's less R code that can go wrong.

Yes, I did build all R objects in another R script, which I open at the same time as the manuscript in R studio. However, when I compose pdfs, it would show an error: "Error: at manuscript.Rnw:221, object 'evergreenspost_fitC_fall_mat' not found Execution halted"

Do I need to link this R script where I build R objects to the .Rnw file somehow?

lizzieinvancouver commented 1 year ago

@alinazeng Yes, you need to source the R file that makes the object inside your Sweave file somewhere for Sweave to see it... if you're still having issues after you try that I would suggest showing it to Deirdre in person as it may be easier to troubleshoot.

alinazeng commented 1 year ago

Hi Lizzie, @lizzieinvancouver

Thanks! I have added in /Sexpr{} (for some, just as an experiment)

Could you take a look at this paragraph "In contrast, fall events (budset, leaf senescence, leaf abcission) advanced strongly with provenance latitude and mean annual temperature" (line 218) and see if the pdf compiles on your end?

I am not sure how to format the UI (Or if I calculated it correctly).

This relationship, however, was observed mostly in North America where fall events advance \Sexpr{format(abs(median(northamericanpost_fitC_fall_lat)),digits = 2)} (UI = \Sexpr{format(abs(quantile(northamericanpost_fitC_fall_lat, 0.9)), digits = 2)}) days per degree we move north.

Please let me know.

Cheers.

lizzieinvancouver commented 1 year ago

@alinazeng I tried to compile, but there are a number of absolute paths so this does not seem a good idea (and it does not compile because of this). If you want me to trouble-shoot Sweave code you'd need to make a small test.rnw file with the minimal code to recreate the error, then I can try to debug.

Your code for the 90% quantile looks fine to me, but you should check it makes sense compared to your plots and other information about the model.

alinazeng commented 1 year ago

@lizzieinvancouver

Hi Lizzie. Actually, I used relative paths in the .RNW so you should be able to compile if you have the necessary files saved on your local disk in the same fashion as I did, unless there's a better way to approach this?

On my end I could compile successfully and there's no error to troubleshoot.

Could you please check you have the following files/folders saved on your local disk?

localadaptclim/Analyses/source_fall_models_for_sweave.R # sourced just this one script localadaptclim/Output/model_fit/ # all .RDS in this /model_fit/ folder would be good to have localadaptclim/Docs

when we are in the ms.RNW, the working directory is the folder containing the .RNW, which is /Docs, therefore I set up the relative paths as ".../Analyses/source_fall_models_for_sweave.R", since /Analyses and /Docs are within the same parent directory /localadaptclim, and so is /Output

Let me know if this does not work out. Thank you.

lizzieinvancouver commented 1 year ago

@alinazeng I did try to compile it before I wrote you and it snagged on the concordence file. Sweave files can be hard to get to run from afar on two machines (especially with any differences in OS) which is annoying (I know), so I don't think me trying to compile your full file is a good approach.

I am happy to try to compile one short file you make that recreates the issue you're having.

alinazeng commented 1 year ago

@lizzieinvancouver Hi lizzie,

Thanks.

No, I had no issues compiling pdfs on my side, so all good. I just wanted to ask you last time if this is how you envision the text to look.

I put (UI = XX) in parentheses right after the mean. Does the format look okay? Or should I not say "UI ="?

image

lizzieinvancouver commented 1 year ago

@alinazeng Ah! Okay. If we're doing 90% intervals you need to request the 0.05 and 0.95 quantiles (that gives 90% total) and then report as something like: 'fall events advance 4.2 days (3.3-4.8) days per degree northward' -- I made up 3.3 and 4.8 but that would be the range of the 0.05 to 0.95 values. You do not need to say UI. We just need to add to the methods section that 'we report all values as posteriors means (I saw you were using medians ... I would use mean) with 90% uncertainty intervals given parenthetically.

alinazeng commented 1 year ago

Hi Lizzie, Thank you very much for the clarification. Something unforeseen has come up and is keeping me occupied. I will try my best to find time to work on this in a week's time but there might be a delay.

Cheers,Alina On Wednesday, March 15, 2023 at 03:49:38 AM PDT, lizzieinvancouver @.***> wrote:

@alinazeng Ah! Okay. If we're doing 90% intervals you need to request the 0.05 and 0.95 quantiles (that gives 90% total) and then report as something like: 'fall events advance 4.2 days (3.3-4.8) days per degree northward' -- I made up 3.3 and 4.8 but that would be the range of the 0.05 to 0.95 values. You do not need to say UI. We just need to add to the methods section that 'we report all values as posteriors means (I saw you were using medians ... I would use mean) with 90% uncertainty intervals given parenthetically.

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: @.***>

lizzieinvancouver commented 1 year ago

I will try my best to find time to work on this in a week's time but there might be a delay.

@alinazeng No worries! Happy to pick this up as your schedule allows. (Of course, let me know though if there is anything I can do to help.)

alinazeng commented 1 year ago

@lizzieinvancouver Hi Lizzie,

Thank you!

I added /Sexpr where we needed them. Since the file paths can be messy to link on your end, you can take a look at the ms.pdf I saved in /Docs. Let me know what needs to be fixed and whether you think it's better to display "4.2 days (3.3-4.8) " or "4.2 +-0.75 days" (made up numbers here)

I know the figures in the ms right now are messy. I will work on them once we are sure where in the ms to include what figures. I will be working on the tables in the supp soon too.

(p.s. I got into Oxford! just received offers from them two days ago. Thank you so much for the recommendation letter!!)

Cheers, Alina

lizzieinvancouver commented 1 year ago

@alinazeng Congrats! That is wonderful news about Oxford. I am really thrilled for you and want to see pictures and hear how it is once you arrive. Sorry I went MIA, traveling this week so have been a mix of busy and then spotty internet on trains.

Nice work on /Sxpr! Since posteriors can be slightly non-normal we don't us the +/- (that's more for Fisherian stats, not Bayesian) and we avoid = since the values are never exact. So, where you have "mostly in North America where fall events advance 4.2 (3.9-4.5) days per degree northward" is great! But you should change "spring events were not related to provenance latitude (mean slope = 0.1 ± 0.16)" to something more like, "spring events were not related to provenance latitude (0.1 [-0.01 - 0.26] days per degree northward)."

If you're uncertainty intervals overlap zero I would also adjust the phrasing to acknowledge this. So:

Spring events advanced slightly with provenance latitude (mean slope = 0.1 ± 0.33) and MAT (mean slope = -0.2 ± 0.39) in Europe, meaning spring events are slightly earlier where provenance MAT is lower (i.e., higher, more northern latitudes)(Fig. 1 & Supplement Table. 1).

Would be:

Spring events advanced slightly with provenance latitude (mean slope = 0.1 ± 0.33) and MAT (mean slope = -0.2 ± 0.39) in Europe, meaning spring events are slightly earlier where provenance MAT is lower (i.e., higher, more northern latitudes)(Fig. 1 & Supplement Table. 1), though uncertainty intervals strongly overlapped zero.

alinazeng commented 1 year ago

@lizzieinvancouver Hi Lizzie,

Thank you! Are you still in Europe? When are you coming back to Canada?

I have updated the expressions as per your guidance: image

At the moment I am fighting the flu and juggling Capstone projects, so turnarounds might be slower than usual.

lizzieinvancouver commented 1 year ago

@alinazeng

Feel better! And good luck on Capstones.

That looks good. Two things... (1) Can you add units and white space around the dash? So the first one would be (0.1 days/degree [-0.06 - 0.26]). (2) Keep the decimals the same -- so either two or one decimal places. Two is probably best (even though we really don't have accuracy to that level IMHO).

I am not back until early August. I am sorry as I think I may miss seeing you again in person in Vancouver, but I am very hopeful our paths will cross again.

alinazeng commented 1 year ago

@lizzieinvancouver

Hi Lizzie.

Thank you for your swift response and good thoughts! I have updated the numbers. Is there anything else you'd like me to add to the Results section? At the moment I am putting down ideas for Introduction in another document. Do you have an outline in mind and topics I need to cover?

I might still be around till the end of August. If not, I am hopeful that I will get to see you after my Masters!

lizzieinvancouver commented 1 year ago

@alinazeng

Okay, I should check the full results again soon so we can aim to finalize that section. Have you finalized figures mostly?

For the intro, you probably should open with the broad relevance of plasticity versus adaptation of phenology and climate change. I would check out this paper by Chevin "Adaptation, Plasticity, and Extinction in a Changing Environment: Towards a Predictive Theory" and look at the list of papers that have cited this since then and a review a couple of those. You might also check out this older on (Range Shifts and Adaptive Responses to Quaternary Climate Change). Then you should cite the review papers on this already, as well as give examples of some specific studies finding a trend and not finding a trend. You probably need to establish WHICH way we expect trends to go with latitude.

There's a couple good papers the lab published recently to look at what we said about latitude the refs we used when discussing latitude: Ettinger et al. 2020 (Nature Climate Change) Ettinger et al. 2021 (New Phyt) and this one.

Eventually we need to review the co-gradient and counter-gradient literature on phenology ... but I am not sure if we need that in the intro.

That's great that you might still be around! I hope to be back around the end of the first week of August (but have not booked tickets).

alinazeng commented 1 year ago

@lizzieinvancouver Thank you Lizzie. I was busy with the Capstone presentation, which is now done. I have two more finals and then farewell undergrad!

Thanks for these suggestions! They are great, and I will incorporate them in the introduction. I expect to send you the draft introduction in one week.

The figures in the MS are updated and they reflect our latest data and models, so you can refer to them while reading the results. I haven't edited font size, and point size, etc. for better readability for publication, and might have missed moving a couple plots into the supp (I took note regarding which ones you wanted me to move to supp, so no need to reiterate. I will move them after I finish the intro and we finalize the result.

Cheers

lizzieinvancouver commented 1 year ago

@alinazeng I just pushed some edits. Look for %EMWApr20. There's no rush but once you have updated the figures as requested -- including with captions, please ask me to review and then when you have also updated the text I can review the full results.

Good luck with all your finals and wrapping up!

alinazeng commented 1 year ago

@lizzieinvancouver

Hi Lizzie,

Thank you so much for the feedback. I need a few days to update the figures and texts for Result and Supplement (am in the middle of moving so things are a bit hectic). I expect to have everything updated in a week.

In the meantime, could you take a look at my rough draft for the Introduction here? I am sharing it to you via google drive because I am yet to update the results in the ms.RNW and don't want to mess up the push and pull.

Please excuse the messy intext-citations. I will clean them up when I put the Intro into the ms.RNW. Let me know what you'd like me to add and fix.

Cheers

alinazeng commented 1 year ago

@lizzieinvancouver Hi lizzie. please excuse the delay. I have uploaded the results and the main figures (all critiques are welcome). reorganized the Supp too (there are several figures dangling at the end of the supp cuz I am unsure if we want them, but the ones we definitely want to include are at the front, and tables).

I also added the draft introduction to the ms if you could take some time to suggest further edits. ms supp

Thank you so much and look forward to hear what you think.

lizzieinvancouver commented 1 year ago

@alinazeng Nice! I am in the field this week so will aim to look at this after most likely.