# set working directory to your local R directory
setwd("C:/Users/Paul/Documents/Rwork")
# load vegan package - vegan does a wide range of community analyses
library(lattice)
library(permute)
library(vegan)
## This is vegan 2.4-4
# read in data files into data frames
bfly_dat<-read.table("butterflies.txt", header=TRUE)
hab_dat<-read.table("habitat.txt", header=TRUE)
# use structure command to get a summary of the data frames
str(bfly_dat)
## 'data.frame': 26 obs. of 51 variables:
## $ patch : Factor w/ 26 levels "p10CH","p11TK",..: 11 19 20 21 22 23 24 25 26 1 ...
## $ cluster: Factor w/ 6 levels "BR","CH","CL",..: 1 1 1 1 4 4 4 2 2 2 ...
## $ EUMA : int 0 0 0 0 0 0 0 1 0 1 ...
## $ PAPO : int 0 0 0 0 0 0 0 1 1 0 ...
## $ PACR : int 0 0 0 1 0 0 0 0 0 0 ...
## $ PAGL : int 2 2 1 7 2 2 5 8 4 0 ...
## $ PATR : int 4 3 0 2 1 3 2 4 3 2 ...
## $ BAPH : int 3 1 0 5 0 0 5 1 0 1 ...
## $ COPH : int 6 2 3 3 2 0 1 2 2 2 ...
## $ COEU : int 7 0 0 5 1 0 4 0 5 3 ...
## $ PHSE : int 0 0 0 2 0 0 0 0 1 0 ...
## $ EULI : int 0 0 0 0 0 0 0 0 0 0 ...
## $ PIRA : int 1 0 0 0 0 0 0 0 1 0 ...
## $ DAPL : int 0 1 0 0 0 0 0 0 0 0 ...
## $ ENAN : int 0 0 0 0 0 0 1 0 0 0 ...
## $ CYGE : int 1 0 2 1 0 1 1 0 5 0 ...
## $ MECY : int 0 6 2 2 3 0 3 0 12 1 ...
## $ HESO : int 0 0 0 0 1 0 0 0 0 0 ...
## $ CEPE : int 0 4 0 0 0 0 6 1 4 2 ...
## $ ASCE : int 0 0 0 0 0 0 0 0 0 0 ...
## $ JUCO : int 0 0 0 0 0 0 0 0 0 0 ...
## $ LIAR : int 0 0 0 0 0 0 0 0 0 0 ...
## $ POIN : int 0 0 0 0 0 0 0 0 0 0 ...
## $ POCO : int 0 1 2 0 0 0 1 1 1 1 ...
## $ PHTH : int 23 4 6 10 9 8 15 6 15 2 ...
## $ SPCY : int 2 0 2 1 1 1 10 3 0 1 ...
## $ SPAP : int 0 0 0 0 0 0 0 0 0 0 ...
## $ CHNY : int 0 0 0 1 0 0 0 0 0 0 ...
## $ SACA : int 0 0 0 0 0 0 0 0 0 0 ...
## $ CACE : int 2 1 0 1 1 2 3 0 1 0 ...
## $ STME : int 1 1 0 0 0 0 2 0 1 0 ...
## $ EVCO : int 4 1 1 4 0 0 1 3 1 2 ...
## $ CELA : int 0 1 3 4 1 0 0 0 0 0 ...
## $ CABO : int 0 0 0 0 0 0 0 0 0 0 ...
## $ ANNU : int 0 0 1 0 0 0 0 0 0 0 ...
## $ POPE : int 0 1 0 0 2 0 1 0 1 0 ...
## $ POTH : int 0 0 0 0 0 0 0 0 0 0 ...
## $ POOR : int 0 0 0 3 0 1 0 0 0 0 ...
## $ WAEG : int 1 0 1 0 0 0 2 2 0 0 ...
## $ POVE : int 0 0 0 0 0 0 0 1 1 0 ...
## $ ATCA : int 0 0 0 0 0 0 0 1 0 0 ...
## $ ATLO : int 2 0 0 0 0 1 0 1 0 0 ...
## $ POZA : int 4 0 1 0 0 0 0 0 3 2 ...
## $ EUVE : int 4 0 0 2 0 0 0 3 1 1 ...
## $ NAIL : int 0 0 0 0 0 0 0 1 0 0 ...
## $ AMVI : int 0 0 1 0 0 0 0 1 1 1 ...
## $ EPCL : int 0 0 1 0 3 0 0 0 2 0 ...
## $ THPY : int 3 1 0 0 1 1 0 0 0 0 ...
## $ THBA : int 0 0 0 1 2 0 0 0 0 0 ...
## $ ERSP : int 0 0 0 0 0 1 0 0 0 0 ...
## $ ERHO : int 3 0 0 0 0 0 1 0 0 0 ...
str(hab_dat)
## 'data.frame': 26 obs. of 4 variables:
## $ area : num 1.838 0.087 0.063 0.331 0.026 ...
## $ connect: num 1.49 2.66 2.45 1.87 1.77 ...
## $ flowers: int 218 178 160 269 42 52 249 230 143 142 ...
## $ hosts : int 17 13 15 18 16 15 19 19 22 19 ...
# extract some useful information from data frame
ncol<-length(names(bfly_dat))
totspec<-ncol-2
nsamp<-length(rownames(bfly_dat))
# create a species x sample matrix for community analysis
specsamp<-as.matrix(bfly_dat[,3:ncol])
rownames(specsamp)<-bfly_dat$patch
# calculate summary statistics for abundance and richness by patch
richcount<-function(x) length(x[x>0])
sampabun<-apply(specsamp,MARGIN=1,sum)
samprich<-apply(specsamp,MARGIN=1,richcount)
specabun<-apply(specsamp,MARGIN=2,sum)
hist(sampabun)
hist(samprich)
# calculate summary statistics for abundance and richness by cluster
clustabun<-tapply(sampabun,bfly_dat$cluster,sum)
clustabun
## BR CH CL HP LX TK
## 185 129 364 115 441 100
clustmat<-rowsum(specsamp,bfly_dat$cluster)
clustrich<-apply(clustmat,MARGIN=1,richcount)
clustrich
## BR CH CL HP LX TK
## 32 28 40 26 36 28
# calculate alpha and beta diversity for patches and cluster
alpha_patch<-mean(samprich)
alpha_cluster<-mean(clustrich)
beta_add_patch<-alpha_cluster-alpha_patch
beta_add_cluster<-totspec-alpha_cluster
beta_mult_patch<-alpha_cluster/alpha_patch
beta_mult_cluster<-totspec/alpha_cluster
alpha_patch
## [1] 16.69231
alpha_cluster
## [1] 31.66667
beta_add_patch
## [1] 14.97436
beta_add_cluster
## [1] 17.33333
beta_mult_patch
## [1] 1.897081
beta_mult_cluster
## [1] 1.547368
# create a rank-abundance plot
specrank<-sort(specabun,decreasing=TRUE)
specrank<-specrank[1:(length(specrank))]
totabun<-sum(specsamp)
nrank<-1:totspec
spec_ra<-specrank/totabun
lognorm<-rad.lognormal(specabun)
lognorm<-lognorm$fitted.values/totabun
preempt<-rad.preempt(specabun)
preempt<-preempt$fitted.values/totabun
plot(spec_ra~nrank,log="y",xlab="Species Rank",ylab="Relative Abundance")
lines(nrank,lognorm)
lines(nrank,preempt,lty=2)
legend(30,0.1,c("Log Normal","Niche Pre-Emption"),lty=c(1,2))
# create a species rarefaction curve and estimate "true" species richness
rc<-specaccum(specsamp,method="rarefaction")
cum_ind<-rc$individuals
cum_rich<-rc$richness
rich_sd<-rc$sd
rich_low<-cum_rich-1.96*rich_sd
rich_hi<-cum_rich+1.96*rich_sd
plot(cum_rich~cum_ind,type="o",xlab="Number of individuals",ylab="Species richness",
xlim=c(0,1400),ylim=c(10,55))
segments(cum_ind,rich_low,cum_ind,rich_hi)
rich_est<-specpool(specsamp)
abline(a=rich_est$chao,b=0,lty=2)
abline(a=rich_est$jack1,b=0,lty=3)
text(1100,52,"Chao Estimate")
text(1100,55,"Jackknife Estimate")
# Try a few regressions of habitat variables on species richness
# First a simple species-area relationship
m1<-lm(log(samprich)~log(hab_dat$area))
summary(m1)
##
## Call:
## lm(formula = log(samprich) ~ log(hab_dat$area))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.47179 -0.08134 0.04584 0.13562 0.36152
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.01453 0.08370 36.018 < 2e-16 ***
## log(hab_dat$area) 0.11912 0.03733 3.191 0.00393 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2175 on 24 degrees of freedom
## Multiple R-squared: 0.2979, Adjusted R-squared: 0.2686
## F-statistic: 10.18 on 1 and 24 DF, p-value: 0.003926
fitted.m1<-fitted(m1)
plot(samprich~hab_dat$area,log="xy",xlab="Area(ha)",ylab="Species richness")
lines(hab_dat$area,exp(fitted.m1))
m2<-lm(log(samprich)~hab_dat$connect)
summary(m2)
##
## Call:
## lm(formula = log(samprich) ~ hab_dat$connect)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.60983 -0.12903 0.06389 0.11728 0.52714
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.86027 0.21595 13.24 1.58e-12 ***
## hab_dat$connect -0.03713 0.10324 -0.36 0.722
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2589 on 24 degrees of freedom
## Multiple R-squared: 0.005362, Adjusted R-squared: -0.03608
## F-statistic: 0.1294 on 1 and 24 DF, p-value: 0.7222
m3<-lm(log(samprich)~log(hab_dat$flowers))
summary(m3)
##
## Call:
## lm(formula = log(samprich) ~ log(hab_dat$flowers))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.43294 -0.15134 -0.00064 0.13908 0.43787
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.90996 0.34515 5.534 1.08e-05 ***
## log(hab_dat$flowers) 0.17910 0.07006 2.556 0.0173 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2302 on 24 degrees of freedom
## Multiple R-squared: 0.214, Adjusted R-squared: 0.1813
## F-statistic: 6.536 on 1 and 24 DF, p-value: 0.01732
m4<-lm(log(samprich)~log(hab_dat$hosts))
summary(m4)
##
## Call:
## lm(formula = log(samprich) ~ log(hab_dat$hosts))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52955 -0.15731 0.07279 0.14124 0.40373
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.3034 0.5974 2.182 0.0392 *
## log(hab_dat$hosts) 0.5024 0.2020 2.487 0.0202 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2315 on 24 degrees of freedom
## Multiple R-squared: 0.2049, Adjusted R-squared: 0.1718
## F-statistic: 6.184 on 1 and 24 DF, p-value: 0.02024
m5<-lm(log(samprich)~log(hab_dat$area)+log(hab_dat$hosts))
summary(m5)
##
## Call:
## lm(formula = log(samprich) ~ log(hab_dat$area) + log(hab_dat$hosts))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.45879 -0.09194 0.05048 0.10523 0.38646
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.03814 0.62339 3.269 0.00337 **
## log(hab_dat$area) 0.09516 0.03927 2.423 0.02366 *
## log(hab_dat$hosts) 0.31545 0.19969 1.580 0.12783
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2111 on 23 degrees of freedom
## Multiple R-squared: 0.3666, Adjusted R-squared: 0.3115
## F-statistic: 6.656 on 2 and 23 DF, p-value: 0.00524
# Create a pairwise dissimilarity matrix
# Try it with method="bray" and method="jaccard", and with bray and a 0.25 power transformation
bfly_dist<-vegdist(specsamp^0.25,method="bray")
bfly_dist
## p1BR p2BR p3BR p4BR p5HP p6HP
## p2BR 0.5048577
## p3BR 0.5430266 0.5283207
## p4BR 0.4008339 0.4876173 0.5132401
## p5HP 0.5429288 0.3912076 0.5046936 0.4165500
## p6HP 0.4660400 0.5840756 0.6659329 0.5257201 0.5002420
## p7HP 0.3375371 0.3239224 0.4883728 0.4042651 0.4755070 0.5771989
## p8CH 0.4554558 0.5215481 0.5158059 0.5378349 0.6796713 0.6248694
## p9CH 0.4140405 0.3876513 0.4649403 0.4507994 0.4918254 0.6658811
## p10CH 0.4763538 0.4617694 0.4664344 0.4707424 0.5959493 0.7481389
## p11TK 0.4691921 0.4939783 0.5102992 0.4479906 0.4711939 0.5821432
## p12TK 0.7735074 0.6517823 0.6397265 0.7543109 0.4683509 0.6451050
## p13TK 0.5934577 0.6091540 0.4987260 0.5226617 0.4748482 0.7227104
## p14LX 0.5978263 0.5252858 0.7352103 0.5697935 0.4429608 0.6252413
## p15LX 0.4555030 0.4523365 0.5507462 0.4521491 0.4438905 0.6799990
## p16LX 0.4170527 0.4927008 0.3502513 0.4781222 0.4771381 0.5417380
## p17LX 0.5483589 0.4199869 0.3873936 0.5640114 0.4445733 0.7029833
## p18LX 0.6233016 0.5184423 0.4517906 0.6100267 0.5193037 0.7217959
## p19LX 0.4944103 0.4157984 0.5799965 0.5092367 0.4130345 0.6295178
## p20LX 0.5166476 0.2634856 0.4644675 0.4402261 0.4231453 0.6282188
## p21LX 0.6575202 0.3984862 0.4628214 0.4444201 0.4776182 0.5987655
## p22CL 0.3803065 0.4688162 0.4707839 0.4431650 0.4684459 0.6513361
## p23CL 0.5634834 0.6230281 0.4571839 0.6121269 0.5761919 0.7017578
## p24CL 0.3644667 0.4295793 0.3310057 0.4181743 0.4681710 0.6666618
## p25CL 0.4882161 0.5797231 0.5662920 0.4537289 0.5299536 0.6961945
## p26CL 0.5801448 0.5092081 0.5885525 0.5602313 0.5346331 0.6640765
## p7HP p8CH p9CH p10CH p11TK p12TK
## p2BR
## p3BR
## p4BR
## p5HP
## p6HP
## p7HP
## p8CH 0.4645141
## p9CH 0.3529208 0.4647166
## p10CH 0.4326585 0.3547126 0.4034518
## p11TK 0.3988694 0.4730277 0.4201489 0.4374305
## p12TK 0.7565433 0.8247145 0.7091536 0.7390172 0.5662788
## p13TK 0.5386992 0.5123302 0.3944259 0.4303293 0.4757348 0.6862818
## p14LX 0.4909726 0.4988437 0.4942159 0.5160519 0.3536024 0.5354286
## p15LX 0.3358194 0.5547760 0.4324000 0.4732619 0.4165278 0.6580244
## p16LX 0.3501555 0.4715769 0.4703100 0.4364500 0.3724521 0.6117193
## p17LX 0.3789078 0.4248073 0.4398020 0.5015961 0.4768713 0.6195947
## p18LX 0.4856413 0.5203776 0.5879469 0.4168852 0.5224900 0.5640589
## p19LX 0.3921752 0.5843804 0.5292918 0.5136322 0.5249559 0.6163634
## p20LX 0.3176426 0.4343996 0.4321506 0.4744351 0.4842206 0.6372500
## p21LX 0.4586152 0.4837956 0.5326561 0.4367365 0.4961768 0.5814688
## p22CL 0.3266891 0.5009465 0.4159172 0.4803003 0.4272161 0.6849786
## p23CL 0.5191652 0.4821693 0.6200347 0.5696708 0.5327747 0.6581208
## p24CL 0.3544160 0.4178620 0.3436166 0.3428263 0.3634316 0.5980692
## p25CL 0.4802387 0.4228279 0.6494804 0.4637781 0.5733720 0.7745236
## p26CL 0.4694395 0.6463139 0.5490857 0.5194183 0.3519226 0.5848980
## p13TK p14LX p15LX p16LX p17LX p18LX
## p2BR
## p3BR
## p4BR
## p5HP
## p6HP
## p7HP
## p8CH
## p9CH
## p10CH
## p11TK
## p12TK
## p13TK
## p14LX 0.5504063
## p15LX 0.4493190 0.4287170
## p16LX 0.5541809 0.5157374 0.4216641
## p17LX 0.5438798 0.5008861 0.4104789 0.4323427
## p18LX 0.6769924 0.4380328 0.4546411 0.3107559 0.3366977
## p19LX 0.6368964 0.3552087 0.3928882 0.4320647 0.4217489 0.3759954
## p20LX 0.6160277 0.4588400 0.3870462 0.3790062 0.3239093 0.3668108
## p21LX 0.6287398 0.4450896 0.4303527 0.4258669 0.2719066 0.2515546
## p22CL 0.4958121 0.5146273 0.3968295 0.4020976 0.4099299 0.4717395
## p23CL 0.5124556 0.5441246 0.4977945 0.4094096 0.4448257 0.4306863
## p24CL 0.4566686 0.4723662 0.3872896 0.3867203 0.3797128 0.4068510
## p25CL 0.7297495 0.5469830 0.5308854 0.3845409 0.4072291 0.3760107
## p26CL 0.6165401 0.4117374 0.4946048 0.4169585 0.4943768 0.4448968
## p19LX p20LX p21LX p22CL p23CL p24CL
## p2BR
## p3BR
## p4BR
## p5HP
## p6HP
## p7HP
## p8CH
## p9CH
## p10CH
## p11TK
## p12TK
## p13TK
## p14LX
## p15LX
## p16LX
## p17LX
## p18LX
## p19LX
## p20LX 0.3570982
## p21LX 0.3267641 0.3122489
## p22CL 0.3118888 0.4325231 0.4785683
## p23CL 0.6224766 0.5330998 0.5260812 0.5951137
## p24CL 0.3670338 0.4688510 0.4079470 0.2718700 0.5492458
## p25CL 0.3687000 0.4346370 0.4014865 0.4462090 0.6538235 0.4289148
## p26CL 0.4079548 0.4692783 0.4698984 0.4393239 0.5737891 0.4609106
## p25CL
## p2BR
## p3BR
## p4BR
## p5HP
## p6HP
## p7HP
## p8CH
## p9CH
## p10CH
## p11TK
## p12TK
## p13TK
## p14LX
## p15LX
## p16LX
## p17LX
## p18LX
## p19LX
## p20LX
## p21LX
## p22CL
## p23CL
## p24CL
## p25CL
## p26CL 0.5679179
mean(bfly_dist)
## [1] 0.4950415
# Do a simple ordination using non-metric multidimensional scaling
bfly_nmds<-metaMDS(specsamp^0.25,distance="bray",k=2,trmax=20,autotransform=FALSE)
## Run 0 stress 0.2278945
## Run 1 stress 0.2862206
## Run 2 stress 0.2484412
## Run 3 stress 0.2908836
## Run 4 stress 0.2278925
## ... New best solution
## ... Procrustes: rmse 0.0009565999 max resid 0.003490566
## ... Similar to previous best
## Run 5 stress 0.2338432
## Run 6 stress 0.2278925
## ... Procrustes: rmse 5.384086e-05 max resid 0.0002287994
## ... Similar to previous best
## Run 7 stress 0.3869948
## Run 8 stress 0.2816801
## Run 9 stress 0.2450963
## Run 10 stress 0.2581451
## Run 11 stress 0.2278925
## ... Procrustes: rmse 7.465985e-05 max resid 0.0002330791
## ... Similar to previous best
## Run 12 stress 0.2611744
## Run 13 stress 0.2278925
## ... Procrustes: rmse 0.0001436274 max resid 0.0004715374
## ... Similar to previous best
## Run 14 stress 0.2454978
## Run 15 stress 0.2487352
## Run 16 stress 0.2480336
## Run 17 stress 0.2445437
## Run 18 stress 0.2599537
## Run 19 stress 0.2278932
## ... Procrustes: rmse 0.0005156531 max resid 0.001925015
## ... Similar to previous best
## Run 20 stress 0.2445792
## *** Solution reached
bfly_nmds
##
## Call:
## metaMDS(comm = specsamp^0.25, distance = "bray", k = 2, autotransform = FALSE, trmax = 20)
##
## global Multidimensional Scaling using monoMDS
##
## Data: specsamp^0.25
## Distance: bray
##
## Dimensions: 2
## Stress: 0.2278925
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'specsamp^0.25'
plot(bfly_nmds,type="text",display="sites",cex=1.1)
# Do a constrained ordination distance-based redundancy analysis using environmental variables
bfly_dbrda<-capscale(specsamp^0.25~log(hab_dat$area)+hab_dat$connect,dist="bray")
anova(bfly_dbrda,by="term")
## Permutation test for capscale under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: capscale(formula = specsamp^0.25 ~ log(hab_dat$area) + hab_dat$connect, distance = "bray")
## Df SumOfSqs F Pr(>F)
## log(hab_dat$area) 1 0.32701 2.8880 0.001 ***
## hab_dat$connect 1 0.26687 2.3568 0.002 **
## Residual 23 2.60430
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bfly_dbrda
## Call: capscale(formula = specsamp^0.25 ~ log(hab_dat$area) +
## hab_dat$connect, distance = "bray")
##
## Inertia Proportion Eigenvals Rank
## Total 3.1982 1.0000 3.5009
## Constrained 0.5939 0.1857 0.6051 2
## Unconstrained 2.6043 0.8143 2.8958 18
## Imaginary -0.3027 7
## Inertia is squared Bray distance
##
## Eigenvalues for constrained axes:
## CAP1 CAP2
## 0.3333 0.2718
##
## Eigenvalues for unconstrained axes:
## MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
## 0.4649 0.4225 0.3641 0.3259 0.2614 0.2225 0.1931 0.1694
## (Showed only 8 of all 18 unconstrained eigenvalues)
plot(bfly_dbrda,display=c("sites","bp"),cex=1.1)
