Nothing
fit_t_comp<-function(phylo,data,error=NULL, model=c("MC","DDexp","DDlin"),pars=NULL,geography.object=NULL, regime.map=NULL){
#check to make sure data are univariate, with names matching phylo object
if(length(data)!=length(phylo$tip.label)){stop("length of data does not match length of tree")}
if(is.null(names(data))){stop("data missing taxa names")}
if(!is.null(dim(data))){stop("data needs to be a single trait")}
is_tip <- phylo$edge[,2] <= length(phylo$tip.label)
if(sum(diff(phylo$edge[is_tip, 2])<0)>0){ stop('fit_t_comp cannot be used with ladderized phylogenies')}
if(length(unique(ape::branching.times(phylo)))<length(ape::branching.times(phylo))){stop("fit_t_comp requires phylogenies where no nodes occur at precisely the same time [see ape::branching.times(phylo)]")}
if(is.null(error)){
if(is.null(geography.object) & is.null(regime.map)){ #single-slope version for sympatric clades
if(is.null(pars)){pars<-c(log(sqrt(var(data)/max(nodeHeights(phylo)))),0)}
params0<-c(0,pars)
if(model=="MC"){
mc.ob<-.createModel_MC(phylo)
opt<-fitTipData(mc.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
S<-opt$inferredParams[3]
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*3 - 2*(-opt$value)), aicc = (2*3 - 2*(-opt$value))+((2*3*(3+1))/(length(phylo$tip.label)-3-1)), free.parameters = 3, sig2 = as.numeric(sig2), S = as.numeric(S), z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
if(model=="DDexp"){
ddexp.ob<-.createModel_DDexp(phylo)
opt<-fitTipData(ddexp.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
r<-opt$inferredParams[3]
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*3 - 2*(-opt$value)), aicc = (2*3 - 2*(-opt$value))+((2*3*(3+1))/(length(phylo$tip.label)-3-1)), free.parameters = 3, sig2 = as.numeric(sig2), r = as.numeric(r), z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
if(model=="DDlin"){
ddlin.ob<-.createModel_DDlin(phylo)
opt<-fitTipData(ddlin.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
b<-opt$inferredParams[3]
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*3 - 2*(-opt$value)), aicc = (2*3 - 2*(-opt$value))+((2*3*(3+1))/(length(phylo$tip.label)-3-1)), free.parameters = 3, sig2 = as.numeric(sig2), b = as.numeric(b), z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
}
if(!is.null(geography.object) & is.null(regime.map)){ #single-slope version with biogeography
if(is.null(pars)){pars<-c(log(sqrt(var(data)/max(nodeHeights(phylo)))),0)}
params0<-c(0,pars)
if(length(geography.object$geography.object)<phylo$Nnode){stop("geography object cannot have fewer components than internode intervals in phylo")}
sgeo<-.resortGeoObject(phylo,geography.object) #resorts geo.object to match tip label order in Marc code
if(model=="MC"){
mc.ob<-.createModel_MC_geo(phylo,sgeo)
opt<-fitTipData(mc.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
S<-opt$inferredParams[3]
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*3 - 2*(-opt$value)), aicc = (2*3 - 2*(-opt$value))+((2*3*(3+1))/(length(phylo$tip.label)-3-1)), free.parameters = 3, sig2 = as.numeric(sig2), S = as.numeric(S), z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
if(model=="DDexp"){
ddexp.ob<-.createModel_DDexp_geo(phylo,sgeo)
opt<-fitTipData(ddexp.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
r<-opt$inferredParams[3]
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*3 - 2*(-opt$value)), aicc = (2*3 - 2*(-opt$value))+((2*3*(3+1))/(length(phylo$tip.label)-3-1)), free.parameters = 3, sig2 = as.numeric(sig2), r = as.numeric(r), z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
if(model=="DDlin"){
ddlin.ob<-.createModel_DDlin_geo(phylo,sgeo)
opt<-fitTipData(ddlin.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
b<-opt$inferredParams[3]
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*3 - 2*(-opt$value)), aicc = (2*3 - 2*(-opt$value))+((2*3*(3+1))/(length(phylo$tip.label)-3-1)), free.parameters = 3, sig2 = as.numeric(sig2), b = as.numeric(b), z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
}
if(is.null(geography.object) & !is.null(regime.map)){ #two-slope version for sympatric clades
class.object<-try(CreateClassObject(regime.map))
if(inherits(class.object, "try-error")){
class.object<-CreateClassObject(regime.map,rnd=6)
}
SMatrix<-.CreateSMatrix(class.object)
smat<-.resortSMatrix(phylo, SMatrix)
if(model=="MC"){
if(is.null(pars)){pars<-c(log(sqrt(var(data)/max(nodeHeights(phylo)))),-0.1,-0.1)}
params0<-c(0,pars)
mc.ob<-.createModel_MC_twoS(phylo,S.object=smat)
opt<-fitTipData(mc.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
S1<-opt$inferredParams[3]
S2<-opt$inferredParams[4]
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*4 - 2*(-opt$value)),", aicc = ",(2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)),", free.parameters = 4, sig2 = ",as.numeric(sig2),", S1_",SMatrix$S1," = ",as.numeric(S1),", S2_",SMatrix$S2," = ",as.numeric(S2),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
if(model=="DDexp"){
if(is.null(pars)){pars<-c(log(sqrt(var(data)/max(nodeHeights(phylo)))),0,0)}
params0<-c(0,pars)
ddexp.ob<-.createModel_DDexp_multi(phylo,r.object=smat)
opt<-fitTipData(ddexp.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
r1<-opt$inferredParams[3]
r2<-opt$inferredParams[4]
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*4 - 2*(-opt$value)),", aicc = ",(2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)),", free.parameters = 4, sig2 = ",as.numeric(sig2),", r1_",SMatrix$S1," = ",as.numeric(r1),", r2_",SMatrix$S2," = ",as.numeric(r2),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
if(model=="DDlin"){
if(is.null(pars)){pars<-c(log(sqrt(var(data)/max(nodeHeights(phylo)))),0,0)}
params0<-c(0,pars)
ddlin.ob<-.createModel_DDlin_multi(phylo,r.object=smat)
opt<-fitTipData(ddlin.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
b1<-opt$inferredParams[3]
b2<-opt$inferredParams[4]
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*4 - 2*(-opt$value)),", aicc = ",(2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)),", free.parameters = 4, sig2 = ",as.numeric(sig2),", b1_",SMatrix$S1," = ",as.numeric(b1),", b2_",SMatrix$S2," = ",as.numeric(b2),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
}
if(!is.null(geography.object) & !is.null(regime.map)){ #two-slope version with biogeography
if(length(geography.object$geography.object)<phylo$Nnode){stop("geography object cannot have fewer components than internode intervals in phylo")}
sgeo0<-.resortGeoObject(phylo,geography.object) #resorts geo.object to match tip label order in code
class.object<-try(CreateClassObject(regime.map))
if(inherits(class.object, "try-error")){
class.object<-CreateClassObject(regime.map,rnd=6)
}
SMatrix<-.CreateSMatrix(class.object)
smat0<-.resortSMatrix(phylo, SMatrix)
int<-try(.ReconcileGeoObjectSMatrix(geo.object=sgeo0,S.matrix=smat0))
#some catches in case there are small rounding issues (happens when events anagenetic in biogeography or regimes happen at a very similar time)
if(inherits(int, "try-error")){
int<-try(.ReconcileGeoObjectSMatrix(geo.object=sgeo,S.matrix=smat,rnd=6))
}
if(inherits(int, "try-error")){
int<-try(.ReconcileGeoObjectSMatrix(geo.object=sgeo,S.matrix=smat,rnd=7))
}
if(inherits(int, "try-error")){
int<-try(.ReconcileGeoObjectSMatrix(geo.object=sgeo,S.matrix=smat,rnd=4))
}
sgeo<-int$geo.object
smat<-int$S.matrix
if(model=="MC"){
if(is.null(pars)){pars<-c(log(sqrt(var(data)/max(nodeHeights(phylo)))),-0.1,-0.1)}
params0<-c(0,pars)
mc.ob<-.createModel_MC_twoS_geo(phylo,geo.object=sgeo,S.object=smat)
opt<-fitTipData(mc.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
S1<-opt$inferredParams[3]
S2<-opt$inferredParams[4]
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*4 - 2*(-opt$value)),", aicc = ",(2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)),", free.parameters = 4, sig2 = ",as.numeric(sig2),", S1_",SMatrix$S1," = ",as.numeric(S1),", S2_",SMatrix$S2," = ",as.numeric(S2),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
if(model=="DDexp"){
if(is.null(pars)){pars<-c(log(sqrt(var(data)/max(nodeHeights(phylo)))),0,0)}
params0<-c(0,pars)
ddexp.ob<-.createModel_DDexp_multi_geo(phylo,geo.object=sgeo,r.object=smat)
opt<-fitTipData(ddexp.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
r1<-opt$inferredParams[3]
r2<-opt$inferredParams[4]
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*4 - 2*(-opt$value)),", aicc = ",(2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)),", free.parameters = 4, sig2 = ",as.numeric(sig2),", r1_",SMatrix$S1," = ",as.numeric(r1),", r2_",SMatrix$S2," = ",as.numeric(r2),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
if(model=="DDlin"){
if(is.null(pars)){pars<-c(log(sqrt(var(data)/max(nodeHeights(phylo)))),0,0)}
params0<-c(0,pars)
ddlin.ob<-.createModel_DDlin_multi_geo(phylo,geo.object=sgeo,r.object=smat)
opt<-fitTipData(ddlin.ob,data,error=NULL,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
b1<-opt$inferredParams[3]
b2<-opt$inferredParams[4]
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*4 - 2*(-opt$value)),", aicc = ",(2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)),", free.parameters = 4, sig2 = ",as.numeric(sig2),", b1_",SMatrix$S1," = ",as.numeric(b1),", b2_",SMatrix$S2," = ",as.numeric(b2),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
}
} else {
if(is.null(geography.object) & is.null(regime.map)){ #single-slope version for sympatric clades
if(model=="MC"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),-0.1,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
mc.ob<-.createModel_MC_ME(phylo)
opt<-fitTipData(mc.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
S<-opt$inferredParams[3]
nuisance<-exp(opt$inferredParams[4])
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*4 - 2*(-opt$value)), aicc = (2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)), free.parameters = 4, sig2 = as.numeric(sig2), S = as.numeric(S), nuisance=as.numeric(nuisance),z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
if(model=="DDexp"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),0,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
ddexp.ob<-.createModel_DDexp_ME(phylo)
opt<-fitTipData(ddexp.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
r<-opt$inferredParams[3]
nuisance<-exp(opt$inferredParams[4])
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*3 - 2*(-opt$value)), aicc = (2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)), free.parameters = 4, sig2 = as.numeric(sig2), r = as.numeric(r), nuisance=as.numeric(nuisance), z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
if(model=="DDlin"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),0,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
ddlin.ob<-.createModel_DDlin_ME(phylo)
opt<-fitTipData(ddlin.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
b<-opt$inferredParams[3]
nuisance<-exp(opt$inferredParams[4])
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*3 - 2*(-opt$value)), aicc = (2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)), free.parameters = 4, sig2 = as.numeric(sig2), b = as.numeric(b), nuisance=as.numeric(nuisance), z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
}
if(!is.null(geography.object) & is.null(regime.map)){ #single-slope version with biogeography
if(length(geography.object$geography.object)<phylo$Nnode){stop("geography object cannot have fewer components than internode intervals in phylo")}
sgeo<-.resortGeoObject(phylo,geography.object) #resorts geo.object to match tip label order in code
if(model=="MC"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),-0.1,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
mc.ob<-.createModel_MC_geo_ME(phylo,geo.object=sgeo)
opt<-fitTipData(mc.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
S<-opt$inferredParams[3]
nuisance<-exp(opt$inferredParams[4])
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*3 - 2*(-opt$value)), aicc =(2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)), free.parameters = 4, sig2 = as.numeric(sig2), S = as.numeric(S), nuisance=as.numeric(nuisance), z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
if(model=="DDexp"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),0,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
ddexp.ob<-.createModel_DDexp_geo_ME(phylo,geo.object=sgeo)
opt<-fitTipData(ddexp.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
r<-opt$inferredParams[3]
nuisance<-exp(opt$inferredParams[4])
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*3 - 2*(-opt$value)), aicc = (2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)), free.parameters = 4, sig2 = as.numeric(sig2), r = as.numeric(r), nuisance=as.numeric(nuisance), z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
if(model=="DDlin"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),0,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
ddlin.ob<-.createModel_DDlin_geo_ME(phylo,geo.object=sgeo)
opt<-fitTipData(ddlin.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
b<-opt$inferredParams[3]
nuisance<-exp(opt$inferredParams[4])
z0<-opt$inferredParams[1]
results<-list(LH = -opt$value, aic = (2*3 - 2*(-opt$value)), aicc = (2*4 - 2*(-opt$value))+((2*4*(4+1))/(length(phylo$tip.label)-4-1)), free.parameters = 4, sig2 = as.numeric(sig2), b = as.numeric(b), nuisance=as.numeric(nuisance), z0 = as.numeric(z0), convergence = opt$convergence)
return(results)
}
}
if(is.null(geography.object) & !is.null(regime.map)){ #multi-slope version for sympatric clades (i.e., no biogeography)
class.object<-try(CreateClassObject(regime.map))
if(inherits(class.object, "try-error")){
class.object<-CreateClassObject(regime.map,rnd=6)
}
SMatrix<-.CreateSMatrix(class.object)
smat<-.resortSMatrix(phylo, SMatrix)
if(model=="MC"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),-0.1,-0.1,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
mc.ob<-.createModel_MC_twoS_ME(phylo,S.object=smat)
opt<-fitTipData(mc.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
S1<-opt$inferredParams[3]
S2<-opt$inferredParams[4]
nuisance<-exp(opt$inferredParams[5])
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*5 - 2*(-opt$value)),", aicc = ",(2*5 - 2*(-opt$value))+((2*5*(5+1))/(length(phylo$tip.label)-5-1)),", free.parameters = 5, sig2 = ",as.numeric(sig2),", S1_",SMatrix$S1," = ",as.numeric(S1),", S2_",SMatrix$S2," = ",as.numeric(S2),", nuisance = ",as.numeric(nuisance),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
if(model=="DDexp"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),0,0,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
ddexp.ob<-.createModel_DDexp_multi_ME(phylo,r.object=smat)
opt<-fitTipData(ddexp.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
r1<-opt$inferredParams[3]
r2<-opt$inferredParams[4]
nuisance<-exp(opt$inferredParams[5])
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*5 - 2*(-opt$value)),", aicc = ",(2*5 - 2*(-opt$value))+((2*5*(5+1))/(length(phylo$tip.label)-5-1)),", free.parameters = 5, sig2 = ",as.numeric(sig2),", r1_",SMatrix$S1," = ",as.numeric(r1),", r2_",SMatrix$S2," = ",as.numeric(r2),", nuisance = ",as.numeric(nuisance),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
if(model=="DDlin"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),0,0,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
ddlin.ob<-.createModel_DDlin_multi_ME(phylo,r.object=smat)
opt<-fitTipData(ddlin.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
b1<-opt$inferredParams[3]
b2<-opt$inferredParams[4]
nuisance<-exp(opt$inferredParams[5])
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*5 - 2*(-opt$value)),", aicc = ",(2*5 - 2*(-opt$value))+((2*5*(5+1))/(length(phylo$tip.label)-5-1)),", free.parameters = 5, sig2 = ",as.numeric(sig2),", b1_",SMatrix$S1," = ",as.numeric(b1),", b2_",SMatrix$S2," = ",as.numeric(b2),", nuisance = ",as.numeric(nuisance),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
}
if(!is.null(geography.object) & !is.null(regime.map)){ #multi-slope version with biogeography
if(length(geography.object$geography.object)<phylo$Nnode){stop("geography object cannot have fewer components than internode intervals in phylo")}
sgeo0<-.resortGeoObject(phylo,geography.object) #resorts geo.object to match tip label order in code
class.object<-try(CreateClassObject(regime.map))
if(inherits(class.object, "try-error")){
class.object<-CreateClassObject(regime.map,rnd=6)
}
SMatrix<-.CreateSMatrix(class.object)
smat0<-.resortSMatrix(phylo, SMatrix)
int<-try(.ReconcileGeoObjectSMatrix(geo.object=sgeo0,S.matrix=smat0))
#some catches in case there are small rounding issues (happens when events anagenetic in biogeography or regimes happen at a very similar time)
if(inherits(int, "try-error")){
int<-try(.ReconcileGeoObjectSMatrix(geo.object=sgeo,S.matrix=smat,rnd=6))
}
if(inherits(int, "try-error")){
int<-try(.ReconcileGeoObjectSMatrix(geo.object=sgeo,S.matrix=smat,rnd=7))
}
if(inherits(int, "try-error")){
int<-try(.ReconcileGeoObjectSMatrix(geo.object=sgeo,S.matrix=smat,rnd=4))
}
sgeo<-int$geo.object
smat<-int$S.matrix
if(model=="MC"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),-0.1,-0.1,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
mc.ob<-.createModel_MC_twoS_geo_ME(phylo,geo.object=sgeo,S.object=smat)
opt<-fitTipData(mc.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
S1<-opt$inferredParams[3]
S2<-opt$inferredParams[4]
nuisance<-exp(opt$inferredParams[5])
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*5 - 2*(-opt$value)),", aicc = ",(2*5 - 2*(-opt$value))+((2*5*(5+1))/(length(phylo$tip.label)-5-1)),", free.parameters = 5, sig2 = ",as.numeric(sig2),", S1_",SMatrix$S1," = ",as.numeric(S1),", S2_",SMatrix$S2," = ",as.numeric(S2),", nuisance = ",as.numeric(nuisance),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
if(model=="DDexp"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),0,0,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
ddexp.ob<-.createModel_DDexp_multi_geo_ME(phylo,geo.object=sgeo,r.object=smat)
opt<-fitTipData(ddexp.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
r1<-opt$inferredParams[3]
r2<-opt$inferredParams[4]
nuisance<-exp(opt$inferredParams[5])
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*5 - 2*(-opt$value)),", aicc = ",(2*5 - 2*(-opt$value))+((2*5*(5+1))/(length(phylo$tip.label)-5-1)),", free.parameters = 5, sig2 = ",as.numeric(sig2),", r1_",SMatrix$S1," = ",as.numeric(r1),", r2_",SMatrix$S2," = ",as.numeric(r2),", nuisance = ",as.numeric(nuisance),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
if(model=="DDlin"){
if(is.null(pars)){pars<-c(log(0.95*sqrt(var(data)/max(nodeHeights(phylo)))),0,0,log(0.05*sqrt(var(data)/max(nodeHeights(phylo)))))}
params0<-c(0,pars)
ddlin.ob<-.createModel_DDlin_multi_geo_ME(phylo,geo.object=sgeo,r.object=smat)
opt<-fitTipData(ddlin.ob,data,error,params0=params0,GLSstyle=TRUE)
sig2<-(exp(opt$inferredParams[2]))^2
b1<-opt$inferredParams[3]
b2<-opt$inferredParams[4]
nuisance<-exp(opt$inferredParams[5])
z0<-opt$inferredParams[1]
eval(parse(text=paste0("results<-list(LH = ",-opt$value,", aic = ",(2*5 - 2*(-opt$value)),", aicc = ",(2*5 - 2*(-opt$value))+((2*5*(5+1))/(length(phylo$tip.label)-5-1)),", free.parameters = 5, sig2 = ",as.numeric(sig2),", b1_",SMatrix$S1," = ",as.numeric(b1),", b2_",SMatrix$S2," = ",as.numeric(b2),", nuisance = ",as.numeric(nuisance),", z0 = ",as.numeric(z0),", convergence = ",opt$convergence,")")))
return(results)
}
}
}
}
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.