Nothing
################################################################################
## ##
## mvMORPH: mvSHIFT ##
## ##
## Created by Julien Clavel - 16-07-2013 ##
## (julien.clavel@hotmail.fr/ julien.clavel@univ-lyon1.fr) ##
## require: phytools, ape, corpcor, subplex ##
## ##
################################################################################
mvSHIFT<-function(tree,data,error=NULL,param=list(age=NULL,sigma=NULL,alpha=NULL,sig=NULL,beta=NULL),model=c("ER","RR","EC","RC","SR","EBOU","OUEB","EBBM","BMEB"), method=c("rpf","sparse","inverse","pseudoinverse"),scale.height=FALSE, optimization=c("L-BFGS-B","Nelder-Mead","subplex"),control=list(maxit=20000),precalc=NULL,diagnostic=TRUE,echo=TRUE){
if(missing(tree)) stop("The tree object is missing!")
if(missing(data)) stop("You must provide a dataset along with your tree!")
#set data as a matrix if a vector is provided instead
if(!is.matrix(data)){data<-as.matrix(data)}
# number of species (tip)
n<-dim(data)[1]
# number of variables
p<-dim(data)[2]
# Check the order of the dataset and the phylogeny
if(!is.null(rownames(data))) {
if(any(tree$tip.label==rownames(data))){
data<-data[tree$tip.label,]
}else if(echo==TRUE){
cat("row names of the data matrix must match tip names of your phylogeny!","\n")
}
}else if(echo==TRUE){
cat("species in the matrix are assumed to be in the same order as in the phylogeny, otherwise specify rownames of 'data'","\n")
}
# choose method for the likelihood computation
method=method[1]
if(method!="rpf" & method!="inverse" & method!="pseudoinverse" & method!="sparse"){
stop("Please choose either the \"rpf\", \"inverse\", \"pseudoinverse\", or \"sparse\" method")
}
# Check if there is missing cases
NA_val<-FALSE
Indice_NA<-NULL
if(any(is.na(data))){
if(method!="pic" & method!="sparse"){
NA_val<-TRUE
}else{
stop("NA values are allowed only with the \"rpf\",\"inverse\" or \"pseudoinverse\" methods")
}
Indice_NA<-which(is.na(as.vector(data)))
}
# choose model
model=model[1]
# root (current default value)
root=FALSE
# choose a method for the optimizer
optimization=optimization[1]
if(is.null(param[["up"]])){up<-param$up<-2}else{up<-param$up}
if(is.null(param[["low"]])){low<-param$low<- -2}else{low<-param$low}
#set age shift with make.era from phytools if the age is provided
if(!is.null(param[["age"]])){
tot<-max(nodeHeights(tree))
limit=tot-param$age
tree<-make.era.map(tree, c(0,limit))
}
#Shift rate model
if(model=="SR"){
mvBM(tree,data,model="BMM",param=param,scale.height=scale.height,diagnostic=diagnostic,method=method,optimization=optimization,echo=echo,control=control,error=error)
}else{
##---------------------------Precalc-parameters-----------------------------##
if(is.null(precalc)==FALSE & inherits(precalc, "mvmorph.precalc")){
tree<-precalc$tree
if(is.null(tree[["mapped.edge"]])){
stop("The tree is not in SIMMAP format with mapped time slices")
}
if(!is.null(param[["age"]])){
if(echo==TRUE) cat("The tree specified in precalc is assumed to represent a shift at :",param$age)
}
# multiple vcv list
vcvList<-precalc$C2
# weight matrix
W<-precalc$D
# Yale sparse format
JAr<-precalc$JAr
IAr<-precalc$IAr
ch<-precalc$ch
precalcMat<-precalc$V
}else{
# scale tree
if(scale.height==TRUE){
maxHeight<-max(nodeHeights(tree))
tree$edge.length<-tree$edge.length/maxHeight
tree$mapped.edge<-tree$mapped.edge/maxHeight
}
# compute the design matrix
W <-multD(tree,p,n,smean=TRUE)
# compute the variance-covariances matrices for each map
multi.tre<-list()
class(multi.tre)<-"multiPhylo"
vcvList<-vcvSplit(tree)
if(method=="sparse"){
## sparse precalcs
#require(spam)
# Temporary multivariate VCV
V<-kronecker((matrix(1,p,p)+diag(p)), vcv.phylo(tree))
# spam object
precalcMat<-as.spam(V);
# precal the cholesky
if(is.null(param[["pivot"]])){pivot<-"MMD"}else{pivot<-param$pivot}
ch<-chol(precalcMat,pivot=pivot)
# Yale Sparse Format indices
JAr<-precalcMat@colindices-1
IAr<-precalcMat@rowpointers-1
# list precalc
}
}
##---------------------------Define initial values--------------------------##
# Pull (alpha) matrix decomposition
if(is.null(param[["decomp"]])){
decomp<-param$decomp<-"cholesky"
}else{
decomp<-param$decomp[1]
}
# scatter (sigma) matrix decomposition
if(is.null(param[["decompSigma"]])){
decompSigma<-param$decompSigma<-"cholesky"
}else{
decompSigma<-param$decompSigma[1]
}
# initial sigma matrix if not provided
if(is.null(param[["sigma"]])){sigma<-param$sigma<-startParamSigma(p, decompSigma, tree, data)}else{sigma<-param$sigma}
# alpha matrix
if(is.null(param[["alpha"]])){alpha<-param$alpha<-startParam(p,decomp,tree)}
if(!is.numeric(param$alpha[[1]])){
if(length(param$alpha)==2){
if(is.numeric(param$alpha[[2]])){
alpha=param$alpha[[2]]}}else{ alpha=rep(0.1,p)}
decomp<-"diagonal"
}else{
alpha<-param$alpha
}
# second sigma matrix
if(is.null(param[["sig"]])){sig<-param$sig<-startParamSigma(p, decompSigma, tree, data)}else{sig<-param$sig}
# beta matrix
if(is.null(param[["beta"]])){beta=0.01}else{beta=param$beta} # 1 beta cf. mvEB
nalpha<-length(alpha)
nsigma<-length(sigma)
#set parameters for ecological constraint model:
if(model=="EC"|| model=="RC" || model=="BMOU" || model=="BMOUi"){
before<-2
after<-1
if(model=="EC"||model=="BMOU"){
model<-"ER"
nmod<-"Ecological Constraint model"
nsig<-nalpha
}else if(model=="RC" || model=="BMOUi"){
model<-"RR"
nmod<-"Radiation and Ecological Constraint model"
nsig<-nalpha+nsigma
}
}else if(model=="ER"|| model=="RR" || model=="OUBM" || model=="OUBMi"){
before<-1
after<-2
if(model=="RR" || model=="OUBMi"){
model<-"RR"
nmod<-"Release and Radiate model"
nsig<-nalpha+nsigma
}else if(model=="ER" || model=="OUBM"){
model<-"ER"
nmod<-"Ecological Release model"
nsig<-nalpha
}
}else if(model=="EBOU" || model=="EBOUi"){# OU-ACDC models
before<-2
after<-1
if(model=="EBOU"){
model<-"OV"
nmod<-"ACDC to OU process with the same rate"
nsig<-nalpha
nbeta<-nsig+nsigma
}else if(model=="EBOUi"){
model<-"OVG"
nmod<-"ACDC to OU process with independent rate"
nsig<-nalpha+nsigma
nbeta<-nsig+nsigma
}
}else if(model=="OUEB" || model=="OUEBi"){
before<-1
after<-2
if(model=="OUEB"){
model<-"OV"
nmod<-"OU to ACDC process with the same rate"
nsig<-nalpha
nbeta<-nsig+nsigma
}else if(model=="OUEBi"){
model<-"OVG"
nmod<-"OU to ACDC process with independent rate"
nsig<-nalpha+nsigma
nbeta<-nsig+nsigma
}
}else if(model=="EBBM" || model=="EBBMi"){ # BM-ACDC models
before<-2
after<-1
if(model=="EBBM"){
model<-"CV"
nmod<-"ACDC to BM process with the same rate"
nalpha<-1
nbeta<-1
nsig<-nbeta
}else if(model=="EBBMi"){
model<-"CVG"
nmod<-"ACDC to BM process with independent rate"
nalpha<-1
nbeta<-1
nsig<-nbeta+nsigma
}
}else if(model=="BMEB" || model=="BMEBi"){
before<-1
after<-2
if(model=="BMEB"){
model<-"CV"
nmod<-"BM to ACDC process with the same rate"
nalpha<-1
nbeta<-1
nsig<-nbeta
}else if(model=="BMEBi"){
model<-"CVG"
nmod<-"BM to ACDC process with independent rate"
nalpha<-1
nbeta<-1
nsig<-nbeta+nsigma
}
}
# function for alpha parameterization
buildA<-function(x,p){matrixParam(x,p,decomp)}
# Matrix parameterization used for the A matrix
decompfun <- function(par){buildA(par,p)}
# function for sigma parameterization
sigmafun <- function(par) {symPar(par, decomp=decompSigma, p=p)}
# bind error to a vector
if(!is.null(error)){
error<-as.vector(error)
error[is.na(error)] <- 0
}
# bind data to a vector
if(is.matrix(data)){ dat<-as.vector(data) }else{ dat<-as.vector(as.matrix(data))}
# number of columns of the design matrix
sizeD<-ncol(W)
##------------------------Likelihood function---------------------------------##
# ou-bm
lik.shift_oubm<-function(alpha,sigma,sig,dat,error,vcvList,p,theta_mle=TRUE,theta=NULL){
if(p==1){
Vou<-.Call(mvmorph_covar_ou_fixed,A=vcvList[[before]],alpha=alpha$values, sigma=sigma)
}else{
Vou<-.Call(mvmorph_covar_mat,nterm=as.integer(n),bt=vcvList[[before]],lambda=alpha$values,S=alpha$vectors,sigmasq=sigma, S1=alpha$invectors)
}
V<-.Call(kronecker_shift, R=sig, C=vcvList[[after]], Rrows=as.integer(p), Crows=as.integer(n), V=Vou)
loglik<-loglik_mvmorph(dat,V,W,n,p,error,precalc,method,sizeD=sizeD,NA_val=NA_val,Indice_NA=Indice_NA,theta_mle=theta_mle,theta=theta)
list(loglik=-loglik$logl, theta=loglik$anc)
}
# eb-ou
lik.shift_ebou<-function(alpha,sigma,sig,beta,dat,error,vcvList,p,theta_mle=TRUE,theta=NULL){
if(p==1){
Vou<-.Call(mvmorph_covar_ou_fixed,A=vcvList[[before]],alpha=alpha$values, sigma=sigma)
}else{
Vou<-.Call(mvmorph_covar_mat,nterm=as.integer(n),bt=vcvList[[before]],lambda=alpha$values,S=alpha$vectors,sigmasq=sigma, S1=alpha$invectors)
}
V<-.Call(kronecker_shiftEB_OU, R=sig, C=vcvList[[after]], beta=matrix(beta,p,p), Rrows=as.integer(p), Crows=as.integer(n), V=Vou)
loglik<-loglik_mvmorph(dat,V,W,n,p,error,precalc,method,sizeD=sizeD,NA_val=NA_val,Indice_NA=Indice_NA,theta_mle=theta_mle,theta=theta)
list(loglik=-loglik$logl, theta=loglik$anc)
}
# eb-bm
lik.shift_ebbm<-function(beta,sigma,sig,dat,error,vcvList,p,theta_mle=TRUE,theta=NULL){
# C2 is the EB matrix
V<-.Call(kronecker_shiftEB_BM, R1=sig, R2=sigma, C1=vcvList[[before]], C2=vcvList[[after]], beta=matrix(beta,p,p), Rrows=as.integer(p), Crows=as.integer(n))
loglik<-loglik_mvmorph(dat,V,W,n,p,error,precalc,method,sizeD=sizeD,NA_val=NA_val,Indice_NA=Indice_NA,theta_mle=theta_mle,theta=theta)
list(loglik=-loglik$logl, theta=loglik$anc)
}
# sparse matrix ou-bm
lik.shift_oubm_sparse<-function(alpha,sigma,sig,dat,error,vcvList,p,theta_mle=TRUE,theta=NULL){
Vou<-.Call(mvmorph_covar_mat,nterm=as.integer(n),bt=vcvList[[before]],lambda=alpha$values,S=alpha$vectors,sigmasq=sigma, S1=alpha$invectors)
V<-.Call(kroneckerSpar_shift, R=sig, C=vcvList[[after]], Rrows=as.integer(p), Crows=as.integer(n), V=Vou, IA=as.integer(IAr), JA=as.integer(JAr), A=as.double(precalcMat@entries))
loglik<-loglik_mvmorph(dat,V,W,n,p,error,precalc,method="sparse",ch=ch,precalcMat=precalcMat,sizeD=sizeD,NA_val=NA_val,Indice_NA=Indice_NA,theta_mle=theta_mle,theta=theta)
list(loglik=-loglik$logl, theta=loglik$anc)
}
# sparse matrix eb-ou
#lik.shift_ebou_sparse<-function(alpha,sigma,sig,beta,dat,error,vcvList,p){
# eig<-eigen(alpha)
# Vou<-.Call("mvmorph_covar_ou_sparse", A=as.double(precalcMat@entries), JA=as.integer(JAr), IA=as.integer(IAr), as.integer(n*p), bt=vcvList[[before]], lambda=eig$values, S=eig$vectors, sigmasq=sigma)
# V<-.Call("kroneckerSpar_shift_OU_EB", R=sig, C=vcvList[[after]], beta=beta, Rrows=as.integer(p), Crows=as.integer(n), IA=as.integer(IAr), JA=as.integer(JAr), A=as.double(precalcMat@entries))
# loglik<-loglik_mvmorph(dat,V,W,n,p,error,precalc,method="sparse")
# list(loglik=-loglik$logl, theta=loglik$anc)
#}
# sparse matrix eb-ou
lik.shift_ebou_sparse<-function(alpha,sigma,sig,beta,dat,error,vcvList,p,theta_mle=TRUE,theta=NULL){
Vou<-.Call(mvmorph_covar_mat,nterm=as.integer(n),bt=vcvList[[before]],lambda=alpha$values,S=alpha$vectors,sigmasq=sigma, S1=alpha$invectors)
V<-.Call(kroneckerSpar_shift_OU_EB, R=sig, C=vcvList[[after]], beta=matrix(beta,p,p), Rrows=as.integer(p), Crows=as.integer(n), V=Vou, IA=as.integer(IAr), JA=as.integer(JAr), A=as.double(precalcMat@entries))
loglik<-loglik_mvmorph(dat,V,W,n,p,error,precalc,method="sparse", ch=ch, precalcMat=precalcMat,sizeD=sizeD,NA_val=NA_val,Indice_NA=Indice_NA,theta_mle=theta_mle,theta=theta)
list(loglik=-loglik$logl, theta=loglik$anc)
}
# sparse matrix eb-bm
lik.shift_ebbm_sparse<-function(beta,sigma,sig,dat,error,vcvList,p,theta_mle=TRUE,theta=NULL){
V<-.Call(kroneckerSpar_shift_EB_BM, R1=sig, R2=sigma, C1=vcvList[[before]], C2=vcvList[[after]], beta=matrix(beta,p,p), Rrows=as.integer(p), Crows=as.integer(n), IA=as.integer(IAr), JA=as.integer(JAr), A=as.double(precalcMat@entries))
loglik<-loglik_mvmorph(dat,V,W,n,p,error,precalc, method="sparse",ch=ch, precalcMat=precalcMat,sizeD=sizeD,NA_val=NA_val,Indice_NA=Indice_NA,theta_mle=theta_mle,theta=theta)
list(loglik=-loglik$logl, theta=loglik$anc)
}
##----------------------Optimization------------------------------------------##
# select the loglik functions & starting values
switch(model,
"RR"={
# Brownian and OU models with different rates
if(method=="sparse"){
lik.shift<-lik.shift_oubm_sparse
}else{
lik.shift<-lik.shift_oubm
}
starting<-c(alpha,sigma,sig)
},
"ER"={
# Brownian and OU models with the same rates
if(method=="sparse"){
lik.shift<-lik.shift_oubm_sparse
}else{
lik.shift<-lik.shift_oubm
}
starting<-c(alpha,sigma)
},
"CV"={
# Brownian & ACDC models with the same rates
if(method=="sparse"){
lik.shift<-lik.shift_ebbm_sparse
}else{
lik.shift<-lik.shift_ebbm
}
starting<-c(beta,sigma)
},
"CVG"={
# Brownian & ACDC models with different rates
if(method=="sparse"){
lik.shift<-lik.shift_ebbm_sparse
}else{
lik.shift<-lik.shift_ebbm
}
starting<-c(beta,sigma,sig)
},
"OV"={
# OU & ACDC models with the same rates
if(method=="sparse"){
lik.shift<-lik.shift_ebou_sparse
}else{
lik.shift<-lik.shift_ebou
}
starting<-c(alpha,sigma,beta)
},
"OVG"={
# OU & ACDC models with independent rates
if(method=="sparse"){
lik.shift<-lik.shift_ebou_sparse
}else{
lik.shift<-lik.shift_ebou
}
starting<-c(alpha,sigma,sig,beta)
})
# optimizer
if(model=="OVG" || model =="OV"){
##-------------------Function log-likelihood----------------------------------##
llfun <- function(par,...){
args <- list(...)
if(is.null(args[["root.mle"]])) args$root.mle <- TRUE
if(is.null(args[["theta"]])) args$theta <- FALSE
if(args$root.mle==TRUE){
result <- -as.numeric(lik.shift(alpha=buildA(par[seq_len(nalpha)],p),sigma=sym.par(par[nalpha+seq_len(nsigma)]), sig=sym.par(par[nsig+seq_len(nsigma)]), beta=ratevalue(up,low,par[nbeta+1]), error=error, dat=dat, vcvList=vcvList,p=p)$loglik)
if(args$theta==TRUE){
theta <- as.numeric(lik.shift(alpha=buildA(par[seq_len(nalpha)],p),sigma=sym.par(par[nalpha+seq_len(nsigma)]), sig=sym.par(par[nsig+seq_len(nsigma)]), beta=ratevalue(up,low,par[nbeta+1]), error=error, dat=dat, vcvList=vcvList,p=p)$theta)
result<-list(logl=result, theta=theta)
}
}else{
result <- -as.numeric(lik.shift(alpha=buildA(par[seq_len(nalpha)],p),sigma=sym.par(par[nalpha+seq_len(nsigma)]), sig=sym.par(par[nsig+seq_len(nsigma)]), beta=ratevalue(up,low,par[nbeta+1]), error=error, dat=dat, vcvList=vcvList,p=p,theta_mle=FALSE,theta=par[nsig+nsigma+seq_len(p)])$loglik)
}
return(result)
}
# attributes to the loglik function
attr(llfun, "model")<-model
attr(llfun, "alpha")<-nalpha
attr(llfun, "sigma")<-nsigma
attr(llfun, "sig")<-nsig
attr(llfun, "beta")<-1
attr(llfun, "theta")<-p
class(llfun) = c("mvmorph.llik")
## Check first if we want to return the log-likelihood function only
if(optimization=="fixed"){
message("No optimization performed, only the Log-likelihood function is returned with default parameters.")
param$sigmafun<-sigmafun
param$alphafun<-decompfun
param$nbspecies<-n
param$ntraits<-p
param$nregimes<-2
param$method<-method
param$optimization<-optimization
param$traits<-colnames(data)
param$names_regimes<-colnames(tree$mapped.edge)
param$model<-model
param$root<-root
param$decomp<-decomp
param$decompSigma<-decompSigma
param$before<-before
param$after<-after
theta.mat<-rep(0,p)
results<-list(llik=llfun, theta=theta.mat, alpha=buildA(alpha,p), sigma=sigmafun(sigma), sig=sigmafun(sig), beta=beta, param=param)
class(results)<-c("mvmorph")
invisible(results)
return(results)
}
if(optimization!="subplex"){
estim <- optim(par=starting,fn = function (par) {lik.shift(alpha=buildA(par[seq_len(nalpha)],p),sigma=sigmafun(par[nalpha+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), beta=ratevalue(up,low,par[nbeta+1]), error=error, dat=dat, vcvList=vcvList,p=p)$loglik },gr=NULL,hessian=TRUE,method=optimization,control=control)
}else{
estim <- subplex(par=starting,fn = function (par) {lik.shift(alpha=buildA(par[seq_len(nalpha)],p),sigma=sigmafun(par[nalpha+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), beta=par[nbeta+1], error=error, dat=dat, vcvList=vcvList,p=p)$loglik},hessian=TRUE,control=control)
}
}else if(model=="CV" || model=="CVG"){
##-------------------Function log-likelihood----------------------------------##
llfun <- function(par,...){
args <- list(...)
if(is.null(args[["root.mle"]])) args$root.mle <- TRUE
if(is.null(args[["theta"]])) args$theta <- FALSE
if(args$root.mle==TRUE){
result <- -as.numeric(lik.shift(beta=ratevalue(up,low,par[nbeta]),sigma=sigmafun(par[nbeta+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p)$loglik)
if(args$theta==TRUE){
theta <- as.numeric(lik.shift(beta=ratevalue(up,low,par[nbeta]),sigma=sigmafun(par[nbeta+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p)$theta)
result<-list(logl=result, theta=theta)
}
}else{
result <- -as.numeric(lik.shift(beta=ratevalue(up,low,par[nbeta]),sigma=sigmafun(par[nbeta+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p,theta_mle=FALSE,theta=par[nsig+nsigma+seq_len(p)])$loglik)
}
return(result)
}
# attributes to the loglik function
attr(llfun, "model")<-model
attr(llfun, "beta")<-nbeta
attr(llfun, "sigma")<-nsigma
attr(llfun, "sig")<-nsig
attr(llfun, "theta")<-p
class(llfun) = c("mvmorph.llik")
## Check first if we want to return the log-likelihood function only
if(optimization=="fixed"){
message("No optimization performed, only the Log-likelihood function is returned with default parameters.")
param$sigmafun<-sigmafun
param$alphafun<-decompfun
param$nbspecies<-n
param$ntraits<-p
param$nregimes<-2
param$method<-method
param$optimization<-optimization
param$traits<-colnames(data)
param$names_regimes<-colnames(tree$mapped.edge)
param$model<-model
param$root<-root
param$decomp<-decomp
param$decompSigma<-decompSigma
param$before<-before
param$after<-after
theta.mat<-rep(0,p)
results<-list(llik=llfun, theta=theta.mat, sigma=sigmafun(sigma), sig=sigmafun(sig), beta=beta, param=param)
class(results)<-c("mvmorph")
invisible(results)
return(results)
}
if(optimization!="subplex"){
estim <- optim(par=starting,fn = function (par) {lik.shift(beta=ratevalue(up,low,par[nbeta]),sigma=sigmafun(par[nbeta+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p)$loglik},gr=NULL,hessian=TRUE,method=optimization,control=control)
}else{
estim <- subplex(par=starting,fn = function (par) {lik.shift(beta=ratevalue(up,low,par[nbeta]),sigma=sigmafun(par[nbeta+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p)$loglik},hessian=TRUE,control=control)
}
}else{
##-------------------Function log-likelihood----------------------------------##
llfun <- function(par,...){
args <- list(...)
if(is.null(args[["root.mle"]])) args$root.mle <- TRUE
if(is.null(args[["theta"]])) args$theta <- FALSE
if(args$root.mle==TRUE){
result <- -as.numeric(lik.shift(alpha=buildA(par[seq_len(nalpha)],p),sigma=sigmafun(par[nalpha+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p)$loglik)
if(args$theta==TRUE){
theta <- as.numeric(lik.shift(alpha=buildA(par[seq_len(nalpha)],p),sigma=sigmafun(par[nalpha+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p)$theta)
result<-list(logl=result, theta=theta)
}
}else{
result <- -as.numeric(lik.shift(alpha=buildA(par[seq_len(nalpha)],p),sigma=sigmafun(par[nalpha+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p,theta_mle=FALSE,theta=par[nsig+nsigma+seq_len(p)])$loglik)
}
return(result)
}
# attributes to the loglik function
attr(llfun, "model")<-model
attr(llfun, "alpha")<-nalpha
attr(llfun, "sigma")<-nsigma
attr(llfun, "sig")<-nsig
attr(llfun, "theta")<-p
class(llfun) = c("mvmorph.llik")
## Check first if we want to return the log-likelihood function only
if(optimization=="fixed"){
message("No optimization performed, only the Log-likelihood function is returned with default parameters.")
param$sigmafun<-sigmafun
param$alphafun<-decompfun
param$nbspecies<-n
param$ntraits<-p
param$nregimes<-2
param$method<-method
param$optimization<-optimization
param$traits<-colnames(data)
param$names_regimes<-colnames(tree$mapped.edge)
param$model<-model
param$root<-root
param$decomp<-decomp
param$decompSigma<-decompSigma
param$before<-before
param$after<-after
theta.mat<-rep(0,p)
results<-list(llik=llfun,theta=theta.mat, alpha=buildA(alpha,p), sigma=sigmafun(sigma), sig=sigmafun(sig), param=param)
class(results)<-c("mvmorph")
invisible(results)
return(results)
}
if(optimization!="subplex"){
estim <- optim(par=starting,fn = function (par) {lik.shift(alpha=buildA(par[seq_len(nalpha)],p),sigma=sigmafun(par[nalpha+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p)$loglik},gr=NULL,hessian=TRUE,method=optimization,control=control)
}else{
estim <- subplex(par=starting,fn = function (par) {lik.shift(alpha=buildA(par[seq_len(nalpha)],p),sigma=sigmafun(par[nalpha+seq_len(nsigma)]), sig=sigmafun(par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p)$loglik},hessian=TRUE,control=control)
}
}
##---------------------theta estimation---------------------------------------##
if(model=="OVG" || model=="OV"){
res.theta<-lik.shift(alpha=buildA(estim$par[seq_len(nalpha)],p),sigma=sigmafun(estim$par[nalpha+seq_len(nsigma)]), sig=sigmafun(estim$par[nsig+seq_len(nsigma)]), beta=ratevalue(up,low,estim$par[nbeta+1]), error=error, dat=dat, vcvList=vcvList,p=p)$theta
}else if(model=="CV" || model=="CVG"){
res.theta<-lik.shift(beta=ratevalue(up,low,estim$par[nbeta]),sigma=sigmafun(estim$par[nbeta+seq_len(nsigma)]), sig=sigmafun(estim$par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p)$theta
}else{
res.theta<-lik.shift(alpha=buildA(estim$par[seq_len(nalpha)],p),sigma=sigmafun(estim$par[nalpha+seq_len(nsigma)]), sig=sigmafun(estim$par[nsig+seq_len(nsigma)]), error=error, dat=dat, vcvList=vcvList,p=p)$theta
}
##---------------------Diagnostics--------------------------------------------##
hess<-eigen(estim$hessian)$values
if(estim$convergence==0 & diagnostic==TRUE){
cat("successful convergence of the optimizer","\n")
}else if(estim$convergence==1 & diagnostic==TRUE){
cat("maximum limit iteration has been reached, please consider increase maxit","\n")
}else if(diagnostic==TRUE){
cat("convergence of the optimizer has not been reached, try simpler model","\n")
}
if(any(hess<0)){
hess.val<-1
if(diagnostic==TRUE){
cat("unreliable solution has been reached, check hessian eigenvectors or try simpler model","\n")}
}else{
hess.val<-0
if(diagnostic==TRUE){
cat("a reliable solution has been reached","\n")}
}
##-------------------Summarize Results----------------------------------------##
# Loglikelihood estimate
LL<- -estim$value
# free parameters
nparam=p+length(estim$par)
# maximum likelihood estimates of alpha and sigma
estim.alpha<-estim$par[seq_len(nalpha)]
estim.sigma<-estim$par[nalpha+seq_len(nsigma)]
if(model=="RR" || model=="CVG" || model=="OVG"){
estim.sig<-estim$par[nsig+seq_len(nsigma)]
}
# names for the plot
if(is.null(colnames(data))){
names_data_matrix<-rep("",p)
names_data<-list(names_data_matrix,names_data_matrix)
}else{
names_data_matrix<-colnames(data)
names_data<-list(names_data_matrix,names_data_matrix)
}
# alpha matrix
alpha.mat<-buildA(estim.alpha,p)$A
dimnames(alpha.mat)<-names_data
# sigma matrix
sigma.mat<-sigmafun(estim.sigma)
dimnames(sigma.mat)<-names_data
# sig matrix
if(model=="RR" || model=="CVG" || model=="OVG"){
sig.mat<-sigmafun(estim.sig)
dimnames(sig.mat)<-names_data
}
# beta matrix
if(model=="OVG" || model=="OV"){
estim.beta<-ratevalue(up,low,estim$par[nbeta+1])
beta.mat<-c(estim.beta)
}
if(model=="CV" || model=="CVG"){
estim.beta<-ratevalue(up,low,estim$par[nbeta])
beta.mat<-c(estim.beta)
}
# AIC
AIC<- -2*LL+2*nparam
# AIC corrected
nobs <- length(which(!is.na(data)))
AICc<-AIC+((2*nparam*(nparam+1))/(nobs-nparam-1)) #Hurvich et Tsai, 1989
# matrix of estimated theta values
theta.mat<-matrix(res.theta,1)
rownames(theta.mat)<-"theta"
colnames(theta.mat)<-names_data_matrix
##-------------------Print results--------------------------------------------##
if(echo==TRUE){
cat("\n")
cat("-- Summary results for the",nmod," --","\n")
cat("LogLikelihood:","\t",LL,"\n")
cat("AIC:","\t",AIC,"\n")
cat("AICc:","\t",AICc,"\n")
cat(nparam,"parameters")
cat("\n")
cat("Estimated theta values","\n")
cat("______________________","\n")
print(theta.mat)
cat("\n")
if(model=="CV" || model=="CVG"){
cat("ML beta values","\n")
cat("______________________","\n")
print(beta.mat)
}else{
cat("ML alpha values","\n")
cat("______________________","\n")
print(alpha.mat)
}
cat("\n")
cat("ML sigma values","\n")
cat("______________________","\n")
print(sigma.mat)
if(model=="RR" || model=="radiate"){
cat("\n")
cat("ML sigma radiation values","\n")
cat("______________________","\n")
print(sig.mat)
}
if(model=="CVG" || model=="OVG"){
cat("\n")
cat("ML sigma values ( recent slice:",colnames(tree$mapped.edge)[2],")","\n")
cat("______________________","\n")
print(sig.mat)
}
if(model=="OVG" || model=="OV"){
cat("\n")
cat("ML beta values","\n")
cat("______________________","\n")
print(beta.mat)
}
}
##-------------------Save infos in parameters---------------------------------##
param$nparam<-nparam
param$nbspecies<-n
param$ntraits<-p
param$nregimes<-2
param$method<-method
param$optimization<-optimization
param$traits<-colnames(data)
param$names_regimes<-colnames(tree$mapped.edge)
param$model<-c(model,nmod)
param$before<-before
param$after<-after
param$sigmafun<-sigmafun
param$alphafun<-decompfun
param$opt<-estim
##------------------List results----------------------------------------------##
if(model=="ER"){
results<-list(LogLik=LL, AIC=AIC, AICc=AICc, theta=theta.mat, alpha=alpha.mat, sigma=sigma.mat, convergence=estim$convergence, hess.values=hess.val, param=param, llik=llfun)
}else if(model=="RR"){
results<-list(LogLik=LL, AIC=AIC, AICc=AICc, theta=theta.mat, alpha=alpha.mat, sigma=sigma.mat, sig=sig.mat, convergence=estim$convergence, hess.values=hess.val, param=param, llik=llfun)
}else if(model=="CV"){
results<-list(LogLik=LL, AIC=AIC, AICc=AICc, theta=theta.mat, beta=beta.mat, sigma=sigma.mat, convergence=estim$convergence, hess.values=hess.val, param=param, llik=llfun)
}else if(model=="CVG"){
results<-list(LogLik=LL, AIC=AIC, AICc=AICc, theta=theta.mat, beta=beta.mat, sigma=sigma.mat, sig=sig.mat, convergence=estim$convergence, hess.values=hess.val, param=param, llik=llfun)
}else if(model=="OV"){
results<-list(LogLik=LL, AIC=AIC, AICc=AICc, theta=theta.mat, alpha=alpha.mat, sigma=sigma.mat, beta=beta.mat, convergence=estim$convergence, hess.values=hess.val, param=param, llik=llfun)
}else if(model=="OVG"){
results<-list(LogLik=LL, AIC=AIC, AICc=AICc, theta=theta.mat, alpha=alpha.mat, sigma=sigma.mat, sig=sig.mat, beta=beta.mat, convergence=estim$convergence, hess.values=hess.val, param=param, llik=llfun)
}
class(results)<-c("mvmorph","mvmorph.shift")
invisible(results)
}
} #fin de la fonction
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.