rm(list=ls())
load("~/Dropbox/FCU/Teaching/Mentoring/2016Spring-2021Fall_AssociateProfessor/2020Fall/ChiYuWu/rcodeCYW/lizardyoboot.rda")
library(gridExtra)
library(phytools)
## Loading required package: ape
## Loading required package: maps
library(ape)
library(MASS)
library(ggplot2)
library(geepack)
library(geeM)
## Loading required package: Matrix
library(gee)
library(ggtree)
## ggtree v3.8.0 For help: https://yulab-smu.top/treedata-book/
##
## If you use the ggtree package suite in published research, please cite
## the appropriate paper(s):
##
## Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam.
## ggtree: an R package for visualization and annotation of phylogenetic
## trees with their covariates and other associated data. Methods in
## Ecology and Evolution. 2017, 8(1):28-36. doi:10.1111/2041-210X.12628
##
## S Xu, Z Dai, P Guo, X Fu, S Liu, L Zhou, W Tang, T Feng, M Chen, L
## Zhan, T Wu, E Hu, Y Jiang, X Bo, G Yu. ggtreeExtra: Compact
## visualization of richly annotated phylogenetic data. Molecular Biology
## and Evolution. 2021, 38(9):4039-4042. doi: 10.1093/molbev/msab166
##
## Guangchuang Yu. Using ggtree to visualize data on tree-like structures.
## Current Protocols in Bioinformatics. 2020, 69:e96. doi:10.1002/cpbi.96
##
## Attaching package: 'ggtree'
## The following object is masked from 'package:Matrix':
##
## expand
## The following object is masked from 'package:ape':
##
## rotate
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)
## % latex table generated in R 4.3.1 by xtable 1.8-4 package
## % Mon Aug 14 16:04:08 2023
## \begin{table}[ht]
## \centering
## \begin{tabular}{rllll}
## \hline
## & 1 & 2 & 3 & 4 \\
## \hline
## 1 & 3.397(0.237) & 3.411(0.344) & 3.88(0.172) & 3.753(0.176) \\
## 2 & -1.188(0.735) & -1.258(1.065) & -3.302(0.566) & -2.831(0.569) \\
## \hline
## \end{tabular}
## \end{table}
boot.nb2.mean
## b0 b1
## 3.411 -1.258
boot.nb2.sd
## b0 b1
## 0.344 1.065
boot.phypoi.mean
## b0 b1
## 3.880 -3.302
boot.phypoi.sd
## b0 b1
## 0.172 0.566
boot.phynb2.mean
## b0 b1
## 3.753 -2.831
boot.phynb2.sd
## b0 b1
## 0.176 0.569
summary(fitPOIS)
##
## Call:
## glm(formula = eggs_per_year ~ eggs_mass, family = poisson(link = "log"),
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.3968 0.2334 14.553 <2e-16 ***
## eggs_mass -1.1792 0.7276 -1.621 0.105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 42.962 on 16 degrees of freedom
## Residual deviance: 40.280 on 15 degrees of freedom
## AIC: 126.01
##
## Number of Fisher Scoring iterations: 4
summary(fitNB2)
##
## Call:
## glm.nb(formula = eggs_per_year ~ eggs_mass, data = data, init.theta = 14.97916062,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.4256 0.3583 9.560 <2e-16 ***
## eggs_mass -1.2707 1.1045 -1.151 0.25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(14.9792) family taken to be 1)
##
## Null deviance: 18.319 on 16 degrees of freedom
## Residual deviance: 17.088 on 15 degrees of freedom
## AIC: 119.28
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 14.98
## Std. Err.: 8.89
##
## 2 x log-likelihood: -113.278
phyfitPOIS
## Call: compar.gee(formula = eggs_per_year ~ eggs_mass, family = "poisson",
## phy = tree, scale.fix = TRUE)
## Number of observations: 17
## Model:
## Link: log
## Variance to Mean Relation: poisson
##
## QIC: -1393.149
##
## Summary of Residuals:
## Min 1Q Median 3Q Max
## -12.829250 3.055108 4.414483 6.225469 17.902768
##
##
## Coefficients:
## Estimate S.E. t Pr(T > |t|)
## (Intercept) 3.877618 0.1661241 23.341700 4.529787e-06
## eggs_mass -3.259023 0.5492056 -5.934068 2.335199e-03
##
## Estimated Scale Parameter: 1
## "Phylogenetic" df (dfP): 6.733333
phyfitNB2
## geem(formula = eggs_per_year ~ eggs_mass, id = c(1, 1, 1, 1,
## 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), data = <environment>,
## family = list(family = "Negative Binomial(14.9792)", link = "log",
## linkfun = function (mu)
## log(mu), linkinv = function (eta)
## pmax(exp(eta), .Machine$double.eps), variance = function (mu)
## mu + mu^2/.Theta, dev.resids = function (y, mu, wt)
## 2 * wt * (y * log(pmax(1, y)/mu) - (y + .Theta) * log((y +
## .Theta)/(mu + .Theta))), aic = function (y, n, mu,
## wt, dev)
## {
## term <- (y + .Theta) * log(mu + .Theta) - y * log(mu) +
## lgamma(y + 1) - .Theta * log(.Theta) + lgamma(.Theta) -
## lgamma(.Theta + y)
## 2 * sum(term * wt)
## }, mu.eta = function (eta)
## pmax(exp(eta), .Machine$double.eps), initialize = expression(
## {
## if (any(y < 0))
## stop("negative values not allowed for the negative binomial family")
## n <- rep(1, nobs)
## mustart <- y + (y == 0)/6
## }), validmu = function (mu)
## all(mu > 0), valideta = function (eta)
## TRUE, simulate = function (object, nsim)
## {
## ftd <- fitted(object)
## rnegbin(nsim * length(ftd), ftd, .Theta)
## }), corstr = "fixed", corr.mat = c(1, 0.933333333, 0.866666667,
## 0.666666667, 0.666666667, 0.666666667, 0.6, 0.4, 0.4, 0.4,
## 0.4, 0.266666667, 0.266666667, 0.133333333, 0.133333333,
## 0.066666667, 0, 0.933333333, 1, 0.866666667, 0.666666667,
## 0.666666667, 0.666666667, 0.6, 0.4, 0.4, 0.4, 0.4, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.866666667,
## 0.866666667, 1, 0.666666667, 0.666666667, 0.666666667, 0.6,
## 0.4, 0.4, 0.4, 0.4, 0.266666667, 0.266666667, 0.133333333,
## 0.133333333, 0.066666667, 0, 0.666666667, 0.666666667, 0.666666667,
## 1, 0.933333333, 0.866666667, 0.6, 0.4, 0.4, 0.4, 0.4, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.666666667,
## 0.666666667, 0.666666667, 0.933333333, 1, 0.866666667, 0.6,
## 0.4, 0.4, 0.4, 0.4, 0.266666667, 0.266666667, 0.133333333,
## 0.133333333, 0.066666667, 0, 0.666666667, 0.666666667, 0.666666667,
## 0.866666667, 0.866666667, 1, 0.6, 0.4, 0.4, 0.4, 0.4, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.6,
## 0.6, 0.6, 0.6, 0.6, 0.6, 1, 0.4, 0.4, 0.4, 0.4, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.4,
## 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 1, 0.933333333, 0.8, 0.8, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.4,
## 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.933333333, 1, 0.8, 0.8, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.4,
## 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.8, 0.8, 1, 0.933333333, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.4,
## 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.8, 0.8, 0.933333333, 1, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.266666667,
## 0.266666667, 0.266666667, 0.266666667, 0.266666667, 0.266666667,
## 0.266666667, 0.266666667, 0.266666667, 0.266666667, 0.266666667,
## 1, 0.933333333, 0.133333333, 0.133333333, 0.066666667, 0,
## 0.266666667, 0.266666667, 0.266666667, 0.266666667, 0.266666667,
## 0.266666667, 0.266666667, 0.266666667, 0.266666667, 0.266666667,
## 0.266666667, 0.933333333, 1, 0.133333333, 0.133333333, 0.066666667,
## 0, 0.133333333, 0.133333333, 0.133333333, 0.133333333, 0.133333333,
## 0.133333333, 0.133333333, 0.133333333, 0.133333333, 0.133333333,
## 0.133333333, 0.133333333, 0.133333333, 1, 0.933333333, 0.066666667,
## 0, 0.133333333, 0.133333333, 0.133333333, 0.133333333, 0.133333333,
## 0.133333333, 0.133333333, 0.133333333, 0.133333333, 0.133333333,
## 0.133333333, 0.133333333, 0.133333333, 0.933333333, 1, 0.066666667,
## 0, 0.066666667, 0.066666667, 0.066666667, 0.066666667, 0.066666667,
## 0.066666667, 0.066666667, 0.066666667, 0.066666667, 0.066666667,
## 0.066666667, 0.066666667, 0.066666667, 0.066666667, 0.066666667,
## 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1),
## maxit = 25, tol = 1e-04)
##
## Coefficients:
## (Intercept) eggs_mass
## 3.747 -2.784
##
## Scale Parameter: 1.566
##
## Correlation Model: fixed
## Working Correlation[1:4,1:4]:
## GA OH AL NJ
## GA 1.0000 0.9333 0.8667 0.6667
## OH 0.9333 1.0000 0.8667 0.6667
## AL 0.8667 0.8667 1.0000 0.6667
## NJ 0.6667 0.6667 0.6667 1.0000
##
## Number of clusters: 1 Maximum cluster size: 17
## Number of observations with nonzero weight: 17
# 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)
## [1] 2.776969 3.530315
quantile(bt.X,probs = c(0.025,0.975))
## 2.5% 97.5%
## 2.776617 3.526464
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
## geem(formula = eggs_per_year ~ eggs_mass, id = c(1, 1, 1, 1,
## 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), data = <environment>,
## family = list(family = "Negative Binomial(14.9792)", link = "log",
## linkfun = function (mu)
## log(mu), linkinv = function (eta)
## pmax(exp(eta), .Machine$double.eps), variance = function (mu)
## mu + mu^2/.Theta, dev.resids = function (y, mu, wt)
## 2 * wt * (y * log(pmax(1, y)/mu) - (y + .Theta) * log((y +
## .Theta)/(mu + .Theta))), aic = function (y, n, mu,
## wt, dev)
## {
## term <- (y + .Theta) * log(mu + .Theta) - y * log(mu) +
## lgamma(y + 1) - .Theta * log(.Theta) + lgamma(.Theta) -
## lgamma(.Theta + y)
## 2 * sum(term * wt)
## }, mu.eta = function (eta)
## pmax(exp(eta), .Machine$double.eps), initialize = expression(
## {
## if (any(y < 0))
## stop("negative values not allowed for the negative binomial family")
## n <- rep(1, nobs)
## mustart <- y + (y == 0)/6
## }), validmu = function (mu)
## all(mu > 0), valideta = function (eta)
## TRUE, simulate = function (object, nsim)
## {
## ftd <- fitted(object)
## rnegbin(nsim * length(ftd), ftd, .Theta)
## }), corstr = "fixed", corr.mat = c(1, 0.933333333, 0.866666667,
## 0.666666667, 0.666666667, 0.666666667, 0.6, 0.4, 0.4, 0.4,
## 0.4, 0.266666667, 0.266666667, 0.133333333, 0.133333333,
## 0.066666667, 0, 0.933333333, 1, 0.866666667, 0.666666667,
## 0.666666667, 0.666666667, 0.6, 0.4, 0.4, 0.4, 0.4, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.866666667,
## 0.866666667, 1, 0.666666667, 0.666666667, 0.666666667, 0.6,
## 0.4, 0.4, 0.4, 0.4, 0.266666667, 0.266666667, 0.133333333,
## 0.133333333, 0.066666667, 0, 0.666666667, 0.666666667, 0.666666667,
## 1, 0.933333333, 0.866666667, 0.6, 0.4, 0.4, 0.4, 0.4, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.666666667,
## 0.666666667, 0.666666667, 0.933333333, 1, 0.866666667, 0.6,
## 0.4, 0.4, 0.4, 0.4, 0.266666667, 0.266666667, 0.133333333,
## 0.133333333, 0.066666667, 0, 0.666666667, 0.666666667, 0.666666667,
## 0.866666667, 0.866666667, 1, 0.6, 0.4, 0.4, 0.4, 0.4, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.6,
## 0.6, 0.6, 0.6, 0.6, 0.6, 1, 0.4, 0.4, 0.4, 0.4, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.4,
## 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 1, 0.933333333, 0.8, 0.8, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.4,
## 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.933333333, 1, 0.8, 0.8, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.4,
## 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.8, 0.8, 1, 0.933333333, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.4,
## 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.8, 0.8, 0.933333333, 1, 0.266666667,
## 0.266666667, 0.133333333, 0.133333333, 0.066666667, 0, 0.266666667,
## 0.266666667, 0.266666667, 0.266666667, 0.266666667, 0.266666667,
## 0.266666667, 0.266666667, 0.266666667, 0.266666667, 0.266666667,
## 1, 0.933333333, 0.133333333, 0.133333333, 0.066666667, 0,
## 0.266666667, 0.266666667, 0.266666667, 0.266666667, 0.266666667,
## 0.266666667, 0.266666667, 0.266666667, 0.266666667, 0.266666667,
## 0.266666667, 0.933333333, 1, 0.133333333, 0.133333333, 0.066666667,
## 0, 0.133333333, 0.133333333, 0.133333333, 0.133333333, 0.133333333,
## 0.133333333, 0.133333333, 0.133333333, 0.133333333, 0.133333333,
## 0.133333333, 0.133333333, 0.133333333, 1, 0.933333333, 0.066666667,
## 0, 0.133333333, 0.133333333, 0.133333333, 0.133333333, 0.133333333,
## 0.133333333, 0.133333333, 0.133333333, 0.133333333, 0.133333333,
## 0.133333333, 0.133333333, 0.133333333, 0.933333333, 1, 0.066666667,
## 0, 0.066666667, 0.066666667, 0.066666667, 0.066666667, 0.066666667,
## 0.066666667, 0.066666667, 0.066666667, 0.066666667, 0.066666667,
## 0.066666667, 0.066666667, 0.066666667, 0.066666667, 0.066666667,
## 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1),
## maxit = 25, tol = 1e-04)
##
## Coefficients:
## (Intercept) eggs_mass
## 3.747 -2.784
##
## Scale Parameter: 1.566
##
## Correlation Model: fixed
## Working Correlation[1:4,1:4]:
## GA OH AL NJ
## GA 1.0000 0.9333 0.8667 0.6667
## OH 0.9333 1.0000 0.8667 0.6667
## AL 0.8667 0.8667 1.0000 0.6667
## NJ 0.6667 0.6667 0.6667 1.0000
##
## Number of clusters: 1 Maximum cluster size: 17
## Number of observations with nonzero weight: 17
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)
