Schefferville - generalists/specialists; pairwise co-occurrence; modelling ## ; gen-spec based on jaccard rankings
Includes OLS analyses with 34 focal species (including T. alpinum)
Author: Tammy Elliott
Date: March 4, 2016
R version 3.1.0
GAM for absolute difference of generalism/specialism (jaccard)
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.1580 0.1153 -1.371 0.171
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## te(phy.dist.sqrt,gen.spec) 4.488 5.186 7.954 2.18e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0655 Deviance explained = 7.3%
## GCV = 7.53 Scale est. = 7.4563 n = 561
gam.check(cooc.gam.sqrt.te)
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 15 iterations.
## The RMS GCV score gradiant at convergence was 8.284997e-07 .
## The Hessian was positive definite.
## The estimated model rank was 25 (maximum possible: 25)
## Model rank = 25 / 25
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## te(phy.dist.sqrt,gen.spec) 24.00 4.49 1.00 0.54
GAM for mean of generalism/specialism (jaccard)
mean.phy.dist.sqrt<-sqrt(mean.data.cooc$mean.phy.dist)
mean.genspec.te<-gam(mean.cooc.z~te(mean.phy.dist.sqrt, mean.gen.spec), dat=mean.data.cooc)
summary(mean.genspec.te)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## mean.cooc.z ~ te(mean.phy.dist.sqrt, mean.gen.spec)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1580 0.1145 -1.38 0.168
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## te(mean.phy.dist.sqrt,mean.gen.spec) 8.167 8.723 5.897 1.17e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0786 Deviance explained = 9.2%
## GCV = 7.4745 Scale est. = 7.3523 n = 561
gam.check(mean.genspec.te)
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 13 iterations.
## The RMS GCV score gradiant at convergence was 4.048159e-07 .
## The Hessian was positive definite.
## The estimated model rank was 25 (maximum possible: 25)
## Model rank = 25 / 25
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## te(mean.phy.dist.sqrt,mean.gen.spec) 24.00 8.17 1.01 0.66
GAM for absolute difference of generalism and specialism < 30 my (jaccard co-occurrence)
cooc.gam.30.te<-gam(cooc.z.30~te(phy.dist.30.sqrt,gen.spec.30), dat=data.cooc.30)
gam.check(cooc.gam.30.te, k.rep = 1000)
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 15 iterations.
## The RMS GCV score gradiant at convergence was 5.098094e-07 .
## The Hessian was positive definite.
## The estimated model rank was 25 (maximum possible: 25)
## Model rank = 25 / 25
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## te(phy.dist.30.sqrt,gen.spec.30) 24.00 3.80 1.03 0.6
summary(cooc.gam.30.te)
##
## 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.6962 0.2265 -3.073 0.00263 **
## ---
## 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.799 4.293 2.329 0.0559 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0568 Deviance explained = 8.64%
## GCV = 6.5176 Scale est. = 6.2612 n = 122
Gam for mean gen.spec < 30 my
mean.cooc.gam.30<-gam(mean.cooc.z.30~te(mean.phy.dist.30.sqrt,mean.gen.spec.30), dat=mean.data.cooc.30)
gam.check(mean.cooc.gam.30, k.rep = 1000)
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 15 iterations.
## The RMS GCV score gradiant at convergence was 5.70918e-07 .
## The Hessian was positive definite.
## The estimated model rank was 25 (maximum possible: 25)
## Model rank = 25 / 25
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## te(mean.phy.dist.30.sqrt,mean.gen.spec.30) 24.00 4.59 1.04 0.62
summary(mean.cooc.gam.30)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## mean.cooc.z.30 ~ te(mean.phy.dist.30.sqrt, mean.gen.spec.30)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6962 0.2246 -3.1 0.00243 **
## ---
## 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,mean.gen.spec.30) 4.595 5.389 2.149 0.0598 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0729 Deviance explained = 10.8%
## GCV = 6.4502 Scale est. = 6.1544 n = 122
Plot four GAMs
#dev.new(width=4.9, height=6)
#layout(matrix(1:4,ncol=2), width = c(1,1),height = c(1,1))
par(mfrow=c(2,2))
par(mar=c(3.5,3,2.75,0.75), mgp=c(1.5,0.3,0), las=0)
my.vis.gam(cooc.gam.sqrt.te, axes=FALSE, type="response", plot.type = "contour", color="gray", main="", mgp=c(1.75,0.5,0),xlab="Phylogenetic distance (sqrt)",cex.main=0.1, lwd=2,cex.lab=1.15, cex=2, ylab="Difference in niche width", font.lab=1)
axis(1, tck=0.025, cex.lab=1.25, cex.axis=1)
axis(2, tck=0.025, cex.lab=1.25, cex.axis=1)
box(bty="o", lwd=2)
my.vis.gam(cooc.gam.30.te, axes=FALSE, type="response", plot.type = "contour", color="gray", main="", mgp=c(1.75,0.5,0),xlab="Phylogenetic distance (sqrt)",cex.main=0.1, lwd=2,cex.lab=1.15, cex=2, ylab="Difference in niche width", font.lab=1)
axis(1, tck=0.025, cex.lab=1.25, cex.axis=1)
axis(2, tck=0.025, cex.lab=1.25, cex.axis=1)
box(bty="o", lwd=2)
my.vis.gam(mean.genspec.te, axes=FALSE, type="response", plot.type = "contour", color="gray", main="", mgp=c(1.75,0.5,0),xlab="Phylogenetic distance (sqrt)",cex.main=0.1, lwd=2,cex.lab=1.15, cex=2, ylab="Mean niche width", font.lab=1)
axis(1, tck=0.025, cex.lab=1.25, cex.axis=1)
axis(2, tck=0.025, cex.lab=1.25, cex.axis=1)
box(bty="o", lwd=2)
my.vis.gam(mean.cooc.gam.30, axes=FALSE, type="response", plot.type = "contour", color="gray", main="", mgp=c(1.75,0.5,0),xlab="Phylogenetic distance (sqrt)",cex.main=0.1, lwd=2,cex.lab=1.15, cex=2, ylab="Mean niche width", font.lab=1)
axis(1, tck=0.025, cex.lab=1.5, cex.axis=1)
axis(2, tck=0.025, cex.lab=1.5, cex.axis=1)
box(bty="o", lwd=2)
#textClick("(a)", cex=1.25)
#textClick("(b)", cex=1.25)
#textClick("(c)", cex=1.25)
#textClick("(d)", cex=1.25)
Plot pairwise relationships
#dev.new(width=6.8, height=6.8)
#layout(matrix(1:4,ncol=2), width = c(1,1),height = c(1,1))
par(mfrow=c(2,3))
par(mar=c(3.5,3,2.75,0.75), mgp=c(1.5,0.3,0), las=0)
plot(data.cooc$gen.spec,data.cooc$cooc.z, ylab="Co-occurrence", xlab="Difference in niche width",lwd=2, pch=16, col="black", cex.axis=1.25, cex.lab=1.5,bty="n", cex=1,tck=0.03)
#abline(cooc.z.gen.spec.lm, lwd=2)
box( lwd=3)
plot(mean.data.cooc$mean.gen.spec,mean.data.cooc$mean.cooc.z, ylab="Co-occurrence", xlab="Mean niche width",lwd=2, pch=16, col="black", cex.axis=1.25, cex.lab=1.5,bty="n", cex=1,tck=0.03)
#abline(cooc.z.gen.spec.lm, lwd=2)
box( lwd=3)
plot(sqrt(data.cooc$phy.dist),data.cooc$cooc.z, ylab="Co-occurrence", xlab="Phylogenetic distance (sqrt)",lwd=2, pch=16, col="black", cex.axis=1.25, cex.lab=1.5,bty="n", cex=1,tck=0.03)
#abline(cooc.z.phy.dist.lm, lwd=2)
box( lwd=3)
plot((data.cooc.30$gen.spec),data.cooc.30$cooc.z, ylab="Co-occurrence",xlab="Difference in niche width", lwd=2, pch=16, col="black", cex.axis=1.25, cex.lab=1.5,bty="n", cex=1,tck=0.03)
#abline(cooc.z.gen.spec.lm.30, lwd=2)
box(lwd=3)
plot((mean.data.cooc.30$mean.gen.spec),mean.data.cooc.30$mean.cooc.z, ylab="Co-occurrence",xlab="Mean niche width", lwd=2, pch=16, col="black", cex.axis=1.25, cex.lab=1.5,bty="n", cex=1,tck=0.03)
#abline(cooc.z.gen.spec.lm.30, lwd=2)
box(lwd=3)
plot(sqrt(data.cooc.30$phy.dist),data.cooc.30$cooc.z, ylab="Co-occurrence", xlab="Phylogenetic distance (sqrt)",lwd=2, pch=16, col="black", cex.axis=1.25, cex.lab=1.5,bty="n", cex=1,tck=0.03)
box(lwd=3)
#abline(cooc.z.phy.dist.lm.30, lwd=2)
#textClick("(a)",cex=c(size=1.5, font=2) )
#textClick("(b)", cex=c(size=1.5, font=2))
#textClick("(c)", cex=c(size=1.5, font=2))
#textClick("(d)", cex=c(size=1.5, font=2))
#textClick("(e)", cex=c(size=1.5, font=2))
#textClick("(f)", cex=c(size=1.5, font=2))
Mantel and partial Mantel tests for mean niches (jaccard co-occurrence)
# Partial mantel
sedges.mean.partial.mantel<-mantel.partial(phy.dist.cyp, mean.jaccard, z.mat, method = "pearson", permutations = 999)
sedges.mean.partial.mantel
##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = phy.dist.cyp, ydis = mean.jaccard, zdis = z.mat, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.1901
## Significance: 0.058
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.154 0.206 0.250 0.297
## Permutation: free
## Number of permutations: 999
# Mantel (phylogenetic distance and co-occurrence)
sedges.mean.phy.dist.mantel<-mantel(phy.dist.cyp, z.mat, method = "pearson", permutations = 999)
sedges.mean.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.0957
## Significance: 0.004
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0569 0.0651 0.0728 0.0803
## Permutation: free
## Number of permutations: 999
# Mantel (specialism/generalism and co-occurrence)
sedges.mean.gen.spec.mantel<-mantel(mean.jaccard, z.mat,method = "pearson", permutations = 999)
sedges.mean.gen.spec.mantel
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = mean.jaccard, ydis = z.mat, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.09204
## Significance: 0.033
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0655 0.0818 0.0962 0.1093
## Permutation: free
## Number of permutations: 999
Mantel and partial Mantel tests for absolute differences in niches
# Partial mantel
sedges.partial.mantel<-mantel.partial(phy.dist.cyp,final.jaccard, z.mat, method = "pearson", permutations = 999)
sedges.partial.mantel
##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = phy.dist.cyp, ydis = final.jaccard, zdis = z.mat, method = "pearson", permutations = 999)
##
## Mantel statistic r: 0.01593
## Significance: 0.429
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.130 0.171 0.209 0.237
## Permutation: free
## Number of permutations: 999
# Mantel (phylogenetic distance and co-occurrence)
sedges.phy.dist.mantel<-mantel(phy.dist.cyp, z.mat, method = "pearson", permutations = 999)
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.0957
## Significance: 0.008
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0578 0.0708 0.0784 0.0908
## Permutation: free
## Number of permutations: 999
#Mantel (specialism/generalism and co-occurrence)
sedges.gen.spec.mantel<-mantel(final.jaccard, z.mat,method = "pearson", permutations = 999)
sedges.gen.spec.mantel
##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = final.jaccard, ydis = z.mat, method = "pearson", permutations = 999)
##
## Mantel statistic r: -0.2006
## Significance: 1
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0584 0.0712 0.0844 0.0949
## Permutation: free
## Number of permutations: 999