hmsc: Hierarchical modelling of species community

In HMSC: Hierarchical modelling of species community

Description

Hierarchical modelling of species community. hmscH is designed to model a community with habitat characteristics only whereas hmscHT also uses species traits to construct the model.

Usage

## S3 method for class 'formula'
hmscH(formula, data, start, priors = NULL, niter = 12000, nburn = 2000, nrot = 100, thin=1, likelihoodtype = "logit", details = TRUE, verbose=TRUE, ...)
## S3 method for class 'HMSCdata'
hmscH(data, start, priors = NULL, niter = 12000, nburn = 2000, nrot = 100, thin=1, likelihoodtype = "logit",spAsso = FALSE, details = TRUE, verbose=TRUE, ...)

## S3 method for class 'HMSCdata'
hmscHT(data, start, priors = NULL, niter = 12000, nburn = 2000, nrot = 100, thin=1, likelihoodtype = "logit", details = TRUE, verbose=TRUE, ...)

Arguments

formula	Model formula. The left side presents the response matrix and the right side gives the descriptors with their parameters. (See details)
data	Object of class HMSCdata.
start	Object of class HMSCparam.
priors	Object of class HMSCprior.
niter	A positive integer defining the number of iterations to be carried out. Default is 12000.
nburn	Numeric. A positive integer defining the number of iterations to be used in the burning (first) phase of the algorithm. Default is 2000.
nrot	Numeric. A positive integer defining the number of iterations that needs to be carried out before eigenvector rotation is performed.
thin	Numeric. A positive integer defining thinning. Default is 1, that is all values are considered. (see details)
likelihoodtype	Character string defining how the likelihood should be calculated. Any unambiguous variation of wording used in this argument is accepted. (See details)
spAsso	Logical. Whether species association should be estimated (TRUE) or not (FALSE).
details	Logical. Whether detailed information about the burning, the accept ratio, and the standard deviation should be given. Default is TRUE.
verbose	Logical. Whether a counter printing on the screen how many iterations have been performed so far is used or not. Default is TRUE.
...	Parameters passed to other functions.

Details

For the user the main difference between hmscH and hmscHT is within the data and start arguments.

The formula was designed so that it is necessary to include the parameters and the descriptors. This function is thus capable of carrying out non-linear analyses (although the current implement is far from being very efficient).

When building a model for the first time, the parameters provided in start can be randomly generated using rnorm[stats], as it is shown in the examples. However, nothing prevents the user to predefine the starting values so that the function can more quickly converge to optimal set of parameters.

When thinning, the value given to thin means that all but the k^th values will be considered.

Currently, the following likelihood have been implemented:

- logit: Binomial likelihood using a logit link function

- poisson: Poisson likelihood using a log link function

Value

model	List with the following items: paramX A 3-dimensional array where, for each species (first dimension), the regression parameters of X (second dimension) are stored for all iterations (third dimension) after burning. paramTr A 3-dimensional array where, for each regression parameter (first dimension), the species traits (second dimension) are stored for all iterations (third dimension) after burning. means Matrix of means measuring the response of a typical response variables to the descriptors. Each row represent an iteration result after burning. sigma A 3-dimensional array where covariance matrices measuring how the species vary in their responses to the descriptors (first and second dimension) are stored for all iterations (third dimension) after burning. sdvalue Matrix of standard deviation used when sampling for the regression parameters (columns) of each response variables (row) after burning. acceptRatio Matrix of accept ratio after burning for each species (rows) when sampling for the regression parameters of X (columns). These values should be close to 0.44. If they are too far from 0.44, the model may not be fitted well. The number of burning iterations should be increased. paramRandom A 3-dimensional array where, for each species (first dimension), the regression parameters for each levels of the random effect (second dimension) are stored for all iterations (third dimension) after burning. sdvalueRandom Matrix of standard deviation used when sampling for the regression parameters for each levels of the random effect (columns) of each response variables (row) after burning. acceptRatioRandom Matrix of accept ratio after burning for each species (rows) when sampling for the regression parameters for each levels of the random effect (columns). These values should be close to 0.44. If they are too far from 0.44, the model may not be fitted well. The number of burning iterations should be increased. RandomVar Scalar defining the within group variance in Random after burning. sdvalueRandomCov Standard deviation used when sampling for RandomCov after burning. acceptRatioRandomCov accept ratio after burning for RandomCov after burning. This values should be close to 0.44. If it is too far from 0.44, the model may not be fitted well. The number of burning iterations should be increased.
burning	List with the following items : paramX A 3-dimensional array where, for each species (first dimension), the regression parameters (second dimension) are stored for all iterations (third dimension) during burning. paramTr A 3-dimensional array where, for each species traits (second dimension), the regression parameter (first dimension) are stored for all iterations (third dimension) during burning. means Matrix of means measuring the response of a typical response variables to the descriptors. Each row represent an iteration result during burning. sigma A 3-dimensional array where covariance matrices measuring how the species vary in their responses to the descriptors (first and second dimension) are stored for all iterations (third dimension) during burning. paramRandom A 3-dimensional array where, for each species (first dimension), the regression parameters for each levels of the random effect (second dimension) are stored for all iterations (third dimension) during burning. RandomVar A vector defining the within group variance in Random during burning.
details	This part of the result is only available when the argument details is TRUE List with the following items : sdvalue A 3-dimensional array where the first two dimensions represent the standard deviation used when sampling for the regression parameters (columns) of each response variables (row) for the different burning iterations (third dimension). acceptRatio A 3-dimensional array where the first two dimensions represent the accept ratio calculated during the burning phase of the model construction when sampling for the regression parameters (columns) of each response variables (row). The accept ratio for all burning iterations are presented as the third dimension. sdvalueRandom A 3-dimensional array where the first two dimensions represent the standard deviation used when sampling for the regression parameters of each levels of the random effect (columns) for each response variables (row) for the different burning iterations (third dimension). acceptRatioRandom A 3-dimensional array where the first two dimensions represent the accept ratio calculated during the burning phase of the model construction when sampling for the regression parameters of each levels of the random effect (columns) for each response variables (row). The accept ratio for all burning iterations are presented as the third dimension. sdvalueRandomCov A vector showing how the standard deviation varies when sampling for RandomCov for the different burning iterations. acceptRatioRandomCov A vector showing how the accept ratio varies when sampling for RandomCov for the different burning iterations.
data	An object of class HMSCdata that includes the data used to construct the model. Note that if traits and/or a random effect are considered in the model, data could also have a HMSCdataTrait and/or a HMSCdataRandom class.

Note

These functions are essentially glorified for loops using scamH or scamHT and/or scamRandom, repetitively.

Author(s)

F. Guillaume Blanchet

References

Ovaskainen, O. and Soininen, J. (2011) Making more out of sparse data: hierarchical modeling of species communities. Ecology 92, 289-295.

See Also

scamH

Examples

#==================
### Simulating data
#==================
desc<-cbind(1,scale(1:50),scale(1:50)^2)
nspecies<-20
dataBase<-communitySimulH(desc,nsp=nspecies)

#=================================
### Formatting data and parameters
#=================================
formdata<-as.HMSCdata(dataBase$data$Y,desc,Ypattern="sp",interceptX=FALSE,interceptTr=FALSE)

#========================
### Formatting parameters
#========================
startParamX<-matrix(rnorm(nspecies*ncol(desc)),nrow=nspecies,ncol=ncol(desc))
startMean<-rnorm(ncol(desc))
startSigma<-rWishart(1,ncol(desc)+1,diag(ncol(desc)))[,,1]

formparam<-as.HMSCparam(formdata,paramX=startParamX,means=startMean,sigma=startSigma)
formpriors<-as.HMSCprior(formparam,rep(0,length(formparam$means)),0,length(formparam$means),diag(length(formparam$means)))

#=============================================================================================================
### Building model using a linear modelling approach (currently 2000 iterations takes about 50 seconds to run)
#=============================================================================================================
system.time(
modelLinear<-hmscH(formdata,formparam,formpriors,niter=200,nburn=100)
)

#=============================================================================================================
### Building model using a formula modelling approach (currently 2000 iterations takes about 2 minutes to run)
#=============================================================================================================
formule<-sp~p1*x1+p2*x2+p3*x3
system.time(
modelFormula<-hmscH(formule,formdata,formparam,formpriors,niter=200,nburn=100)
)

#===============
### Some results
#===============
### Linear model
summaryLin<-summary(modelLinear)
summaryLin[1:5,1:5,]

### Formula model
summaryForm<-summary(modelFormula)
summaryForm[1:5,1:5,]

#------------------------------
### Draw model from modelLinear
#------------------------------

par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	plot(summaryLin[,i,1],type="l",ylim=c(0,1),las=1,lwd=3,xaxt="n",col="red")
	title(paste("Species",i))
	lines(summaryLin[,i,3],col="red",lwd=1)
	lines(summaryLin[,i,5],col="red",lwd=1)
	lines(dataBase$probMat[,i],col="black",lwd=3)
}

#-------------------------------
### Draw model from modelFormula
#-------------------------------
dev.new()
par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	plot(summaryForm[,i,1],type="l",ylim=c(0,1),las=1,lwd=3,xaxt="n",col="red")
	title(paste("Species",i))
	lines(summaryForm[,i,3],col="red",lwd=1)
	lines(summaryForm[,i,5],col="red",lwd=1)
	lines(dataBase$probMat[,i],col="black",lwd=3)
}

#------------------
### Mixing - Linear
#------------------
#___________________
#### Community Sigma
#___________________
par(mfrow=c(3,3),mar=c(1,4,3,1))
for(i in 1:3){
	for(j in 1:3){
		mixing(modelLinear,toplot="sigma",param=c(i,j),col="red",xlab="Sites",ylab="Sigma",ylim=range(c(modelLinear$model$sigma[i,j,],dataBase$param$sigma[i,j])),las=1)
		abline(h=dataBase$param$sigma[i,j],col="red",lwd=3)
		title(paste("Community Sigma",i,"-",j))
	}
}

#__________________
### Community means
#__________________
par(mfrow=c(1,1),mar=c(5,4,4,2))
mixing(modelLinear,toplot="means",param=1,ylim=range(rbind(modelLinear$model$means,dataBase$param$means)),xlab="Sites",ylab="Means",las=1)
title("Community means")
mixing(modelLinear,toplot="means",param=2,add=TRUE,col="blue",xlab="Sites",ylab="Means",las=1)
mixing(modelLinear,toplot="means",param=3,add=TRUE,col="red",xlab="Sites",ylab="Means",las=1)
abline(h=dataBase$param$means,col=c("black","blue","red"),lwd=3)

#_____________________________
### Parameter for each species 
#_____________________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelLinear,toplot="paramX",param=c(i,1),ylim=range(cbind(modelLinear$model$paramX[i,,],dataBase$param$paramX[i,])),xlab="Sites",ylab="paramX",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelLinear,toplot="paramX",param=c(i,2),add=TRUE,col="blue")
	mixing(modelLinear,toplot="paramX",param=c(i,3),add=TRUE,col="red")
	abline(h=dataBase$param$paramX[i,],col=c("black","blue","red"),lwd=3)
}

#------------------------------------------
### Check the distribution of the posterior
#------------------------------------------
#__________________
### Community means
#__________________
pairs(modelLinear$model$means,asp=1,pch=19,cex=0.2)

#-------------------
### Mixing - Formula
#-------------------
#___________________
#### Community Sigma
#___________________
par(mfrow=c(3,3),mar=c(1,3,3,1))
for(i in 1:3){
	for(j in 1:3){
		mixing(modelFormula,toplot="sigma",param=c(i,j),col="red",xlab="Sites",ylab="Sigma",ylim=range(c(modelFormula$model$sigma[i,j,],dataBase$param$sigma[i,j])),las=1)
		abline(h=dataBase$param$sigma[i,j],col="red",lwd=3)
		title(paste("Community Sigma",i,"-",j))
	}
}

#__________________
### Community means
#__________________
par(mfrow=c(1,1),mar=c(5,4,4,2))
mixing(modelFormula,toplot="means",param=1,ylim=range(rbind(modelFormula$model$means,dataBase$param$means)),xlab="Sites",ylab="Means",las=1)
title("Community means")
mixing(modelFormula,toplot="means",param=2,add=TRUE,col="blue",xlab="Sites",ylab="Means",las=1)
mixing(modelFormula,toplot="means",param=3,add=TRUE,col="red",xlab="Sites",ylab="Means",las=1)
abline(h=dataBase$param$means,col=c("black","blue","red"),lwd=3)

#_____________________________
### Parameter for each species 
#_____________________________
par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	mixing(modelFormula,toplot="paramX",param=c(i,1),ylim=range(cbind(modelFormula$model$paramX[i,,],dataBase$param$param[i,])),xlab="Sites",ylab="param",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelFormula,toplot="paramX",param=c(i,2),add=TRUE,col="blue")
	mixing(modelFormula,toplot="paramX",param=c(i,3),add=TRUE,col="red")
	abline(h=dataBase$param$param[i,],col=c("black","blue","red"),lwd=3)
}

#------------------------------------------
### Check the distribution of the posterior
#------------------------------------------
#__________________
### Community means
#__________________
pairs(modelFormula$model$means,asp=1,pch=19,cex=0.2)

########################
### Example using Traits
########################
#==============================
### Simulating data with traits
#==============================
desc <- cbind(1,scale(1:50),scale(1:50)^2)
trait <- rbind(1,1:20/5)
dataTrait<-communitySimulHT(desc,trait)

#=================================
### Formatting data and parameters
#=================================
formdata<-as.HMSCdata(dataTrait$data$Y,desc,trait,Ypattern="sp",interceptX=FALSE,interceptTr=FALSE)

#========================
### Formatting parameters
#========================
startParamX<-matrix(rnorm(nrow(trait)*ncol(desc)),nrow=ncol(dataTrait$data$Y),ncol=ncol(desc))
startParamTr<-matrix(rnorm(ncol(trait)*ncol(desc)),nrow=ncol(desc),ncol=nrow(trait))
startSigma<-rWishart(1,ncol(desc)+1,diag(ncol(desc)))[,,1]

formparam<-as.HMSCparam(formdata,paramX=startParamX,paramTr=startParamTr,sigma=startSigma)
formpriors<-as.HMSCprior(formparam,kappa0=0,nu0=length(formparam$means),Lambda0=diag(nrow(formparam$sigma)))

#=============================================================================================================
### Building model using a linear modelling approach (currently 2000 iterations takes about 50 seconds to run)
#=============================================================================================================
system.time(
modelTraits<-hmscHT(formdata,formparam,formpriors,niter=200,nburn=100)
)

### Traits model
summaryTraits<-summary(modelTraits)
summaryTraits[1:5,1:5,]

#------------------
### Mixing - Traits
#------------------
#___________________
#### Community Sigma
#___________________
par(mfrow=c(3,3),mar=c(1,3,3,1))
for(i in 1:3){
	for(j in 1:3){
		mixing(modelTraits,toplot="sigma",param=c(i,j),col="red",xlab="Sites",ylab="Sigma",ylim=range(c(modelTraits$model$sigma[i,j,],dataTrait$param$sigma[i,j])),las=1)
		abline(h=dataTrait$param$sigma[i,j],col="red",lwd=3)
		title(paste("Community Sigma",i,"-",j))
	}
}

#__________________
### Community means
#__________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelTraits,toplot="means",param=c(1,i),ylim=range(cbind(modelTraits$model$means[,i,],dataTrait$param$means[,i])),xlab="Sites",ylab="means",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelTraits,toplot="means",param=c(2,i),add=TRUE,col="blue")
	mixing(modelTraits,toplot="means",param=c(3,i),add=TRUE,col="red")
	abline(h=dataTrait$param$means[,i],col=c("black","blue","red"),lwd=3)
}

#_____________________________
### Parameter for each species 
#_____________________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelTraits,toplot="paramX",param=c(i,1),ylim=range(cbind(modelTraits$model$paramX[i,,],dataTrait$param$paramX[i,])),xlab="Sites",ylab="paramX",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelTraits,toplot="paramX",param=c(i,2),add=TRUE,col="blue")
	mixing(modelTraits,toplot="paramX",param=c(i,3),add=TRUE,col="red")
	abline(h=dataTrait$param$paramX[i,],col=c("black","blue","red"),lwd=3)
}

#____________________________
### Parameter for each traits 
#____________________________
par(mfrow=c(3,1),mar=c(1,4,2,1))
for(i in 1:3){
	mixing(modelTraits,toplot="paramTr",param=c(i,1),ylim=range(cbind(modelTraits$model$paramTr[i,,],dataTrait$param$paramTr[i,])),xlab="Sites",ylab="paramTr",las=1,xaxt="n")
	title(paste("Habitat",i))
	mixing(modelTraits,toplot="paramTr",param=c(i,2),add=TRUE,col="blue")
	abline(h=dataTrait$param$paramTr[i,],col=c("black","blue"),lwd=3)
}

#------------------------------
### Draw model from modelTraits
#------------------------------
par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	plot(summaryTraits[,i,1],type="l",ylim=c(0,1),las=1,lwd=3,xaxt="n",col="red")
	title(paste("Species",i))
	lines(summaryTraits[,i,3],col="red",lwd=1)
	lines(summaryTraits[,i,5],col="red",lwd=1)
	lines(dataTrait$probMat[,i],col="black",lwd=3)
}

#################################
### Example using a Random effect
#################################
#==============================
### Simulating data with traits
#==============================
desc <- cbind(1,scale(1:50),scale(1:50)^2)
arbi<-as.factor(c(1,rep(1:3,each=16),3))
dataRandom<-communitySimulH(desc,Random=arbi,nsp=20)

#=================================
### Formatting data and parameters
#=================================
formdata<-as.HMSCdata(dataRandom$data$Y,desc,Random=arbi,Ypattern="sp",interceptX=FALSE,scaleX=FALSE)

#========================
### Formatting parameters
#========================
startParamX<-matrix(rnorm(ncol(dataRandom$data$Y)*ncol(desc)),nrow=ncol(dataRandom$data$Y),ncol=ncol(desc))
startSigma<-rWishart(1,ncol(desc)+1,diag(ncol(desc)))[,,1]
startRandom<-matrix(rnorm(ncol(dataRandom$data$Y)*nlevels(arbi)),nrow=ncol(dataRandom$data$Y),ncol=nlevels(arbi))
startMean<-rnorm(ncol(desc))

formparam<-as.HMSCparam(formdata,paramX=startParamX,paramRandom=startRandom,means=startMean,sigma=startSigma,RandomVar=dataRandom$param$RandomCommSp+1,RandomCommSp=dataRandom$param$RandomCommSp)

#============================================================================
### Building model (currently 2000 iterations takes about 2.5 minutes to run)
#============================================================================
system.time(
modelRandom<-hmscH(formdata,formparam,niter=200,nburn=100)
)

### Traits model
summaryRandom<-summary(modelRandom)
summaryRandom[1:5,1:5,]

#------------------
### Mixing - Random
#------------------
#___________________
#### Community Sigma
#___________________
par(mfrow=c(3,3),mar=c(1,3,3,1))
for(i in 1:3){
	for(j in 1:3){
		mixing(modelRandom,toplot="sigma",param=c(i,j),col="red",xlab="Sites",ylab="Sigma",ylim=range(c(modelRandom$model$sigma[i,j,],dataRandom$param$sigma[i,j])),las=1)
		abline(h=dataRandom$param$sigma[i,j],col="red",lwd=3)
		title(paste("Community Sigma",i,"-",j))
	}
}

#__________________
### Community means
#__________________
par(mfrow=c(3,1),mar=c(1,4,2,1))
for(i in 1:3){
	mixing(modelRandom,toplot="means",param=i,col="red",ylim=range(cbind(modelRandom$model$means[,i],dataRandom$param$means[i])),xlab="Sites",ylab="means",las=1,xaxt="n")
	abline(h=dataRandom$param$means[i],col="red",lwd=3)
}

#__________________________________
### Parameter of X for each species 
#__________________________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelRandom,toplot="paramX",param=c(i,1),ylim=range(cbind(modelRandom$model$paramX[i,,],dataRandom$param$paramX[i,])),xlab="Sites",ylab="paramX",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelRandom,toplot="paramX",param=c(i,2),add=TRUE,col="blue")
	mixing(modelRandom,toplot="paramX",param=c(i,3),add=TRUE,col="red")
	abline(h=dataRandom$param$paramX[i,],col=c("black","blue","red"),lwd=3)
}

#_________________________________________________
### Parameter for each levels of the Random effect
#_________________________________________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelRandom,toplot="paramRandom",param=c(i,1),ylim=range(cbind(modelRandom$model$paramRandom[i,,],dataRandom$param$paramRandom[i,])),xlab="Sites",ylab="paramRandom",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelRandom,toplot="paramRandom",param=c(i,2),add=TRUE,col="blue")
	mixing(modelRandom,toplot="paramRandom",param=c(i,3),add=TRUE,col="red")
	abline(h=dataRandom$param$paramRandom[i,],col=c("black","blue","red"),lwd=3)
}

#________________________________
### Variance of the Random effect
#________________________________
par(mfrow=c(1,1),mar=c(1,4,2,1))
mixing(modelRandom,toplot="RandomVar",param=i,col="red",ylim=range(c(modelRandom$model$RandomVar,dataRandom$param$RandomVar)),xlab="Sites",ylab="RandomVar",las=1,xaxt="n")
abline(h=dataRandom$param$RandomVar,col="red",lwd=3)

#------------------------------
### Draw model from modelRandom
#------------------------------
par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	plot(summaryRandom[,i,1],type="l",ylim=c(0,1),las=1,lwd=3,xaxt="n",col="red")
	title(paste("Species",i))
	lines(summaryRandom[,i,3],col="red",lwd=1)
	lines(summaryRandom[,i,5],col="red",lwd=1)
	lines(dataRandom$probMat[,i],col="black",lwd=3)
}

##########################################
### Example using Traits and Random effect
##########################################
#==============================
### Simulating data with traits
#==============================
desc <- cbind(1,scale(1:50),scale(1:50)^2)
arbi<-as.factor(c(1,rep(1:3,each=16),3))
trait <- rbind(1,1:20/5)
dataTraitsRandom<-communitySimulHT(desc,trait,Random=arbi)

#=================================
### Formatting data and parameters
#=================================
formdata<-as.HMSCdata(dataTraitsRandom$data$Y,desc,trait,Random=arbi,Ypattern="sp",interceptX=FALSE,interceptTr=FALSE)

#========================
### Formatting parameters
#========================
startParamX<-matrix(rnorm(nrow(trait)*ncol(desc)),nrow=ncol(dataTraitsRandom$data$Y),ncol=ncol(desc))
startParamTr<-matrix(rnorm(ncol(trait)*ncol(desc)),nrow=ncol(desc),ncol=nrow(trait))
startSigma<-rWishart(1,ncol(desc)+1,diag(ncol(desc)))[,,1]
startRandom<-matrix(rnorm(ncol(dataTraitsRandom$data$Y)*nlevels(arbi)),nrow=ncol(dataTraitsRandom$data$Y),ncol=nlevels(arbi))

formparam<-as.HMSCparam(formdata,paramX=startParamX,paramRandom=startRandom,paramTr=startParamTr,sigma=startSigma,RandomVar=dataTraitsRandom$param$RandomVar,RandomCommSp=dataTraitsRandom$param$RandomCommSp)

#=============================================================================================================
### Building model using a linear modelling approach (currently 2000 iterations takes about 50 seconds to run)
#=============================================================================================================
system.time(
modelTraitsRandom<-hmscHT(formdata,formparam,niter=200,nburn=100)
)

### Traits model
summaryTraitsRandom<-summary(modelTraitsRandom)
summaryTraitsRandom[1:5,1:5,]

#-------------------------
### Mixing - Traits Random
#-------------------------
#___________________
#### Community Sigma
#___________________
par(mfrow=c(3,3),mar=c(1,3,3,1))
for(i in 1:3){
	for(j in 1:3){
		mixing(modelTraitsRandom,toplot="sigma",param=c(i,j),col="red",xlab="Sites",ylab="Sigma",ylim=range(c(modelTraitsRandom$model$sigma[i,j,],dataTraitsRandom$param$sigma[i,j])),las=1)
		abline(h=dataTraitsRandom$param$sigma[i,j],col="red",lwd=3)
		title(paste("Community Sigma",i,"-",j))
	}
}

#__________________
### Community means
#__________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelTraitsRandom,toplot="means",param=c(1,i),ylim=range(cbind(modelTraitsRandom$model$means[,i,],dataTraitsRandom$param$means[,i])),xlab="Sites",ylab="means",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelTraitsRandom,toplot="means",param=c(2,i),add=TRUE,col="blue")
	mixing(modelTraitsRandom,toplot="means",param=c(3,i),add=TRUE,col="red")
	abline(h=dataTraitsRandom$param$means[,i],col=c("black","blue","red"),lwd=3)
}

#__________________________________
### Parameter of X for each species 
#__________________________________
par(mfrow=c(5,4),mar=c(1,4,2,1))
for(i in 1:20){
	mixing(modelTraitsRandom,toplot="paramX",param=c(i,1),ylim=range(cbind(modelTraitsRandom$model$paramX[i,,],dataTraitsRandom$param$paramX[i,])),xlab="Sites",ylab="paramX",las=1,xaxt="n")
	title(paste("Species",i))
	mixing(modelTraitsRandom,toplot="paramX",param=c(i,2),add=TRUE,col="blue")
	mixing(modelTraitsRandom,toplot="paramX",param=c(i,3),add=TRUE,col="red")
	abline(h=dataTraitsRandom$param$paramX[i,],col=c("black","blue","red"),lwd=3)
}

#____________________________
### Parameter for each traits 
#____________________________
par(mfrow=c(3,1),mar=c(1,4,2,1))
for(i in 1:3){
	mixing(modelTraitsRandom,toplot="paramTr",param=c(i,1),ylim=range(cbind(modelTraitsRandom$model$paramTr[i,,],dataTraitsRandom$param$paramTr[i,])),xlab="Sites",ylab="paramTr",las=1,xaxt="n")
	title(paste("Habitat",i))
	mixing(modelTraitsRandom,toplot="paramTr",param=c(i,2),add=TRUE,col="blue")
	abline(h=dataTraitsRandom$param$paramTr[i,],col=c("black","blue"),lwd=3)
}

#------------------------------------
### Draw model from modelTraitsRandom
#------------------------------------
par(mfrow=c(5,4),mar=c(1,3,2,1))
for(i in 1:20){
	plot(summaryTraitsRandom[,i,1],type="l",ylim=c(0,1),las=1,lwd=3,xaxt="n",col="red")
	title(paste("Species",i))
	lines(summaryTraitsRandom[,i,3],col="red",lwd=1)
	lines(summaryTraitsRandom[,i,5],col="red",lwd=1)
	lines(dataTraitsRandom$probMat[,i],col="black",lwd=3)
}

HMSC documentation built on May 2, 2019, 6:53 p.m.