SHORT_TITLE: Population Changes END_SHORT_TITLE

FLUCTUATION IN ABUNDANCE OF TREE SPECIES

R. CONDIT

Date: July 11, 2012.

DESCRIPTION: The question addressed here is how much abundances of all the tree species in a plot fluctuate from census to census. At one extreme, abundances are stabilized, and populations do not change. In contrast, there may be environmental drivers favoring some species and not others, so that some species show marked (ie, more than neutral) increases and others decreases. In between the two extremes is a community in which no external drivers affect abundances but populations fluctuate due to demographic stochasticity: a neutral forest.

The test of these models hinges on a histogram of the rate of population change of every species in the forest. Define λ = N1
N0, the ratio of a future population N1 to the current N0, and r = little.r = logλ. In a highly stabilized community, λ 1 for every species, and r = 0, whereas in a forest with environmental fluctionas, λ may be distributed widely but ought to be centered close to 1. The neutral community would fall in between, with λ varying a small amount due to chance alone.

The distribution of λ and r under neutrality can be predicted mathematically, assuming deaths are binomially distributed across individuals and a Poisson distribution of births whose mean matches the death rate. Because of this statistical model of neutral fluctuations, any environmentally caused variation can be estimated precisely. In a conceptual sense, it means subtracting the predicted neutral variance of λ from the observed variance of λ; what’s left is the environmental variances. The neutral prediction depends on the mortality rate, which is also fitted by the model. The following demonstrates an R program for running this model with forest plot abundances.

The model allows various functional forms for the distribution of population growth. The example below is based on a power probability distribution, where p(x) x + 1k (k < 1). The exponent k describes how broad the distribution, with more extreme k (ie, more negative) meaning a narrower distribution. Indeed, frac1k is close to the standard deviation of the distribution of r. The power distribution, though, allows k on the negative size of x to differ from k for x positive. It is thus an asymmetric power distribution. END_DESCRIPTION

1. Formatting data

The calculations start with the full plot dataframes described at http://ctfs.arnarb.harvard.edu/Public/CTFSRPackage/help/RoutputFull.pdf. The BCI data can be downloaded from http://ctfs.arnarb.harvard.edu/bci/RAnalyticalTables after filling out the request form http://ctfs.arnarb.harvard.edu/webatlas/datasets/bci/. Assume for the following example that dataframes for three censues are loaded, bci.full1, bci.full2, and bci.full3. To make it easier to follow this tutorial, the three dataframes are copied to three new dataframes with the generic names census1, census2, and census3.

  [1] 5 4 3

  > head(bci.full1)

    treeID stemID    tag StemTag     sp quadrat    gx    gy MeasureID CensusID dbh  pom hom  ExactDate DFstatus codes nostems     date status          agb
  1      1      1 -05599         swars1    4007 800.2 152.2         1        1  90  1.3 1.3 1981-06-12    alive  <NA>       1 7833.000      A 0.0426371077
  2      2      1 -22851         hybapr    0718 151.5 378.8         2        1  35  1.3 1.3 1982-07-01    alive     L       1 8217.000      A 0.0028894861
  3      3      1 -24362         aegipa    0417  95.2 357.5         3        1  10  1.3 1.3 1982-07-02    alive  <NA>       1 8218.000      A 0.0001466688
  4      4     NA -26589    <NA> beilpe    0007  11.7 151.1        NA       NA  NA <NA>  NA       <NA>    prior  <NA>      NA 8237.215      P 0.0000000000
  5      5     NA -26590    <NA> faraoc    0004   7.7  96.2        NA       NA  NA <NA>  NA       <NA>    prior  <NA>      NA 8233.467      P 0.0000000000
  6      6     NA -26703    <NA> hybapr    0210  50.1 215.4        NA       NA  NA <NA>  NA       <NA>    prior  <NA>      NA 8234.400      P 0.0000000000

  > head(bci.full2)

    treeID stemID    tag StemTag     sp quadrat    gx    gy MeasureID CensusID dbh  pom hom  ExactDate DFstatus codes nostems date status          agb
  1      1      1 -05599         swars1    4007 800.2 152.2         1        2  94  1.3 1.3 1985-06-28    alive     Z       1 9310      A 0.0478479558
  2      2      1 -22851         hybapr    0718 151.5 378.8         2        2  NA <NA>  NA 1985-03-06     dead  D,MC       1 9196      D 0.0000000000
  3      3      1 -24362         aegipa    0417  95.2 357.5         3        2  15  1.3 1.3 1985-02-06    alive  <NA>       1 9168      A 0.0003618661
  4      4      1 -26589         beilpe    0007  11.7 151.1         4        2  10  1.3 1.3 1985-01-11    alive     Z       1 9142      A 0.0001132402
  5      5      1 -26590         faraoc    0004   7.7  96.2         5        2  10  1.3 1.3 1985-01-08    alive     Z       1 9139      A 0.0001313070
  6      6      1 -26703         hybapr    0210  50.1 215.4         6        2  10  1.3 1.3 1985-01-23    alive     Z       1 9154      A 0.0001496468

  > head(bci.full3)

    treeID stemID    tag StemTag     sp quadrat    gx    gy MeasureID CensusID dbh  pom hom  ExactDate DFstatus codes nostems  date status agb
  1      1      1 -05599         swars1    4007 800.2 152.2         1        3  NA <NA>  NA 1990-10-20     dead   D,N       1 11250      D   0
  2      2     NA -22851    <NA> hybapr    0718 151.5 378.8        NA       NA  NA <NA>  NA       <NA>     dead  <NA>      NA 11135      D   0
  3      3      1 -24362         aegipa    0417  95.2 357.5         2        3  NA <NA>  NA 1990-04-06     dead   D,N       1 11053      D   0
  4      4      1 -26589         beilpe    0007  11.7 151.1         3        3  NA <NA>  NA 1990-02-23     dead   D,N       1 11011      D   0
  5      5      1 -26590         faraoc    0004   7.7  96.2         4        3  NA <NA>  NA 1990-02-16     dead   D,N       1 11004      D   0
  6      6      1 -26703         hybapr    0210  50.1 215.4         5        3  NA <NA>  NA 1990-03-07     dead   D,N       1 11023      D   0

2. Running the model of environmental variance

A function in the CTFSRPackage (topic population change) called model.littleR.Gibbs runs the model to estimate leftover environmental variance over two census intervals. To run with your own data, simply replace census1 and census2 with your own R Analytical Tables for a CTFS plot. The fitting procedure requires a Bayesian parameter-search, which must run at least 1000 steps. I recommend here running the model once for 1200 steps, as in the example below in order to explore results. In order to get really solid parameter estimates, a run of 10,000 steps should be set. Ddepending on your computer, 1000 steps should take about 3 minutes and 10,000 steps 30 minutes. Note that mindbh needs to be set. The example below assumes dbh is in mm, and is set to 10.

  > mod=model.littleR.Gibbs(cns1=bci.full1,cns2=bci.full2,modeltype='asympower',mindbh=10,start.param=c(-3,.8,.01,-.5),
  +                         bad.modelparam=bad.asympower.param,steps=1200,burn=200,showstep=25)

3. Model results

The result, saved object as mod, is a Bayesian model object including the estimated parameters and their confidence, plus the complete parameter search. The element of the list mean has the fitted parameters. The first two describe the estimated distribution of mortality rates:

The next three parameters describe the distribution of population growth rates, as defined by r = logλ:

Upper and lower confidence limits for the five parameters are returned in the components upper and lower.

  > mod$means

       hyperMu      hyperSD       hyperR   hyperSDlow    hyperSDup
  -3.143127188  0.799571923  0.007376671 -0.033373853 -0.021036037

  > mod$upper

     hyperMu.97.5%    hyperSD.97.5%     hyperR.97.5% hyperSDlow.97.5%  hyperSDup.97.5%
       -3.04806004       0.86449199       0.01205397      -0.02831436      -0.01603845

  > mod$lower

     hyperMu.2.5%    hyperSD.2.5%     hyperR.2.5% hyperSDlow.2.5%  hyperSDup.2.5%
      -3.23394167      0.69632734      0.00236593     -0.03942421     -0.02513421

There is another function in the same topic of the CTFSRPackage called graph.abundmodel that creates a graphical view of the fitted model. Pass it the name of the model object and nothing more and it will show a histogram of observed rates of population change r (open circles for all species, blue circles for abundant species), along with the fit to the histogram (green line).

The green curve is the estimated variation in in abundance due to the environment: what is left over after stochastic demography is accounted for. The points, which give the observed histogram, include stochastic plus environmental, thus should be broader than the green curve.

The x-limits of the graph can be adjusted using the argument xrange. In addition, confidence limits can be added by setting conf=50, or any other positive number (the number of lines to draw). The confidence lines are based on all the steps in the Bayesian sampling and reflect confidence in the fitted model.

  > graph.abundmodel(fit=mod)

  $Fastest_increases
           N1   N2    S     time    mortrate   little.r    fitmort   lowermean  uppermean       fitr  lowermean.1 uppermean.1
  paligu  376  661  247 3.053886 0.137595454 0.18473667 0.13580836 0.114566514 0.16348333 0.17947287  0.160145553  0.19976713
  clidde    8   14    7 2.796704 0.047745989 0.20009834 0.05042199 0.007695854 0.12725222 0.09222749  0.012910178  0.20616605
  ochrpy    4    8    3 2.691700 0.106877479 0.25751282 0.07209790 0.011949876 0.22434263 0.07789825 -0.003086927  0.19329611
  micoar  528  677  389 3.195909 0.095596265 0.07777913 0.09400349 0.079693081 0.10928199 0.07352084  0.054718502  0.09178355
  ingath   61   78   58 2.649370 0.019035038 0.09278998 0.02570207 0.010488673 0.05038648 0.07164201  0.037096014  0.11302327
  cha2sc  194  239  190 2.895034 0.007196491 0.07205628 0.01102246 0.004923829 0.02073393 0.06783491  0.052186529  0.08665592
  cuparu   56   73   49 3.273421 0.040792615 0.08098799 0.04214098 0.018211934 0.07011654 0.06306315  0.025184741  0.10037320
  eugega  963 1162  917 3.345195 0.014631716 0.05615354 0.01487109 0.010590150 0.01983953 0.05491912  0.046786853  0.06450339
  pipeco 3142 3706 2548 3.019620 0.069396421 0.05467364 0.06896105 0.064359946 0.07414464 0.05397923  0.046916452  0.06017186
  ingas1  211  256  204 3.477802 0.009700995 0.05558663 0.01276202 0.006178632 0.02195466 0.05284752  0.036660953  0.07092626
  
  $Biggest_losses
          N1  N2   S     time  mortrate    little.r   fitmort  lowermean uppermean        fitr lowermean.1  uppermean.1
  bactc1 241  84  84 3.500939 0.3010564 -0.30105644 0.2919958 0.25117819 0.3444808 -0.28853423  -0.3320738 -0.239333266
  bactba 111  55  55 3.388399 0.2072357 -0.20723565 0.1993966 0.15314841 0.2600171 -0.18436620  -0.2363322 -0.119488754
  pipecu 120  65  62 3.012689 0.2191920 -0.20350736 0.2081870 0.15105210 0.2717600 -0.17911681  -0.2410535 -0.109033896
  bactc2  38  17  17 3.290464 0.2444557 -0.24445571 0.2093493 0.11905326 0.3218468 -0.15663695  -0.2823522 -0.019543977
  hampap  76  49  41 3.217818 0.1917950 -0.13640085 0.1745544 0.12478946 0.2437310 -0.10887454  -0.1927807 -0.022790842
  sennda 201 135 114 3.363347 0.1686137 -0.11834345 0.1628252 0.13144848 0.2036193 -0.10781092  -0.1461227 -0.073600848
  cecrob  61  38  27 3.123925 0.2609016 -0.15150417 0.2447049 0.16448578 0.3436948 -0.10137663  -0.1867371 -0.009922695
  conoci 389 279 232 3.555180 0.1453772 -0.09348825 0.1419213 0.11942632 0.1659425 -0.09128813  -0.1148889 -0.066206010
  solaha 125  89  63 3.055432 0.2242495 -0.11117163 0.2181227 0.17405402 0.2763442 -0.09093615  -0.1412767 -0.039738815
  bactma 469 355 341 2.993622 0.1064664 -0.09302609 0.1056323 0.08544943 0.1261681 -0.09020806  -0.1066020 -0.071652058

PIC

  > graph.abundmodel(fit=mod,xrange=c(-.1,.1))

  $Fastest_increases
           N1   N2    S     time    mortrate   little.r    fitmort   lowermean  uppermean       fitr  lowermean.1 uppermean.1
  paligu  376  661  247 3.053886 0.137595454 0.18473667 0.13580836 0.114566514 0.16348333 0.17947287  0.160145553  0.19976713
  clidde    8   14    7 2.796704 0.047745989 0.20009834 0.05042199 0.007695854 0.12725222 0.09222749  0.012910178  0.20616605
  ochrpy    4    8    3 2.691700 0.106877479 0.25751282 0.07209790 0.011949876 0.22434263 0.07789825 -0.003086927  0.19329611
  micoar  528  677  389 3.195909 0.095596265 0.07777913 0.09400349 0.079693081 0.10928199 0.07352084  0.054718502  0.09178355
  ingath   61   78   58 2.649370 0.019035038 0.09278998 0.02570207 0.010488673 0.05038648 0.07164201  0.037096014  0.11302327
  cha2sc  194  239  190 2.895034 0.007196491 0.07205628 0.01102246 0.004923829 0.02073393 0.06783491  0.052186529  0.08665592
  cuparu   56   73   49 3.273421 0.040792615 0.08098799 0.04214098 0.018211934 0.07011654 0.06306315  0.025184741  0.10037320
  eugega  963 1162  917 3.345195 0.014631716 0.05615354 0.01487109 0.010590150 0.01983953 0.05491912  0.046786853  0.06450339
  pipeco 3142 3706 2548 3.019620 0.069396421 0.05467364 0.06896105 0.064359946 0.07414464 0.05397923  0.046916452  0.06017186
  ingas1  211  256  204 3.477802 0.009700995 0.05558663 0.01276202 0.006178632 0.02195466 0.05284752  0.036660953  0.07092626
  
  $Biggest_losses
          N1  N2   S     time  mortrate    little.r   fitmort  lowermean uppermean        fitr lowermean.1  uppermean.1
  bactc1 241  84  84 3.500939 0.3010564 -0.30105644 0.2919958 0.25117819 0.3444808 -0.28853423  -0.3320738 -0.239333266
  bactba 111  55  55 3.388399 0.2072357 -0.20723565 0.1993966 0.15314841 0.2600171 -0.18436620  -0.2363322 -0.119488754
  pipecu 120  65  62 3.012689 0.2191920 -0.20350736 0.2081870 0.15105210 0.2717600 -0.17911681  -0.2410535 -0.109033896
  bactc2  38  17  17 3.290464 0.2444557 -0.24445571 0.2093493 0.11905326 0.3218468 -0.15663695  -0.2823522 -0.019543977
  hampap  76  49  41 3.217818 0.1917950 -0.13640085 0.1745544 0.12478946 0.2437310 -0.10887454  -0.1927807 -0.022790842
  sennda 201 135 114 3.363347 0.1686137 -0.11834345 0.1628252 0.13144848 0.2036193 -0.10781092  -0.1461227 -0.073600848
  cecrob  61  38  27 3.123925 0.2609016 -0.15150417 0.2447049 0.16448578 0.3436948 -0.10137663  -0.1867371 -0.009922695
  conoci 389 279 232 3.555180 0.1453772 -0.09348825 0.1419213 0.11942632 0.1659425 -0.09128813  -0.1148889 -0.066206010
  solaha 125  89  63 3.055432 0.2242495 -0.11117163 0.2181227 0.17405402 0.2763442 -0.09093615  -0.1412767 -0.039738815
  bactma 469 355 341 2.993622 0.1064664 -0.09302609 0.1056323 0.08544943 0.1261681 -0.09020806  -0.1066020 -0.071652058

PIC

  > graph.abundmodel(fit=mod,xrange=c(-.1,.1),conf=50)

  $Fastest_increases
           N1   N2    S     time    mortrate   little.r    fitmort   lowermean  uppermean       fitr  lowermean.1 uppermean.1
  paligu  376  661  247 3.053886 0.137595454 0.18473667 0.13580836 0.114566514 0.16348333 0.17947287  0.160145553  0.19976713
  clidde    8   14    7 2.796704 0.047745989 0.20009834 0.05042199 0.007695854 0.12725222 0.09222749  0.012910178  0.20616605
  ochrpy    4    8    3 2.691700 0.106877479 0.25751282 0.07209790 0.011949876 0.22434263 0.07789825 -0.003086927  0.19329611
  micoar  528  677  389 3.195909 0.095596265 0.07777913 0.09400349 0.079693081 0.10928199 0.07352084  0.054718502  0.09178355
  ingath   61   78   58 2.649370 0.019035038 0.09278998 0.02570207 0.010488673 0.05038648 0.07164201  0.037096014  0.11302327
  cha2sc  194  239  190 2.895034 0.007196491 0.07205628 0.01102246 0.004923829 0.02073393 0.06783491  0.052186529  0.08665592
  cuparu   56   73   49 3.273421 0.040792615 0.08098799 0.04214098 0.018211934 0.07011654 0.06306315  0.025184741  0.10037320
  eugega  963 1162  917 3.345195 0.014631716 0.05615354 0.01487109 0.010590150 0.01983953 0.05491912  0.046786853  0.06450339
  pipeco 3142 3706 2548 3.019620 0.069396421 0.05467364 0.06896105 0.064359946 0.07414464 0.05397923  0.046916452  0.06017186
  ingas1  211  256  204 3.477802 0.009700995 0.05558663 0.01276202 0.006178632 0.02195466 0.05284752  0.036660953  0.07092626
  
  $Biggest_losses
          N1  N2   S     time  mortrate    little.r   fitmort  lowermean uppermean        fitr lowermean.1  uppermean.1
  bactc1 241  84  84 3.500939 0.3010564 -0.30105644 0.2919958 0.25117819 0.3444808 -0.28853423  -0.3320738 -0.239333266
  bactba 111  55  55 3.388399 0.2072357 -0.20723565 0.1993966 0.15314841 0.2600171 -0.18436620  -0.2363322 -0.119488754
  pipecu 120  65  62 3.012689 0.2191920 -0.20350736 0.2081870 0.15105210 0.2717600 -0.17911681  -0.2410535 -0.109033896
  bactc2  38  17  17 3.290464 0.2444557 -0.24445571 0.2093493 0.11905326 0.3218468 -0.15663695  -0.2823522 -0.019543977
  hampap  76  49  41 3.217818 0.1917950 -0.13640085 0.1745544 0.12478946 0.2437310 -0.10887454  -0.1927807 -0.022790842
  sennda 201 135 114 3.363347 0.1686137 -0.11834345 0.1628252 0.13144848 0.2036193 -0.10781092  -0.1461227 -0.073600848
  cecrob  61  38  27 3.123925 0.2609016 -0.15150417 0.2447049 0.16448578 0.3436948 -0.10137663  -0.1867371 -0.009922695
  conoci 389 279 232 3.555180 0.1453772 -0.09348825 0.1419213 0.11942632 0.1659425 -0.09128813  -0.1148889 -0.066206010
  solaha 125  89  63 3.055432 0.2242495 -0.11117163 0.2181227 0.17405402 0.2763442 -0.09093615  -0.1412767 -0.039738815
  bactma 469 355 341 2.993622 0.1064664 -0.09302609 0.1056323 0.08544943 0.1261681 -0.09020806  -0.1066020 -0.071652058

PIC

4. Running the model for several censuses

There is a wrapper function runSeveralAbundFit that can handle any number of census datasets. It simply fits the abundance model for each pair of censuses, then for the first and last census, saving all into one big list, plus it repeats all the models for larger trees, setting mindbh to 10 times the submitted value. It’s only purpose is to finish a series of models from one plot in a single model run. Since it takes about 30 minutes per census, you might leave it running overnight. The following code saves the results to an R data object in the folder abunddist.

  [1] 2

  > modall=fitSeveralAbundModel(allcns=list(bci.full1,bci.full2,bci.full3),mindbh=10,steps=1200,burn=200,show=25,bad.modelparam=bad.asympower.param)
  > save(modall,file='abunddist/abundFitModel.rdata')