tests/tests.r

In mvMORPH: Multivariate Comparative Tools for Fitting Evolutionary Models to Morphometric Data

#library(testthat)
#library(mvMORPH)
# Comment test_check in order to not run tests with R CMD check
# test_check("mvMORPH")

# test-suite for the functions mvBM/mvOU/mvEB/mvSHIFT/mvLL/mvSIM/mvOUTS/mvRWTS are under preparation

tests <- FALSE

# if(tests==TRUE){
#     require(mvMORPH)
#     # General tree to use in fit and simulations
#     set.seed(1)
#     
#     # Generating a random tree
#     #tree<-phylo<-pbtree(n=50)
#     # Setting the regime states of tip species
#     #sta<-as.vector(c(rep("Forest",20),rep("Savannah",30))); names(sta)<-tree$tip.label
#     # Making the simmap tree with mapped states
#     #tree<-make.simmap(tree,sta , model="ER", nsim=1)
#     
#     # for reproducibility
#     tree <- read.simmap(text="((t37:{Forest,0.17911043},t38:{Forest,0.17911043}):{Forest,4.19021273},(((t39:{Forest,0.1406569},t40:{Forest,0.1406569}):{Forest,0.07210538},t36:{Forest,0.21276228}):{Forest,3.77896996},((((((t21:{Forest,0.474973},t22:{Forest,0.474973}):{Forest,0.11994112},(t28:{Forest,0.3176896},(t41:{Forest,0.11335236},t42:{Forest,0.11335236}):{Forest,0.20433723}):{Forest,0.27722453}):{Forest,0.34242237},t8:{Forest,0.9373365}):{Forest,0.87640728},(((t12:{Forest,0.73431592},t13:{Forest,0.73431592}):{Forest,0.69159925},(t2:{Forest,1.31658647},(t18:{Forest,0.57994406},t19:{Forest,0.57994406}):{Forest,0.73664241}):{Forest,0.1093287}):{Forest,0.23284738},(t1:{Forest,1.54484389},((t16:{Forest,0.59318045},t17:{Forest,0.59318045}):{Forest,0.33458062},(t34:{Forest,0.25021113},t35:{Savannah,0.17781996:Forest,0.07239117}):{Forest,0.67754994}):{Forest,0.61708282}):{Forest,0.11391867}):{Forest,0.15498123}):{Forest,0.33873822},(t4:{Savannah,1.234117},t5:{Savannah,1.234117}):{Savannah,0.0748023:Forest,0.8435627}):{Forest,1.17266697},(((t3:{Savannah,1.25848874},((t14:{Savannah,0.60935948},t15:{Savannah,0.60935948}):{Savannah,0.08390608},(t20:{Savannah,0.57975625},((t24:{Savannah,0.32641206},(t47:{Savannah,0.01395592},t48:{Savannah,0.01395592}):{Savannah,0.31245613}):{Savannah,0.19000652},(t29:{Savannah,0.30478906},t30:{Savannah,0.30478906}):{Savannah,0.21162951}):{Savannah,0.06333768}):{Savannah,0.1135093}):{Savannah,0.56522318}):{Savannah,0.47279525},((((t31:{Savannah,0.28631849},(t45:{Savannah,0.01659224},t46:{Savannah,0.01659224}):{Savannah,0.26972625}):{Savannah,0.15599405},t23:{Savannah,0.44231253}):{Savannah,0.33983399},t11:{Savannah,0.78214652}):{Savannah,0.11039434},(t25:{Savannah,0.32607817},(t49:{Savannah,0.00334081},t50:{Savannah,0.00334081}):{Savannah,0.32273736}):{Savannah,0.5664627}):{Savannah,0.83874312}):{Savannah,0.87012285},(((t26:{Savannah,0.32102325},t27:{Savannah,0.32102325}):{Savannah,0.59890242},t9:{Savannah,0.91992568}):{Savannah,1.43882358},(((t6:{Savannah,1.06510523},t7:{Savannah,1.06510523}):{Savannah,0.02627092},(t43:{Savannah,0.02516062},t44:{Savannah,0.02516062}):{Savannah,1.06621553}):{Savannah,0.79310082},((t32:{Savannah,0.27743011},t33:{Savannah,0.27743011}):{Savannah,0.5706055},t10:{Savannah,0.84803561}):{Savannah,1.03644137}):{Savannah,0.47427228}):{Savannah,0.24265758}):{Savannah,0.19050932:Forest,0.53323281}):{Forest,0.66658327}):{Forest,0.37759092});")
#     
#     col<-c("blue","orange"); names(col)<-c("Forest","Savannah")
#     # Plot of the phylogeny for illustration
#     plotSimmap(tree,col,fsize=0.6,node.numbers=FALSE,lwd=3, pts=FALSE)
#     
#     # for reproducibility without seed
#     #trait <-
#    
#     # General time-series
#     timeseries <- 0:49
#     ## tests on mvSIM
#     
#     ## Bivariate case
#     # generate sigma
#     sigma <- matrix(m <- runif(4), 2, 2) %*% t(matrix(m, 2, 2))
#     sigbmm <- list(matrix(m <- runif(4), 2, 2) %*% t(matrix(m, 2, 2)),
#                     matrix(m <- runif(4), 2, 2) %*% t(matrix(m, 2, 2)))
#     # generate alpha
#     alpha <- matrix(m <- runif(4), 2, 2) %*% t(matrix(m, 2, 2))
#     # generate theta
#     ou1 <- bm <- c(0, 2)
#     oum <- bmD <- matrix(c(0,1,0.5,2),2,2) # traits are on columns
#     # generate beta
#     beta <- -0.05
#     # generate trend
#     trend <- c(0.02, 0.02)
#     
#     # param BM1
#     parbm1 <- list(sigma=sigma, theta=bm)
#     # param BM1 trend
#     parbm1t <- list(sigma=sigma, theta=bm, trend=trend)
#     # param BM1 multiple means
#     parbm1D <- list(sigma=sigma, theta=bmD, smean=FALSE)
#     # param BMM
#     parbmm <- list(sigma=sigbmm, theta=bm)
#     # param BMM trend
#     parbmmt <- list(sigma=sigbmm, theta=bm, trend=trend)
#     # param BMM multiple means
#     parbmmD <- list(sigma=sigbmm, theta=bmD, smean=FALSE)
#     
#     # param OU1
#     parou1 <- list(sigma=sigma, alpha=alpha, theta=ou1)
#     # param OUM
#     paroum <- list(sigma=sigma, alpha=alpha, theta=oum)
#     # param OUM
#     paroumR <- list(sigma=sigma, alpha=alpha, theta=ou1, root=FALSE)
#     
#     # param EB
#     pareb <- list(sigma=sigma, theta=bm, beta=beta)
#     
#     # Simulate datasets
#     dat1 <- mvSIM(tree, model="BM1", nsim=1, param=parbm1)
#      dat2 <- mvSIM(tree, model="BM1", nsim=1, param=parbm1t)
#       dat3 <- mvSIM(tree, model="BM1", nsim=1, param=parbm1D)
#        dat4 <- mvSIM(tree, model="BMM", nsim=1, param=parbmm)
#         dat5 <- mvSIM(tree, model="BMM", nsim=1, param=parbmmt)
#          dat6 <- mvSIM(tree, model="BMM", nsim=1, param=parbmmD)
#           dat7 <- mvSIM(tree, model="OU1", nsim=1, param=parou1)
#            dat8 <- mvSIM(tree, model="OUM", nsim=1, param=paroum)
#             dat9 <- mvSIM(timeseries, model="OUTS", nsim=1, param=paroum)
#              dat10 <- mvSIM(timeseries, model="OUTS", nsim=1, param=paroumR)
#               dat11 <- mvSIM(tree, model="EB", nsim=1, param=pareb)
#              
#     # Intentional errors
#     test1 <- mvSIM(tree, model="OUTS", nsim=1, param=paroum) # tree instead of ts
#      test2 <- mvSIM(timeseries, model="OU1", nsim=1, param=parou1) # ts instead of tree
#       test3 <- mvSIM(timeseries, model="OUTS", nsim=1, param=parou1) # wrong number of theta
#        test4 <- mvSIM(model="OUTS", nsim=1, param=paroum) # missing ts/tree
#        
#        
#        ## Model fit
#        
#          fit1 <- mvBM(tree, dat1, model="BM1")
#            fit2 <- mvBM(tree, dat2, model="BM1", param=list(trend=T)) # must be fitted on a n-ultrametric tree
#              fit3 <- mvBM(tree, dat3, model="BM1", param=list(smean=F))
#                fit4 <- mvBM(tree, dat4, model="BMM")
#                  fit5 <- mvBM(tree, dat5, model="BMM", param=list(trend=T)) # must be fitted on a n-ultrametric tree
#                    fit6 <- mvBM(tree, dat6, model="BMM", param=list(smean=F))
#                      fit7 <- mvOU(tree, dat7, model="OU1")
#                      fit8 <- mvOU(tree, dat8, model="OUM", param=list(vcv="randomRoot")) # avoid large values sometimes; known problem to fix?
#                          fit9 <- mvOUTS(timeseries, dat9)
#                            fit10 <- mvOUTS(timeseries, dat10, param=list(root=F))
#                              fit11 <- mvRWTS(timeseries, dat10)
#                                fit11 <- mvEB(tree, dat11)
#                                
#         # Tests for errors
#         test5 <- mvBM(tree, dat2, model="BM1", param=list(trend=c(1,1,1,1), smean=F))
#         
#         ## Model fit with NA values and defaults methods
#         
#         dat1[8,2]<- dat1[25,1] <- NA
#           dat2[8,2]<- dat2[25,1] <- NA
#             dat3[8,2]<- dat3[25,1] <- NA
#               dat4[8,2]<- dat4[25,1] <- NA
#                 dat5[8,2]<- dat5[25,1] <- NA
#                   dat6[8,2]<- dat6[25,1] <- NA
#                     dat7[8,2]<- dat7[25,1] <- NA
#                       dat8[8,2]<- dat8[25,1] <- NA
#                         dat9[8,2]<- dat9[25,1] <- NA
#                           dat10[8,2]<- dat10[25,1] <- NA
#                           # then fit
#                           fit1 <- mvBM(tree, dat1, model="BM1")
#                           fit2 <- mvBM(tree, dat2, model="BM1", param=list(trend=T))
#                           fit3 <- mvBM(tree, dat3, model="BM1", param=list(smean=F))
#                           fit4 <- mvBM(tree, dat4, model="BMM")
#                           fit5 <- mvBM(tree, dat5, model="BMM", param=list(trend=T))
#                           fit6 <- mvBM(tree, dat6, model="BMM", param=list(smean=F))
#                           fit7 <- mvOU(tree, dat7, model="OU1")
#                           fit8 <- mvOU(tree, dat8, model="OUM")
#                           fit9 <- mvOUTS(timeseries, dat9)
#                           fit10 <- mvOUTS(timeseries, dat10, param=list(root=F))
#                           fit11 <- mvRWTS(timeseries, dat10)
#                           fit11 <- mvEB(tree, dat11)
# 
# 
# # Expected likelihood
# 
# # for reproducibility (dat1 converted to vector)
# traits <- matrix(c(0.759912758269539, 0.614545055791663, 1.99523967616929, 2.06050645201129, 1.72847166301803, -2.74185937878331, -1.74984543950876, -2.06239818030884, -2.06041993850937, -1.79080629892104, -1.20081330679844, -0.54349593713857, -0.0607971864579239, -0.00233277135380217, -0.400693887343343, -1.65772161485507, -0.182544756090116, -0.654109529770913, -0.737499486430531, -1.72515976761732, -1.72562879593648, -0.25330550757395, 0.729897907633583, -0.246081522346491, 0.151078889032263, 0.0138479096177918, -0.848009457276881, -0.581243223593088, -0.667057834272238, -0.66847007470311, 0.0981507462333929, 0.172861363132646, -0.353171766681764, -0.451300105973155, -0.368139106479026, -0.0627361956502809, -0.709892073985014, -0.875683848248537, -0.570848855977929, -0.532265981711379, -0.369679020727696, 0.0776255560504741, 0.138142231533873, -0.833365792498594, -0.309928572800132, -1.35885805922583, -1.14218922343049, 1.48087206844568, 1.09136693791236, -0.0180692260844691, 3.23281192282276, 3.00136528638948, 5.32074612474669, 5.41402685056098, 4.913298892713, -2.23542682300963, -0.726407146095314, -1.1868584417652, -1.24452289327373, -0.784150890488637, 0.146694082761811, 1.32180581471685, 2.05153383820109, 2.03959665464931, 1.48623453081428, -0.518344824420335, 1.70143196101035, 1.04697657897734, 0.892292938362765, -0.630853524174634, -0.629376986734026, 1.6045525143076, 3.11635168256392, 1.62949197480133, 2.32593600511521, 2.02205029936445, 0.748804111684117, 1.16117939166309, 1.05667448503384, 1.05133732067517, 2.22897875073374, 2.34707599071464, 1.45209863662888, 1.34341733329564, 1.48063807048236, 1.95499546327074, 1.0100600946259, 0.727548824795377, 1.15112488908316, 1.20921909380616, 1.36592576262338, 2.05376926747503, 2.14226651915817, 0.712214997646626, 1.47055230647709, -0.116941097964031, 0.212904236741522, 4.41536257330546, 3.8195839768853, 2.09983331019748), ncol=2)
# 
# # univariate:
# test6 <- mvLL(tree, traits[,1], method="pic")
# 
#     #$logl
#     #[1] -37.28764
#     #
#     #$theta
#     #[1] 0.3443453
#     #
#     #$sigma
#     #[,1]
#     #[1,] 0.3028476
#     #
#     #attr(,"class")
#     #[1] "mvmorph" "loglik"
# 
#  test7 <- mvLL(tree, traits[,1], method="rpf")
#   test8 <- mvLL(tree, traits[,1], method="sparse")
#    test9 <- mvLL(tree, traits[,1], method="inverse")
#     test10 <- mvLL(tree, traits[,1], method="pseudoinverse")
#     
#     all.equal(test6$loglik, test7$loglik)
#     all.equal(test6$loglik, test8$loglik)
#     all.equal(test6$loglik, test9$loglik)
#     all.equal(test6$loglik, test10$loglik)
#     
# # multivariate:
#     test11 <- mvLL(tree, traits, method="pic")
#     
#     #$logl
#     #[1] 56.40704
#     #
#     #$theta
#     #[1] 0.3443453 2.6059401
#     #
#     #$sigma
#     #[,1]      [,2]
#     #[1,] 0.3028476 0.4721665
#     #[2,] 0.4721665 0.7377560
#     #
#     #attr(,"class")
#     #[1] "mvmorph" "loglik"
#     
#     test12 <- mvLL(tree, traits, method="pic", param=list(estim=FALSE, mu=c(1,1)))
#     
#     #$logl
#     #[1] -1983.862
#     #
#     #$theta
#     #[1] 1 1
#     #
#     #attr(,"class")
#     #[1] "mvmorph" "loglik"
#     
#     test13 <- mvLL(tree, traits, method="pic", param=list(estim=FALSE, mu=c(1,1), sigma=matrix(c(2,1,1,1.5),2,2)))
#     
#     #$logl
#     #[1] -115.824
#     #
#     #$theta
#     #[1] 1 1
#     #
#     #attr(,"class")
#     #[1] "mvmorph" "loglik"
# 
# ## ----------------- 3 traits ------------------- ##
# 
# ## Bivariate case
# # generate sigma
# sigma <- matrix(m <- runif(9), 3, 3) %*% t(matrix(m, 3, 3))
# sigbmm <- list(matrix(m <- runif(9), 3, 3) %*% t(matrix(m, 3, 3)),
# matrix(m <- runif(9), 3, 3) %*% t(matrix(m, 3, 3)))
# # generate alpha
# alpha <- matrix(m <- runif(9), 3, 3) %*% t(matrix(m, 3, 3))
# # generate theta
# ou1 <- bm <- c(0, 2, 1)
# oum <- bmD <- matrix(c(0,1,0,1.5,0.5,2),ncol=3, nrow=2) # traits are on columns
# # generate beta
# beta <- -0.05
# # generate trend
# trend <- c(0.02, 0.02, 0.02)
# 
# # param BM1
# parbm1 <- list(sigma=sigma, theta=bm)
# # param BM1 trend
# parbm1t <- list(sigma=sigma, theta=bm, trend=trend)
# # param BM1 multiple means
# parbm1D <- list(sigma=sigma, theta=bmD, smean=FALSE)
# # param BMM
# parbmm <- list(sigma=sigbmm, theta=bm)
# # param BMM trend
# parbmmt <- list(sigma=sigbmm, theta=bm, trend=trend)
# # param BMM multiple means
# parbmmD <- list(sigma=sigbmm, theta=bmD, smean=FALSE)
# 
# # param OU1
# parou1 <- list(sigma=sigma, alpha=alpha, theta=ou1)
# # param OUM
# paroum <- list(sigma=sigma, alpha=alpha, theta=oum)
# # param OUM
# paroumR <- list(sigma=sigma, alpha=alpha, theta=ou1, root=FALSE)
# 
# # param EB
# pareb <- list(sigma=sigma, theta=bm, beta=beta)
# 
# # Simulate datasets
# dat1 <- mvSIM(tree, model="BM1", nsim=1, param=parbm1)
#  dat2 <- mvSIM(tree, model="BM1", nsim=1, param=parbm1t)
#   dat3 <- mvSIM(tree, model="BM1", nsim=1, param=parbm1D)
#    dat4 <- mvSIM(tree, model="BMM", nsim=1, param=parbmm)
#     dat5 <- mvSIM(tree, model="BMM", nsim=1, param=parbmmt)
#      dat6 <- mvSIM(tree, model="BMM", nsim=1, param=parbmmD)
#       dat7 <- mvSIM(tree, model="OU1", nsim=1, param=parou1)
#        dat8 <- mvSIM(tree, model="OUM", nsim=1, param=paroum)
#         dat9 <- mvSIM(timeseries, model="OUTS", nsim=1, param=paroum)
#          dat10 <- mvSIM(timeseries, model="OUTS", nsim=1, param=paroumR)
#           dat11 <- mvSIM(tree, model="EB", nsim=1, param=pareb)
# 
# # Intentional errors
# test1 <- mvSIM(tree, model="OUTS", nsim=1, param=paroum) # tree instead of ts
#     test2 <- mvSIM(timeseries, model="OU1", nsim=1, param=parou1) # ts instead of tree
#         test3 <- mvSIM(timeseries, model="OUTS", nsim=1, param=parou1) # wrong number of theta
#             test4 <- mvSIM(model="OUTS", nsim=1, param=paroum) # missing ts/tree
# 
# 
# ## Model fit
# 
# fit1 <- mvBM(tree, dat1, model="BM1")
#  fit2 <- mvBM(tree, dat2, model="BM1", param=list(trend=T))
#   fit3 <- mvBM(tree, dat3, model="BM1", param=list(smean=F))
#    fit4 <- mvBM(tree, dat4, model="BMM")
#     fit5 <- mvBM(tree, dat5, model="BMM", param=list(trend=T))
#      fit6 <- mvBM(tree, dat6, model="BMM", param=list(smean=F))
#       fit7 <- mvOU(tree, dat7, model="OU1")
#        fit8 <- mvOU(tree, dat8, model="OUM")
#         fit9 <- mvOUTS(timeseries, dat9)
#          fit10 <- mvOUTS(timeseries, dat10, param=list(root=F))
#           fit11 <- mvRWTS(timeseries, dat10)
#            fit11 <- mvEB(tree, dat11)
#            
#            fit12 <- mvOU(tree, dat8, model="OUM", param=list(vcv="randomRoot"))
#            fit13 <- mvOU(tree, dat8, model="OUM", param=list(vcv="randomRoot", root=TRUE))
#            fit14 <- mvOU(tree, dat8, model="OUM", param=list(vcv="randomRoot", alpha=alpha)) # error
#            fit15 <- mvOU(tree, dat8, model="OUM", param=list(vcv="randomRoot", alpha=alpha, decomp="qr")) # error
#            fit16 <- mvOU(tree, dat8, model="OUM", param=list(vcv="randomRoot", sigma=sigma)) # error
#            fit17 <- mvBM(tree, dat1, model="BM1", param=list(sigma=sigma))
#            fit18 <- mvBM(tree, dat1, model="BMM", param=list(sigma=sigbmm))
#            fit19 <- mvOU(tree, dat8, model="OUM", param=list(vcv="randomRoot", alpha=runif(6), decomp="choles")) # error for the name
#            fit20 <- mvOUTS(timeseries, dat9, param=list(vcv="randomRoot", alpha=alpha, decomp="qr"))
#            fit21 <- mvOUTS(timeseries, dat9, param=list(vcv="randomRoot", decomp="qr", alpha=alpha, decompSigma="eigen+"))
# 
# ## Model fit with NA values and defaults methods
# 
# dat1[8,2]<- dat1[25,1] <- NA
# dat2[8,2]<- dat2[25,1] <- NA
# dat3[8,2]<- dat3[25,1] <- NA
# dat4[8,2]<- dat4[25,1] <- NA
# dat5[8,2]<- dat5[25,1] <- NA
# dat6[8,2]<- dat6[25,1] <- NA
# dat7[8,2]<- dat7[25,1] <- NA
# dat8[8,2]<- dat8[25,1] <- NA
# dat9[8,2]<- dat9[25,1] <- NA
# dat10[8,2]<- dat10[25,1] <- NA
# # then fit
# fit1 <- mvBM(tree, dat1, model="BM1")
#     fit2 <- mvBM(tree, dat2, model="BM1", param=list(trend=T))
#         fit3 <- mvBM(tree, dat3, model="BM1", param=list(smean=F))
#             fit4 <- mvBM(tree, dat4, model="BMM")
#                 fit5 <- mvBM(tree, dat5, model="BMM", param=list(trend=T))
#                     fit6 <- mvBM(tree, dat6, model="BMM", param=list(smean=F))
#                         fit7 <- mvOU(tree, dat7, model="OU1")
#                             fit8 <- mvOU(tree, dat8, model="OUM")
#                                 fit9 <- mvOUTS(timeseries, dat9)
#                                     fit10 <- mvOUTS(timeseries, dat10, param=list(root=F))
#                                         fit11 <- mvRWTS(timeseries, dat10)
#                                             fit12 <- mvEB(tree, dat11)
#                                               fit13 <- mvRWTS(timeseries, dat10, param=list(trend=TRUE))
# 
# 
# # simulate generic
# d1 <- simulate(fit1,nsim=1, tree=tree)
#     d2 <- simulate(fit2,nsim=1, tree=tree)
#         d3 <- simulate(fit3,nsim=1, tree=tree)
#             d4 <- simulate(fit4,nsim=1, tree=tree)
#                 d5 <- simulate(fit5,nsim=1, tree=tree)
#                     d6 <- simulate(fit6,nsim=1, tree=tree)
#                         d7 <- simulate(fit7,nsim=1, tree=tree)
#                             d8 <- simulate(fit8,nsim=1, tree=tree)
#                                 d9 <- simulate(fit9,nsim=1, tree=timeseries)
#                                     d10 <- simulate(fit10,nsim=1, tree=timeseries)
#                                         d11 <- simulate(fit11,nsim=1, tree=timeseries)
#                                             d12 <- simulate(fit12,nsim=1, tree=tree)
#                                                 d13 <- simulate(fit13,nsim=1, tree=timeseries)
# 
# } # end of tests