Calculate GAM models
cooc.gam.sqrt.te<-gam(cooc.z~te(phy.dist.sqrt,gen.spec), dat=data.cooc)
summary(cooc.gam.sqrt.te)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## cooc.z ~ te(phy.dist.sqrt, gen.spec)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1695 0.1006 -1.685 0.0926 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## te(phy.dist.sqrt,gen.spec) 4.964 5.656 5.192 5.72e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0441 Deviance explained = 5.26%
## GCV = 5.7365 Scale est. = 5.6755 n = 561
#only those species pairs <30my
cooc.gam.30<-gam(cooc.z.30~te(phy.dist.30.sqrt,gen.spec.30), dat=data.cooc.30)
summary(cooc.gam.30)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## cooc.z.30 ~ te(phy.dist.30.sqrt, gen.spec.30)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6138 0.1840 -3.336 0.00114 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## te(phy.dist.30.sqrt,gen.spec.30) 3.045 3.088 3.657 0.0138 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0642 Deviance explained = 8.78%
## GCV = 4.2729 Scale est. = 4.1312 n = 122
#Plot two GAMs side-by-side
#dev.new(width=5.95, height=4)
layout(matrix(1:2,ncol=2), width = c(1,1),height = c(1,1))
par(mar=c(4,2.5,2.75,0.75), mgp=c(1.5,0.5,0), las=0)
#par(mgp=c(axis.title.position, axis.label.position, axis.line.position))
my.vis.gam(cooc.gam.sqrt.te, type="response", plot.type = "contour", color="gray", ylab="Difference in niche width", main="",
xlab="", cex.lab=1.1, cex.axis=0.9, cex.main=0.1, lwd=2, cex=2)
box(bty="o", lwd=2)
my.vis.gam(cooc.gam.30,type="response", plot.type = "contour", color="gray", ylab="Difference in niche width", main="",
xlab="", cex.lab=1.1, cex.axis=0.9, cex.main=0.1, lwd=2, cex=2)
box(bty="o", lwd=2)
Set up mantels
Partial mantel
sedges.02.partial.mantel<-mantel.partial(phy.dist.cyp, final.h.inv, z.mat, method = "pearson", permutations = 999)
sedges.02.partial.mantel
##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = phy.dist.cyp, ydis = final.h.inv, zdis = z.mat, method = "pearson", permutations = 999)
##
## Mantel statistic r: -0.09168
## Significance: 0.866
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.104 0.149 0.187 0.230
## Permutation: free
## Number of permutations: 999
#Correlation between phy.dist and cooccurrence, conditional on generalism and specialism
sedges.02.partial.mantel.phy.dist<-mantel.partial(phy.dist.cyp,z.mat, final.h.inv,method = "pearson", permutations = 999)
sedges.02.partial.mantel.phy.dist
##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = phy.dist.cyp, ydis = z.mat, zdis = final.h.inv, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.09147
## Significance: 0.011
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0582 0.0721 0.0815 0.0914
## Permutation: free
## Number of permutations: 999
#Correlation between gen.spec and cooccurrence, conditional on phy.dist
sedges.02.partial.mantel.gen.spec<-mantel.partial(final.h.inv,z.mat,phy.dist.cyp ,method = "pearson", permutations = 999)
sedges.02.partial.mantel.gen.spec
##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = final.h.inv, ydis = z.mat, zdis = phy.dist.cyp, method = "pearson", permutations = 999)
##
## Mantel statistic r: -0.07963
## Significance: 0.957
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0550 0.0686 0.0798 0.0944
## Permutation: free
## Number of permutations: 999
Mantel
sedges.02.phy.dist.mantel<-mantel(phy.dist.cyp, z.mat, method = "pearson", permutations = 999)
sedges.02.phy.dist.mantel
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = phy.dist.cyp, ydis = z.mat, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.09951
## Significance: 0.006
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0582 0.0716 0.0856 0.0942
## Permutation: free
## Number of permutations: 999
sedges.02.gen.spec.mantel<-mantel(final.h.inv, z.mat,method = "pearson", permutations = 999)
sedges.02.gen.spec.mantel
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = final.h.inv, ydis = z.mat, method = "pearson", permutations = 999)
##
## Mantel statistic r: -0.08876
## Significance: 0.958
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0586 0.0734 0.0876 0.1016
## Permutation: free
## Number of permutations: 999
#I don't think that this is possible for <30 my because not a symmetrical matrix
Generalism and specialism - mean
GAM models for mean of generalism and specialism
# final GAM models
mean.cooc.gam.sqrt.te<-gam(cooc.z~te(mean.phy.dist.sqrt,gen.spec), dat=mean.data.cooc)
summary(mean.cooc.gam.sqrt.te)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## cooc.z ~ te(mean.phy.dist.sqrt, gen.spec)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1695 0.1007 -1.682 0.0931 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## te(mean.phy.dist.sqrt,gen.spec) 7.595 9.646 2.668 0.00397 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0411 Deviance explained = 5.41%
## GCV = 5.7821 Scale est. = 5.6935 n = 561
#only those species pairs <30my
mean.cooc.gam.30<-gam(cooc.z~te(mean.phy.dist.30.sqrt,gen.spec), dat=mean.data.cooc.30)
summary(mean.cooc.gam.30)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## cooc.z ~ te(mean.phy.dist.30.sqrt, gen.spec)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6138 0.1711 -3.587 0.00051 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## te(mean.phy.dist.30.sqrt,gen.spec) 16.05 18.13 2.079 0.011 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.191 Deviance explained = 29.8%
## GCV = 4.1527 Scale est. = 3.5725 n = 122
#Plot two GAMs side-by-side
#dev.new(width=5.95, height=4)
layout(matrix(1:2,ncol=2), width = c(1,1),height = c(1,1))
par(mar=c(4,2.5,2.75,0.75), mgp=c(1.5,0.5,0), las=0)
#par(mgp=c(axis.title.position, axis.label.position, axis.line.position))
my.vis.gam(mean.cooc.gam.sqrt.te, type="response", plot.type = "contour", color="gray", ylab="Average niche width", main="",
xlab="", cex.lab=1.1, cex.axis=0.9, cex.main=0.1, lwd=2, cex=2)
box(bty="o", lwd=2)
my.vis.gam(mean.cooc.gam.30,type="response", plot.type = "contour", color="gray", ylab="Average niche width", main="",
xlab="", cex.lab=1.1, cex.axis=0.9, cex.main=0.1, lwd=2, cex=2)
box(bty="o", lwd=2)
#textClick("(a)", cex=c(size=1, font=2))
#textClick("(b)", cex=1)
#textClick("Phylogenetic distance", cex=1)
#textClick("Phylogenetic distance", cex=1)
Partial mantel for average of gen.spec
# Partial mantel
mean.sedges.partial.mantel<-mantel.partial(phy.dist.cyp, mean.final.h.inv, z.mat, method = "pearson", permutations = 999)
mean.sedges.partial.mantel
##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = phy.dist.cyp, ydis = mean.final.h.inv, zdis = z.mat, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.1186
## Significance: 0.181
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.162 0.194 0.229 0.259
## Permutation: free
## Number of permutations: 999
#Correlation between phy.dist and cooccurrence, conditional on gen.spec
mean.sedges.partial.mantel.phy.dist<-mantel.partial(phy.dist.cyp,z.mat, mean.final.h.inv,method = "pearson", permutations = 999)
mean.sedges.partial.mantel.phy.dist
##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = phy.dist.cyp, ydis = z.mat, zdis = mean.final.h.inv, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.1038
## Significance: 0.005
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0582 0.0734 0.0844 0.0990
## Permutation: free
## Number of permutations: 999
#Correlation between gen.spec and cooccurrence, conditional on phy.dist
mean.sedges.partial.mantel.gen.spec<-mantel.partial(mean.final.h.inv,z.mat,phy.dist.cyp ,method = "pearson", permutations = 999)
mean.sedges.partial.mantel.gen.spec
##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = mean.final.h.inv, ydis = z.mat, zdis = phy.dist.cyp, method = "pearson", permutations = 999)
##
## Mantel statistic r: -0.0432
## Significance: 0.793
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0656 0.0803 0.0969 0.1140
## Permutation: free
## Number of permutations: 999
Mantel for average of gen.spec
mean.sedges.phy.dist.mantel<-mantel(phy.dist.cyp, z.mat, method = "pearson", permutations = 999)
mean.sedges.phy.dist.mantel
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = phy.dist.cyp, ydis = z.mat, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.09951
## Significance: 0.006
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0597 0.0744 0.0856 0.0967
## Permutation: free
## Number of permutations: 999
mean.sedges.gen.spec.mantel<-mantel(mean.final.h.inv, z.mat,method = "pearson", permutations = 999)
mean.sedges.gen.spec.mantel
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = mean.final.h.inv, ydis = z.mat, method = "pearson", permutations = 999)
##
## Mantel statistic r: -0.03127
## Significance: 0.714
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0675 0.0923 0.1061 0.1233
## Permutation: free
## Number of permutations: 999
Plot of average differences between specialism and generalism
layout(matrix(1:2,nrow=1),widths=c(0.8,0.2))
colfunc <- colorRampPalette(c("grey90", "black"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(x=x,y=y, col = cols, pch=16,cex=1.25, xlab="Phylogenetic Distance", ylab="Average Generalism/specialism")
xl <- 1
yb <- 1
xr <- 1.5
yt <- 2
par(mar=c(5.1,0.5,4.1,0.5))
plot(NA,type="n",ann=FALSE,xlim=c(1,2),ylim=c(1,2),xaxt="n",yaxt="n",bty="n")
rect(
xl,
head(seq(yb,yt,(yt-yb)/20),-1),
xr,
tail(seq(yb,yt,(yt-yb)/20),-1),
col=colfunc(20)
)
text(1.45,0.975,"Co-occurrence",cex=1)
text(1.7,1.02, "-6.5", cex=1)
text(1.7, 2.0, "12", cex=1 )
layout(matrix(1:2,nrow=1),widths=c(0.8,0.2))
colfunc <- colorRampPalette(c("grey90", "black"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(x=x.30,y=y.30, col = cols, pch=16,cex=1.25, xlab="Phylogenetic Distance", ylab="Average Generalism/specialism")
xl <- 1
yb <- 1
xr <- 1.5
yt <- 2
par(mar=c(5.1,0.5,4.1,0.5))
plot(NA,type="n",ann=FALSE,xlim=c(1,2),ylim=c(1,2),xaxt="n",yaxt="n",bty="n")
rect(
xl,
head(seq(yb,yt,(yt-yb)/20),-1),
xr,
tail(seq(yb,yt,(yt-yb)/20),-1),
col=colfunc(20)
)
text(1.45,0.975,"Co-occurrence",cex=1)
text(1.7,1.02, "-4.5", cex=1)
text(1.7, 2.0, "9.5", cex=1 )
