forecasting figure in phenophylo?

[lizzieinvancouver]

@dbuona mentioned that one way to make the ms stronger might be to add some forecasting to the paper. 'Show me the forecasting!' he suggested. Jonathan liked this idea and quickly suggested that we ...

1. Grab all the GBIF points for every species in the model, then calculate chilling, forcing and photoperiod for each location in the past and under a warming scenario and forecast for both the phylo and lambda=0 model. I was horrified as this seems like a TON of work to do, mostly because chilling is SO painful.
2. One alternative could be similar but do a forcing only model and use mean monthly temperatures.
3. Dan suggest we forecast for every PEP725 location (or John O'Keefe HF data), I was again not into this as I don't think our model is designed to do that sort of forecasting ... it's a first big step.
4. I suggested perhaps a more conceptual + real data figure using a few species that have little data (and so would be pooled strongly in the lambda=0 model) and some well sampled and show the mu from the lambda=0 model under a 'now' and '+2 scenario' ... though we still need to get the data and decide how to do it ....
5. We could do something without a warming scenario using just the species where we already have chilling from the OSPREE ranges?
6. Riffing on Dan's O'Keefe idea we could just fit DOY ~ spring temp for the O'Keefe data ... different data and model predictors but perhaps shows what we want? We run the PMM and the lambda=0 model and then do compare versus observed to see which fits better (imagine the HF species as a phylogeny ... and then show on each tip the raw data species mean versus estimate from the lambda=0 and the PMM model for that species).

I am aiming to try 6 in the next couple of days and see if it looks hopeful. It's not critical to have this figure in the paper, but may help pull it together? We definitely don't want to slow progress down as getting this ms wrapped up soon would be WONDERFUL.

@MoralesCastilla Let us know what you think.

[lizzieinvancouver]

@MoralesCastilla Just looking at the output this morning and it lines up pretty well with the what we'd expect for non-z scored predictors from our previous OSPREE work:

> summary(fitsum)$summary
> fitsum
                         mean     se_mean        sd         2.5%         25%        50%
a_z                54.4918680 0.078407741 5.4014429 43.710714930 50.93773566 54.6195510
sigma_interceptsa  19.9266289 0.028208155 1.8075433 16.723107529 18.66767717 19.7934645
b_zf               -1.0362409 0.004809582 0.2878906 -1.584971186 -1.22419262 -1.0449172
sigma_interceptsbf  0.7386676 0.003583073 0.1279565  0.519193922  0.64877186  0.7288803
lam_interceptsbf    0.7652868 0.007855429 0.1797287  0.331135629  0.65641293  0.8050351
b_zc               -2.3587098 0.008153032 0.5615012 -3.441373167 -2.73418043 -2.3714403
sigma_interceptsbc  1.6333045 0.004897133 0.2153155  1.253219260  1.48477142  1.6173415
lam_interceptsbc    0.6349412 0.004725026 0.1665610  0.272974382  0.52310672  0.6455514
b_zp               -0.2066890 0.002532047 0.1309802 -0.465674067 -0.29123893 -0.2069776
sigma_interceptsbp  0.4455723 0.002910595 0.0772849  0.302027030  0.39335716  0.4419757
lam_interceptsbp    0.2717079 0.008439277 0.2166749  0.008975798  0.09317178  0.2192407
sigma_y            12.7329327 0.001933436 0.1788492 12.388030678 12.61245445 12.7283894
                          75%       97.5%     n_eff      Rhat
a_z                58.1442041 64.63319313 4745.7155 1.0006533
sigma_interceptsa  21.0322524 23.78187000 4106.0858 1.0007966
b_zf               -0.8568560 -0.44795642 3582.9466 1.0011502
sigma_interceptsbf  0.8174979  1.01800814 1275.3024 1.0025560
lam_interceptsbf    0.9086135  0.99074085  523.4737 1.0141756
b_zc               -1.9950135 -1.22758008 4743.1088 1.0002117
sigma_interceptsbc  1.7678268  2.09301271 1933.1553 1.0002740
lam_interceptsbc    0.7596506  0.91777167 1242.6201 1.0009346
b_zp               -0.1252764  0.05423467 2675.8862 1.0009594
sigma_interceptsbp  0.4945935  0.60578574  705.0593 1.0084032
lam_interceptsbp    0.4066513  0.78333721  659.1840 1.0071314
sigma_y            12.8508395 13.09521046 8556.8590 0.9998985

The prior on a is a little low so I will push updated code and models later today (I hope), but I doubt it will change anything as we have a lot of data. I only did a quick check so someday you should check more closely that nothing is too close to the prior.... if we end up using these models. No need to do it before then. Whatsapp me if you need help anytime.

[lizzieinvancouver]

Also, this is what I was looking at for comparison values ...

[lizzieinvancouver]

OLD UPDATE that I should have included here ... on Harvard Forest...

In early January Dan, Jonathan and I discussed somehow adding more forecasting to the PMM paper. Our best idea (we thought) was to build a model of the Harvard Forest long-term phenology data (collected by John O'Keefe), and compare predictions with PMM versus no phylogeny. With Deirdre's excellent help I got the PMM running and it runs well (the lambda is really high -- 0.9 or so for slope; for intercept it's 0.5) but I struggled on the regular model. Eventually I just used rstanarm and moved along, but I realized I am not sure what is a good comparison here (We have no 'true' value to compare to ...) so I am giving up. Also, the estimates are pretty similar (see https://github.com/lizzieinvancouver/ospree/blob/master/analyses/phylogeny/okeefePMM/compareslopes_HFPMM.pdf -- black line is 1:1 and blue line is mean across regular model slopes). I left all my work in the repo, in case it is useful.

[lizzieinvancouver]

@MoralesCastilla Okay, I ran the models with slightly tweaked priors and it looks like a was being held back a little, but everything else basically unchanged in the PMM model:

> fitsum
                         mean     se_mean         sd         2.5%         25%        50%
a_z                62.4397212 0.068383473 5.02823518 52.486132260 59.05058807 62.4026158
sigma_interceptsa  19.4449613 0.025787697 1.72169675 16.460308845 18.25208781 19.3222967
b_zf               -1.1151750 0.005028246 0.27608431 -1.648296374 -1.29630518 -1.1182845
sigma_interceptsbf  0.7196323 0.003354999 0.12344080  0.506734154  0.63082021  0.7102837
lam_interceptsbf    0.7376357 0.008150304 0.19789317  0.260238120  0.61865925  0.7770631
b_zc               -2.4388003 0.006654259 0.55048063 -3.508540631 -2.80833486 -2.4426954
sigma_interceptsbc  1.6148279 0.004070717 0.20738649  1.250891134  1.46906480  1.5998184
lam_interceptsbc    0.6155167 0.005424936 0.17190489  0.248357480  0.50285310  0.6286604
b_zp               -0.2189654 0.002528910 0.12958238 -0.483621148 -0.29851189 -0.2192625
sigma_interceptsbp  0.4440849 0.003147445 0.07587849  0.302366244  0.39281014  0.4405405
lam_interceptsbp    0.2535540 0.008264571 0.20692461  0.008324172  0.08897465  0.1990209
sigma_y            12.7350584 0.001792660 0.18041250 12.386719439 12.61153145 12.7332077
                          75%       97.5%      n_eff      Rhat
a_z                65.8471731 72.31010696  5406.6574 1.0005235
sigma_interceptsa  20.5156880 23.11342686  4457.4676 1.0000638
b_zf               -0.9412298 -0.55078760  3014.7438 1.0020230
sigma_interceptsbf  0.7980943  0.98527555  1353.7318 1.0015636
lam_interceptsbf    0.8958557  0.99178927   589.5409 1.0079618
b_zc               -2.0817382 -1.34881518  6843.6006 0.9997975
sigma_interceptsbc  1.7465037  2.05937668  2595.4881 1.0023419
lam_interceptsbc    0.7422882  0.91597968  1004.1236 1.0044686
b_zp               -0.1373734  0.03240224  2625.5795 1.0002217
sigma_interceptsbp  0.4927916  0.59801077   581.1937 1.0047886
lam_interceptsbp    0.3668395  0.76520012   626.8789 1.0036336
sigma_y            12.8579271 13.08631893 10128.3238 0.9997123

I have not looked at the lambda=0 output, but I ran both and they're now updated on the drive. I will push changes on git now.

[MoralesCastilla]

@lizzieinvancouver

Here's a draft of the forecasting figure