nlm = NonlinearModelFit[data, Exp[a x/(b + c x)], {a, b, c}, x];
{bands90[x_], bands95[x_], bands99[x_], bands999[x_]} =
Table[nlm["MeanPredictionBands",
ConfidenceLevel -> cl], {cl, {.9, .95, .99, .999}}];
Show[ListPlot[data],
Plot[{nlm[x], bands90[x], bands95[x], bands99[x], bands999[x]}, {x,
1, 15}, Filling -> {2 -> {1}, 3 -> {2}, 4 -> {3}, 5 -> {4}}],
AspectRatio -> 4/5]