Since MOSFET has finite output resistance in saturation/active mode, the slope of unsignificanlty rising drain current is defined by Ua and slope parameter as lambda:
This parameter (as I know) is not given in any MOSFET datasheet.
Question: Is there any other way to get slope parameter out of the equation? For example with transconductance it can be done but I don't really know how to get it out of it (so this parameter can be defined with other already given variable of parameter).
Although not typically listed on datasheets, the MOSFET parameter \$\lambda\$ can be sometimes be found in the SPICE model provided by the manufacturer. Consider for example the N-channel MOSFET 2N7002. An old SPICE model from Zetex defines \$\lambda\$ of 2N7002 as \$267\cdot 10^{-6} \textrm{ V}^{-1}\$. For more information about how \$\lambda\$ is defined in SPICE models, see pages 128-129 of the HSPICE manual.
However, it is important to keep in mind that the parameter \$\lambda\$ may vary significantly from device to device, so it would be unwise to design a circuit that is sensitive to this parameter. In addition, channel length modulation, which \$\lambda\$ models, is only part of the story for determining MOSFET output resistance (\$r_o\$). Other effects such as drain-induced barrier lowering and substrate current induced body effect (SCBE) may also be important depending on the bias point of the device.
In summary, you may be able to get a very rough idea of the output impedance of a MOSFET by looking at the manufacturer's SPICE model. But the reality is that \$r_o\$ varies from device to device and is a complex function of bias point.
The resistance is a function of the current flow in the channel, which is actually driven by the Poisson's charge transport. You can think of the resistance as the components of the vertical field, the horizontal field and collision; however, you can just work back from the current equations. Starting from the EKV model, I got this:
$$ I_{on}=\frac{W}{L}\frac{\mu C_{cox}}{2 \kappa}\left( 2\kappa\left(V_g-V_{T0}\right)\left(V_{ds}\right) + \left(V_{d}-V_{ds}\right)^2 - V_{d}^2 \right). $$ which I then can substitute with \$V=IR\$ to get:
$$ R_{on} =\frac{V_{ds}}{I_{on}}= \frac{\frac{L}{W}\frac{2 \kappa}{\mu_g C_{cox}}} {\left( 2\kappa\left(V_g-V_{T0}\right) -V_{d} -V_s. \right)}. $$
I dug through my data to find a sweep of a nFET, which is not in saturation to show the difference between real data, BSIM 3.3 and what I did above. The resistance graph looks as
where I fixed 25mV across the device to keep it out of saturation. This is different from what you wanted, but it's what I had as far as data. You can drop the drain terms in the equations. The BSIM models from a commercial process happened to not be all that great for this type of simulation as it's just one of those places where BSIM doesn't work well. Due to the threshold voltage, you can see the resistance go up as we approach subthreshold.
