#================================================== #================================================== #================================================== # Libraries library( cvTools ) library( lars ) library( gregmisc ) library( hierNet ) library( doMC ) registerDoMC( 8 ) #================================================== #================================================== #================================================== # Functions #================================================== # Simple linear model function (applying cross-validation similar to Edwards et al.) linear.model.crossvalidation <- function( data, repeats = 1000 ) { #================================================== # For each fold (1000 cases) split the data into a training- and testset (1/2) (randomly) n <- nrow( data ) folds <- cvFolds( n, K = 2, R = repeats ) #================================================== # Loop through each fold and perform analysis R2s <- NULL MSEs <- NULL coefficients <- NULL for( i in 1:repeats ) { #================================================== # Get training and test data from folds if( ( n %% 2 ) == 0 ) { training.rows <- folds$subsets[ seq( 1, n - 1, 2 ), i ] test.rows <- folds$subsets[ seq( 2, n, 2 ), i ] } else { training.rows <- folds$subsets[ seq( 1, n, 2 ), i ] test.rows <- folds$subsets[ seq( 2, n - 1, 2 ), i ] } #================================================== # Training and test data training.set <- data[ training.rows, ] test.set <- data[ test.rows, ] #================================================== # Perform linear modelling with training data for 3 different types of linear models # 1. Simple model 'Y ~ X1' formula <- 'Y ~ X1' training.set.simple <- training.set[ , c( 1, 4 ) ] test.set.simple <- test.set[ , c( 1, 4 ) ] # Perform linear modelling training.lm.simple <- lm( formula, training.set ) R2.simple <- round( summary( training.lm.simple )$adj.r.squared, 4 ) # Test/Predict linear model with test data test.predict <- as.vector( predict( training.lm.simple, test.set.simple, interval = 'prediction' )[ , 1 ] ) n.test <- nrow( test.set.simple ) MSE.simple <- sum( ( test.predict - test.set.simple$Y )^2 ) / n.test # 2. Additive model 'Y ~ X1 + X2' formula <- 'Y ~ X1 + X2' training.set.additive <- training.set[ , c( 1:2, 4 ) ] test.set.additive <- test.set[ , c( 1:2, 4 ) ] # Perform linear modelling training.lm.additive <- lm( formula, training.set.additive ) R2.additive <- round( summary( training.lm.additive )$adj.r.squared, 4 ) # Test/Predict linear model with test data test.predict <- as.vector( predict( training.lm.additive, test.set.additive, interval = 'prediction' )[ , 1 ] ) n.test <- nrow( test.set.additive ) MSE.additive <- sum( ( test.predict - test.set.additive$Y )^2 ) / n.test # 3. interaction model 'Y ~ X1 + X2 + X1X2' formula <- 'Y ~ X1 + X2 + X1X2' training.set.interaction <- training.set test.set.interaction <- test.set # Perform linear modelling training.lm.interaction <- lm( formula, training.set.interaction ) R2.interaction <- round( summary( training.lm.interaction )$adj.r.squared, 4 ) # Test/Predict linear model with test data test.predict <- as.vector( predict( training.lm.interaction, test.set.interaction, interval = 'prediction' )[ , 1 ] ) n.test <- nrow( test.set.interaction ) MSE.interaction <- sum( ( test.predict - test.set.interaction$Y )^2 ) / n.test #================================================== # Save R-squared and MSEs R2s <- rbind( R2s, c( R2.simple, R2.additive, R2.interaction ) ) MSEs <- rbind( MSEs, c( MSE.simple, MSE.additive, MSE.interaction ) ) } #================================================== # Output return( list( R2s = R2s, MSEs = MSEs ) ) } #================================================== # LASSO and BIC model crossvalidation optimisation LASSO.BIC.model.crossvalidation <- function( data, repeats = 1000 ) { #================================================== # For each fold (1000 cases) split the data into a training- and testset (1/2) (randomly) n <- nrow( data ) folds <- cvFolds( n, K = 2, R = repeats ) #================================================== # Loop through each fold and perform analysis R2s <- NULL predictors <- NULL min.BIC.steps <- NULL for( i in 1:repeats ) { #================================================== # Get training and test data from folds if( ( n %% 2 ) == 0 ) { training.rows <- folds$subsets[ seq( 1, n - 1, 2 ), i ] test.rows <- folds$subsets[ seq( 2, n, 2 ), i ] } else { training.rows <- folds$subsets[ seq( 1, n, 2 ), i ] test.rows <- folds$subsets[ seq( 2, n - 1, 2 ), i ] } #================================================== # Training and test data training.set <- data[ training.rows, ] test.set <- data[ test.rows, ] #================================================== # Apply improved 'lasso.adapt.bic2' function with training data lasso.bic <- LASSO.BIC( training.set, test.set ) #================================================== # Save output values R2s <- c( R2s, lasso.bic$R2 ) predictors <- rbind( predictors, lasso.bic$predictor.occurance ) min.BIC.steps <- c( min.BIC.steps, lasso.bic$min.bic2.step ) } return( list( R2s = R2s, predictors = predictors, min.BIC.steps = min.BIC.steps ) ) } #================================================== # Single LASSO and BIC optimisation LASSO.BIC <- function( training.set, test.set ) { #================================================== # Prepare data x <- training.set[ , 1:3 ] y <- training.set$Y ok<-complete.cases(x,y) #================================================== # This if statement insures that if only one predictor is available the function does not crash if (length( dim( x )) == 0 ) { x <- x[ ok ] y <- y[ ok ] m <- 1 n <- length( x ) } else { x <- x[ ok, ] # get rid of na's y <- y[ ok ] # since regsubsets can't handle na's m <- ncol( x ) n <- nrow( x ) } x <- as.matrix( x ) # in case x is not a matrix #================================================== # standardize variables like lars does one <- rep( 1, n ) meanx <- drop( one %*% x ) / n xc <- scale( x, meanx, FALSE ) # first subtracts mean normx <- sqrt( drop( one %*% ( xc^2 ) ) ) names( normx ) <- NULL xs <- scale( xc, FALSE, normx ) # now rescales with norm (not sd) #================================================== # Perform Lasso regression out.ls <- lm( y ~ xs ) # ols fit on standardized beta.ols <- out.ls$coeff[ 2:( m + 1 ) ] # ols except for intercept w <- abs( beta.ols ) # weights for adaptive lasso xs <- scale( xs, center = FALSE, scale = 1 / w ) # xs times the weights object <- lars( xs, y, type = "lasso", normalize = FALSE ) #================================================== # Get min BIC # bic=log(n)*object$df+n*log(as.vector(object$RSS)/n) # rss/n version sig2f <- summary( out.ls )$sigma^2 # full model mse bic2 <- log( n )*object$df + as.vector( object$RSS ) / sig2f # Cp version step.bic2 <- which.min( bic2 ) # step with min BIC #================================================== # Minimum BIC step min.bic2.step <- step.bic2 - 1 #================================================== # Get the importance importance <- names( unlist( object$actions ) ) # Importance order if( min.bic2.step == 0 ) { importance.bic2.min <- '' } else { importance.bic2.min <- importance[ 1:min.bic2.step ] } predictor.occurance <- ( c( 'X1', 'X2', 'X1X2' ) %in% importance.bic2.min ) * 1 if( length( importance.bic2.min ) < 3 ) { importance.bic2.min.tuned <- rep( 'NA', 3 ) importance.bic2.min.tuned[1:length(importance.bic2.min)] <- importance.bic2.min } else { importance.bic2.min.tuned <- importance.bic2.min } #================================================== # Calculate R2 (R2 from training data) R2.min.BIC <- as.numeric( round( object$R2[ step.bic2 ], 4 ) ) n <- dim( training.set )[ 1 ] p <- length( importance.bic2.min ) R2.adj.min.BIC <- round( 1 - ( 1 - R2.min.BIC ) * ( ( n - 1 ) / ( n - p - 1 ) ), 4 ) #================================================== # Output return( list( R2 = R2.adj.min.BIC, predictor.occurance = predictor.occurance, min.bic2.step = as.numeric( min.bic2.step ) ) ) } #================================================== # HierNet optimisation hierNet.crossvalidation <- function( deletion, strong.option = FALSE, repeats = 1000 ) { #================================================== # Load data data.X1X2 <- data.frame( read.table( 'X1X2.data', header=FALSE ) ) colnames( data.X1X2 ) <- c( 'X1', 'X2' ) data.Y.raw <- data.frame( read.table( paste( deletion, 'Y.data', sep='_' ), header=FALSE ) ) colnames( data.Y.raw ) <- paste( 'Y', 1:ncol(data.Y.raw), sep='' ) # Box-Cox transformation of Y data data.Y <- apply( data.Y.raw, 2, function( y, lambda ) { return ( unlist( lapply( y, function( y.elem, lambda ) { return( ( y.elem^lambda - 1 ) / lambda ) } , lambda=lambda ) ) ) } , lambda=0.25 ) # , lambda=0.5 ) data.Y <- data.frame( data.Y ) # Input data X1 <- c( rep( 1, ncol( data.Y) * 2 ), rep( 0, ncol( data.Y) * 2 ) ) X2 <- c( rep( 1, ncol( data.Y) ), rep( 0, ncol( data.Y) ), rep( 1, ncol( data.Y) ), rep( 0, ncol( data.Y) ) ) Y <- c( as.vector( t( as.matrix( data.Y[ 1:2, ] ) ) ), as.vector( t( as.matrix( data.Y[ 3:4, ] ) ) ) ) Y.na <- which( is.na(Y), arr.ind=TRUE ) Y.index <- 1:length( Y ) Y <- Y[ !( Y.index %in% Y.na ) ] # Remove NA Y <- Y / sd( Y ) # standardize data so that coefficients are interpretable X1 <- X1[ !( Y.index %in% Y.na ) ] # Remove NA X2 <- X2[ !( Y.index %in% Y.na ) ] # Remove NA X1X2 <- X1 * X2 data.additive <- data.frame( cbind( X1, X2, Y ) ) data <- data.additive #================================================== # For each fold (1000 cases) split the data into a training- and testset (1/2) (randomly) n <- nrow( data ) folds <- cvFolds( n, K = 2, R = repeats ) #================================================== # Loop through each fold and perform analysis coefficients.list <- NULL for( i in 1:repeats ) { # Get training and test data from folds if( ( n %% 2 ) == 0 ) { training.rows <- folds$subsets[ seq( 1, n - 1, 2 ), i ] test.rows <- folds$subsets[ seq( 2, n, 2 ), i ] } else { training.rows <- folds$subsets[ seq( 1, n, 2 ), i ] test.rows <- folds$subsets[ seq( 2, n - 1, 2 ), i ] } # Training and test data training.set <- data[ training.rows, ] test.set <- data[ test.rows, ] # Apply hierNet path fit.path <- hierNet.path( as.matrix( training.set[ , 1:2 ] ), as.vector( training.set[ , 3 ] ), nlam = 100, maxlam = 100, flmin = 0.0001, strong = strong.option ) # Apply hierNet crossvalidation fit.cv <- hierNet.cv( fit.path, as.matrix( training.set[ , 1:2 ] ), as.vector( training.set[ , 3 ] ) ) # Get best lamda hat according to the ... best.lamhat <- round( fit.cv$lamhat.1se, 2 ) mat <- cbind(round( fit.cv$lamlist, 2), fit.cv$nonzero, round( fit.cv$cv.err, 2 ), round( fit.cv$cv.se, 2 ) ) dimnames( mat ) <- list(NULL, c( "Lambda", "Number of nonzero", "Mean CV error", "SE" ) ) mat <- data.frame( mat ) best.lamhat.index <- which( mat$Lambda == best.lamhat )[ 1 ] # Get coefficients that have been included in the analysis coefficients <- c( as.numeric( fit.path$bp[ , best.lamhat.index ] - fit.path$bn[ , best.lamhat.index ] ), fit.path$th[ , , best.lamhat.index ][ 1, 2 ] ) # Save output coefficients.list <- rbind( coefficients.list, coefficients ) } coefficients.list.df <- data.frame( coefficients.list != 0 ) predictors.hierNet <- apply( coefficients.list.df, 2, sum ) predictors.hierNet <- data.frame( t ( predictors.hierNet ) ) colnames( predictors.hierNet ) <- c( 'X1', 'X2', 'X1X2' ) best.model.hierNet <- table( apply( coefficients.list.df, 1, function( vec ) { return( length( which( vec == TRUE ) ) ) } ) + 1 ) best.model.numPredictors.hierNet <- rep( 0, 4 ) best.model.numPredictors.hierNet[ as.numeric( names( best.model.hierNet ) ) ] <- as.numeric( best.model.hierNet ) best.model.numPredictors.hierNet <- data.frame( t ( best.model.numPredictors.hierNet ) ) # colnames( best.model.numPredictors.hierNet ) <- paste( 0:3, 'Predictors' ) # Save results write.table( predictors.hierNet, file = paste( deletion, '_predictors_hierNet_', strong.option, '.txt', sep = '' ), row.names = FALSE ) write.table( best.model.numPredictors.hierNet, file = paste( deletion, '_numPredictors_hierNet_', strong.option, '.txt', sep = '' ), row.names = FALSE, col.names = FALSE ) } #================================================== #================================================== #================================================== # Main #================================================== # start time measurement start <- proc.time() #================================================== #================================================== # Set working directory setwd('~/Dropbox/projects/qcri/fGWAS/data/images_reproduction/linear_lasso_hierNet/') #================================================== # Variables deletions <- c( 'MCM22', 'PEP12', 'PEP7', 'PHO88', 'SKI8', 'VPS16', 'PHO88_MITO_KILL', 'PHO88_MITO_NO-KILL', 'SKI8_MITO_KILL', 'SKI8_MITO_NO-KILL', 'SKI7' ) R2s.linear <- NULL MSEs.linear <- NULL R2s.lasso.bic <- NULL predictors.lasso.bic <- NULL min.BIC.steps.lasso.bic <- NULL #================================================== # Loop through different deletion experiments and perform analysis for( deletion in deletions ) { print( '==================================================' ) print( deletion ) #================================================== # Load data data.X1X2 <- data.frame( read.table( 'X1X2.data', header=FALSE ) ) colnames( data.X1X2 ) <- c( 'X1', 'X2' ) data.Y.raw <- data.frame( read.table( paste( deletion, 'Y.data', sep='_' ), header=FALSE ) ) colnames( data.Y.raw ) <- paste( 'Y', 1:ncol(data.Y.raw), sep='' ) #================================================== # Box-Cox transformation of Y data data.Y <- apply( data.Y.raw, 2, function( y, lambda ) { return ( unlist( lapply( y, function( y.elem, lambda ) { return( ( y.elem^lambda - 1 ) / lambda ) } , lambda=lambda ) ) ) } , lambda=0.25 ) data.Y <- data.frame( data.Y ) #================================================== # Form input data matrix/dataframe X1 <- c( rep( 1, ncol( data.Y) * 2 ), rep( 0, ncol( data.Y) * 2 ) ) X2 <- c( rep( 1, ncol( data.Y) ), rep( 0, ncol( data.Y) ), rep( 1, ncol( data.Y) ), rep( 0, ncol( data.Y) ) ) Y <- c( as.vector( t( as.matrix( data.Y[ 1:2, ] ) ) ), as.vector( t( as.matrix( data.Y[ 3:4, ] ) ) ) ) Y.na <- which( is.na(Y), arr.ind=TRUE ) Y.index <- 1:length( Y ) Y <- Y[ !( Y.index %in% Y.na ) ] # Remove NA Y <- Y / sd( Y ) # standardize data so that coefficients are interpretable X1 <- X1[ !( Y.index %in% Y.na ) ] # Remove NA X2 <- X2[ !( Y.index %in% Y.na ) ] # Remove NA X1X2 <- X1 * X2 data.interaction <- data.frame( cbind( X1, X2, X1X2, Y ) ) #================================================== # Run linear model analyse lm.model <- linear.model.crossvalidation( data.interaction ) lm.model.R2s <- lm.model$R2s R2s.linear <- rbind( R2s.linear, c( apply( lm.model.R2s, 2, mean ), apply( lm.model.R2s, 2, sd ) ) ) lm.model.MSEs <- lm.model$MSEs MSEs.min.table <- table( apply( lm.model.MSEs, 1, function( vec ){ which( vec == min( vec ) ) } ) ) MSEs.min.freq <- rep( 0, 3 ) for( name in names( MSEs.min.table ) ) { MSEs.min.freq[ as.numeric( name ) ] <- as.numeric( MSEs.min.table[ names( MSEs.min.table ) == name ] ) } MSEs.linear <- rbind( MSEs.linear, MSEs.min.freq ) # Save output write.table( lm.model.R2s, file = paste( deletion, '_linear_R2s.txt', sep = '' ), row.names = FALSE, col.names = FALSE ) write.table( lm.model.MSEs, file = paste( deletion, '_linear_MSEs.txt', sep = '' ), row.names = FALSE, col.names = FALSE ) #================================================== # Run LASSO and BIC analysis lasso.bic.model <- LASSO.BIC.model.crossvalidation( data.interaction ) lasso.bic.predictors <- data.frame( lasso.bic.model$predictors ) colnames( lasso.bic.predictors ) <- c( 'X1', 'X2', 'X1X2' ) R2s.lasso.bic <- rbind( R2s.lasso.bic, c( mean( lasso.bic.model$R2s ), sd( lasso.bic.model$R2s ) ) ) predictors.lasso.bic <- rbind( predictors.lasso.bic, apply( lasso.bic.predictors, 2, sum ) ) table.min.BIC.steps <- table( lasso.bic.model$min.BIC.steps ) min.BIC.steps <- rep( 0, 4 ) min.BIC.steps[ as.numeric( names( table.min.BIC.steps ) ) + 1 ] <- as.numeric( table.min.BIC.steps ) min.BIC.steps.lasso.bic <- rbind( min.BIC.steps.lasso.bic, min.BIC.steps ) # Save output write( lasso.bic.model$R2s, file = paste( deletion, '_LASSO-BIC_R2s.txt', sep = '' ), ncolumns = 1 ) write.table( lasso.bic.predictors, file = paste( deletion, '_LASSO-BIC_predictors.txt', sep = '' ), row.names = FALSE ) write( lasso.bic.model$min.BIC.steps, file = paste( deletion, '_LASSO-BIC_minBIC-steps.txt', sep = '' ), ncolumns = 1 ) } #================================================== # Handle linear model output R2s.linear <- data.frame( R2s.linear ) rownames( R2s.linear ) <- deletions colnames( R2s.linear ) <- c( paste( 'R2mean', c( 'simple', 'additive', 'interaction'), sep = '_' ), paste( 'R2sd', c( 'simple', 'additive', 'interaction'), sep = '_' ) ) write.table( R2s.linear, file = 'LinearModel_R2s.txt' ) MSEs.linear <- data.frame( MSEs.linear ) rownames( MSEs.linear ) <- deletions colnames( MSEs.linear ) <- c( 'X1', 'X1+X2', 'X1+X2+X1X2' ) write.table( MSEs.linear, file = 'LinearModel_MSEs.txt' ) #================================================== # Handle LASSO+BIC model output R2s.lasso.bic <- data.frame( R2s.lasso.bic ) rownames( R2s.lasso.bic ) <- deletions colnames( R2s.lasso.bic ) <- c( 'R2mean', 'R2sd' ) write.table( R2s.lasso.bic, file = 'LASSO-BIC_R2s.txt' ) predictors.lasso.bic <- data.frame( predictors.lasso.bic ) rownames( predictors.lasso.bic ) <- deletions write.table( predictors.lasso.bic, file = 'LASSO-BIC_predictors.txt' ) min.BIC.steps.lasso.bic <- data.frame( min.BIC.steps.lasso.bic ) rownames( min.BIC.steps.lasso.bic ) <- deletions colnames( min.BIC.steps.lasso.bic ) <- paste( 0:3, 'Predictors' ) write.table( min.BIC.steps.lasso.bic, file = 'LASSO-BIC_minSteps.txt' ) #================================================== #================================================== # Generate plots #================================================== # 1. Plot...Compare simple, additive and interaction linear models according to R2 ci.l.matrix <- as.matrix( R2s.linear[ , 1:3] - R2s.linear[ , 4:6] ) negative.sd <- which( ci.l.matrix < 0, arr.ind = TRUE ) # if( dim( negative.sd ) == NULL ) # { # for( row in 1:dim( negative.sd )[ 1 ] ) # { # print( negative.sd[ row, ] ) # ci.l.matrix[ negative.sd[ row, 1 ], negative.sd[ row, 2 ] ] <- 0 # } # } else # { # ci.l.matrix[ negative.sd[ 1 ], negative.sd[ 2 ] ] <- 0 # } pdf( 'RecoveredHeritability_linearModels_R2.pdf', width = 12, height = 8 ) bp <- barplot2( t( R2s.linear[ , 1:3] ), beside = TRUE, horiz = FALSE, ylim = c( 0, 1 ), ylab='Fraction of phenotypic variance explained', axisnames = FALSE, plot.grid = TRUE, ci.lwd = 1.5, cex.axis = 1.2, cex.lab= 1.5, ci.u = t( as.matrix( R2s.linear[ , 1:3] + R2s.linear[ , 4:6] ) ), ci.l = t( ci.l.matrix ), plot.ci = TRUE ) legend( "topleft", legend = c( 'Simple', 'Additive', 'Interaction' ), fill = c( 'red', 'orange', 'yellow'), text.font = 2 ) axis(1, labels = FALSE, tick = FALSE, at=seq(2.5,38.5,4) ) labels <- c( deletions[ 1:6 ], paste( 'PHO88', c( '(killer)', '(non-killer)' ), sep='' ), paste( 'SKI8', c( '(killer)', '(non-killer)' ), sep='' ), deletions[11] ) text( bp[ 1, ] + 1, par("usr")[ 3 ] - 0.005, labels = labels, srt = 40, offset = 1, adj = 1, xpd = TRUE, font = 4 ) dev.off() #================================================== #================================================== # 2. Plot...Compare linear with LASSO+BIC model according to R2 ci.l.matrix <- as.matrix( cbind( R2s.linear[ , 3] - R2s.linear[ , 6], R2s.lasso.bic[ , 1 ] - R2s.lasso.bic[ , 2 ] ) ) negative.sd <- which( ci.l.matrix < 0, arr.ind = TRUE ) # for( row in 1:dim( negative.sd )[ 1 ] ) # { # print( negative.sd[ row, ] ) # ci.l.matrix[ negative.sd[ row, 1 ], negative.sd[ row, 2 ] ] <- 0 # } pdf( 'RecoveredHeritability_linearModels-LASSO-BIC_R2.pdf', width = 12, height = 8 ) bp <- barplot2( t( cbind( R2s.linear[ , 3], R2s.lasso.bic[ , 1 ] ) ), beside = TRUE, horiz = FALSE, ylim = c( 0, 1 ), ylab='Fraction of phenotypic variance explained', axisnames = FALSE, plot.grid = TRUE, ci.lwd = 1.5, cex.axis = 1.2, cex.lab= 1.5, ci.u = t( as.matrix( cbind( R2s.linear[ , 3] + R2s.linear[ , 6], R2s.lasso.bic[ , 1 ] + R2s.lasso.bic[ , 2 ] ) ) ), ci.l = t( ci.l.matrix ), plot.ci = TRUE, col = c( 'yellow', 'blue' ) ) legend( "topleft", legend = c( 'Linear', 'LASSO+BIC' ), fill = c( 'yellow', 'blue' ), text.font = 2 ) axis(1, labels = FALSE, tick = FALSE, at=seq(2.5,38.5,4) ) labels <- c( deletions[ 1:6 ], paste( 'PHO88', c( '(killer)', '(non-killer)' ), sep='' ), paste( 'SKI8', c( '(killer)', '(non-killer)' ), sep='' ), deletions[11] ) text( bp[ 1, ] + 1, par("usr")[ 3 ] - 0.005, labels = labels, srt = 40, offset = 1, adj = 1, xpd = TRUE, font = 4 ) dev.off() #================================================== #================================================== #================================================== # Apply hierNet Crossvalidation (1000 repeats) on multiple processes (hierNet strong option set to FALSE) strong.option <- FALSE # foreach(i = 1:length( deletions ) ) %dopar% # { # hierNet.crossvalidation( deletions[ i ], strong.option = strong.option ) # } hierNet.crossvalidation( deletions[ 11 ], strong.option = strong.option ) # Handle output num.predictors.hierNet <- NULL for( deletion in deletions ) { num.predictors.hierNet <- rbind( num.predictors.hierNet, read.table( paste( deletion, '_numPredictors_hierNet_', strong.option, '.txt', sep = '' ), header = FALSE ) ) } num.predictors.hierNet <- data.frame( num.predictors.hierNet ) rownames( num.predictors.hierNet ) <- deletions colnames( num.predictors.hierNet ) <- paste( 0:3, 'Predictors' ) write.table( num.predictors.hierNet, file = paste( 'HierNet_minSteps_', strong.option, '.txt', sep = '' ) ) predictors.hierNet <- NULL for( deletion in deletions ) { predictors.hierNet <- rbind( predictors.hierNet, read.table( paste( deletion, '_predictors_hierNet_', strong.option, '.txt', sep = '' ), header = TRUE ) ) } predictors.hierNet <- data.frame( predictors.hierNet ) rownames( predictors.hierNet ) <- deletions write.table( predictors.hierNet, file = paste( 'HierNet_predictors_', strong.option, '.txt', sep = '' ) ) #================================================== #================================================== #================================================== #================================================== # Apply hierNet Crossvalidation (1000 repeats) on multiple processes (hierNet strong option set to TRUE) strong.option <- TRUE # foreach(i = 1:length( deletions ) ) %dopar% # { # hierNet.crossvalidation( deletions[ i ], strong.option = strong.option ) # } hierNet.crossvalidation( deletions[ 11 ], strong.option = strong.option ) # Handle output num.predictors.hierNet <- NULL for( deletion in deletions ) { num.predictors.hierNet <- rbind( num.predictors.hierNet, read.table( paste( deletion, '_numPredictors_hierNet_', strong.option, '.txt', sep = '' ), header = FALSE ) ) } num.predictors.hierNet <- data.frame( num.predictors.hierNet ) rownames( num.predictors.hierNet ) <- deletions colnames( num.predictors.hierNet ) <- paste( 0:3, 'Predictors' ) write.table( num.predictors.hierNet, file = paste( 'HierNet_minSteps_', strong.option, '.txt', sep = '' ) ) predictors.hierNet <- NULL for( deletion in deletions ) { predictors.hierNet <- rbind( predictors.hierNet, read.table( paste( deletion, '_predictors_hierNet_', strong.option, '.txt', sep = '' ), header = TRUE ) ) } predictors.hierNet <- data.frame( predictors.hierNet ) rownames( predictors.hierNet ) <- deletions write.table( predictors.hierNet, file = paste( 'HierNet_predictors_', strong.option, '.txt', sep = '' ) ) #================================================== #================================================== #================================================== # Measure and print out time needed for the script end <- proc.time() duration <- end-start print( paste( "Script duration:", round( duration[3] / 60, 2 ), "min") ) #================================================== # Main end #================================================== #================================================== #==================================================