problems returning negative slopes

lizzieinvancouver commented 2 years ago

I ran into this issue about 6 weeks ago when trying to check we can recover the parameters in @MoralesCastilla OSPREE results, which TJ Davies asked about. I was editing a version of simulate_fit_uber_threeslopeintercept.R and found issues with negative slopes. @FaithJones and @willpearse asked if it was a priors issue.

Here's where I have gotten:

My original code (ubermini_2.R) has no problem with negative slopes.
simulate_fit_uber_threeslopeintercept.R seems to do better with 1-2 negative slopes than all three negative and I think perhaps better when you set the slope sigmas to 20 instead of 10, but they were centered at 0.
I think there are some similar issues here in getting back the slopes, lamba and phylogeny sigma as when you ask for the mu, sigma in a partial pooling model: that is, if you have a relatively small number of species you may not randomly draw species that match well to your original parameters. So it's hard to test this fully without doing more species, and this model scales REALLY poorly.
I run fitContinuous(spetree, intercepts, model="lambda") to try to check what the lamba might be and often it is more similar to the model output, but I am not sure how good to feel about using different model output to check my model output.

Below is the general code I have been running:

nspecies = 40
nind = 10

# Now set up the trait parameters
param <- list(a_z = 30, # root value intercept
              lam_interceptsa = 0.4, # lambda intercept
              sigma_interceptsa = 16, # rate of evolution intercept
              b_zf = -6, # root value trait1
              lam_interceptsbf = 0.7, # lambda trait1
              sigma_interceptsbf = 6, # rate of evolution trait1
              b_zc = -7, # root value trait2
              lam_interceptsbc = 0.5, # lambda trait2
              sigma_interceptsbc = 7, # rate of evolution trait2
              b_zp = -1.4, # root value trait3
              lam_interceptsbp = 0.3, # lambda trait3
              sigma_interceptsbp = 2, # rate of evolution trait3
              sigma_y = 13 # overall sigma
              )
# Set priors
phypriors <- list(
    a_z_prior_mu = 30, 
    a_z_prior_sigma = 10,
    lam_interceptsa_prior_alpha = 6, 
    lam_interceptsa_prior_beta = 6,  
    sigma_interceptsa_prior_mu = 10, 
    sigma_interceptsa_prior_sigma = 10,
    b_zf_prior_mu = 0, 
    b_zf_prior_sigma = 20, # was 10
    lam_interceptsbf_prior_alpha = 6, 
    lam_interceptsbf_prior_beta = 6,  
    sigma_interceptsbf_prior_mu = 10, 
    sigma_interceptsbf_prior_sigma = 10,
    b_zc_prior_mu = 0, 
    b_zc_prior_sigma = 20,
    lam_interceptsbc_prior_alpha = 6, 
    lam_interceptsbc_prior_beta = 6,  
    sigma_interceptsbc_prior_mu = 10,
    sigma_interceptsbc_prior_sigma = 10,
    b_zp_prior_mu = 0,
    b_zp_prior_sigma = 20,
    lam_interceptsbp_prior_alpha = 6, 
    lam_interceptsbp_prior_beta = 6,  
    sigma_interceptsbp_prior_mu = 10, 
    sigma_interceptsbp_prior_sigma = 10,
    sigma_y_mu_prior = 10, 
    sigma_y_mu_sigma = 10
)

The output is not bad (run yourself to check, it does not take long) but the 'bz_f' is around -3, instead of -6 at the mean.

A few later and adding the output to show (this is using the above priors):

> summary(testme)$summary[names(param), c("2.5%", "25%", "mean", "75%", "97.5%")]
                         2.5%        25%       mean        75%       97.5%
a_z                19.3528176 25.3524746 28.5706692 31.6546694 38.21952816
lam_interceptsa     0.2177077  0.3709254  0.4573911  0.5450373  0.70031401
sigma_interceptsa  11.2528407 13.4342992 14.9789751 16.2963273 19.67653284
b_zf               -6.9771484 -4.6332610 -3.5498141 -2.4537150  0.01957988
lam_interceptsbf    0.1990238  0.3482895  0.4337008  0.5162417  0.68039231
sigma_interceptsbf  3.7275138  4.3569528  4.8567989  5.2678574  6.34250309
b_zc               -5.3317917 -2.6123743 -1.0504175  0.4447112  3.44478011
lam_interceptsbc    0.2423342  0.4053982  0.4967730  0.5899430  0.74072775
sigma_interceptsbc  4.5742528  5.3599002  5.8998710  6.3551530  7.69787851
b_zp               -3.4470700 -2.1792929 -1.5644058 -0.9672777  0.33200885
lam_interceptsbp    0.2387825  0.4007225  0.4906196  0.5811677  0.73848148
sigma_interceptsbp  1.7179999  2.1325269  2.4139126  2.6724261  3.33259666
sigma_y            12.9084348 13.6967776 14.1266715 14.5447753 15.46559918
> t(param)
     a_z lam_interceptsa sigma_interceptsa b_zf lam_interceptsbf sigma_interceptsbf b_zc lam_interceptsbc
[1,] 30  0.4             16                -6   0.7              6                  -7   0.5             
     sigma_interceptsbc b_zp lam_interceptsbp sigma_interceptsbp sigma_y
[1,] 7                  -1.4 0.3              2                  13

lizzieinvancouver commented 2 years ago

[x] Right now I am running the above code but with 180 species and 20 ind each which I think is more similar to Nacho's data (did I mention that my model is slow? It is .... )
[ ] Next step is to run the code using Nacho's priors.
[ ] Also, TJD requested "another approach would be to use a delta transform of branch lengths (versus lambda) – the delta transform is a more direct measure of early vs late burst evolution. I think you only want to do this if it can be implemented by a simple substitution in your Stan code (if not, forget I ever mentioned it)."

lizzieinvancouver commented 2 years ago

Above code with 180 species and 20 ind:

> summary(testme)$summary[names(param), c("2.5%", "25%", "mean", "75%", "97.5%")]
                         2.5%        25%       mean        75%       97.5%
a_z                25.6356160 31.0807784 33.6526045 36.2896353 41.31497921
lam_interceptsa     0.2136373  0.3280266  0.3951826  0.4622820  0.58097259
sigma_interceptsa  14.6138985 15.8084821 16.5755666 17.2753153 18.96086368
b_zf               -7.4069107 -5.1098833 -3.9471401 -2.7631003 -0.43034857
lam_interceptsbf    0.5259026  0.6312682  0.6775191  0.7290883  0.80322904
sigma_interceptsbf  4.7166738  5.1446615  5.4179780  5.6687920  6.20620375
b_zc               -8.5517207 -6.0484456 -4.7931694 -3.4787639 -1.15332992
lam_interceptsbc    0.2128007  0.3252506  0.3888928  0.4500825  0.57154220
sigma_interceptsbc  6.3172444  6.8138248  7.1657860  7.4803275  8.18800558
b_zp               -2.1562684 -1.4609121 -1.1190482 -0.7877977 -0.09188855
lam_interceptsbp    0.2556062  0.3775644  0.4443258  0.5099231  0.63283385
sigma_interceptsbp  1.6375313  1.7953956  1.8921524  1.9803540  2.19612227
sigma_y            12.7356985 12.9369959 13.0556637 13.1692029 13.41790472
> t(param)
     a_z lam_interceptsa sigma_interceptsa b_zf lam_interceptsbf sigma_interceptsbf b_zc lam_interceptsbc sigma_interceptsbc b_zp lam_interceptsbp
[1,] 30  0.4             16                -6   0.7              6                  -7   0.5              7                  -1.4 0.3             
     sigma_interceptsbp sigma_y
[1,] 2  

> fitContinuous(spetree, intercepts, model="lambda")
GEIGER-fitted comparative model of continuous data
 fitted ‘lambda’ model parameters:
    lambda = 0.280980

> fitContinuous(spetree, slopes_bf, model="lambda")
GEIGER-fitted comparative model of continuous data
 fitted ‘lambda’ model parameters:
    lambda = 0.713862

> fitContinuous(spetree, slopes_bc, model="lambda")
GEIGER-fitted comparative model of continuous data
 fitted ‘lambda’ model parameters:
    lambda = 0.277327

> fitContinuous(spetree, slopes_bp, model="lambda")
GEIGER-fitted comparative model of continuous data
 fitted ‘lambda’ model parameters:
    lambda = 0.264521

lizzieinvancouver commented 2 years ago

@MoralesCastilla Gave me some Stan code so I tried to guess at the priors, I will push that code. I presume they are not the current priors as they are a bit too informative ... here's the output you get with them for 40 species and 10 reps:

> summary(testme)$summary[names(param), c("2.5%", "25%", "mean", "75%", "97.5%")]
                          2.5%        25%       mean        75%      97.5%
a_z                19.51776127 25.0751158 28.0261827 30.7739547 37.6448647
lam_interceptsa     0.01764048  0.1868317  0.3739923  0.5461470  0.8193567
sigma_interceptsa  10.41653533 12.9021608 14.7569487 16.2593202 21.0160887
b_zf               -6.64033604 -4.6751918 -3.7827552 -2.9707643 -0.7419553
lam_interceptsbf    0.01560908  0.1341071  0.2924646  0.4259815  0.7032562
sigma_interceptsbf  3.44549992  4.0623032  4.5736742  4.9805250  6.1986628
b_zc               -6.32305088 -3.2747549 -2.0883295 -0.7854263  1.4240409
lam_interceptsbc    0.02652788  0.2382849  0.4283001  0.6133275  0.8476398
sigma_interceptsbc  4.31036003  5.1008849  5.7993734  6.3609889  8.0261278
b_zp               -3.30534305 -2.1783033 -1.6632799 -1.1726396  0.1043179
lam_interceptsbp    0.03463622  0.2526582  0.4587370  0.6565015  0.9139612
sigma_interceptsbp  1.50638397  1.9452331  2.2650324  2.5389862  3.2943143
sigma_y            15.04382998 15.8405956 16.3121778 16.7664784 17.7234642
> t(param)
     a_z lam_interceptsa sigma_interceptsa b_zf lam_interceptsbf sigma_interceptsbf b_zc lam_interceptsbc
[1,] 30  0.4             16                -6   0.7              6                  -7   0.5             
     sigma_interceptsbc b_zp lam_interceptsbp sigma_interceptsbp sigma_y
[1,] 7                  -1.4 0.3              2                  13

And with 180 spp and 20 reps (on quick glance I think it looks better, but need to look more closely):


> summary(testme)$summary[names(param), c("2.5%", "25%", "mean", "75%", "97.5%")]
                          2.5%        25%       mean        75%      97.5%
a_z                26.16788397 31.3460049 33.6749868 35.9931946 40.9447214
lam_interceptsa     0.09471524  0.2387768  0.3275394  0.4136576  0.5765438
sigma_interceptsa  14.28478736 15.4613583 16.3121416 17.0517536 18.9718984
b_zf               -7.47987922 -5.1066596 -3.9563395 -2.8043435 -0.4454207
lam_interceptsbf    0.56754361  0.6819431  0.7276882  0.7810711  0.8506121
sigma_interceptsbf  4.85441833  5.3100113  5.6109722  5.8874111  6.5047324
b_zc               -8.43974922 -6.2264660 -5.1973400 -4.1534252 -2.1338918
lam_interceptsbc    0.12708536  0.2383223  0.3194505  0.3941891  0.5463607
sigma_interceptsbc  6.15913589  6.6442033  6.9841236  7.2589029  8.0955819
b_zp               -2.12458174 -1.4542851 -1.1306831 -0.8014007 -0.1530833
lam_interceptsbp    0.15850374  0.3158795  0.4061253  0.4972156  0.6585695
sigma_interceptsbp  1.61490859  1.7654698  1.8723188  1.9613890  2.2091046
sigma_y            12.92811206 13.1500490 13.2643711 13.3785709 13.6069277
> t(param)
     a_z lam_interceptsa sigma_interceptsa b_zf lam_interceptsbf sigma_interceptsbf b_zc lam_interceptsbc
[1,] 30  0.4             16                -6   0.7              6                  -7   0.5             
     sigma_interceptsbc b_zp lam_interceptsbp sigma_interceptsbp sigma_y
[1,] 7                  -1.4 0.3              2                  13

lizzieinvancouver / pmm

problems returning negative slopes #20