I am a beginner with Julia and tried it out by fitting data.
I got unexpected results for vcov(fit) and had a look in the source code.

There:

I noticed 2 issues:

1. if one has no weights (errors of the values to be fit), the Covariance matrix is calculated using QR decomposition of the Jacobian but when one has weights no QR decomposition is used.
   This is strange because both methods deliver the same results but the QR decomposition might be more stable numerically
2. the MSE is ignored when one has weights. I searched in literature and cannot find why one should omit the MSE. Also Mathlab uses the MSE:
   https://www.mathworks.com/help/stats/nlinfit.html#btk7ign-CovB

The second point is my main problem because the errors of the fit model parameters were larger than expected.
I think it is a bug to omit the MSE. When there are weights, the MSE includes them already, meaning the MSE does not become 1 because of the presence of weights. And therefore it cannot just be omitted.

I the tutorial I see there correctly:
https://julianlsolvers.github.io/LsqFit.jl/latest/tutorial/#Weighted-Least-Squares-1
Cov(γ∗)=σ²[J′WJ]^−1
and σ² is the MSE
there you also write

the covariance of the estimated parameter is calculated as covar = inv(J'*fit.wt*J)

But this is incorrect because the dimensions won't fit and I get the error:
DimensionMismatch and looking at the source code the calculations is
covar = covar = inv(J' * J)
so the weights are ignored.

I use Julia 1.11.1 and LsqFit v0.15.0

p.s. as a newbie I might report the error at the wrong place, please point me then to the correct place.

p.s. the LsqFit docs:
https://julianlsolvers.github.io/LsqFit.jl/latest/tutorial
I see the command ''estimate_var'' but this is marked deprecated.