try our model with early and late burst parameters ...

lizzieinvancouver commented 2 years ago

Started this in commit f997e3acfe3f6a8bf33168b40c9d44d29c229a78 ... and here's the output with that code (late burst on FORCE, early burst on PHOTO):

> summary(testme)$summary[stan.param, c("mean", "2.5%", "25%", "75%", "97.5%")]
                         mean       2.5%        25%        75%      97.5%
a_z                27.6951331 22.0462882 25.7502902 29.6261824 33.5134532
lam_interceptsa     0.9397507  0.8710453  0.9242030  0.9599822  0.9794277
sigma_interceptsa  10.2218437  8.8453391  9.6919297 10.7118913 11.7902140
b_zf                2.4472631 -6.6214711 -0.5585790  5.5477252 11.4018477
lam_interceptsbf    0.8940128  0.7938451  0.8710322  0.9251940  0.9571256
sigma_interceptsbf 13.3056983 11.4626112 12.5730236 13.9537248 15.6357609
b_zc                4.6628782 -2.7014232  2.1705943  7.2143950 11.9426805
lam_interceptsbc    0.8848155  0.7539248  0.8560962  0.9233084  0.9611836
sigma_interceptsbc 11.0329314  9.3965430 10.3970163 11.6127176 12.9502789
b_zp                2.4455405 -3.5528481  0.4482955  4.5096052  8.1738182
lam_interceptsbp    0.9621981  0.9191936  0.9529100  0.9753914  0.9878321
sigma_interceptsbp  8.5642518  7.4698526  8.1331164  8.9512171  9.8495968
sigma_y             0.5307887  0.5028968  0.5207836  0.5403287  0.5605405
> t(param)
     a_z lam_interceptsa sigma_interceptsa b_zf delta_interceptsbf sigma_interceptsbf b_zc
[1,] 30  1               1                 5    2                  1                  5   
     lam_interceptsbc sigma_interceptsbc b_zp delta_interceptsbp sigma_interceptsbp sigma_y
[1,] 1

lizzieinvancouver commented 2 years ago

Setting the two deltas to 1 yields more similar sigmas (as predicted), but we're not sure why the sigmas are way off ... could be the seed for what we picked? .... err, later on we tried a simulated tree and we think this insane high sigma is because of our trees polytomies. So? So, interpret the sigmas relatively -- do not interpret absolute sigmas.


> # Compare to true values
> summary(testme)$summary[stan.param, c("mean", "2.5%", "25%", "75%", "97.5%")]
                         mean       2.5%        25%        75%      97.5%
a_z                27.7104189 22.1013351 25.7679418 29.6813925 33.2954755
lam_interceptsa     0.9390172  0.8705272  0.9239161  0.9593582  0.9795226
sigma_interceptsa  10.2278795  8.8930707  9.7131538 10.7042667 11.7572576
b_zf                2.9120608 -4.1962589  0.4205901  5.3754163 10.0868883
lam_interceptsbf    0.9295525  0.8606502  0.9125605  0.9515868  0.9740501
sigma_interceptsbf 11.1345690  9.6213810 10.5494977 11.6732470 12.9224146
b_zc                4.6380239 -2.9644258  2.2326259  7.1109585 12.0556003
lam_interceptsbc    0.8869743  0.7633047  0.8590656  0.9244432  0.9607681
sigma_interceptsbc 11.0417058  9.4137503 10.4105732 11.6213860 12.9757376
b_zp                3.9325644 -3.4471438  1.2900425  6.5177994 11.3374772
lam_interceptsbp    0.9512807  0.8947846  0.9390435  0.9677564  0.9836953
sigma_interceptsbp 10.8223000  9.4420828 10.2826130 11.3138922 12.4970926
sigma_y             0.5317054  0.5037541  0.5215293  0.5415428  0.5611063
> t(param)
     a_z lam_interceptsa sigma_interceptsa b_zf delta_interceptsbf sigma_interceptsbf b_zc
[1,] 30  1               1                 5    1                  1                  5   
     lam_interceptsbc sigma_interceptsbc b_zp delta_interceptsbp sigma_interceptsbp sigma_y
[1,] 1                1                  5    1                  1                  0.5

lizzieinvancouver commented 2 years ago

So, take home from today is that it does seem like sigma goes up in our model when we generate a late-burst model of evolution, compared to when we generate an early-burst ... but it's still a little wonky. Tomorrow we should:

[x] Re-run this code on a larger tree ... and maybe better on a simulated tree of 200 spp. with deltas as early/late
[x] Re-run this code on a larger tree ... and maybe better on a simulated tree of 200 spp. with deltas as 1

lizzieinvancouver commented 2 years ago

Here's the output with 100 spp on a simulated tree with early and late burst parameter (commit 7d8447f33812f259ecbbdf49d822a98a68303cf1):

> summary(testme)$summary[stan.param, c("mean", "2.5%", "25%", "75%", "97.5%")]
                         mean       2.5%        25%        75%      97.5%
a_z                30.1225791 29.4832299 29.8997455 30.3405035 30.7627357
lam_interceptsa     0.8568671  0.7142513  0.8238300  0.9003062  0.9490959
sigma_interceptsa   0.7991273  0.6699118  0.7485277  0.8453812  0.9545045
b_zf                5.3679132  4.4538769  5.0670506  5.6734923  6.2744453
lam_interceptsbf    0.8576148  0.7008668  0.8163461  0.9092798  0.9599088
sigma_interceptsbf  1.1578423  0.9760485  1.0858617  1.2221945  1.3756715
b_zc                4.8960033  4.2100518  4.6649864  5.1285638  5.5804659
lam_interceptsbc    0.8769102  0.7322242  0.8443720  0.9206030  0.9619812
sigma_interceptsbc  0.8512979  0.7238783  0.8016513  0.8984148  0.9957862
b_zp                5.3386106  4.8620858  5.1786941  5.4986526  5.8229542
lam_interceptsbp    0.8330571  0.6688969  0.7910936  0.8865465  0.9431528
sigma_interceptsbp  0.6467007  0.5423435  0.6051488  0.6838663  0.7755461
sigma_y             0.4997369  0.4729627  0.4897495  0.5095549  0.5285383
> t(param)
     a_z lam_interceptsa sigma_interceptsa b_zf delta_interceptsbf sigma_interceptsbf b_zc
[1,] 30  1               1                 5    2                  1                  5   
     lam_interceptsbc sigma_interceptsbc b_zp delta_interceptsbp sigma_interceptsbp sigma_y
[1,] 1                1                  5    0.5                1                  0.5    
>

lizzieinvancouver commented 2 years ago

Delta=2 for forcing, delta=0.5 for photoperiod with 200 species ... input:

if(!useOSPREEtree){
nspecies = 200
spetree <- pbtree(n=nspecies, nsim=1, b=1, complete=FALSE,scale=1)
}

spetree$tip.label <- paste("s", 1:nspecies, sep="")
nind = 10

# Now set up the trait parameters
# But let's making photoperiod (delta=0.5) as early-burst and forcing as 'late-burst' (delta=2)
# ... and make the other lambda values 1
param <- list(a_z = 30, # root value intercept
              lam_interceptsa = 1, # lambda intercept
              sigma_interceptsa = 1, # rate of evolution intercept
              b_zf = 5, # root value trait1
              delta_interceptsbf = 2, # DELTA trait1
              sigma_interceptsbf = 1, # rate of evolution trait1
              b_zc = 5, # root value trait2
              lam_interceptsbc = 1, # lambda trait2
              sigma_interceptsbc = 1, # rate of evolution trait2
              b_zp = 5, # root value trait3
              delta_interceptsbp = 0.5, # DELTA trait3
              sigma_interceptsbp = 1, # rate of evolution trait3
              sigma_y = 0.5 # overall sigma, is much higher in real data (~13)
              )

output:

> summary(testme)$summary[stan.param, c("mean", "2.5%", "25%", "75%", "97.5%")]
                         mean       2.5%        25%        75%      97.5%
a_z                30.1844053 29.4622663 29.9385782 30.4314817 30.8878616
lam_interceptsa     0.9495341  0.8969649  0.9374892  0.9657216  0.9820712
sigma_interceptsa   0.9603797  0.8499757  0.9193581  0.9978730  1.0840033
b_zf                5.2192975  4.4119651  4.9406629  5.4957518  6.0434671
lam_interceptsbf    0.9360963  0.8693480  0.9212255  0.9561286  0.9769717
sigma_interceptsbf  1.0609016  0.9338990  1.0143295  1.1057135  1.1956891
b_zc                4.6329808  3.9739297  4.4114667  4.8498414  5.2798110
lam_interceptsbc    0.9616244  0.9158482  0.9516720  0.9751486  0.9876680
sigma_interceptsbc  0.8601785  0.7653740  0.8231191  0.8946711  0.9656902
b_zp                4.4645270  3.8609830  4.2650992  4.6687802  5.0588140
lam_interceptsbp    0.9791349  0.9558645  0.9739592  0.9860182  0.9929168
sigma_interceptsbp  0.7540692  0.6803173  0.7263490  0.7800521  0.8388815
sigma_y             0.5083600  0.4887897  0.5013336  0.5152053  0.5287668
> t(param)
     a_z lam_interceptsa sigma_interceptsa b_zf delta_interceptsbf sigma_interceptsbf b_zc
[1,] 30  1               1                 5    2                  1                  5   
     lam_interceptsbc sigma_interceptsbc b_zp delta_interceptsbp sigma_interceptsbp sigma_y
[1,] 1                1                  5    0.5                1                  0.5

lizzieinvancouver commented 2 years ago

Here's the output from the two deltas set to 1....

> summary(testme)$summary[stan.param, c("mean", "2.5%", "25%", "75%", "97.5%")]
                         mean       2.5%        25%        75%      97.5%
a_z                30.1898759 29.4678953 29.9407168 30.4432525 30.8974613
lam_interceptsa     0.9492131  0.8978151  0.9368854  0.9650256  0.9821730
sigma_interceptsa   0.9599533  0.8520185  0.9176748  0.9992404  1.0843214
b_zf                5.0604587  4.4137636  4.8384969  5.2841681  5.7030853
lam_interceptsbf    0.9565439  0.9081314  0.9463136  0.9704030  0.9843410
sigma_interceptsbf  0.8555322  0.7615738  0.8203219  0.8878487  0.9648854
b_zc                4.6383152  3.9710540  4.4085760  4.8633481  5.3077488
lam_interceptsbc    0.9619735  0.9180933  0.9528188  0.9749864  0.9875800
sigma_interceptsbc  0.8606549  0.7659279  0.8257314  0.8949230  0.9627126
b_zp                4.6692980  3.9643547  4.4170912  4.9212168  5.4002904
lam_interceptsbp    0.9743676  0.9446363  0.9681480  0.9830964  0.9914985
sigma_interceptsbp  0.9543205  0.8582755  0.9175998  0.9894868  1.0613455
sigma_y             0.5079991  0.4881305  0.5007217  0.5150257  0.5287013
> t(param)
     a_z lam_interceptsa sigma_interceptsa b_zf delta_interceptsbf sigma_interceptsbf b_zc
[1,] 30  1               1                 5    1                  1                  5   
     lam_interceptsbc sigma_interceptsbc b_zp delta_interceptsbp sigma_interceptsbp sigma_y
[1,] 1                1                  5    1                  1                  0.5

And the Geiger output suggests we pulled a low sigma on the bf:

> fitContinuous(spetree, intercepts, model="lambda")
GEIGER-fitted comparative model of continuous data
 fitted ‘lambda’ model parameters:
    lambda = 0.996476
    sigsq = 1.026503
    z0 = 30.187936

 model summary:
    log-likelihood = -146.251277
    AIC = 298.502555
    AICc = 298.625004
    free parameters = 3

Convergence diagnostics:
    optimization iterations = 100
    failed iterations = 0
    number of iterations with same best fit = 12
    frequency of best fit = 0.12

 object summary:
    'lik' -- likelihood function
    'bnd' -- bounds for likelihood search
    'res' -- optimization iteration summary
    'opt' -- maximum likelihood parameter estimates
> fitContinuous(spetree, slopes_bf, model="lambda")
GEIGER-fitted comparative model of continuous data
 fitted ‘lambda’ model parameters:
    lambda = 0.999434
    sigsq = 0.894939
    z0 = 5.070060

 model summary:
    log-likelihood = -128.680724
    AIC = 263.361447
    AICc = 263.483896
    free parameters = 3

Convergence diagnostics:
    optimization iterations = 100
    failed iterations = 0
    number of iterations with same best fit = 1
    frequency of best fit = 0.01

 object summary:
    'lik' -- likelihood function
    'bnd' -- bounds for likelihood search
    'res' -- optimization iteration summary
    'opt' -- maximum likelihood parameter estimates
> fitContinuous(spetree, slopes_bc, model="lambda")
GEIGER-fitted comparative model of continuous data
 fitted ‘lambda’ model parameters:
    lambda = 0.999825
    sigsq = 0.873764
    z0 = 4.634840

 model summary:
    log-likelihood = -125.675042
    AIC = 257.350083
    AICc = 257.472532
    free parameters = 3

Convergence diagnostics:
    optimization iterations = 100
    failed iterations = 0
    number of iterations with same best fit = 1
    frequency of best fit = 0.01

 object summary:
    'lik' -- likelihood function
    'bnd' -- bounds for likelihood search
    'res' -- optimization iteration summary
    'opt' -- maximum likelihood parameter estimates
> fitContinuous(spetree, slopes_bp, model="lambda")
GEIGER-fitted comparative model of continuous data
 fitted ‘lambda’ model parameters:
    lambda = 1.000000
    sigsq = 1.002001
    z0 = 4.652357

 model summary:
    log-likelihood = -139.085129
    AIC = 284.170257
    AICc = 284.292706
    free parameters = 3

lizzieinvancouver commented 2 years ago

Meeting today the next step is to build a figure like Fig 3b -- showing sigmas for forcing, chilling, photoperiod.... use a pb (simulated tree) with forcing at late burst, photo is early burst, chilling is 1 ... make sigmas higher (5) and set other traits similar to what we already in sims.

[x] @lizzieinvancouver does the code
[ ] @MoralesCastilla does the figure!

lizzieinvancouver commented 2 years ago

@MoralesCastilla I did my tasks and put in some holder figures ... see commit 30a29ec471c10a50059ebb9bac16c2891feeedee

lizzieinvancouver commented 1 year ago

I moved the relevant files over to ospree repo in commit # a157206bb032059d50a90e171e784056b198efa7. Closing this as we'll finalize the figure there.

lizzieinvancouver / pmm

try our model with early and late burst parameters ... #22