As @ user20650 noted, the argument softmaxis different. Inside nnet.formulathere is a section:
if (length(lev) == 2L) {
y <- as.vector(unclass(y)) - 1
res <- nnet.default(x, y, w, entropy = TRUE, ...)
res$lev <- lev
}
else {
y <- class.ind(y)
res <- nnet.default(x, y, w, softmax = TRUE, ...)
res$lev <- lev
}
Here softmaxset to TRUE. Putting it in a call nnetfixes the problem and they now match.
fit <- nnet(x, y, size = 3, decay = .1, softmax = TRUE)
pred <- predict(fit, iris[,1:4])
rowSums(head(pred))