forked from leonelhalsina/lemad
-
Notifications
You must be signed in to change notification settings - Fork 0
/
test_lemad_geosse
80 lines (68 loc) · 2.59 KB
/
test_lemad_geosse
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
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
lemad_test_geosse <- function(){
#test geosse
pars <- c(1.5, 0.5, 1.0, 0.7, 0.7, 2.5, 0.5)
names(pars) <- c("sA", "sB", "sAB", "xA", "xB", "dA", "dB")
#set.seed(5)
#phy <- diversitree::tree.geosse(pars, max.t=4, x0=0)
phy <- NULL; rm(phy);
utils::data('example_phy_GeoSSE', package = 'lemad');
traits <- as.numeric(phy$tip.state)
testit::assert(!is.null(phy))
lik.g <- diversitree::make.geosse(phy, phy$tip.state)
pars.g <- c(1.5, 0.5, 1.0, 0.7, 0.7, 1.4, 1.3)
names(pars.g) <- diversitree::argnames(lik.g)
lik.c <- diversitree::make.classe(phy, phy$tip.state+1, 3)
pars.c <- 0 * diversitree::starting.point.classe(phy, 3)
pars.c['lambda222'] <- pars.c['lambda112'] <- pars.g['sA']
pars.c['lambda333'] <- pars.c['lambda113'] <- pars.g['sB']
pars.c['lambda123'] <- pars.g['sAB']
pars.c['mu2'] <- pars.c['q13'] <- pars.g['xA']
pars.c['mu3'] <- pars.c['q12'] <- pars.g['xB']
pars.c['q21'] <- pars.g['dA']
pars.c['q31'] <- pars.g['dB']
lik.g(pars.g) # -175.7685
classe_diversitree_LL <- lik.c(pars.c) # -175.7685
## lemad part
traits_coded_letters <- NULL
for(ii in 1:length(traits)){
if(traits[ii] == 0){
traits_coded_letters <- c(traits_coded_letters,"A")
}
if(traits[ii] == 1){
traits_coded_letters <- c(traits_coded_letters,"B")
}
if(traits[ii] == 2){
traits_coded_letters <- c(traits_coded_letters,"AB")
}
}
lambdas<-list()
lambdas[[1]] <- matrix(0,ncol = 3,nrow = 3,byrow = TRUE)
#lambdas[[1]][1,1] <- 1.5
lambdas[[1]][2,1] <- 1.5
lambdas[[1]][3,1] <- 0.5
lambdas[[1]][3,2] <- 1
lambdas[[2]]<-matrix(0,ncol = 3,nrow = 3,byrow = TRUE)
lambdas[[2]][2,2] <- 1.5
#lambdas[[2]][2,2] <- 1.1
lambdas[[3]] <- matrix(0,ncol=3,nrow=3,byrow=TRUE)
lambdas[[3]][3,3] <- 0.5
#lambdas[[3]][3,3] <- 1
mus<-c(0,0.7,0.7)
q<-matrix(0,ncol = 3,nrow = 3,byrow = TRUE)
q[2,1] <- 1.4
q[3,1] <- 1.3
q[1,2] <- 0.7
q[1,3] <- 0.7
parameter <- list()
parameter[[1]] <- lambdas
parameter[[2]] <- mus
parameter[[3]] <- q
lemad_LL <- cla_lemad_loglik(parameter, phy, traits_coded_letters,
num_max_multiregion = 2,
use_fortran = TRUE, methode = "ode45", cond = "maddison_cond",
root_state_weight = "maddison_weights", sampling_fraction = c(1,1,1),
run_parallel = FALSE, setting_calculation = NULL,
setting_parallel = NULL, see_ancestral_states = TRUE,
loglik_penalty = 0)
testthat::expect_equal(classe_diversitree_LL,lemad_LL$LL)
}