R Under development (unstable) (2024-08-26 r87056 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #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
> 
> 
> proc.time()
   user  system elapsed 
   0.10    0.03    0.14