R/TrueTrait.R

In WGCNA: Weighted Correlation Network Analysis

TrueTrait=function(datX, y,datXtest=NULL, corFnc = "bicor", corOptions = "use = 'pairwise.complete.obs'", LeaveOneOut.CV=FALSE,skipMissingVariables=TRUE,addLinearModel=FALSE){
datX=as.matrix(datX)
no.variables=dim(as.matrix(datX))[[2]]
datVariableInfo=data.frame(matrix(NA, nrow=no.variables, ncol=4))
names(datVariableInfo)=c("Variable","center", "scale", "weights.y.true2")
if ( is.null(colnames(datX))){datVariableInfo$Variable=1:no.variables} else { datVariableInfo$Variable=colnames(datX)}
no.observations=dim(as.matrix(datX))[[1]]
if (no.observations !=  length(y) ) {stop("The number of rows of datX does not correspond to the length of y. Consider transposing your input matrix or use a different input.") }
if (no.observations==1 ) {
warning("Only 1 observations, i.e. the length of y is 1. The function cannot be used. For your convenience, the estimated true values will be set to the input value of y.")
y.true1=y; y.true2=y;  y.true3=y}
if (no.observations>1 ) {
y.true1=rep(NA, length(y) )
y.true2=rep(NA, length(y) )
y.true3=rep(NA, length(y) )
y.lm=rep(NA, length(y) )
restNonMissingY= !is.na(y)
r.characteristic=NA
SD.ytrue2=NA
SD.ytrue3=NA
SsquaredBE=NA

if (sum(restNonMissingY,na.rm=TRUE) >3 ){
corX = parse(text = paste(corFnc, "(datX,y ",prepComma(corOptions), ")"))
rVector= as.numeric(eval(corX))
datCoef=t(coef(lm(datX~ y,na.action="na.exclude")))
datVariableInfo$center=datCoef[,1]	#intercept
datVariableInfo$scale=datCoef[,2]	#slope
datXscaled=scale(datX,center=datCoef[,1],scale=datCoef[,2] )
weights0=rVector^2/((1-rVector^2)*var(y,na.rm=TRUE)) # Steve, this is where I made the one change
weights=weights0/sum(weights0)
datVariableInfo$weights.y.true2=weights
y.true1=as.numeric(apply(as.matrix(datXscaled),1,mean))
y.true2=as.numeric(as.matrix(datXscaled)%*%weights)

if (skipMissingVariables ) {
y.true1= as.numeric(apply(as.matrix(datXscaled),1,mean,na.rm=TRUE))
weightsMatrix=matrix(weights,byrow=TRUE,nrow=dim(as.matrix(datXscaled))[[1]],ncol=length(weights) )
weightsMatrix[is.na(datXscaled)]=0
rowsum.weightsMatrix=apply(as.matrix(weightsMatrix),1,sum)
weightsMatrix=t(scale(t(as.matrix(weightsMatrix)),center=F,scale= rowsum.weightsMatrix))
datXscaledweighted= as.matrix(datXscaled* weightsMatrix)
# this corresponds to formula 25 in Klemera et al 2006
y.true2=as.numeric(apply(datXscaledweighted,1,sum,na.rm=TRUE))
 } #end of if (skipMissingVariables )


# the following is different from Klemera in that it has an absolute value
r.characteristic=sum(rVector^2/sqrt(1-rVector^2) )/sum(abs(rVector)/sqrt(1-rVector^2) )
no.missing=sum(apply(  as.matrix( is.na(datX)),1,sum))
if (sum(no.missing)>0) {warning("The input datX contains missing values.\n
I recommend you impute missing values in datX before running this function. ") 
} # end of if (sum(no.missing)>0)
# formula 37 from Klemera
SsquaredBE=var( y.true2-y,na.rm=TRUE) -(1- r.characteristic^2)/r.characteristic^2*var(y,na.rm=TRUE)/no.variables
# this corresponds to formula 34 in Klemera
y.true3=(as.numeric( as.matrix(datXscaled)%*% weights0)+y/SsquaredBE )/( sum(weights0)+ 1/SsquaredBE)
y.true3[is.na(y.true3) ]=y[is.na(y.true3)]
} # end of if (no.observations>1 ) 
SD.ytrue2=sqrt(1-r.characteristic^2)/r.characteristic*sqrt(var(y,na.rm=TRUE)/no.variables)
# now formula 42
SD.ytrue3=SD.ytrue2/sqrt(1+SD.ytrue2^2/SsquaredBE )
} # end of if (sum(restNonMissingY,na.rm=TRUE) >3 )
datEstimates=data.frame(y, y.true1,y.true2,y.true3)

if (!is.null(datXtest)){
datXtest=as.matrix(datXtest)
no.variablestest=dim(as.matrix(datXtest))[[2]]
if (no.variablestest != no.variables) {stop("the number of variables in the test data is not the same as in the training data")}

y.true1test=rep(NA, length(y) )
y.true2test=rep(NA, length(y) )
y.true3test=rep(NA, length(y) )
restNonMissingY= !is.na(y)
if (sum(restNonMissingY,na.rm=TRUE) >3 ){
datXtestscaled=scale(datXtest,center=datCoef[,1],scale=datCoef[,2] )
y.true1test=as.numeric(apply(as.matrix(datXtestscaled),1,mean) )
y.true2test=as.numeric( as.matrix(datXtestscaled)%*%weights)
if (skipMissingVariables ) {
y.true1test= as.numeric(apply(as.matrix(datXtestscaled),1,mean,na.rm=TRUE))
weightsMatrixtest=matrix(weights,byrow=TRUE,nrow=dim(as.matrix(datXtestscaled))[[1]],ncol=length(weights) )
weightsMatrixtest[is.na(datXtestscaled)]=0
rowsum.weightsMatrixtest=apply(as.matrix(weightsMatrixtest),1,sum)
weightsMatrixtest=t(scale(t(as.matrix(weightsMatrixtest)),center=F,scale= rowsum.weightsMatrixtest))
datXscaledweightedtest= as.matrix(datXtestscaled* weightsMatrixtest)
# this corresponds to formula 25 in Klemera et al 2006
y.true2test=as.numeric(apply(datXscaledweightedtest,1,sum,na.rm=TRUE))
 } #end of if (skipMissingVariables )
} # end of if (sum(restNonMissingY,na.rm=TRUE) >3 )


datEstimatestest=data.frame(y.true1= y.true1test,y.true2= y.true2test)

} # end of if (!is.null(datXtest))

if ( LeaveOneOut.CV  ) {
y.true1test.LOO=rep(NA,no.observations)
y.true2test.LOO=rep(NA,no.observations)
y.lmLOO= rep(NA,no.observations)
for ( i in 1:no.observations ){
rm(datCoef); rm(corX);

datX.LOO=datX[-i,]
datXtest.LOO= matrix(datX[i,],nrow=1)
y.LOO=y[-i]
no.variables=dim(as.matrix(datX.LOO))[[2]]
no.observations=dim(as.matrix(datX.LOO))[[1]]

if (no.observations==1 ) {
warning("When dealing with leave one out cross validation, there is only 1 observations in the training data")}

if (no.observations>1 ) {

if (addLinearModel) {
lmLOO=lm(y.LOO~., data=data.frame(datX.LOO),na.action=na.exclude)
y.lmLOO[i]= sum(datXtest.LOO*lmLOO$coeff[-1])+lmLOO$coeff[[1]]
}

corX = parse(text = paste(corFnc, "(datX.LOO,y.LOO ", 
                prepComma(corOptions), ")"))
rVector= as.numeric(eval(corX))
datCoef=t(coef(lm(datX.LOO~ y.LOO,na.action="na.exclude")))
datX.LOOscaled=scale(datX.LOO,center=datCoef[,1],scale=datCoef[,2] )
weights0=rVector^2/(1-rVector^2)
weights=weights0/sum(weights0)
datXtest.LOOscaled=(datXtest.LOO-datCoef[,1])/datCoef[,2] 
y.true1test.LOO[i]= mean(datXtest.LOOscaled)
y.true2test.LOO[i]=sum(datXtest.LOOscaled*weights)
if (skipMissingVariables ) {
y.true1test.LOO[i]= mean(datXtest.LOOscaled,na.rm=TRUE)
weightsMatrixLOO= weights
weightsMatrixLOO[is.na(datXtest.LOOscaled)]=0
rowsum.weightsMatrixLOO=sum(weightsMatrixLOO)
weightsMatrixLOO=weightsMatrixLOO/rowsum.weightsMatrixLOO
datXscaledweightedLOO= datXtest.LOOscaled* weightsMatrixLOO
y.true2test.LOO[i]=sum(datXscaledweightedLOO,na.rm=TRUE)
 } #end of if (skipMissingVariables )
} # end of for loop
} # end of if (no.observations>1 ) 
datEstimates.LeaveOneOut.CV=data.frame(y.true1= y.true1test.LOO,
y.true2= y.true2test.LOO)
} # end of if ( LeaveOneOut.CV  ) 

if (addLinearModel) {
y.lmTest=rep(NA, dim(as.matrix(datXtest))[[1]] )
restNonMissingY= !is.na(y)
if (sum(restNonMissingY,na.rm=TRUE) >3 ){
lm1=lm(y~., data=data.frame(datX),na.action=na.exclude)
y.lmTraining=predict(lm1)
y.lmTraining=predict(lm1)
if( !is.null(datXtest)) {y.lmTest=predict(lm1,newdata=data.frame(datXtest))}
}
}
if ( !is.null(datXtest) & LeaveOneOut.CV & !addLinearModel ) {out=  list( datEstimates=datEstimates, datEstimatestest= datEstimatestest,
datEstimates.LeaveOneOut.CV= datEstimates.LeaveOneOut.CV, SD.ytrue2=SD.ytrue2, SD.ytrue3=SD.ytrue3, datVariableInfo= datVariableInfo ) }
if ( !is.null(datXtest) & !LeaveOneOut.CV & !addLinearModel) {out=  list( datEstimates=datEstimates, datEstimatestest= datEstimatestest, SD.ytrue2=SD.ytrue2, SD.ytrue3=SD.ytrue3, datVariableInfo= datVariableInfo ) }
if ( is.null(datXtest) & LeaveOneOut.CV & !addLinearModel) {out=  list( datEstimates=datEstimates, datEstimates.LeaveOneOut.CV= datEstimates.LeaveOneOut.CV, SD.ytrue2=SD.ytrue2, SD.ytrue3=SD.ytrue3, datVariableInfo= datVariableInfo ) }
if ( is.null(datXtest) & !LeaveOneOut.CV & !addLinearModel) {out=  list( datEstimates=datEstimates, SD.ytrue2=SD.ytrue2, SD.ytrue3=SD.ytrue3, datVariableInfo= datVariableInfo ) }

if ( !is.null(datXtest) & LeaveOneOut.CV & addLinearModel ) {out=  list( datEstimates=data.frame(datEstimates, y.lm=y.lmTraining), datEstimatestest= data.frame(datEstimatestest, y.lm=y.lmTest),
datEstimates.LeaveOneOut.CV= data.frame(datEstimates.LeaveOneOut.CV,y.lm=y.lmLOO), SD.ytrue2=SD.ytrue2, SD.ytrue3=SD.ytrue3, datVariableInfo= datVariableInfo ) }

if ( !is.null(datXtest) & !LeaveOneOut.CV & addLinearModel) {out=  list( datEstimates=data.frame(datEstimates,y.lm=y.lmTraining), datEstimatestest= data.frame(datEstimatestest,y.lm=y.lmTest), SD.ytrue2=SD.ytrue2, SD.ytrue3=SD.ytrue3, datVariableInfo= datVariableInfo ) }

if ( is.null(datXtest) & LeaveOneOut.CV & addLinearModel) {out=  list( datEstimates=data.frame(datEstimates,y.lm=y.lmTraining), datEstimates.LeaveOneOut.CV= data.frame(datEstimates.LeaveOneOut.CV,y.lm=y.lmLOO),
 SD.ytrue2=SD.ytrue2, SD.ytrue3=SD.ytrue3, datVariableInfo= datVariableInfo ) }
if ( is.null(datXtest) & !LeaveOneOut.CV & addLinearModel) {out=  list( datEstimates=data.frame(datEstimates,y.lm=y.lmTraining), SD.ytrue2=SD.ytrue2, SD.ytrue3=SD.ytrue3, datVariableInfo= datVariableInfo ) }
out
} # end of function