How to include site

[DeirdreLoughnan]

@lizzieinvancouver I have been thinking about how to include site in the joint model. In our last meeting we discussed including latitude in the model as a continuous variable instead of having sites as dummy variables.

I rewrote the model to include latitude as a continuous variable, and now the model generates a slope estimate (b) for each species' traits for every degree latitude, the model code is:

...
transformed parameters{
  // Traits
  vector[N] y_hat; 
  vector[n_spec] mu_grand_sp;
  // Phenology
  real betaForceSp[n_spec];     //species level beta forcing
  real betaPhotoSp[n_spec];     //species level beta photoperiod
  real betaChillSp[n_spec];     //species level beta chilling
  // Traits
  for(i in 1:n_spec){
    mu_grand_sp[i] = mu_grand + muSp[i];
  }
  for (i in 1:N){
    y_hat[i] = mu_grand + muSp[trait_species[i]] + b[trait_species[i]]*lat[i] ;
  }
  // Phenology
  for (isp in 1:n_spec){
    betaForceSp[isp] = alphaForceSp[isp] + betaTraitxForce * (mu_grand_sp[isp]);
  }
  for (isp in 1:n_spec){
    betaPhotoSp[isp] = alphaPhotoSp[isp] + betaTraitxPhoto * (mu_grand_sp[isp]);
  }  
  for (isp in 1:n_spec){
    betaChillSp[isp] = alphaChillSp[isp] + betaTraitxChill * (mu_grand_sp[isp]);
  }
}

Most species have very small changes in traits with latitude and the UI cross zero:

[lizzieinvancouver]

@DeirdreLoughnan Sounds super promising!

[DeirdreLoughnan]

@lizzieinvancouver here is my progress on including site as a continuous variable. For the trait only part of the model, the test data I am using is here and the stan model is here.

This model does a pretty good job of returning the main parameter estimates:

      Parameter Test.data.values  Estiamte        X2.5     X97.5
1       mu.tran                2  2.185294 -0.62549202  9.847331
2      mu.grand               10 10.962462  9.11307772 13.382674
3    sigma.tran                1  1.513281  0.05128542  8.202570
4 sigma.species                5  5.986751  5.98675109  7.490044
5       trt.var                1  1.004396  1.00439634  1.022306

But you see a similar issue as we found with the traitor modeling, where the estimates of muSp are slightly off from simulated values:

[DeirdreLoughnan]

The full model does a fairly good job of returning the model estimates, but with a few issues:

       Parameter Test.data.values    Estiamte         X2.5        X97.5
1         mutran              2.0   2.5638443   0.51764220   4.07779746
2     mu_forcesp            -10.0 -10.1794081 -10.96022586  -9.39852179
3     mu_chillsp            -14.0 -13.4769342 -14.33031321 -12.57029303
4     mu_photosp            -15.0 -15.1733402 -16.19189355 -14.20576673
5     sigma_tran              1.0   0.8203776   0.08279876   3.17578980
6  sigma_forcesp              1.0   0.9499618   0.75386986   1.20518018
7  sigma_chillsp              1.0   1.1246461   0.89149014   1.43157227
8  sigma_photosp              1.0   0.9938081   0.68721399   1.38303431
9   sigma_phenoy              3.0   3.0010979   2.96439322   3.03847445
10         betaF              0.3   0.3469834   0.27130091   0.41905600
11         betaC             -0.4  -0.4331277  -0.51908645  -0.35185634
12         betaP             -0.2  -0.1541842  -0.24529521  -0.05841521

I will continue to work on it and see if I can improve the estimates for mutran, muchillsp in particular.

[DeirdreLoughnan]

@lizzieinvancouver Sadly, running the model with 60 spp instead of 40 did not improve the model estimates, some got a little worse, especially sigmaphotosp:

       Parameter Test.data.values    Estiamte        X2.5         X25           X50         X75       X97.5
1         mutran              2.0   2.8045284   0.1082778   2.2635966   2.7997938   3.3583215   5.3464928
2     mu_forcesp            -10.0 -10.0333392 -10.6280109 -10.2388581  -10.0354260  -9.8293535  -9.4330053
3     mu_chillsp            -14.0 -14.1888138 -14.8266683 -14.4087388  -14.1896654 -13.9771447 -13.5173965
4     mu_photosp            -15.0 -15.0333082 -15.5070913 -15.1951475  -15.0330016 -14.8717375 -14.5617723
5     sigma_tran              1.0   1.5042649   0.4327008   0.8226801   1.2625048   1.9629560   3.8137010
6  sigma_forcesp              1.0   0.9489035   0.7527134   0.8673541   0.9377218   1.0198850   1.2041155
7  sigma_chillsp              1.0   1.0280614   0.8159177   0.9360183   1.0192629   1.1073484   1.2927864
8  sigma_photosp              1.0   0.6972528   0.5504254   0.6370115   0.6926494   0.7493224   0.8729458
9   sigma_phenoy              3.0   2.9976858   2.9613572   2.9854449   2.9973486   3.0099108   3.0347593
10         betaF              0.3   0.3254536   0.2730196   0.3069642   0.3265998   0.3444778   0.3754525
11         betaC             -0.4  -0.3909004  -0.4500457  -0.4092934  -0.3906983  -0.3721440  -0.3347560
12         betaP             -0.2  -0.2031380  -0.2454293  -0.2176397  -0.2033093  -0.1889107  -0.1615634

But on a positive note, removing the grand mean does improve the species estimates!

[lizzieinvancouver]

@DeirdreLoughnan I would check what might be up with the sigma_photosp and sigma_tran -- since it's not getting those.I am actually not sure what sigma_tran is, given that there are only two transects and thus would not be good for partial pooling ... (and I don't see where latitude is, but it would not need a sigma either since it's continuous).

[DeirdreLoughnan]

@lizzieinvancouver if you have time this afternoon, could we chat about this in more detail?

[DeirdreLoughnan]

@lizzieinvancouver I continue to have issues with how to incorporate site into this model. After our last conversation, we talked about the model including transect as a dummy variable and an interaction between the transect dummy variable and latitude:

y_hat[i] =   muSp[species[i]] + b_tranE * tranE[i]  + b_tranlat * (tranE * latitude)

But I seem stuck at getting the test data for the interaction to return the true values.

I took several steps back and was able to return the test data values fairly well for model with muSp and the transect dummy variable:

     Parameter Test.data.values Estiamte      X2.5      X25      X50      X75     X97.5
1       bTranE                4 4.049068 4.0022755 4.032535 4.048875 4.065402  4.097475
2     sigma_sp                5 5.883881 0.3411229 2.837319 5.342898 8.338506 14.544073
3 sigma_traity                1 1.007117 0.9904040 1.001276 1.007029 1.013001  1.024129

The estimate for sigma_sp is a little high, but the estimates for individual species look good:

But when I try adding latitude, either as a simple slope or as an interactions, the estimates are really bad:

     Parameter Test.data.values   Estiamte        X2.5        X25        X50        X75     X97.5
1       bTranE                4 -31.352190 -41.0354092 -34.667673 -31.318554 -28.003578 -21.34691
2     bTranLat                2 108.616257 108.3107336 108.511392 108.617082 108.720196 108.91374
3     sigma_sp                5   5.826951   0.3217277   2.884909   5.327186   8.299395  13.88553
4 sigma_traity                1 287.484793 283.6445453 286.149891 287.472565 288.839452 291.43939

I think part of the issue is the test data, but I don't have a rational reason for changing it and I think I could really benefit from talking about it more and a second set of eyes.

[DeirdreLoughnan]

@lizzieinvancouver I figured out the issue was with the way I was coding the interaction parameter and then adding it to calculate yTraiti.

The model output for the subsets of the data looked much better, so I went ahead and ran the model for the full model. The output looks pretty good, but not perfect:

       Parameter Test.data.values    Estiamte        X2.5         X25          X50         X75       X97.5
1         bTranE              4.0   3.4784021   2.4083588   3.1270973    3.4829672   3.8416681   4.5031540
2       bTranLat              2.0   2.0106743   1.9900568   2.0034633    2.0105538   2.0178052   2.0320489
3     mu_forcesp            -10.0  -9.4105091 -10.4492824  -9.7543079  -9.4324426  -9.0706732  -8.3310481
4     mu_chillsp            -14.0 -14.3780959 -15.4860732 -14.7371280  -14.3678980 -14.0156398 -13.3100039
5     mu_photosp            -15.0 -15.4811907 -16.3143809 -15.7422628  -15.4818317 -15.2066908 -14.6500921
6       sigma_sp              5.0   5.0071388   3.7836999   4.4637372   4.9248207   5.4657316   6.6598418
7   sigma_traity              1.0   1.0059972   0.9844312   0.9983096   1.0057236   1.0137590   1.0281675
8  sigma_forcesp              1.0   1.1932571   0.8867785   1.0577457   1.1702944   1.2988199   1.6320695
9  sigma_chillsp              1.0   1.1726527   0.8765847   1.0428704   1.1516221   1.2725267   1.5931743
10 sigma_photosp              1.0   0.9030784   0.6800047   0.8063770   0.8859213   0.9861722   1.2188055
11  sigma_phenoy              3.0   2.9732855   2.9258574   2.9570208   2.9733017   2.9896074   3.0191808
12         betaF              0.3   0.2706549   0.1699054   0.2398130   0.2712550   0.3010979   0.3683886
13         betaC             -0.4  -0.3536558  -0.4505924  -0.3890603  -0.3544528  -0.3202395  -0.2520712
14         betaP             -0.2  -0.1511724  -0.2292865  -0.1757392  -0.1509959  -0.1263453  -0.0758689

But I am definitely happy with the progress for today. Thanks for you help!

[DeirdreLoughnan]

I ran the model with twice as many spp and reps (40 vs 20 above) and the output does look better, except for bTranE and sigma_sp:

       Parameter Test.data.values    Estiamte        X2.5         X25           X50         X75       X97.5
1         bTranE              4.0   3.6153540   3.0777352   3.4202407   3.6155720   3.8080656   4.1759040
2       bTranLat              2.0   2.0069510   1.9965469   2.0033246   2.0069558   2.0105099   2.0169421
3     mu_forcesp            -10.0  -9.9170983 -10.8544477 -10.2529963  -9.9139155  -9.5949105  -8.9577225
4     mu_chillsp            -14.0 -13.9974622 -14.9387317 -14.3091020  -13.9919772 -13.6831298 -13.0615258
5     mu_photosp            -15.0 -15.3347722 -16.2662611 -15.6430644  -15.3370279 -15.0083532 -14.4325579
6       sigma_sp              5.0   4.4006001   3.5620069   4.0456382   4.3571033   4.7024305   5.4796514
7   sigma_traity              1.0   1.0052632   0.9932787   1.0011191   1.0052708   1.0094949   1.0171449
8  sigma_forcesp              1.0   1.1222524   0.9009913   1.0288660   1.1090117   1.2011217   1.4146107
9  sigma_chillsp              1.0   1.1562696   0.9307797   1.0576007   1.1433923   1.2394430   1.4611381
10 sigma_photosp              1.0   1.1187904   0.8856199   1.0167166   1.1084341   1.2036920   1.4181994
11  sigma_phenoy              3.0   2.9930129   2.9667543   2.9839269   2.9928577   3.0018097   3.0196973
12         betaF              0.3   0.3147594   0.2258027   0.2851414   0.3153203   0.3449079   0.3984852
13         betaC             -0.4  -0.3960073  -0.4812834  -0.4234247  -0.3966884  -0.3681156  -0.3099671
14         betaP             -0.2  -0.1962686  -0.2823979  -0.2251928  -0.1964218  -0.1683275  -0.1099325

But the species level estimates are awesome!

[lizzieinvancouver]

@DeirdreLoughnan Nice work! These results make me think the model is good for our purposes. If you want, you could run the full code again (including the R code so you get slightly different species) and check it's not systematically off, but I would be happy with the results myself and move on.