What is really proportional to the distance from the neutral axis (let's denote it $z$) is the strain.

For pure bending and a symmetric cross-section the equation for strain is given by
$$\epsilon(z) = \frac{z}{\rho}$$
where:
- $\rho$ or $R$ is the radius of curvature of the beam

Because in the linear region the constitutive equation is $\sigma(z) = E\cdot \epsilon(z)$, the development of stress at distance z is proportional to the strain.
Because stress in a small area $\Delta A$ is defined as $\sigma(z) =\frac{\Delta F}{\Delta A}$, in that small area $\Delta A$ which is at distance z you are getting a Force :
$$\Delta F(z) = E\cdot \Delta A\cdot \epsilon(z)$$
Those forces produce bending moments the further away they are from the neutral axis (about which the cross-section rotates), the contribution of the bending moment is greater.
$$\Delta M(z) = \Delta F(z) \cdot z= E\cdot \Delta A\cdot \epsilon(z)\cdot z$$
So at a distance (z) any small area $\Delta A$ will produce a bending moment equal to
$$\Delta M(z) = \frac{E\cdot}{\rho} \Delta A\cdot z^2$$
However, as you can see at the flanges you get more area at the same (more or less distance) (e.g. $z=\pm\frac{A}{2}$ for the image below).

So, although the stresses and strain (in the linear region) follow a linear distribution wrt to the distance from the neutral axis, you can see that the actual bending moment eventually is proportional
- to the square of the distance z ,
- to the area at distance z from the neutral axis.
Therefore, because the flanges have a significant area, further away from the neutral axis they have greater contribution to the bending moment which resists bending.