---
title: "Primate Body Mass--Female"
output: html_document
date: "2024-05-01"
---

```{r}
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$BodyMassFemale_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$BodyMassFemale_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$BodyMassFemale_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)
```
