Schefferville - generalists/specialists; pairwise co-occurrence; modelling

Includes OLS analyses with 34 focal species (including T. alpinum)

Author: Tammy L. Elliott

Date:December 8, 2015

R version 3.1.0

Import phylogenetic tree and make distance matrix

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.158      0.115  -1.374     0.17
## 
## Approximate significance of smooth terms:
##                              edf Ref.df     F  p-value    
## te(phy.dist.sqrt,gen.spec) 5.323  6.052 7.363 1.23e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.07   Deviance explained = 7.89%
## GCV = 7.5051  Scale est. = 7.4205    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.6962     0.2253   -3.09  0.00249 **
## ---
## 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      3 3.904  0.0106 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0672   Deviance explained = 9.03%
## GCV = 6.4022  Scale est. = 6.1922    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

# Partial mantel

sedges.partial.mantel<-mantel.partial(phy.dist.cyp, final.h.inv, 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.h.inv, zdis = z.mat,      method = "pearson", permutations = 999) 
## 
## Mantel statistic r: -0.09717 
##       Significance: 0.882 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.120 0.171 0.222 0.269 
## Permutation: free
## Number of permutations: 999

#Correlation between phy.dist and cooccurrence, conditional on generalism and specialism

sedges.partial.mantel.phy.dist<-mantel.partial(phy.dist.cyp,z.mat, final.h.inv,method = "pearson", permutations = 999)
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 = final.h.inv,      method = "pearson", permutations = 999) 
## 
## Mantel statistic r: 0.08214 
##       Significance: 0.025 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0605 0.0716 0.0815 0.0963 
## Permutation: free
## Number of permutations: 999

#Correlation between gen.spec and cooccurrence, conditional on phy.dist

sedges.partial.mantel.gen.spec<-mantel.partial(final.h.inv,z.mat,phy.dist.cyp ,method = "pearson", permutations = 999)
sedges.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.127 
##       Significance: 0.99 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0522 0.0643 0.0761 0.0883 
## Permutation: free
## Number of permutations: 999

Mantel

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.006 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0515 0.0649 0.0752 0.0885 
## Permutation: free
## Number of permutations: 999

sedges.gen.spec.mantel<-mantel(final.h.inv, z.mat,method = "pearson", permutations = 999)
sedges.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.1361 
##       Significance: 0.993 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0538 0.0690 0.0819 0.0898 
## 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.1580     0.1139  -1.387    0.166
## 
## Approximate significance of smooth terms:
##                                   edf Ref.df    F  p-value    
## te(mean.phy.dist.sqrt,gen.spec) 13.83  16.01 3.84 7.41e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0876   Deviance explained =   11%
## GCV = 7.4778  Scale est. = 7.2801    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.6962     0.2142  -3.251  0.00154 **
## ---
## 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) 14.83   16.9 1.795  0.0379 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.157   Deviance explained =   26%
## GCV =   6.43  Scale est. = 5.5959    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

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.06363 
##       Significance: 0.331 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.159 0.195 0.219 0.247 
## 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.09782 
##       Significance: 0.004 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0528 0.0673 0.0794 0.0886 
## 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.03733 
##       Significance: 0.779 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0582 0.0752 0.0913 0.1037 
## 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.0957 
##       Significance: 0.004 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0550 0.0678 0.0769 0.0891 
## 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.03132 
##       Significance: 0.697 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0598 0.0800 0.0961 0.1135 
## 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, "-7", cex=1)
text(1.7, 2.0, "12", 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 )