########################################################################################## # COMMUNITY ECOLOGY PRACTICAL (Anna Norberg & Nerea Abrego) ########################################################################################## # 1 Data ########################################################################################## # acquire data data("dune") data("dune.env") # information about the data ?dune # trnasforming abundances to presence-absence dune.pa<-(dune>0)*1 dune.env names(dune.env) # plotting sp~environment par(mfrow=c(1,3)) plot(dune.pa[,1]~dune.env[,1],xlab=colnames(dune.env)[1],ylab=colnames(dune)[1]) plot(dune.pa[,1]~jitter(as.numeric(dune.env[,2]),amount=0.1),xlab=colnames(dune.env)[2],ylab=colnames(dune)[1]) plot(dune.pa[,1]~jitter(as.numeric(dune.env[,5]),amount=0.1),xlab=colnames(dune.env)[5],ylab=colnames(dune)[1]) mtext(text=paste(colnames(dune)[1],"~ continuous environmental characteristics"),outer=TRUE,side=3,line=-2) ########################################################################################## # 2 Single species model ########################################################################################## # 2.1 Format data and define priors #=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= dune.env.num<-dune.env[,c(1,2,5)] # let's choose the soil thickness, moisture and manure dune.env.num<-apply(dune.env.num,c(1,2),as.numeric) # convert them all to numeric and continuous dune_f <- as.HMSCdata(Y=as.matrix(dune.pa[,1]), X=dune.env.num, interceptX=TRUE, interceptTr=FALSE, scaleX=TRUE, scaleTr=FALSE) # coverting data to right format, adding intercept and scaling the covariates to mean=0 and sd=1 priors <- as.HMSCprior(HMSCdata=dune_f) # by only inserting the HMSCdata object, you get default priors # 2.2 Estimation #=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= hmsc_fit <- hmsc(data=dune_f, # HMSCdata object priors=priors, # HMSCprior object family="probit", # define the link function niter=1500, # define number of MCMC iterations nburn=200, # define the number of MCMC iterations used in burning thin=1, # define thinning verbose=TRUE) # print the iterations or not # 2.3 Estimates #=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= names(hmsc_fit$estimation) # all the results gained after the burning Betas <- hmsc_fit$estimation$paramX # parameter estimates for the regression coefficients quantifying the effect of the environment ncovar<-dim(Betas)[2] # number of covariates (including the intercept) # mixing plots for the parameter estimates par(mfrow=c(1,ncovar),oma=c(1,1,2,1),mar=c(1,2,2,2)) for (i in 1:ncovar) { plot(Betas[1,i,],type="l",main=dimnames(Betas)[[2]][i]) } mtext(text=paste(colnames(dune)[1],"~ environmental characteristics"),outer=TRUE,side=3,line=-0.5) # posterior means Betas_post_means<-apply(Betas,c(1,2),mean) # 95% credible intervals Betas_CIs<-apply(Betas,c(1,2),quantile, probs=c(0.025,0.975)) # posterior plot plot(x=0,y=0,ylim=c(min(Betas_CIs),max(Betas_CIs)),xlim=c(0.5,ncovar+0.5), type="n",xaxt="n",xlab="",ylab="",main="Posterior means and 95% credible intervals",cex.main=1.5,cex.lab=1.5,cex=1.5,cex.axis=1.5) for (i in 1:ncovar) { points(y=Betas_post_means[i], x=i, main=colnames(Betas_post_means)[i], col="grey50",type="p",pch=16) lines(y=c(Betas_CIs[1,,i],Betas_CIs[2,,i]),x=c(i,i),lwd=1.5) } axis(side=1,at=c(1:ncovar),dimnames(Betas)[[2]],cex.axis=1.5) abline(h=0, lty=2) ##########################################################################################