skyline: Skyline Plot Estimate of Effective Population Size

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/skyline.R

Description

skyline computes the generalized skyline plot estimate of effective population size from an estimated phylogeny. The demographic history is approximated by a step-function. The number of parameters of the skyline plot (i.e. its smoothness) is controlled by a parameter epsilon.

find.skyline.epsilon searches for an optimal value of the epsilon parameter, i.e. the value that maximizes the AICc-corrected log-likelihood (logL.AICc).

Usage

1
2
3
4
5
6
7
8
skyline(x, ...)
## S3 method for class 'phylo'
skyline(x, ...)
## S3 method for class 'coalescentIntervals'
skyline(x, epsilon=0, ...)
## S3 method for class 'collapsedIntervals'
skyline(x, old.style=FALSE, ...)
find.skyline.epsilon(ci, GRID=1000, MINEPS=1e-6, ...)

Arguments

x

Either an ultrametric tree (i.e. an object of class "phylo"), or coalescent intervals (i.e. an object of class "coalescentIntervals"), or collapsed coalescent intervals (i.e. an object of class "collapsedIntervals").

epsilon

collapsing parameter that controls the amount of smoothing (allowed range: from 0 to ci$total.depth, default value: 0). This is the same parameter as in collapsed.intervals.

old.style

Parameter to choose between two slightly different variants of the generalized skyline plot (Strimmer and Pybus, pers. comm.). The default value FALSE is recommended.

ci

coalescent intervals (i.e. an object of class "coalescentIntervals")

GRID

Parameter for the grid search for epsilon in find.skyline.epsilon.

MINEPS

Parameter for the grid search for epsilon in find.skyline.epsilon.

...

Any of the above parameters.

Details

skyline implements the generalized skyline plot introduced in Strimmer and Pybus (2001). For epsilon = 0 the generalized skyline plot degenerates to the classic skyline plot described in Pybus et al. (2000). The latter is in turn directly related to lineage-through-time plots (Nee et al., 1995).

Value

skyline returns an object of class "skyline" with the following entries:

time

A vector with the time at the end of each coalescent interval (i.e. the accumulated interval lengths from the beginning of the first interval to the end of an interval)

interval.length

A vector with the length of each interval.

population.size

A vector with the effective population size of each interval.

parameter.count

Number of free parameters in the skyline plot.

epsilon

The value of the underlying smoothing parameter.

logL

Log-likelihood of skyline plot (see Strimmer and Pybus, 2001).

logL.AICc

AICc corrected log-likelihood (see Strimmer and Pybus, 2001).

find.skyline.epsilon returns the value of the epsilon parameter that maximizes logL.AICc.

Author(s)

Korbinian Strimmer

References

Strimmer, K. and Pybus, O. G. (2001) Exploring the demographic history of DNA sequences using the generalized skyline plot. Molecular Biology and Evolution, 18, 2298–2305.

Pybus, O. G, Rambaut, A. and Harvey, P. H. (2000) An integrated framework for the inference of viral population history from reconstructed genealogies. Genetics, 155, 1429–1437.

Nee, S., Holmes, E. C., Rambaut, A. and Harvey, P. H. (1995) Inferring population history from molecular phylogenies. Philosophical Transactions of the Royal Society of London. Series B. Biological Sciences, 349, 25–31.

See Also

coalescent.intervals, collapsed.intervals, skylineplot, ltt.plot.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
# get tree
data("hivtree.newick") # example tree in NH format
tree.hiv <- read.tree(text = hivtree.newick) # load tree

# corresponding coalescent intervals
ci <- coalescent.intervals(tree.hiv) # from tree

# collapsed intervals
cl1 <- collapsed.intervals(ci,0)
cl2 <- collapsed.intervals(ci,0.0119)

#### classic skyline plot ####
sk1 <- skyline(cl1)        # from collapsed intervals 
sk1 <- skyline(ci)         # from coalescent intervals
sk1 <- skyline(tree.hiv)   # from tree
sk1

plot(skyline(tree.hiv))
skylineplot(tree.hiv) # shortcut

plot(sk1, show.years=TRUE, subst.rate=0.0023, present.year = 1997)

#### generalized skyline plot ####

sk2 <- skyline(cl2)              # from collapsed intervals
sk2 <- skyline(ci, 0.0119)       # from coalescent intervals
sk2 <- skyline(tree.hiv, 0.0119) # from tree
sk2

plot(sk2)

# classic and generalized skyline plot together in one plot
plot(sk1, show.years=TRUE, subst.rate=0.0023, present.year = 1997, col=c(grey(.8),1))
lines(sk2,  show.years=TRUE, subst.rate=0.0023, present.year = 1997)
legend(.15,500, c("classic", "generalized"), col=c(grey(.8),1),lty=1)


# find optimal epsilon parameter using AICc criterion
find.skyline.epsilon(ci)

sk3 <- skyline(ci, -1) # negative epsilon also triggers estimation of epsilon
sk3$epsilon

Example output

$time
  [1] 0.021161 0.049822 0.050058 0.054444 0.057580 0.057620 0.058582 0.059336
  [9] 0.059726 0.061847 0.063012 0.065150 0.065542 0.067767 0.068246 0.068285
 [17] 0.069341 0.069631 0.069725 0.070130 0.070506 0.071398 0.077168 0.077774
 [25] 0.078085 0.078439 0.078488 0.079791 0.079902 0.080093 0.080206 0.080907
 [33] 0.080919 0.080923 0.081266 0.081732 0.082433 0.082434 0.082720 0.083124
 [41] 0.083172 0.083436 0.083818 0.084541 0.084566 0.084629 0.085711 0.086004
 [49] 0.087462 0.087519 0.087787 0.088609 0.088644 0.088785 0.088823 0.089216
 [57] 0.089479 0.089540 0.090237 0.090248 0.090579 0.090777 0.091414 0.091575
 [65] 0.092297 0.092865 0.092881 0.093773 0.094006 0.094130 0.094174 0.094384
 [73] 0.094521 0.094683 0.094749 0.095068 0.095107 0.095333 0.095334 0.095369
 [81] 0.096163 0.096270 0.096432 0.096562 0.097150 0.097278 0.097341 0.097342
 [89] 0.099386 0.099467 0.099582 0.099726 0.099948 0.101084 0.101878 0.101892
 [97] 0.101961 0.101962 0.102399 0.102590 0.103060 0.103714 0.103878 0.103926
[105] 0.104386 0.104484 0.104691 0.105017 0.105018 0.105020 0.105181 0.106054
[113] 0.106195 0.107063 0.107287 0.107629 0.108329 0.108551 0.108702 0.109333
[121] 0.109996 0.111546 0.111869 0.112580 0.112666 0.113317 0.113663 0.113664
[129] 0.113731 0.113986 0.114132 0.114392 0.114561 0.114632 0.114774 0.116003
[137] 0.116109 0.117194 0.117310 0.117514 0.117728 0.118538 0.118650 0.118782
[145] 0.118990 0.119083 0.120265 0.120408 0.120491 0.120494 0.122054 0.122422
[153] 0.122572 0.123109 0.124283 0.125240 0.127902 0.127903 0.127904 0.127999
[161] 0.129007 0.130366 0.130367 0.130470 0.131217 0.131462 0.132291 0.135095
[169] 0.135794 0.135795 0.135914 0.140307 0.143055 0.144389 0.146425 0.148331
[177] 0.150189 0.151362 0.155704 0.164770 0.166301 0.166846 0.176701 0.178108
[185] 0.179167 0.181334 0.192334 0.196999 0.197000 0.204899 0.204900 0.209112

$interval.length
  [1] 0.021161 0.028661 0.000236 0.004386 0.003136 0.000040 0.000962 0.000754
  [9] 0.000390 0.002121 0.001165 0.002138 0.000392 0.002225 0.000479 0.000039
 [17] 0.001056 0.000290 0.000094 0.000405 0.000376 0.000892 0.005770 0.000606
 [25] 0.000311 0.000354 0.000049 0.001303 0.000111 0.000191 0.000113 0.000701
 [33] 0.000012 0.000004 0.000343 0.000466 0.000701 0.000001 0.000286 0.000404
 [41] 0.000048 0.000264 0.000382 0.000723 0.000025 0.000063 0.001082 0.000293
 [49] 0.001458 0.000057 0.000268 0.000822 0.000035 0.000141 0.000038 0.000393
 [57] 0.000263 0.000061 0.000697 0.000011 0.000331 0.000198 0.000637 0.000161
 [65] 0.000722 0.000568 0.000016 0.000892 0.000233 0.000124 0.000044 0.000210
 [73] 0.000137 0.000162 0.000066 0.000319 0.000039 0.000226 0.000001 0.000035
 [81] 0.000794 0.000107 0.000162 0.000130 0.000588 0.000128 0.000063 0.000001
 [89] 0.002044 0.000081 0.000115 0.000144 0.000222 0.001136 0.000794 0.000014
 [97] 0.000069 0.000001 0.000437 0.000191 0.000470 0.000654 0.000164 0.000048
[105] 0.000460 0.000098 0.000207 0.000326 0.000001 0.000002 0.000161 0.000873
[113] 0.000141 0.000868 0.000224 0.000342 0.000700 0.000222 0.000151 0.000631
[121] 0.000663 0.001550 0.000323 0.000711 0.000086 0.000651 0.000346 0.000001
[129] 0.000067 0.000255 0.000146 0.000260 0.000169 0.000071 0.000142 0.001229
[137] 0.000106 0.001085 0.000116 0.000204 0.000214 0.000810 0.000112 0.000132
[145] 0.000208 0.000093 0.001182 0.000143 0.000083 0.000003 0.001560 0.000368
[153] 0.000150 0.000537 0.001174 0.000957 0.002662 0.000001 0.000001 0.000095
[161] 0.001008 0.001359 0.000001 0.000103 0.000747 0.000245 0.000829 0.002804
[169] 0.000699 0.000001 0.000119 0.004393 0.002748 0.001334 0.002036 0.001906
[177] 0.001858 0.001173 0.004342 0.009066 0.001531 0.000545 0.009855 0.001407
[185] 0.001059 0.002167 0.011000 0.004665 0.000001 0.007899 0.000001 0.004212

$population.size
  [1] 392.071008 525.528096   4.282220  78.750630  55.714176   0.703120
  [7]  16.730142  12.972570   6.637800  35.709156  19.400745  35.214998
 [13]   6.385680  35.844750   7.630949   0.614367  16.448256   4.466000
 [19]   1.431150   6.095655   5.594128  13.117752  83.866950   8.705190
 [25]   4.414956   4.965912   0.679189  17.844585   1.501830   2.552906
 [31]   1.491939   9.141741   0.154560   0.050880   4.308423   5.779798
 [37]   8.584446   0.012090   3.413410   4.759524   0.558144   3.029664
 [43]   4.326150   8.079525   0.275650   0.685314  11.610942   3.101405
 [49]  15.221520   0.586872   2.721004   8.229042   0.345450   1.371930
 [55]   0.364458   3.715029   2.450108   0.559980   6.304365   0.098021
 [61]   2.905518   1.711908   5.424055   1.349985   5.960832   4.616704
 [67]   0.128016   7.024500   1.805750   0.945624   0.330132   1.550010
 [73]   0.994620   1.156680   0.463386   2.202057   0.264654   1.507420
 [79]   0.006555   0.225435   5.024432   0.665112   0.989010   0.779350
 [85]   3.460968   0.739584   0.357273   0.005565  11.160240   0.433836
 [91]   0.604095   0.741744   1.121100   5.623200   3.851694   0.066542
 [97]   0.321264   0.004560   1.951205   0.834861   2.010660   2.737644
[103]   0.671580   0.192240   1.801360   0.375144   0.774387   1.191530
[109]   0.003570   0.006972   0.547883   2.899233   0.456840   2.742880
[115]   0.690144   1.027026   2.048200   0.632700   0.419025   1.704331
[121]   1.742364   3.961800   0.802655   1.717065   0.201756   1.482978
[127]   0.765006   0.002145   0.139360   0.514080   0.285138   0.491660
[133]   0.309270   0.125670   0.242962   2.031537   0.169176   1.670900
[139]   0.172260   0.291924   0.294892   1.074060   0.142800   0.161700
[145]   0.244608   0.104904   1.277742   0.148005   0.082170   0.002838
[151]   1.408680   0.316848   0.123000   0.418860   0.869934   0.672771
[157]   1.772892   0.000630   0.000595   0.053295   0.532224   0.674064
[163]   0.000465   0.044805   0.303282   0.092610   0.290979   0.911300
[169]   0.209700   0.000276   0.030107   1.014783   0.577080   0.253460
[175]   0.348156   0.291618   0.252688   0.140760   0.455910   0.825006
[181]   0.119418   0.035970   0.542025   0.063315   0.038124   0.060676
[187]   0.231000   0.069975   0.000010   0.047394   0.000003   0.004212

$parameter.count
[1] 192

$epsilon
[1] 0

$logL
[1] 1408.966

$logL.AICc
[1] NA

attr(,"class")
[1] "skyline"
$time
 [1] 0.021161 0.049822 0.061847 0.077168 0.089216 0.101878 0.113986 0.127902
 [9] 0.140307 0.155704 0.176701 0.192334 0.209112

$interval.length
 [1] 0.021161 0.028661 0.012025 0.015321 0.012048 0.012662 0.012108 0.013916
 [9] 0.012405 0.015397 0.020997 0.015633 0.016778

$population.size
 [1] 392.07100800 525.52809600  26.43747675  18.16241385   4.32071145
 [6]   2.19342354   1.06974257   0.55211856   0.27727433   0.33138171
[11]   0.38060475   0.09827875   0.02431880

$parameter.count
[1] 13

$epsilon
[1] 0.0119

$logL
[1] 1239.478

$logL.AICc
[1] 1225.456

attr(,"class")
[1] "skyline"
Searching for the optimal epsilon... epsilon = 0.01191938 
[1] 0.01191938
Searching for the optimal epsilon... epsilon = 0.01191938 
[1] 0.01191938

ape documentation built on April 25, 2021, 9:06 a.m.