In neural networks, how to tell the feature which contributes the most to the output value? I have a neural network which uses input features $N$, some input layers $L$ which predict a continuous variable $X$. Can we say which features or combination of 2 features of the initial $N$ features contribute the most to the change in $X$?
Example: Suppose I am trying to predict the house prices on basis of number of rooms, area of house, porch area, length, width. 
Can I identify the feature $x$ or combination of 2 features that causes the most change in prices of house?
I was thinking of a brute force approach by varying one feature and keeping rest constant and see varies the prices the most. If there's a standard way, please share some resources to read.
 A: This is more of a statistics question than a specific programming one.  If you have your heart set on using neural nets an example using feature selection with Garson's algorithm is here.  Below I have provided the code that you can try.  Hopefully this can give you something to start with.
But please note this is only one possible answer.  There are many other approaches people have taken as this is an active area of research (neural networks are complex!).  There are very likely other methods that may be more suitable, more efficient, etc.  You may not want to even use neural nets (I don't know your specific reasons).  Depending on your data it may be better to use some alternative feature selection up front before the neural net.  A simple google scholar search for 'neural network feature selection' will return several papers on the matter.  There are many very strong opinions on the subject of neural networks so be warned there is no definitive answer out there.
# code from link noted above (slightly updated)
require(clusterGeneration)
require(nnet)

#define number of variables and observations
set.seed(2)
num.vars<-8
num.obs<-10000

#define correlation matrix for explanatory variables
#define actual parameter values
cov.mat<-genPositiveDefMat(num.vars,covMethod=c("unifcorrmat"))$Sigma
rand.vars<-mvrnorm(num.obs,rep(0,num.vars),Sigma=cov.mat)
parms<-runif(num.vars,-10,10)
y<-rand.vars %*% matrix(parms) + rnorm(num.obs,sd=20)

#prep data and create neural network
y<-data.frame((y-min(y))/(max(y)-min(y)))
names(y)<-'y'
rand.vars<-data.frame(rand.vars)
mod1<-nnet(rand.vars,y,size=8,linout=T)

require(devtools)

#import 'gar.fun' from beckmw's Github - this is Garson's algorithm
source_gist('6206737')

#use the function on the model created above
gar.fun('y',mod1)

Here is the output plot.  You can see there are both positive and negative values.  The negative and positive values reflect negative and positive relationships between the variable and the response variable.  
A: caret package offers variable importance for neural network models along with others. You can refer to this here which says:

Neural Networks: The method used here is based on Gevrey et al (2003), which uses combinations of the absolute values of the weights. For classification models, the class-specific importances will be the same.

It also includes feature selection by Recursive Feature:

Recursive Feature Elimination: Variable importance is computed using the ranking method used for feature selection. For the final subset size, the importances for the models across all resamples are averaged to compute an overall value. 

Here's what you can do :
library(caret)

# varimp
control <- trainControl(method="repeatedcv", number=10, repeats=3)
model <- train(y~., data=train, method="nnet",trControl=control)
imp<-varImp(model)
plot(imp)

# rfe
control <- rfeControl(functions=rfFuncs, method="cv", number=10)
results <- rfe(train[,!colnames %in% c("y")],train$y, sizes=c(1:(ncol(train)-1), rfeControl=control)
 print(results)
# features selected 
predictors(results)
plot(results, type=c("g", "o"))

