Hello! I have some data of the form (x,y,z) which I know it is described by a function of the form: ##z=y(a+bx)##, where a and b are parameters to be fitted for. z and y have error associated to them while x doesn't (x is actually an integer going from 0 to 3 for each value of y). I tried to do the fit in 2 different ways. I first made a linear fit of the form ##z=yA## for A (I used this package which accounts for the error on both z and y: https://docs.scipy.org/doc/scipy/reference/odr.html) for each value of x, then I made a fit of the form ##A=a+bx## for a and b, with the error on A obtained from the first fit. In the end I get a value and error for a and b. A second method I used was to fit directly ##z=y(a+bx)## to the whole data at once (it is not really a linear fit anymore, but it can be easily done in Python, with the same package as mentioned above). Now I get a new set of values and errors for a and b. The values obtained using the 2 methods are consistent with each other (within the errors on a and b), but using the first method gives a smaller error than in the second method. Is there anything I am missing? Shouldn't I get the exactly same result both ways? And in case the answer is no, which method should I use and why? Thank you!