title: Using ggplot2 to plot predicted values with robust standard errors
body: <p>I am trying to use ggplot2 to plot the predicted values of negative binomial regression, one with a binary variable turned on, and another with it turned off. So there will be two two plots that can be compared.</p>
  
  <p>The link <a href="http://www.ats.ucla.edu/stat/r/dae/nbreg.htm" rel="nofollow">here</a> demonstrates how to do it at the bottom of the page, but I want to be able to create shading around the plot of the predicted values using robust standard errors.  I don't see how to get this from the predict() function.  Is there any work around from this code example to get robust standard errors to shade around the plotted lines?</p>
  
  <p>I use the code here from <a href="http://www.ats.ucla.edu/stat/r/dae/poissonreg.htm" rel="nofollow">this site</a> to generate robust standard errors:</p>
  
  <pre><code>require(sandwich)
  cov.nb1 &lt;- vcovHC(nb1, type = "HC0")
  std.err &lt;- sqrt(diag(cov.nb1))
  r.est &lt;- cbind(Estimate = coef(nb1), `Robust SE` = std.err, `Pr(&gt;|z|)` = 2 *
      pnorm(abs(coef(nb1)/std.err), lower.tail = FALSE), LL = coef(nb1) - 1.96 *
      std.err, UL = coef(nb1) + 1.96 * std.err)
  
  r.est
  </code></pre>
  
  <p>the model I am using is this:</p>
  
  <pre><code>nb1 &lt;- glm.nb(citecount ~ expbin*novcr + expbin*I(novcr^2) + disease + length +
  as.factor(year), data = nov4d.dt)
  </code></pre>
  
  <p>And a sample of the data I am using is this:</p>
  
  <pre><code>nov4d.dt  &lt;-
      structure(list(PMID = c(1279136L, 1279186L, 1279186L, 1279187L, 
      1279187L, 1279190L, 1279257L, 1279317L, 1279332L, 1279523L), 
          min = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), max = c(32L, 
          32L, 32L, 32L, 32L, 32L, 32L, 32L, 32L, 32L), mean = c(11L, 
          13L, 13L, 19L, 19L, 16L, 24L, 15L, 8L, 19L), length = c(45L, 
          120L, 120L, 78L, 78L, 136L, 45L, 36L, 171L, 78L), threslength = c(13L, 
          20L, 20L, 7L, 7L, 26L, 4L, 6L, 77L, 14L), novlength = c(5L, 
          6L, 6L, 3L, 3L, 6L, 3L, 3L, 36L, 5L), novind = c("TRUE", 
          "TRUE", "TRUE", "TRUE", "TRUE", "TRUE", "TRUE", "TRUE", "TRUE", 
          "TRUE"), novcr = c(0.111111, 0.05, 0.05, 0.0384615, 0.0384615, 
          0.0441176, 0.0666667, 0.0833333, 0.210526, 0.0641026), novcrt = c(0.288889, 
          0.166667, 0.166667, 0.0897436, 0.0897436, 0.191176, 0.0888889, 
          0.166667, 0.450292, 0.179487), year = c(1991L, 1991L, 1992L, 
          1992L, 1992L, 1992L, 1992L, 1992L, 1991L, 1992L), disease = structure(c(1L, 
          4L, 2L, 4L, 2L, 1L, 4L, 4L, 2L, 4L), .Label = c("alz", "bc", 
          "cl", "lc"), class = "factor"), citecount = c(5L, 8L, 8L, 
          12L, 12L, 0L, 1L, 0L, 92L, 0L), novind2 = c(TRUE, TRUE, TRUE, 
          TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE), rad = c(FALSE, 
          FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE
          ), exp = c(260, 351, 351, 65, 65, 480, 104, 273, 223, 0), 
          novind4 = c(FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, 
          FALSE, TRUE, FALSE), novind5 = c(FALSE, FALSE, FALSE, FALSE, 
          FALSE, FALSE, FALSE, FALSE, TRUE, FALSE), novind6 = c(FALSE, 
          FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE
          ), expbin = c(TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, 
          TRUE, TRUE, FALSE), expbin2 = c(TRUE, TRUE, TRUE, FALSE, 
          FALSE, TRUE, FALSE, TRUE, TRUE, FALSE)), .Names = c("PMID", 
      "min", "max", "mean", "length", "threslength", "novlength", "novind", 
      "novcr", "novcrt", "year", "disease", "citecount", "novind2", 
      "rad", "exp", "novind4", "novind5", "novind6", "expbin", "expbin2"
      ), sorted = "PMID", class = c("data.frame"), row.names = c(NA, 
      -10L))
  </code></pre>