rm(list=ls()) setwd("~/Dropbox/FCU/Teaching/Mentoring/2022Spring-Current_Professor/2022Fall/YourueiMin/rcode_dcj/empricaldataanalysis/matchedtreetraitTony/BodyMass") source("~/Dropbox/FCU/Teaching/Mentoring/2022Spring-Current_Professor/2022Fall/YourueiMin/rcode_dcj/empricaldataanalysis/mainPhyarima.R") ########################### ##### diurnal_subtree ##### ########################### diurnal_subtree<-read.tree("diurnal_subtree.nwk") tiplabel<-diurnal_subtree$tip.label tiplabel length(tiplabel) diurnal_BodyMass<-read.csv("diurnal_BodyMass.csv") diurnal_BodyMass length(diurnal_BodyMass$Species) same<-intersect(tiplabel,diurnal_BodyMass$Species) diurnal_BodyMass$Species%in%same # Get the names of the tips to drop tips_to_drop <- setdiff(tiplabel, diurnal_BodyMass$Species) tips_to_drop<-c(tips_to_drop,"Pithecia_hirsuta","Pithecia_vanzolinii","Propithecus_candidus") # Drop the tips diurnal_subtree_with_bodymass <- drop.tip(diurnal_subtree, tips_to_drop) plot(diurnal_subtree_with_bodymass) diurnal_BodyMass <- diurnal_BodyMass[!(diurnal_BodyMass$Species %in% c("Pithecia_hirsuta", "Pithecia_vanzolinii", "Propithecus_candidus")), ] dim(diurnal_BodyMass) y<-diurnal_BodyMass$BodyMassMale_kg names(y)<-diurnal_BodyMass$Species ####################################### ############ diurnal Analysis ####### ####################################### mfit<-modelfit(trait=y, tree=diurnal_subtree_with_bodymass) diurnal.empmodelout<-empmodel(mfit) diurnal.empmodelout$parammleset diurnal.empmodelout$fitset #empmodelout$parammleset[1,][c("phi","sigsq")] xtable(diurnal.empmodelout$parammleset,digits=3) xtable(diurnal.empmodelout$fitset,digits=3) ## add up test.phi here diurnal_betas.rrphylo<-compute_betas_rrphylo(diurnal_subtree_with_bodymass, y)$betas.rrphylo test.phi(param=diurnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=diurnal_betas.rrphylo,tree=diurnal_subtree_with_bodymass) ########################### ##### nocturnal_subtree ### ########################### nocturnal_subtree<-read.tree("nocturnal_subtree.nwk") tiplabel<-nocturnal_subtree$tip.label tiplabel length(tiplabel) nocturnal_BodyMass<-read.csv("nocturnal_BodyMass.csv") length(nocturnal_BodyMass$Species) same<-intersect(tiplabel,nocturnal_BodyMass$Species) nocturnal_BodyMass$Species%in%same # Get the names of the tips to drop tips_to_drop <- setdiff(tiplabel, nocturnal_BodyMass$Species) tips_to_drop # Drop the tips nocturnal_subtree_with_bodymass <- drop.tip(nocturnal_subtree, tips_to_drop) plot(nocturnal_subtree_with_bodymass) nocturnal_BodyMass nocturnal_subtree_with_bodymass y<-nocturnal_BodyMass$BodyMassMale_kg names(y)<-nocturnal_BodyMass$Species tips_to_drop<-names(y)[is.na(y)] tips_to_drop # Drop the tips nocturnal_subtree_with_bodymass <- drop.tip(nocturnal_subtree_with_bodymass, tips_to_drop) nocturnal_subtree_with_bodymass length(y) y y<-y[!is.na(y)] y length(y) ####################################### ############ nocturnal Analysis ####### ####################################### #compute_betas_rrphylo(nocturnal_subtree_with_bodymass, y)$betas.rrphylo mfit<-modelfit(trait=y, tree=nocturnal_subtree_with_bodymass) nocturnal.empmodelout<-empmodel(mfit) nocturnal.empmodelout$parammleset nocturnal.empmodelout$fitset xtable(nocturnal.empmodelout$parammleset,digits=3) xtable(nocturnal.empmodelout$fitset,digits=3) ## add up test.phi here nocturnal_betas.rrphylo<-compute_betas_rrphylo(nocturnal_subtree_with_bodymass, y)$betas.rrphylo test.phi(param=nocturnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=nocturnal_betas.rrphylo,tree=nocturnal_subtree_with_bodymass) #################################### #########Combined tree############## #################################### tree_with_bodymass<-bind.tree(diurnal_subtree_with_bodymass, nocturnal_subtree_with_bodymass, where = "root", position = 0, interactive = FALSE) tree_with_bodymass<-chronopl(tree_with_bodymass, 0, age.min = 1, age.max = NULL,node = "root", S = 1, tol = 1e-8,CV = FALSE, eval.max = 500, iter.max = 500) plot(tree_with_bodymass, show.tip.label = FALSE) BodyMass<-rbind(diurnal_BodyMass,nocturnal_BodyMass) y<-BodyMass$BodyMassMale_kg names(y)<-BodyMass$Species y length(y) names(y)%in% tree_with_bodymass$tip.label tree_with_bodymass$tip.label %in%names(y) ####################################### ############ Full Analysis ########## ####################################### #compute_betas_rrphylo(tree_with_bodymass, y)$betas.rrphylo mfit<-modelfit(trait=y, tree=tree_with_bodymass) comb.empmodelout<-empmodel(mfit) comb.empmodelout$parammleset comb.empmodelout$fitset xtable(comb.empmodelout$parammleset,digits=3) xtable(comb.empmodelout$fitset,digits=3) ## add up test.phi here comb_betas.rrphylo<-compute_betas_rrphylo(tree_with_bodymass, y)$betas.rrphylo test.phi(param=comb.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=comb_betas.rrphylo,tree=tree_with_bodymass) ############## ####AR1 test## ############## #########Combined tree############## comb_betas.rrphylo<-compute_betas_rrphylo(tree_with_bodymass, y)$betas.rrphylo combAr1test <- test.phi(param=comb.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=comb_betas.rrphylo,tree=tree_with_bodymass) #########diurnal tree############## diurnal_betas.rrphylo<-compute_betas_rrphylo(diurnal_subtree_with_bodymass, y)$betas.rrphylo diurnalAr1test <- test.phi(param=diurnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=diurnal_betas.rrphylo,tree=diurnal_subtree_with_bodymass) #########nocturnal tree############## nocturnal_betas.rrphylo<-compute_betas_rrphylo(nocturnal_subtree_with_bodymass, y)$betas.rrphylo nocturnalAr1test <- test.phi(param=nocturnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=nocturnal_betas.rrphylo,tree=nocturnal_subtree_with_bodymass) treecase<-c("combined tree", "diurnal tree", "nocturnal tree") zval<-c(combAr1test$z,diurnalAr1test$z,nocturnalAr1test$z) pval<-c(combAr1test$pval,diurnalAr1test$pval,nocturnalAr1test$pval) sig<-ifelse(pval<0.05,"Y","N") ar1testtb<-data.frame(zval=zval,pval=pval,sig=sig) rownames(ar1testtb)<-treecase ar1testtb xtable(ar1testtb,digits=3) ################################# #### LRT across models ## ################################# LL_diurnal <- diurnal.empmodelout$fitset LL_nocturnal <- nocturnal.empmodelout$fitset LL_comb <- comb.empmodelout$fitset # Compute the test statistic test_statistic <- 2 * ((LL_diurnal$loglikset + LL_nocturnal$loglikset)-LL_comb$loglikset) test_statistic #names(test_statistic)<-c(paste(LL_comb$modelset,".chi2",sep="")) test_statistic # Perform the chi-square test p_value.array <- 1 - pchisq(test_statistic, df=LL_comb$paranumset) # two AR(1) - one AR(1) #names(p_value.array)<-c(paste(LL_comb$modelset,".pal",sep="")) p_value.array # zero seems make sense as the rate is from the y simulated from tworates<-ifelse(p_value.array<0.05,"Y","N") #names(tworates)<-rep("two rates?",5) tworates LL_diurnal$loglikset LL_nocturnal$loglikset LL_comb$loglikset mset<-LL_comb$modelset diurnallikeset<-LL_diurnal$loglikset nocturnallikeset<-LL_nocturnal$loglikset comblikeset<-LL_comb$loglikset df2rates1rate<-data.frame( diurnallikeset=diurnallikeset, nocturnallikeset=nocturnallikeset, comblikeset=comblikeset, test_statistic=test_statistic, p_value.array=p_value.array,tworates=tworates ) rownames(df2rates1rate)<-LL_diurnal$modelset colnames(df2rates1rate)<-c("log L_di","log L_no","log L_dino","chi2","pval") df2rates1rate xtable(df2rates1rate) # Do it separately, just call from the previous diurnal case and nocturnal case. # testrate2subtrees_output <- testrate2subtrees(tree=tree_with_bodymass, y=BodyMass) # tworatedf<- treeARanova(testrate2subtrees_output) # tworatedf # xtable(tworatedf,digits=3)