title: Combine multiple neural network models
body: <p>I run a loop 200 times within which I:</p>
  
  <ul>
  <li><p>randomly divide my dataset into training and test sets</p></li>
  <li><p>fit a neural networks model with R's <code>nnet()</code> command on the training set</p></li>
  <li><p>assess performance on the test set</p></li>
  </ul>
  
  <p>I save each model to a list.</p>
  
  <p>Now I want to use a combined model to make out-of-sample predictions. I've used the <code>combine</code> function for this purpose on <code>randomForest</code> objects. <strong>Is there a similar combine command for <code>nnet</code> objects?</strong></p>
  
  <p>I cannot upload the dataset but below is the code I'm currently working with. It works as is, except for the last line where I seek a command to combine models.</p>
  
  <pre><code>    n &lt;- 200
      nnet_preds &lt;- matrix(NA,  ncol = 1,  nrow = n)
      nnet_predstp &lt;- matrix(NA,  ncol = 1,  nrow = n)
      nnet_predstn &lt;- matrix(NA,  ncol = 1,  nrow = n)
      nnet_predsfptp &lt;- matrix(NA,  ncol = 1,  nrow = n)
      nnet_predsfntp &lt;- matrix(NA,  ncol = 1,  nrow = n)
      NN_predictions &lt;- matrix(NA,  ncol = 1,  nrow = 10) 
      outcome_2010_NN &lt;- testframe2[, "ytest"] 
      nn_model &lt;- vector("list", n)
  
      data&lt;-testframe2[, c("sctownpop", "sctownhh", "scnum_tribes", "sctownmusl_percap", "scmuslim", "scprop_stranger", "sctownstrg_percap", "scprop_domgroup", "scexcom", "sctownexcm_percap", "sctownretn_percap", "scprop_under30", "scprop_male", "scedulevel", "scprop_noeduc", "scprop_anypeace", "scgroup_prtcptn", "scpubcontr", "scsafecommdum", "scciviccommdum", "scoll_action_add", "scngodependent", "scgovtdependent", "scpolicourtscorr", "screlmarry", "scmslmnolead", "sccrime_scale", "scviol_scale", "sclandconf_scale", "sctownnlnd_percap", "scnolandnofarm", "scfarmocc", "scunemployed", "scwealthindex", "scwealth_inequality", "scviol_experienced", "scviol_part", "scanylndtake", "scdisp_ref", "sfacilities", "sfreq_visits", "sctot_resources", "scmeanprice_above75", "scmosquesdum", "scmnrt_ldrshp", "scany_majorconf", "sstate_presence", "sremote", "scmobilec", "scradio_low")]
  
      data = cbind(outcome_2010_NN, data)
  
      sampleSplit = round(nrow(data)-(nrow(data)/5))
  
      for(i in 1:n) {     
  
  set.seed(06511+i)
      data_randomization &lt;- data[sample(1:nrow(data),  dim(data)[1],  replace=FALSE), ]
  
      train &lt;- data_randomization[1:sampleSplit, ]
      test  &lt;- data_randomization[(sampleSplit+1):nrow(data), ]
  
      nn_model[[i]] &lt;- nnet(outcome_2010_NN ~ sctownpop +           sctownhh+ scnum_tribes+ sctownmusl_percap+ scmuslim+          scprop_stranger+   sctownstrg_percap+ scprop_domgroup+     scexcom+  sctownexcm_percap+   sctownretn_percap+   scprop_under30 +  scprop_male+         scedulevel+          scprop_noeduc+       scprop_anypeace+     scgroup_prtcptn+     scpubcontr+          scsafecommdum+       scciviccommdum+      scoll_action_add+    scngodependent+      scgovtdependent+     scpolicourtscorr+    screlmarry+          scmslmnolead+        sccrime_scale+       scviol_scale+        sclandconf_scale+    sctownnlnd_percap+   scnolandnofarm+      scfarmocc+           scunemployed+        scwealthindex+       scwealth_inequality+ scviol_experienced+  scviol_part+         scanylndtake+        scdisp_ref+          sfacilities+         sfreq_visits+        sctot_resources+     scmeanprice_above75+ scmosquesdum+        scmnrt_ldrshp+       scany_majorconf+     sstate_presence+     sremote+             scmobilec+           scradio_low, 
      data=train,  size = 3,  decay = 0.1)# size=number of units/nodes in the (single_hidden layer); decay=parameter for weight decay. Default 0.
  
      predictions &lt;- predict(nn_model[[i]],  test)
  
      nnpredorder&lt;-rank(predictions)
      nncvpredictionsA50 &lt;- ifelse( nnpredorder &gt; 24,  1,  0 )    # manually optimized
  
      errors &lt;- table(test[, "outcome_2010_NN"],  nncvpredictionsA50)             
  
      accuracy.rate &lt;- (errors[1, 1]+errors[2, 2])/sum(errors)
      true.pos.rate &lt;- (errors[2, 2]/(errors[2, 2]+errors[2, 1]))
      true.neg.rate &lt;- (errors[1, 1]/(errors[1, 1]+errors[1, 2]))
      FPTP &lt;- (errors[1, 2]/errors[2, 2])
      FNTP &lt;- (errors[2, 1]/errors[2, 2])
  
      nnet_preds[i, ] &lt;- accuracy.rate
      nnet_predstp[i, ] &lt;- true.pos.rate
      nnet_predstn[i, ] &lt;- true.neg.rate
      nnet_predsfptp[i, ] &lt;- FPTP
      nnet_predsfntp[i, ] &lt;- FNTP
  }
  
      mean(nnet_preds); sd(nnet_preds)
      mean(nnet_predstp); sd(nnet_predstp)
  
      NN_predictions[1, ] &lt;- mean(nnet_predstp) # TP accuracy rate (sensitivity)
      NN_predictions[2, ] &lt;- sd(nnet_predstp) # TP accuracy rate (sensitivity)
      NN_predictions[3, ] &lt;- mean(nnet_predstn)  # TN accuracy rate (specificity)
      NN_predictions[4, ] &lt;- sd(nnet_predstn) # TN accuracy rate (specificity)
      NN_predictions[5, ] &lt;- mean(nnet_preds)  # Accuracy rate
      NN_predictions[6, ] &lt;- sd(nnet_preds) # Accuracy rate
      NN_predictions[7, ] &lt;- mean(nnet_predsfptp)  # Ratio FP:TP
      NN_predictions[8, ] &lt;- sd(nnet_predsfptp) # Ratio FP:TP
      NN_predictions[9, ] &lt;- mean(nnet_predsfntp)  # Ratio FN:TP
      NN_predictions[10, ] &lt;- sd(nnet_predsfntp) # Ratio FN:TP
  
      print(NN_predictions)
  
  ### Combine NN models  #Where `combine` is the randomForest command     
  aggNNmodel &lt;- do.call(combine, nn_model)
  </code></pre>