epcTools example

epcTools - a package for incorporating environment into parameters of Ornstein-Uhlenbeck models

First, lets install epc tools if you don’t have it already.

Prepare a simulated example

We will need to simulate a phylogeny and environmental vector. We can then use these two to simulate an EPC cache object which will vary with the trait of interest.

The sim.EPC function requires several arguments. First, the tree and environmental values but also including other parameters documented elsewhere such as the type of EPC relationship and what variable(s) are affected by the environment. Here we simulate a cache object with alpha and theta varying with the environment and ten time slices.

#Tree and environment
tree_example<-simTree(100,200,100,1)
env_example<-(c(6,8,11,13,9,8,5,4,8,10))

#Parameters to simulate under
base_sig2=0.5
base_b0_t=1
base_b1_t=2
base_b0_a=0.01
base_b1_a=0.03
base_slice=10


sim0.5f_200b_at_linear<-sim.EPC(tree_example,epc_params=c("alpha","theta"),m_type="linear",X=env_example,n_slice=base_slice,1,b0_a=base_b0_a,b1_a=base_b1_a,sig2=base_sig2,b0_t=base_b0_t,b1_t=base_b1_t)

Visualize how the traits are varying with the environment

#EPC phenogram
simEPC.phenogram(sim0.5f_200b_at_linear[[1]],10)

epcTools maximum likelihood search

Now, lets see if we can recover the parameters we simulated our model under. This requires invoking the epc.max.lik() function. Note that this function can take a vector of starting parameters fed by the user or find one automatically using the start.searcher function.

#EPC phenogram
example_results<-epc.max.lik(sim0.5f_200b_at_linear[[1]])

## [1] "Finding starting parameters..."
##         b0t         b1t     sigma        b0a          b1a      lik
## 36 16.08379 -0.01917567 0.9961479 0.05208001 -0.003969108 462.6231
## [1] "Finding maximum likelihood..."
## [1] "Finding maximum likelihood..."
## [1] "Finding starting parameters..."
##          b0t       b1t    sigma        b0a         b1a      lik
## 164 17.10644 0.1719627 1.058908 0.01154724 0.006784836 441.0758
## [1] "Finding maximum likelihood..."
## [1] "Finding maximum likelihood..."
## [1] "Finding starting parameters..."
##         b0t       b1t    sigma       b0a         b1a      lik
## 30 12.41576 0.7752388 1.032144 0.2057524 0.004276592 425.4583
## [1] "Finding maximum likelihood..."
## [1] "Finding maximum likelihood..."
## [1] "Finding starting parameters..."
##       b0t       b1t    sigma        b0a          b1a      lik
## 1 14.0471 0.2581318 0.981456 0.07293743 -0.004093211 460.6198
## [1] "Finding maximum likelihood..."
## [1] "Finding maximum likelihood..."
## [1] "Calculating intervals at a confidence level of 95%"
## [1] "Done replicate 500"
## [1] "CI of values (the 35 replicates within 5.53524884675818 neglnL of the optimum)"
##        neglnL parameter_1 parameter_2 parameter_3 parameter_4 parameter_5
## [1,] 204.3724    1.034274    1.938498   0.4621982  0.03814634  0.02275434
## [2,] 209.8057    1.565158    2.007553   0.5933609  0.05507022  0.02814532
## [1] "Rough volume of good region is 4.38708128236731e-07"
## [1] "Done replicate 1000"
## [1] "CI of values (the 120 replicates within 5.53524884675818 neglnL of the optimum)"
##        neglnL parameter_1 parameter_2 parameter_3 parameter_4 parameter_5
## [1,] 204.3724    1.034274    1.938498   0.4621982  0.03814634  0.02275434
## [2,] 209.8057    1.565158    2.007553   0.5933609  0.05507022  0.02866196
## [1] "Rough volume of good region is 4.80751423769994e-07"
## [1] "Done replicate 1500"
## [1] "CI of values (the 305 replicates within 5.53524884675818 neglnL of the optimum)"
##        neglnL parameter_1 parameter_2 parameter_3 parameter_4 parameter_5
## [1,] 204.3724    1.034274    1.936741   0.4621982  0.03764241  0.02275434
## [2,] 209.9062    1.565158    2.007553   0.6219078  0.05507022  0.02882251
## [1] "Rough volume of good region is 6.34951859436322e-07"
## [1] "Done replicate 2000"
## [1] "CI of values (the 548 replicates within 5.53524884675818 neglnL of the optimum)"
##        neglnL parameter_1 parameter_2 parameter_3 parameter_4 parameter_5
## [1,] 204.3724    1.034274    1.936741   0.4621982  0.03764241  0.02275434
## [2,] 209.9062    1.565158    2.007553   0.6219078  0.05507022  0.02882251
## [1] "Rough volume of good region is 6.34951859436322e-07"
## [1] "Done replicate 2500"
## [1] "CI of values (the 707 replicates within 5.53524884675818 neglnL of the optimum)"
##        neglnL parameter_1 parameter_2 parameter_3 parameter_4 parameter_5
## [1,] 204.3724    1.034274    1.936741   0.4621982  0.03550704  0.02275434
## [2,] 209.9062    1.565158    2.007553   0.6236128  0.05507022  0.02882251
## [1] "Rough volume of good region is 7.20358867206209e-07"

example_results

## $fit
##   EPC_model   Type      Lik Npar      AIC
## 1     alpha linear 204.3724    5 418.7447
## 2     theta linear 204.3724    5 418.7447
## 
## $dentist_output
## This ran 2500 steps looking for all points within 5.53524884675818 negative log likelihood units of the best parameter values.
## 
## Parameters: 
##                       b0_t     b1_t      sig2       b0_a       b1_a
## best             1.1943484 1.983004 0.5164584 0.04357437 0.02452023
## lower.CI         1.0342736 1.936741 0.4621982 0.03550704 0.02275434
## upper.CI         1.5651577 2.007553 0.6236128 0.05507022 0.02882251
## lowest.examined  0.9178561 1.474651 0.3431824 0.03151233 0.01873364
## highest.examined 1.6535847 2.437712 0.6424223 0.06241387 0.03293040
## 
## $search_results
##                    b0_t       b1_t      sig2        b0_a        b1_a likelihood
## BFGS3         1.1943484  1.9830039 0.5164584  0.04357437  0.02452023   204.3724
## BFGS          1.1911752  1.9832842 0.5162067  0.04304112  0.02458740   204.3728
## Nelder-Mead   0.5601993  2.0399975 0.5324286 -0.04912843  0.03649362   206.4201
## Nelder-Mead3  2.1562557  1.9136223 0.5384032  0.33654694 -0.01248023   219.1869
## Nelder-Mead2 12.8115580  0.7455699 1.0863637 -0.21417352  0.05354339   358.8136
## BFGS2        12.8115580  0.7455699 1.0863637 -0.21417352  0.05354339   358.8136
## BFGS1        22.1763587 -0.1099564 0.8909789 -0.05357720  0.01339436   425.1057
## Nelder-Mead1 22.1763587 -0.1099564 0.8909789 -0.05357720  0.01339436   425.1057
##               xtime      AIC
## BFGS3        15.753 418.7447
## BFGS         13.462 418.7455
## Nelder-Mead  35.253 422.8402
## Nelder-Mead3 35.504 448.3739
## Nelder-Mead2 29.840 727.6273
## BFGS2         2.717 727.6273
## BFGS1         2.628 860.2114
## Nelder-Mead1 31.121 860.2114
## 
## attr(,"class")
## [1] "epc"  "list"

Dentist parameter visualization

#EPC phenogram
plot(example_results)