---
title: "Lizard regression curve"
output: html_document
date: "2023-08-14"
---

```{r}
rm(list=ls())


load("~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/ChiYuWu/rcodeCYW/lizardyoboot.rda")

library(gridExtra)
library(phytools)
library(ape)
library(MASS)
library(ggplot2)
library(geepack)
library(geeM)
library(gee)
library(ggtree)


boot.poi.mean<-round(boot.poi.mean,3)
boot.poi.sd<-round(boot.poi.sd,3)

boot.nb2.mean<-round(boot.nb2.mean,3)
boot.nb2.sd<-round(boot.nb2.sd,3)

boot.phypoi.mean<-round(boot.phypoi.mean,3)
boot.phypoi.sd<-round(boot.phypoi.sd,3)

boot.phynb2.mean<-round(boot.phynb2.mean,3)
boot.phynb2.sd<-round(boot.phynb2.sd,3)


out.meansd.array<-array(NA,c(2,4))

#poi
out.meansd.array[1,1]<-paste(boot.poi.mean[1],"(",boot.poi.sd[1],")",sep="")
out.meansd.array[2,1]<-paste(boot.poi.mean[2],"(",boot.poi.sd[2],")",sep="")

#nb2
out.meansd.array[1,2]<-paste(boot.nb2.mean[1],"(",boot.nb2.sd[1],")",sep="")
out.meansd.array[2,2]<-paste(boot.nb2.mean[2],"(",boot.nb2.sd[2],")",sep="")

#phypoi
out.meansd.array[1,3]<-paste(boot.phypoi.mean[1],"(",boot.phypoi.sd[1],")",sep="")
out.meansd.array[2,3]<-paste(boot.phypoi.mean[2],"(",boot.phypoi.sd[2],")",sep="")

#phynb2
out.meansd.array[1,4]<-paste(boot.phynb2.mean[1],"(",boot.phynb2.sd[1],")",sep="")
out.meansd.array[2,4]<-paste(boot.phynb2.mean[2],"(",boot.phynb2.sd[2],")",sep="")

#install.packages("xtable")
library(xtable)
xtable(out.meansd.array)



boot.nb2.mean
boot.nb2.sd

boot.phypoi.mean
boot.phypoi.sd

boot.phynb2.mean
boot.phynb2.sd




summary(fitPOIS)
summary(fitNB2)
phyfitPOIS
phyfitNB2

# Now resampling


# ordsamplep.nb2(n=sims,r=r,mu=mu,Sigma=C)


# bootstrp normal

X<-rnorm(100,mean=3,sd=2)

m.X<-mean(X)
sd.X<-sd(X)
bt.X<-apply(replicate(1000,rnorm(100,mean=m.X,sd=sd.X) ),2,mean)

c(mean(X) - 1.96*sd(X)/10,mean(X) + 1.96*sd(X)/10)
quantile(bt.X,probs = c(0.025,0.975))



library(ggplot2)

# Poisson regression plot

newdat <- data.frame(eggs_mass = sort(eggs_mass))
predPOISSE<-predict(fitPOIS, newdata = newdat, type = "response", se.fit = TRUE)
newdat$predPOIS <-predPOISSE$fit 
newdat$predPOISse <-predPOISSE$se 
newdat$predPOISucl <-newdat$predPOIS + 1.96*newdat$predPOISse
newdat$predPOISlcl <-newdat$predPOIS - 1.96*newdat$predPOISse




p.poi<-ggplot(data=newdat, mapping=aes(x=eggs_mass, y=eggs_per_year)) + 
  geom_point() +  
  geom_line(mapping=aes(y=predPOIS)) +
  geom_ribbon(data=newdat, aes(ymin = predPOISlcl, ymax = predPOISucl), alpha=0.2)+ggtitle("Poisson  regression")


# NB2 regression plot

newdat$predNB2 <- predict(fitNB2, newdata = newdat, type = "response")
predNB2SE<-predict(fitNB2, newdata = newdat, type = "response", se.fit = TRUE)
newdat$predNB2 <-predNB2SE$fit 
newdat$predNB2se <-predNB2SE$se 
newdat$predNB2ucl <-newdat$predNB2 + 1.96*newdat$predNB2se
newdat$predNB2lcl <-newdat$predNB2 - 1.96*newdat$predNB2se


p.nb2<-ggplot(data=newdat, mapping=aes(x=eggs_mass, y=eggs_per_year)) + 
  geom_point() +  
  geom_line(mapping=aes(y=predNB2)) +
  geom_ribbon(data=newdat, aes(ymin = predNB2lcl, ymax = predNB2ucl), alpha=0.2)+ggtitle("Negative Binomial regression")


# phy poisson regression plot

subfitPOIS<-fitPOIS
subfitPOIS$coefficients<-phyfitPOIS$coefficients
subfitPOIS$fitted.values<-phyfitPOIS$fitted.values
phyfitPOIS<-subfitPOIS
predphyPOISSE<-predict(phyfitPOIS, newdata = newdat, type = "response", se.fit = TRUE)
newdat$predphyPOIS <-predphyPOISSE$fit 
newdat$predphyPOISse <-predphyPOISSE$se 
newdat$predphyPOISucl <-newdat$predphyPOIS + 1.96*newdat$predphyPOISse
newdat$predphyPOISlcl <-newdat$predphyPOIS - 1.96*newdat$predphyPOISse

p.phypoi<-ggplot(data=newdat, mapping=aes(x=eggs_mass, y=eggs_per_year)) + 
  geom_point() +  
  geom_line(mapping=aes(y=predphyPOIS)) +
  geom_ribbon(data=newdat, aes(ymin = predphyPOISlcl, ymax = predphyPOISucl), alpha=0.2)+ggtitle("PhyPoisson regression")


# phy NB2 regression plot

phyfitNB2<-compar.gee.nb2(eggs_per_year ~ eggs_mass,phy=tree,theta=theta)

phyfitNB2

subfitNB2<-fitNB2
subfitNB2$coefficients<-phyfitNB2$coefficients
subfitNB2$fitted.values<-c(sort(fitted(phyfitNB2), decreasing = T))#phyfitNB2$fitted.values
phyfitNB2<-subfitNB2
predphyNB2SE<-predict(phyfitNB2, newdata = newdat, type = "response", se.fit = TRUE)
newdat$predphyNB2 <-predphyNB2SE$fit 
newdat$predphyNB2se <-predphyNB2SE$se 
newdat$predphyNB2ucl <-newdat$predphyNB2 + 1.96*newdat$predphyNB2se
newdat$predphyNB2lcl <-newdat$predphyNB2 - 1.96*newdat$predphyNB2se

p.phynb2<-ggplot(data=newdat, mapping=aes(x=eggs_mass, y=eggs_per_year)) + 
  geom_point() +  
  geom_line(mapping=aes(y=predphyNB2)) +
  geom_ribbon(data=newdat, aes(ymin = predphyNB2lcl, ymax = predphyNB2ucl), alpha=0.2)+ggtitle("PhyNB2 regression")


# Poisson Regression
p.poi <- ggplot(data=newdat, mapping=aes(x=eggs_mass, y=eggs_per_year)) + 
  geom_point(color="blue") +  
  geom_line(mapping=aes(y=predPOIS), color="blue") +
  geom_ribbon(data=newdat, aes(ymin=predPOISlcl, ymax=predPOISucl), alpha=0.2, fill="blue") +
  ggtitle("Poisson Regression")+xlab("Egg Mass") + ylab("Egg Per year")

# Negative Binomial Regression
p.nb2 <- ggplot(data=newdat, mapping=aes(x=eggs_mass, y=eggs_per_year)) + 
  geom_point(color="blue") +  
  geom_line(mapping=aes(y=predNB2), color="blue") +
  geom_ribbon(data=newdat, aes(ymin=predNB2lcl, ymax=predNB2ucl), alpha=0.2, fill="blue") +
  ggtitle("Negative Binomial Regression")+xlab("Egg Mass") + ylab("Egg Per year")

# PhyPoisson Regression
p.phypoi <- ggplot(data=newdat, mapping=aes(x=eggs_mass, y=eggs_per_year)) + 
  geom_point(color="purple") +  
  geom_line(mapping=aes(y=predphyPOIS), color="purple") +
  geom_ribbon(data=newdat, aes(ymin=predphyPOISlcl, ymax=predphyPOISucl), alpha=0.2, fill="purple") +
  ggtitle("Phylogenetic Poisson Regression")+xlab("Egg Mass") + ylab("Egg Per year")

# PhyNB2 Regression
p.phynb2 <- ggplot(data=newdat, mapping=aes(x=eggs_mass, y=eggs_per_year)) + 
  geom_point(color="purple") +  
  geom_line(mapping=aes(y=predphyNB2), color="purple") +
  geom_ribbon(data=newdat, aes(ymin=predphyNB2lcl, ymax=predphyNB2ucl), alpha=0.2, fill="purple") +
  ggtitle("Phylogenetic Negative Binomial Regression")+xlab("Egg Mass") + ylab("Egg Per year")

grid.arrange(p.poi, p.nb2, p.phypoi, p.phynb2, nrow=2)

```
