Matrices are useful for proving statements like
The ratio between the areas of a parallelogram and the quadrilateral formed by joining their midpoints is $2$.
The ratio between the volumes of a parallelepiped and the cuboctahedron formed by joining its midpoints is $\frac{6}{5}$.
, because one can apply some matrix $T^{-1}$ to the vertices of the shapes to turn them into cubes/squares (i.e. via an affine transformation). Then, solve the problem in this now simpler case. As the ratios of areas are unchanged under $T$ acting on a shape, the conclusion holds in the general case.

Now, is there a mathematical object, similar to matrices, that deals with non-affine transformations? For example, could similar method be used to solve the first problem for any arbitrary quadrilateral?
The theory of non-affine transformations is probably going to be much more complex than in the affine case, as the ratios of areas are not preserved, with some areas ballooning out and some squashed (although I conjecture a non-affine transformation acts locally as an affine transformation). However, I ask in hope that there exist theorems that could help simplify the situation. I have heard of the shoelace theorem, but am not sure how generally it could be applied to solve general problems.
Using the shoelace formula
As I stated in my other answer that it is not particularly useful to proove the statements you quoted, I'll here focus on a viable solution for that. The key is the shoelace formula you already mentioned. It states that the area of a polygon is
\begin{multline*} A(p_0,p_1,p_2,\dots,p_{n-1})= \frac12\sum_{i=0}^{n-1} \begin{vmatrix} x_i & x_{(i+1)\bmod n} \\ y_i & y_{(i+1)\bmod n} \\ \end{vmatrix} =\\ \tfrac12\bigl( (x_0y_1-x_1y_0) + (x_1y_2-x_2y_1)+\dots+ %(x_{n-2}y_{n-1}-x_{n-1}y_{n-2})+ (x_{n-1}y_0-x_0y_{n-1})\bigr) \end{multline*}
So in your case of $n=4$ the area will be $$A=\tfrac12\bigl(x_0y_1-x_1y_0+x_1y_2-x_2y_1+x_2y_3-x_3y_2+x_3y_0-x_0y_3\bigr)$$
Your midpoints are \begin{align*} x_0' &= \frac{x_0+x_1}2 & x_1' &= \frac{x_1+x_2}2 & x_2' &= \frac{x_2+x_3}2 & x_3' &= \frac{x_3+x_0}2 \\ y_0' &= \frac{y_0+y_1}2 & y_1' &= \frac{y_1+y_2}2 & y_2' &= \frac{y_2+y_3}2 & y_3' &= \frac{y_3+y_0}2 \end{align*}
Plug these new values into the above formula, simplify things a bit, and you will obtain a result which is just half the one you got before.
But in all of this, the only matrices involved were those $2\times2$ matrices used for the determinants.