It is easy to derive the equations for moments of inertia for masses of different figures, by employing the following generalized integral:
$$I=\int_0^Mr^2 dm$$
Which is based on the fact that the moment of inertia of a point mass is $I=mr^2$. But how do we even know this? If moment of inertia really is the rotational equivalent of mass, would it not just be $I=mr$, similar to how linear velocity is just angular velocity times the radius?
At the end of the day, all these quantities are just definitions. So the only thing one can do is show that these definitions are natural.
First I will motivate why angular momentum is important. Noethers' theorem says that for every symmetry of equations of motions there is some conserved quantity. Newtons' equations of motions are translationally invariant and this turns out be equivalent to momentum being conserved. This makes momentum an important quantity.
Similarly, Newtons' laws of motion are invariant under rotations and this leads to angular momentum, $\vec L=\vec r\times\vec p$, to be conserved. It is also natural to study the time derivative of this quantity, the torque, similar to how we study forces as the time derivative of momentum.
Next, we will see that based on parallels with normal momentum it is natural to define $I=mr^2$. We can split regular momentum into mass and velocity. Velocity is the time derivative of the position $\vec r$, which I will call the 'dynamical variable'. $$\vec p=\overbrace{m}^{\text{inertia}}\underbrace{\vec v}_{\text{dynamical}\\\text{ variable}}$$
Angular momentum is defined as $\vec L=\vec r\times\vec p$. Angular momentum is associated with rotations, so it is natural to look at the angle as the dynamical variable. If we assume circular motion for a second we can take the momentum to be perpendicular to $\vec r$. \begin{align} L&=rp\\ &=mrv \end{align} We can also write the velocity as $v=\omega r$, where $\omega=\dot\theta$. This gives us a similar split in an inertial term and a time derivative of the dynamical variable $\theta$. $$ \vec L=\overbrace{mr^2}^{\text{inertia}}\underbrace{\omega}_{\text{dynamical}\\\text{ variable}}$$
This last equation can be interpreted as follows: for something which has high $mr^2$, it is hard to change the rotation rate. We can argue that an $I=mr$ says the same, but because of the equations of motion $I=mr^2$ is more natural.
We can see the factor $r$ as a conversion factor for going from angles to displacements. Similar to how $s=r\theta$, where $s$ the arc length. In the torque equation $\tau=I\alpha$ we get one factor of $r$ for converting from angular acceleration to acceleration and we get one factor for converting from torque to force. See the following (very) handwavy explanation. \begin{align} \tau&=I\alpha\\ rF&=(mr^2)(\frac{a}{r})\\ F&=ma \end{align}
First things first the equation given in OP is not correct and should rather be something like $$ I = \int_\Omega \rho(\vec r)\, \vec r_{\text{othogonal}}^2\text{d}^3r, $$ where $\Omega\subset \mathbb R ^3$ is the set of points where $\rho \neq 0$ and $\vec r_{\text{othogonal}}$ is the distance vector between the point $\vec r$ and the axis of rotation which is orthogonal on that axis.
Assume we are having some collection of point particles of masses $m_i$ and positions $\vec r_i$. Their relative positions are fixed, for example by massles rigid rods. We want to find the rotational equivalent of inertia. So let us write down the kinetic energy of such a collection of particles around some common axis of rotation $\vec\omega$ \begin{align*} E_\text{kin.} &= \frac 1 2 \sum_{i=1}^N m_i\,\vec{v}_i^2 \\ &= \frac 1 2 \sum_{i=1}^N m_i\,(\vec{r}_i \times \vec{\omega})^2\\ &= \frac 1 2 \sum_{i=1}^N m_i\,\left[(\vec{r}_i \cdot \vec{r}_i) (\vec{\omega}\cdot\vec{\omega}) - (\vec{\omega}\cdot\vec{r}_i) (\vec{r_i}\cdot\vec{\omega})\right] \\ &= \frac 1 2 \sum_{i=1}^N m_i\,\left[\vec{\omega}^T (\vec{r}_i \cdot \vec{r}_i) \vec{\omega} - \vec{\omega}^T\vec{r}_i\vec{r}_i^T\vec{\omega}\right] \\ &= \frac 1 2 \sum_{i=1}^N m_i\,\vec{\omega}^T \left[(\vec{r}_i \cdot \vec{r}_i)\mathbb{I} - \vec{r}_i\vec{r}_i^T\right] \vec{\omega} \\ &= \frac 1 2 \vec{\omega}^T \left\{\sum_{i=1}^N m_i\, \left[(\vec{r}_i \cdot \vec{r}_i)\mathbb{I} - \vec{r}_i\vec{r}_i^T\right]\right\} \vec{\omega}\\ &= \frac 1 2 \vec{\omega}^T \Theta\ \vec{\omega}, \end{align*} where $\Theta$ is the tensor of inertia. This simplifies if we just consider one particle and for example rotation around the $z$-axis. Then by basic algebra we arrive at \begin{align} E_\text{kin.} &= \frac 1 2 m\vec r ^2 \omega^2-\frac 1 2 mr_z^2 \omega^2\\ &=\frac m 2(\vec r^2 -r_z^2) \omega^2 \end{align} Since we are interested in the analogy between translation and rotation we compare this expression with the kinetic energy of a free translating particle of mass $m$ in 1 dimension $$ E_\text{kin.}=\frac m 2 v^2. $$ Since we do have the conceptual correspondence between kinematic variables of rotation and translation like $$ \text{Rotation}\leftrightarrow\text{Translation} $$ $$ \vec\omega\leftrightarrow\vec v $$ we could compare $$m\leftrightarrow m(\vec r² - r_z^2)$$ as quantities defining inertia of translation and rotation. And in fact we loosely find this correspondence again in the general equations of motion of a rotating rigid body and a translating mass.
Your statment is incorrect because the MMOI tensor for a point particle located at $\boldsymbol{r}$ is
$$ {\bf I} = \left( (\boldsymbol{r} \cdot \boldsymbol{r}) {\bf 1} - \boldsymbol{r} \odot \boldsymbol{r} \right) m $$
where $\cdot$ is the vector dot product, $\odot$ is the vector outer product and ${\bf 1}$ is the identity matrix.
In integral form for an extended body the above is
$$ {\bf I} = \int \left( (\boldsymbol{r} \cdot \boldsymbol{r}) {\bf 1} - \boldsymbol{r} \odot \boldsymbol{r} \right) {\rm d}m $$
In terms of components, with $\boldsymbol{r} = \pmatrix{x \\ y \\z}$ the above is
$${\bf I} = \int \begin{vmatrix} y^2+z^2 & -x y & - x z \\ -x & x^2+z^2 & -y z \\ -x z & -y z & x^2+y^2 \end{vmatrix} {\rm d}m$$
In our statment, you are looking only at a single diagonal term of the tensor, such as ${\bf I}_{zz} = \int ( x^2+y^2) {\rm d}m$.
To derive the above take the body rotating about the center of mass with some arbitrary $\boldsymbol{\omega}$ and calculate the angular momentum vector $\boldsymbol{L}$. Mass moment of inertia is defined as the tensor that transforms $\boldsymbol{\omega}$ into $\boldsymbol{L}$.
$$\begin{aligned}\boldsymbol{L} & =\int\boldsymbol{r}\times{\rm d}\boldsymbol{p}\\ & =\int\boldsymbol{r}\times\boldsymbol{v}{\rm d}m\\ & =\int\boldsymbol{r}\times\left(\boldsymbol{\omega}\times\boldsymbol{r}\right){\rm d}m\\ & =\int\left(\boldsymbol{\omega}\left(\boldsymbol{r}\cdot\boldsymbol{r}\right)-\boldsymbol{r}\left(\boldsymbol{r}\cdot\boldsymbol{\omega}\right)\right){\rm d}m\\ & =\left(\int\left(\left(\boldsymbol{r}\cdot\boldsymbol{r}\right){\bf 1}-\boldsymbol{r}\odot\boldsymbol{r}\right){\rm d}m\right)\boldsymbol{\omega} \\ & = {\bf I}\, \boldsymbol{\omega} \end{aligned}$$
using the vector triple product identity $a\times(b \times c) = b(a \cdot c) - c ( a \cdot b)$, and the linear algebra identity $a ( b \cdot c) = a ( b^\intercal c) = (a b^\intercal) c = (a \odot b) c$.
Following @AlmostClueless notation, if we assume that the position $\vec{r}^\prime$ of the point mass is perpendicular to the axis of rotation and that its density is given by $\rho(\vec{r}) = m \delta(\vec{r}-\vec{r}^\prime)$, where $\delta(\vec{r}-\vec{r}^\prime)$ is the Dirac delta, then we can write
$$ I = \int_\Omega m \delta(\vec{r}-\vec{r}^\prime) \rVert \vec{r} \rVert^2 \text{d}^3 r = m \rVert \vec{r}^\prime \rVert^2 = m r^{\prime 2} $$
where we integrate over all possible vectors that are perpendicular of the axis of rotation.
