I assume you want $A=\mathbb Z$ here since $[M:N]\in\mathbb Z\cup\{\infty\}$ while $\det(a_{ij})\in A$. Furthermore, we will need to assume that the $n_j$'s are linearly independent (otherwise $[M:N]=\infty$ while $\det(a_{ij})=0$), i.e. we are assuming that $N$ is free (submodules of free $\mathbb Z$-modules are always free) of rank $n$ over $\mathbb Z$. With these assumptions, my answer is basically an expanded version of what you read in Milne's notes.
Since $N\subset M$ are both free of rank $n$ over $\mathbb Z$, the structure theorem for modules over a PID tells us that there exists a basis $m_1',\dots,m_n'$ of $M$ along with scalers $c_1,\dots,c_n\in\mathbb Z$ such that $c_1m_1',\dots,c_nm_n'\in M$ form a basis for $N\subset M$. Using these bases, the inclusion $N\hookrightarrow M$ is given by the diagonal matrix $\mathrm{diag}(c_1,\dots,c_n)$. Since $(a_{ij})$ represents this inclusion in another basis, we have
$$\prod_{i=1}^nc_i=\det\mathrm{diag}(c_1,\dots,c_n)=\det(a_{ij}).$$
At the same time, with this nice choice of bases for $M,N$, one sees that
$$\frac MN\cong\prod_{i=1}^n\frac{\mathbb Z}{c_i\mathbb Z}$$
so also
$$\left|\frac MN\right|=\prod_{i=1}^nc_i=\det(a_{ij}).$$
All that remains is to explain where we get this choice of bases from. We have $N\subset M$ free $\mathbb Z$-modules of the same finite rank $n$. We want to construct a basis $e_1,\dots,e_n$ of $M$ along with elements $a_1,\dots,a_n\in\mathbb Z$ such that $a_1e_1,\dots,a_ne_n$ is a basis for $N$.
One argues by induction on $n=\mathrm{rank}(M)$. For any homomorphism $\lambda:M\to\mathbb Z$, let $I_\lambda=\lambda(N)\subset\mathbb Z$ be the ideal capaturing the image of $\lambda\vert_N$. Choose $\lambda_1$ such that $I_{\lambda_1}\subset\mathbb Z$ is maximal among $\{I_\lambda\}_{\lambda:M\to\mathbb Z}$, write $I_{\lambda_1}=a_1\mathbb Z$, and choose $x_1\in N$ such that $\lambda_1(x_1)=a_1$. Now, briefly consider any basis $f_1,\dots,f_n$ of $M$ with dual basis $\newcommand\dual[1]{#1^\vee}\dual f_1,\dots,\dual f_n$ of $\dual M=\mathrm{Hom}(M,\mathbb Z)$ (so $\dual f_i(f_j)=\delta_{ij}$); if we write $x_1=\sum_{i=1}^nb_if_i$, then
$$b_i=\dual f_i(x_1)\in(a_1)$$
by maximality of $I_{\lambda_1}$, so $a_1\mid b_i$ for all $i$. Hence, $x_1=a_1e_1$ for some $e_1\in M$ (in particular, $\lambda_1(e_1)=1$). We claim now that
$$M=\mathbb Ze_1\oplus\ker\lambda_1.$$
Indeed, for $m\in M$, one has $m=\lambda_1(m)e_1+(m-\lambda_1(m)e_1)\in\mathbb Ze_1+\ker\lambda_1.$ At the same time, if $m\in\mathbb Ze_1\cap\ker\lambda_1$, then $m=\lambda_1(m)e_1=0$ which shows $M=\mathbb Ze_1\oplus\ker\lambda_1$ as claimed. The inductive hypothesis then gives a basis $e_2,\dots,e_n$ of $\ker\lambda_1$ along with $a_2,\dots,a_n\in\mathbb Z$ such that $a_2e_2,\dots,a_ne_n$ is a basis of $N\cap\ker\lambda_1$, and so we win. Our desired bases are $e_1,e_2,\dots,e_n$ and $a_1e_1,a_2,e_2,\dots,a_ne_n.$
