Order of Predictors in R Neuralnet formula I am having trouble understanding why the order of predictors in the formula of the nn function of R neuralnet package changes the resulting network. Here is a reproducible example with 2 inputs, and a categorical output.
Initialization of the data 
set.seed(1)
density=runif(100,0,1)
residual.sugar=runif(100,0,1)
quality=round(runif(100,0,1))

Constructing the data frame
mydata=data.frame(density, residual.sugar, quality)

First neural network. Notice the order of predictors - density first then residual sugar
set.seed(3)
nn<-neuralnet(quality~density+residual.sugar, linear.output=FALSE, 
hidden=2, data=mydata, act.fct="logistic", rep=1)
plot(nn)

Second neural network. Notice that I changed the order of predictors - now residual sugar comes first followed by density
set.seed(3)
nn<-neuralnet(quality~residual.sugar+density, linear.output=FALSE, 
hidden=2, data=mydata, act.fct="logistic", rep=1)
plot(nn)

I expect the 2 runs of the nn function to produce networks that are equivalent to each other, in other words I expect that although the weights and biases would swap on the network - the error and the number of steps that it takes for the network to converge should be identical. However, as you will see when you run the above code the 2 networks are drastically different. Why?
 A: Fixing the seed means that the random initializations for the weights are fixed between the two models. However, changing the order of the features changes which feature is assigned to which weights. This means that even though the weight matrices are the same, you have two different models, and the update procedures will be different. 
Let's consider an example. Suppose that in both models, a weight matrix is initialized as
$$
\begin{bmatrix}
-0.13 & 0.12 \\
0.29 & -0.37
\end{bmatrix}
$$
For the first model, the output of the affine transformation for that unit is given by
$$
\begin{bmatrix}
-0.13 & 0.12 \\
0.29 & -0.37
\end{bmatrix}
\begin{bmatrix}
\text{density}\\
\text{residual sugar}
\end{bmatrix} + \text{bias}.
$$
However, for the second model the output of the affine transformation is 
$$
\begin{bmatrix}
-0.13 & 0.12 \\
0.29 & -0.37
\end{bmatrix}
\begin{bmatrix}
\text{residual sugar} \\
\text{density}
\end{bmatrix} + \text{bias}.
$$
which is clearly a different result whenever $\text{residual sugar} \neq \text{density}$. This difference in result means that the backpropagation will produce different results between the two models, and therefore the weight updates will be different between the two models. 
Neural networks are well-known to be non-convex optimization problems, so it's not surprising that the optimization does not recover the same solutions between the two models. Because of the non-convexity of the model, there are many weight configurations which will have the same loss value.
