R/as.HMSCdata.R

In HMSC: Hierarchical modelling of species community

as.HMSCdata <-
function(Y,X,Tr=NULL,Random=NULL,Ypattern=NULL,scaleX=TRUE,scaleTr=TRUE,interceptX=TRUE,interceptTr=TRUE){
#### F. Guillaume Blanchet - March 2013, April 2013, July 2013
##########################################################################################
	Ynames<-colnames(Y)
	if(is.null(Ypattern)){
		Ypattern<-"y"
		if(!is.null(Ynames)){
			print("'Ypattern' was defined as 'y'")
		}
	}
	
	### Check for NAs
	if(any(is.na(Y))){
		stop("There is at least one NA in 'Y'")
	}
	if(any(is.na(X))){
		stop("There is at least one NA in 'X'")
	}
	if(!is.null(Tr)){
		if(any(is.na(Tr))){
			stop("There is at least one NA in 'Tr'")
		}
	}
	if(!is.null(Random)){
		if(any(is.na(Random))){
			stop("There is at least one NA in 'Random'")
		}
	}
	
	#### Check format
	if(length(dim(Y))!=2){
		stop("'Y' shoulds be a table")
	}
	if(length(dim(X))!=2){
		stop("'X' shoulds be a table")
	}
	if(!is.null(Tr)){
		if(length(dim(Tr))!=2){
			stop("'Tr' shoulds be a table")
		}
	}
	if(!is.null(Random)){
		if(!is.factor(Random)){
			stop("'Random' should be a factor")
		}
	}
	
	#### Check if dimensions of all tables match
	if(nrow(Y)!=nrow(X)){
		stop("'X' and 'Y' should have the same number of rows")
	}
	
	if(!is.null(Tr)){
		if(ncol(Y)!=ncol(Tr)){
			stop("'Y' and 'Tr' should have the same number of columns")
		}
	}
	if(!is.null(Random)){
		if(length(Random)!=nrow(Y)){
			stop("'Random' should have a length equal to the number of rows of 'Y'")
		}
	}
	
	### Transform Random to a 0-1 table
	if(!is.null(Random)){
		RandomMat<-matrix(0, length(Random), nlevels(Random))
		RandomMat[(1:length(Random))+length(Random) * (unclass(Random) - 1)]<-1
		colnames(RandomMat)<-levels(Random)
	}
	
	#### Check column names
	if(is.null(Ynames)){
		colnames(Y)<-paste(Ypattern,1:ncol(Y),sep="")
		print(paste("column names were added to 'Y' using",Ypattern,"as Ypattern"))
	}
	
	if(is.null(colnames(X))){
		colnames(X)<-paste("x",1:ncol(X),sep="")
		print("column names were added to 'X'")
	}
	
	if(!is.null(Tr)){
		if(is.null(colnames(Tr))){
			colnames(Tr)<-paste(Ypattern,1:ncol(Tr),sep="")
			print("column names were added to 'Tr'")
		}
	}
	
	#### Check row names
	if(is.null(rownames(Y))){
		rownames(Y)<-paste("site",1:nrow(Y),sep="")
		print(paste("row names were added to 'Y'"))
	}
	
	if(is.null(rownames(X))){
		rownames(X)<-paste("site",1:nrow(X),sep="")
		print("row names were added to 'X'")
	}

	if(!is.null(Tr)){
		if(is.null(rownames(Tr))){
			rownames(Tr)<-paste("t",1:nrow(Tr),sep="")
			print("row names were added to 'Tr'")
		}
	}
	if(!is.null(Random)){
		rownames(RandomMat)<-paste("site",1:nrow(RandomMat),sep="")
	}
	
	#### Check regular pattern in names of Y
	nY<-ncol(Y)
	
	if(any(unique(regexpr(Ypattern,Ynames))!=1)){
		colnames(Y)<-paste(Ypattern,colnames(Y),sep="")
		print("The column names of 'Y' were modified, they now all start with 'Ypattern'")
	}
		
	### Add an intercept to X and scale
	if(interceptX){
		if(scaleX){
			### Check if any columns have a null variance
			zeroVar<-which(apply(X,2,sd)==0)
			if(length(zeroVar)!=0){
				X[,-zeroVar]<-scale(X[,-zeroVar])
				warning(paste(colnames(X)[zeroVar],"are explanatory variable(s) with a variance of 0, for this reason no intercept were added, check to make sure this is OK"))
			}else{
				X<-cbind(1,scale(X))
				colnames(X)[1]<-"Intercept"
			}
		}else{
			X<-cbind(1,X)
			colnames(X)[1]<-"Intercept"
		}
	}else{
		if(scaleX){
			### Check if any columns have a null variance
			zeroVar<-which(apply(X,2,sd)==0)
			if(length(zeroVar)!=0){
				X[,-zeroVar]<-scale(X[,-zeroVar])
				warning(paste(colnames(X)[zeroVar],"are explanatory variable(s) with a variance of 0, make sure this is OK"))
			}else{
				X<-scale(X)
			}
		}
	}
	
	### Add an intercept to Tr
	if(!is.null(Tr)){
		if(interceptTr){
			if(scaleTr){
				### Check if any columns have a null variance
				zeroVar<-which(apply(Tr,1,sd)==0)
				if(length(zeroVar)!=0){
					Tr[-zeroVar,]<-scale(Tr[-zeroVar,])
					warning(paste(rownames(Tr)[zeroVar],"are trait(s) with a variance of 0, for this reason no intercept were added, check to make sure this is OK"))
				}else{
					Tr<-rbind(1,t(scale(t(Tr))))
					rownames(Tr)[1]<-"Intercept"
				}
			}else{
				Tr<-rbind(1,Tr)
				rownames(Tr)[1]<-"Intercept"
			}
		}else{
			if(scaleTr){
				### Check if any columns have a null variance
				zeroVar<-which(apply(Tr,1,sd)==0)
				if(length(zeroVar)!=0){
					Tr[-zeroVar,]<-scale(Tr[-zeroVar,])
					warning(paste(rownames(Tr)[zeroVar],"are traits with a variance of 0, make sure this is OK"))
				}
			}
		}
	}
	
	#### Check classes
	if(!is.data.frame(Y)){
		Y<-as.data.frame(Y)
		print("'Y' was converted to a data.frame")
	}
	
	if(!is.data.frame(X)){
		X<-as.data.frame(X)
		print("'X' was converted to a data.frame")
	}

	if(!is.null(Tr)){
		if(!is.data.frame(Tr)){
			Tr<-as.data.frame(Tr)
			print("'Tr' was converted to a data.frame")
		}
	}
	if(!is.null(Random)){
		RandomMat<-as.data.frame(RandomMat)
	}
	
	#### Return results
	if(is.null(Tr) & is.null(Random)){
		res<-list(Y=Y,X=X)
		attributes(res)<-list(names=c("Y","X"),Ypattern=Ypattern)
		class(res)<-"HMSCdata"
	}
	
	if(is.null(Random) & !is.null(Tr)){
		res<-list(Y=Y,X=X,Tr=Tr)
		attributes(res)<-list(names=c("Y","X","Tr"),Ypattern=Ypattern)
		class(res)<-c("HMSCdataTrait","HMSCdata")
	}
	
	if(!is.null(Random) & is.null(Tr)){
		res<-list(Y=Y,X=X,Random=RandomMat)
		attributes(res)<-list(names=c("Y","X","Random"),Ypattern=Ypattern)
		class(res)<-c("HMSCdataRandom","HMSCdata")
	}

	if(!is.null(Tr) & !is.null(Random)){
		res<-list(Y=Y,X=X,Tr=Tr,Random=RandomMat)
		attributes(res)<-list(names=c("Y","X","Tr","Random"),Ypattern=Ypattern)
		class(res)<-c("HMSCdataTrait","HMSCdataRandom","HMSCdata")
	}
	
	return(res)
}