How to fit the complex dielectric function of gold

Question

I want to fit the dielectric function of gold measured by Johnson & Christy:

( *omega, epsilon1+I*epsilon2 *)
data={{6.59912, 0.227056 + 3.04128 I}, {6.4693, 
0.295191 + 3.17592 I}, {6.34905, 0.292524 + 3.28568 I}, {6.21911, 
0.203899 + 3.32766 I}, {6.09918, 
0.138171 + 3.39682 I}, {5.97918, -0.010416 + 
3.3904 I}, {5.84914, -0.1325 + 3.51 I}, {5.7294, -0.233769 + 
3.6062 I}, {5.59925, -0.346329 + 3.7102 I}, {5.47924, -0.4155 + 
3.8252 I}, {5.35939, -0.551009 + 3.8922 I}, {5.22941, -0.616896 + 
4.05504 I}, {5.10938, -0.744529 + 
4.16328 I}, {4.97928, -0.891261 + 4.33846 I}, {4.85946, -1.08044 + 
4.49008 I}, {4.73946, -1.2365 + 4.7223 I}, {4.60942, -1.34641 + 
4.97628 I}, {4.48942, -1.36651 + 5.28242 I}, {4.3595, -1.33226 + 
5.49486 I}, {4.2395, -1.30678 + 5.59644 I}, {4.11949, -1.22742 + 
5.78034 I}, {3.98945, -1.24255 + 5.79258 I}, {3.86955, -1.2308 + 
5.84584 I}, {3.73953, -1.35529 + 5.57368 I}, {3.61955, -1.31024 + 
5.53816 I}, {3.49951, -1.23196 + 5.598 I}, {3.36959, -1.40063 + 
5.6092 I}, {3.24957, -1.60489 + 5.64436 I}, {3.11957, -1.6494 + 
5.73888 I}, {2.99964, -1.70216 + 5.71736 I}, {2.8796, -1.6922 + 
5.6492 I}, {2.74964, -1.759 + 5.28264 I}, {2.62968, -1.7027 + 
4.84438 I}, {2.49968, -2.27829 + 3.81264 I}, {2.37969, -3.94616 + 
2.58044 I}, {2.25972, -5.84213 + 2.1113 I}, {2.12973, -8.11267 + 
1.66054 I}, {2.00973, -10.6619 + 1.37424 I}, {1.87977, -13.6482 + 
1.03516 I}, {1.75976, -16.8177 + 1.06678 I}, {1.63979, -20.6102 + 
1.27176 I}, {1.50981, -25.8113 + 1.62656 I}, {1.38982, -32.0407 + 
1.92542 I}, {1.25984, -40.2741 + 2.794 I}, {1.13985, -51.0496 + 
3.861 I}, {1.01987, -66.2185 + 5.7015 I}, {0.889886, -90.4265 + 
8.18634 I}, {0.769902, -125.351 + 
12.5552 I}, {0.639918, -189.042 + 25.3552 I}}


dataReal=Transpose@{data[[All,1]],Re@data[[All,2]]};
dataImag=Transpose@{data[[All,1]],Im@data[[All,2]]};

with the following model (Drude-Lorentz-Model):

$$ \varepsilon=\varepsilon_\infty-\frac{\omega^2_D}{\omega\left(\omega+i\gamma_D \right)}-\sum_{i=1}^{n}\frac{\Delta\epsilon_i\omega ^2_{L,i}}{\omega^2-\omega^2_{L,i}+i\gamma_{L,i}\omega} $$ for $n=2$ and $n=3$. With $\epsilon_\infty=9.84$, $\omega_D=9$ and $\gamma_D=67/1000$

So far I'm fitting real/imaginary part of the function separately to the real/imaginary part of the data (dataReal, dataImag):

nlmReal=NonlinearModelFit[
dataReal, {ComplexExpand[Re[model]],cons}, {{wL1}, {wL2}, {gL1}, {gL2},{de1},{de2}},\[Omega],
Method -> {NMinimize, 
Method -> {"DifferentialEvolution",
"ScalingFactor" -> 0.9,"CrossProbability" -> 0.1, 
  "PostProcess" -> {FindMinimum, Method -> "QuasiNewton"}}}];


nlmImag=NonlinearModelFit[
dataImag, {ComplexExpand[Re[model]],cons}, {{wL1}, {wL2}, {gL1}, {gL2},{de1},{de2}},\[Omega],
Method -> {NMinimize, 
Method -> {"DifferentialEvolution",
"ScalingFactor" -> 0.9,"CrossProbability" -> 0.1, 
  "PostProcess" -> {FindMinimum, Method -> "Newton"}}}];

with

cons={wL1>0&&wL2>0&&gL1>0&&gL2>0&&de1>0&&de2>0}

No matter if I use Newton, QuasiNewton, LevenbergMarquardt or else and adding cons to the fit or not, the parameter of the fits are always different.

Grid[{{nlmReal[{"ParameterTable", "EstimatedVariance","AdjustedRSquared"}], 
  nlmImag[{"ParameterTable", "EstimatedVariance","AdjustedRSquared"}]}}, Frame -> All]

For $n=3$ it is the same. How can I achieve a better and more accurate fit for the given model? Can the fit be done without separating into real and imaginary part?